# Effect of Ambipolar Diffusion on the Non-linear Evolution of Magnetorotational Instability in Weakly Ionized Disks

###### Abstract

We study the role of ambipolar diffusion (AD) on the non-linear evolution of the MRI in protoplanetary disks using the strong coupling limit, which applies when the electron recombination time is much shorter than the orbital time. The effect of AD in this limit is characterized by the dimensionless number , the frequency of which neutral particles collide with ions normalized to the orbital frequency. We perform three-dimensional unstratified shearing-box simulations of the MRI over a wide range of as well as different magnetic field strengths and geometries. The saturation level of the MRI turbulence depends on the magnetic geometry and increases with the net magnetic flux. There is an upper limit to the net flux for sustained turbulence, corresponding to the requirement that the most unstable vertical wavelength be less than the disk scale height. Correspondingly, at a given , there exists a maximum value of the turbulent stress . For , the largest stress is associated with a field geometry that has both net vertical and toroidal flux. In this case, we confirm the results of linear analyses that show the fastest growing mode has a non-zero radial wave number with growth rate exceeding the pure vertical field case. We find there is a very tight correlation between the turbulent stress and the plasma at the saturated state of the MRI turbulence regardless of field geometry, and rapidly decreases with decreasing . In particular, we find for and for .

###### Subject headings:

magnetohydrodynamics — instabilities — methods: numerical — planetary systems: protoplanetary disks — turbulence## 1. Introduction

One of the most fundamental questions in accretion disk dynamics is how the disk transports angular momentum and accretes to the central object. The magnetorotational instability (MRI, Balbus & Hawley, 1991) is widely considered to be the most likely mechanism for the transport process. The non-linear evolution of the MRI under ideal MHD conditions has been studied extensively using both local (Hawley et al., 1995; Stone et al., 1996; Miller & Stone, 2000) and global (Armitage, 1998; Hawley, 2000, 2001; Fromang & Nelson, 2006) numerical simulations. It has been found that MRI generates vigorous MHD turbulence and produce efficient outward transport of angular momentum whose rate is compatible with observations. However, accretion disks in some systems are only partially ionized, and non-ideal MHD effects need to be taken into account. In particular, most regions of the protoplanetary disks (PPDs) are too cold for sufficient thermal ionization, and effective ionization may be achieved only in the disk surface layer due to external ionization sources such as X-ray radiation from the central star and cosmic ray ionization (Gammie, 1996; Stone et al., 2000). Non-ideal MHD effects reflect the incomplete coupling between the disk material and the magnetic field, and substantially affect the growth and saturation of the MRI.

There are three non-ideal MHD effects as manifested in the generalized Ohms’s law, namely the Ohmic resistivity, Hall effect and ambipolar diffusion (AD), with three different regimes associated with the relative importance of these terms. In general, the Ohmic term dominates at high density and very low ionization, the AD term dominates in the opposite limit, and the Hall term is important in between. So far the majority of studies have been focused on the Ohmic regime. In this case, MRI is damped when the Elsasser number falls below order unity (Jin, 1996; Turner et al., 2007), where is the Alfvén velocity in the vertical direction, is the Ohmic resistivity, is the angular frequency of Keplerian rotation. Another often quoted criterion is the magnetic Reynolds number Re for MRI to be self-sustained (where is the sound speed), which has the advantage of being independent of the magnetic field (Fleming et al., 2000). Ohmic resistivity has been used extensively to model the layered accretion in PPDs, where the surface layer of the disk is sufficiently ionized to couple to the magnetic field and drive the MRI turbulence, while the midplane region is too poorly ionized and “dead” (Fleming & Stone, 2003; Turner et al., 2007; Turner & Sano, 2008; Ilgner & Nelson, 2008; Oishi & Mac Low, 2009).

The importance of Hall and AD terms in PPDs has been studied in a number of theoretical works, but relatively little attention has been paid to numerical simulations of the non-linear regime. Linear analysis of the MRI in the Hall regime have been performed by Wardle (1999) and Balbus & Terquem (2001). The growth rate is strongly affected by the Hall term and depends on the sign of . Nevertheless, numerical simulations including both the Hall and Ohmic terms (where the Ohmic term dominates) showed that the Hall term does not strongly affect the saturation amplitude of MRI (Sano & Stone, 2002a, b). It is yet to study the behavior of MRI in the regime where Hall effect dominates over other terms, and to include the Hall term in the more realistic vertically stratified simulations.

The relative motion between ions and neutrals leads to AD. AD is ideally studied using the two-fluid approach, where the ions and neutrals are treated as separate fluids, coupled by the ion-neutral drag via collisions. Moreover, ion and neutral densities are affected by the ionization and recombination processes. Most analytical studies in the linear regime adopt the Boussinesq approximation where ion and neutral densities are kept constant (Blaes & Balbus, 1994; Kunz & Balbus, 2004; Desch, 2004). These studies show that the growth of MRI is suppressed when the collision frequency of a neutral with the ions falls below the orbital frequency. In the mean time, when both vertical and azimuthal field is present, unstable modes always exist due to the effect of AD, and these unstable modes require non-zero radial wavenumbers (Kunz & Balbus, 2004; Desch, 2004). Blaes & Balbus (1994) also studied the effect of ionization and recombination with compressibility (for vertical propagating waves), and found that the presence of azimuthal and radial field allows the coupling to acoustic and ionization modes, and the azimuthal field tends to stabilize the flow when the recombination time is not too long compared with dynamical time.

The effect AD on the MRI in the non-linear regime was first studied by Mac Low et al. (1995). They implemented and tested AD in the “strong coupling” limit (see below) in the ZEUS code and performed simulations with net vertical flux for various ion-neutral coupling strengths. Their results confirmed the stability analysis of Blaes & Balbus (1994), but their simulations are only two-dimensional and did not follow the evolution much beyond the linear stage. In another study, Brandenburg et al. (1995) included the effect of AD (also in the strong coupling limit) in their three-dimensional simulations of a local, vertically stratified disk. They found that turbulence remains self-sustained in a case where AD time is long compared with orbital time, although reduced in strength, and in another case where the AD time was set comparable to , turbulence decayed.

A systematic study on the non-linear evolution of MRI with AD is done by Hawley & Stone (1998) (hereafter HS) using three-dimensional (3D) numerical simulations. They used the two-fluid approach without considering the ionization-recombination processes, therefore ions and neutrals obey their own continuity equations. Both net-vertical and net-toroidal magnetic configurations were considered. They found that the system behaves like fully-ionized gas when the ion-neutral collision frequency is greater than , while ions and neutrals behave independently when the collision frequency fall below . The amplitude of magnetic field at saturation is proportional to the ion density when it is much smaller than the neutral density. The two-fluid approach adopted by HS is valid when the recombination time scale is long compared with the dynamical time. However, AD in PPDs is in general in the “strong coupling” limit (Shu, 1991). Two conditions must be satisfied in this limit: 1. The ion density is negligible compared with the neutral density . 2. The electron recombination time must be much smaller than the orbital frequency . In this limit, the ion density is purely determined by the ionization-recombination equilibrium with the neutrals, and the two-fluid formulation is simplified into a single-fluid formalism (for the neutrals). In PPDs, condition 1 is always satisfied, and we will show in a companion paper (Bai, 2011) that condition 2 is almost always satisfied.

In this paper, we conduct three-dimensional (3D) local shearing-box simulations to explore the effect of AD on the non-linear evolution of MRI in the strong coupling limit. This is conceptually different from the simulations performed by HS in that the ion density does not obey continuity equation, and is set by the neutral density due to chemical equilibrium. Effectively, this allows the coupling of the MRI with acoustic and ionization modes, which leads to more complicated interactions (Blaes & Balbus, 1994). Moreover, our simulations correspond to the limit where the ion density is negligibly small (i.e., in HS), which is difficult for two-fluid simulations due to the stiffness of the equations. In the strong coupling limit, there is only one controlling parameter, namely, the ion-neutral collision frequency . We perform three sets of simulations with net vertical flux, net toroidal flux and both. In each group of runs, we systematically vary as well as the strength of the net field. Our main goal is to study the conditions under which MRI turbulence can be self-sustained or is suppressed due to AD. In addition, we study the properties of the MRI turbulence in the AD dominated regime.

This paper is structured as follows. In Section 2 we provide the formulation of AD in the strong coupling limit, describe the numerical method and code test problems. A series of numerical simulations on the non-linear evolution of MRI with AD are presented and analyzed in Section 3, from which we discuss the condition under which MRI turbulence is sustained or suppressed as well as the properties of the MRI turbulence in AD dominated regime. We conclude and briefly discuss various implications In Section 5.

## 2. Formulation and Numerical Method

### 2.1. Ambipolar Diffusion in Weakly Ionized Plasma

In weakly ionized plasma, the inertia of the ionized species is negligible. The gas dynamics can thus be described by single-fluid equations for the neutrals, modified by non-ideal MHD effects. The effect of ambipolar diffusion (AD) derives from the relative motion between the ions and the neutrals. Since the ion inertia is negligible, the ion velocity is determined by the balance between the Lorentz force and the ion-neutral collisional drag

(1) |

where , and are the ion mass, number density and velocity respectively, , are the density and velocity of the neutrals, , are the current density and magnetic field vectors, , with begin the rate coefficient for momentum transfer between the ions and the neutrals and the neutral mass. The magnetic field is effectively carried by the ions, thus the induction equation now reads

(2) |

where is the ion mass density, is the Alfvén velocity, is the component of the current density that is perpendicular to the direction of the magnetic field. The second term on the right hand side is the AD term, from which the ambipolar diffusivity can be defined as . The rest of the fluid equations remain the same as the ideal MHD equations. In particular, the momentum equation is unchanged since the neutrals effectively feel the magnetic tension and pressure by the ion-neutral collisions.

The above derivation is strictly valid when electrons and ions are the only charge
carriers, and all ions have the same mass. Nevertheless, the effect of multiple ion
and charged grain species can be combined into an effective AD coefficient^{1}^{1}1A
more generalized derivation of all non-ideal MHD terms in the induction equation is
given and discussed in Bai (2011) (and see also Wardle (1999, 2007)), which
also includes Ohmic resistivity and the Hall term.. AD becomes the dominant
non-ideal MHD effect when the gyro-frequency of both the ions and the electrons are
higher than their collision frequency with the neutrals (i.e., both ions and electrons are
coupled to the magnetic field). In practice, this corresponds to regions with low density
and strong magnetic field.

The effect of AD in rotating disks (with angular frequency ) is characterized by the parameter (Chiang & Murray-Clay, 2007):

(3) |

which is the number of effective collisions for a neutral molecule/atom with the ions in a dynamical time . Physically, we see that , which is equivalent to the Elsasser number for Ohmic resistivity . Therefore, measures the ratio of the AD time scale over the critical wavelength and the dynamical time scale.

In this paper, we study AD in the strong coupling limit (Shu, 1991): the electron (ion) recombination time is much shorter than the dynamical time so that the ion density is determined by the local thermodynamical quantities (density and temperature). The strong coupling limit is widely applicable in protoplanetary disks, and we show in our companion paper (Bai, 2011) that the recombination time is typically at least one order of magnitude smaller than the dynamical time, even in the disk coronal regions. Using the strong coupling limit, we assume the ion density depends on the neutral density in the form of

(4) |

where and are the reference density of the ions and the
neutrals. A simple-minded calculation of the ionization-recombination equilibrium
gives , which yields . In the equation is the
ionization rate, of the order s, is the effective rate coefficient
for recombination, of the order cm s g in the absence
of grains (Blaes & Balbus, 1994). In reality, the recombination process is complicated
by a complex network of gas phase and grain phase chemical reactions, and we
can address the complications by exploring different values of . Being a
two-body process in general, and one expects .^{2}^{2}2For example,
Figure 1a and Figure 3 of Bai & Goodman (2009) illustrate the dependence of
electron abundance on the gas density in various chemistry
models with and without dust grains, and holds in essentially all
circumstances. In fact, our simulations show that the properties of the MRI turbulence
is insensitive to (see Section 3.1.3).

### 2.2. Numerical Method

We use Athena, a higher-order Godunov MHD code with constrained transport technique to enforce the divergence-free constraint on the magnetic field (Gardiner & Stone, 2005, 2008; Stone et al., 2008) for all calculations presented in this paper. Non-ideal MHD terms including Ohmic resistivity (Davis et al., 2010), the Hall term (in progress) and AD (this paper) have been developed for Athena. We consider a local patch of a Keplerian disk using the standard shearing-box formalism (Goldreich & Lynden-Bell, 1965), which adopts a local reference frame at a fiducial radius corotating with the disk at orbital frequency . In this frame, we write the MHD equations with AD in a Cartesian coordinate system as

(5) |

(6) |

(7) |

where is the total stress tensor

(8) |

is the identity tensor, is the gas pressure. are unit vectors pointing to the radial, azimuthal and vertical directions respectively, where is along the direction. Disk vertical gravity is ignored and our simulations are vertically unstratified. We use an isothermal equation of state , where is the isothermal sound speed. Periodic boundary conditions are used in the azimuthal and vertical directions, while the radial boundary conditions are shearing periodic as usual.

An orbital advection scheme (Masset, 2000; Johnson et al., 2008) has been implemented in Athena (Stone & Gardiner, 2010). It splits the dynamical equations into two systems, one of which corresponds to linear advection operator with background flow velocity and the other evolves only velocity fluctuations. The orbital advection scheme not only accelerates the calculation in large box size by admitting larger time steps, but also makes the calculation more accurate.

The last term on the right hand side of equation (7) represents the AD term, where is the effective ion density, approximated by equation (4). Non-ideal MHD terms (e.g., the AD term) are implemented in Athena in an operator-split way. Over one time step, one first solves the induction equation using only the non-ideal MHD terms, and then solves ideal MHD equations. The induction equation with AD term is a diffusion equation and is evolved by a fully explicit forward-Euler method. The divergence free condition is maintained using constrained transport: we calculate the AD electromotive force defined at cell edges by interpolating the magnetic field and current density to such locations. The method is stable, with time step constrained to be proportional to grid size squared. The overall method is first order accurate in time.

Numerical method including the AD term has been frequently implemented in the framework of two-fluid models (e.g., Toth, 1994; Stone, 1997; Smith & Mac Low, 1997; Falle, 2003; Li et al., 2006; O’Sullivan & Downes, 2006, 2007; Tilley & Balsara, 2008), with the main applications on star formation, turbulence and shocks in the interstellar medium. However, single-fluid models for AD in the strong coupling limit relevant to weakly ionized disks (Brandenburg et al., 1995; Mac Low et al., 1995) are relatively less studied numerically. Recently, Choi et al. (2009) described an explicit scheme for incorporating AD in the strong coupling limit in an MHD code that is second order accurate, and use super time-stepping to accelerate the calculation when the AD coefficient is large. As we find MRI is suppressed at moderately large AD coefficients (see Section 3), super time-stepping is not needed for the purpose of this paper.

### 2.3. Code Tests

In this subsection, we conduct two test problems to examine the performance of our implementation of the AD term. For all these tests, we consider non-rotating system and turn off the shearing box source terms on the right hand side of equation (6).

#### 2.3.1 Isothermal C-type Shock Test

The effect of ambipolar diffusion is best manifested in C-type shocks (Draine, 1980), which is a shock with continuous transitions consequent of the AD. For the purpose of the code test, here we consider the isothermal C-type test by Mac Low et al. (1995), which has become a standard test problem for AD.

We consider steady-state () shock and work in the shock frame, with upstream gas density moving at velocity . The upstream gas is threaded by a uniform magnetic field that lies at an angle to the velocity. Let be in the direction, and be in the plane. For a continuous shock, the jump conditions reduce to . Assuming the ion density is constant (), the equations that describe the C-type shock read

(9) |

(10) |

(11) |

(12) |

Note that is constant.

The shock is characterized by three dimensionless parameters: the sonic Mach number , the Alfvén Mach number (where ), and the angle of the magnetic field with the upstream flow. The characteristic length scale of the problem is given by . We further define , and . After some algebra, we arrive at a dimensionless first order differential equation for (Mac Low et al., 1995)

(13) |

One can numerically integrate this ordinary differential equation to obtain the C-type shock profile. In Figure 1 we show a semi-analytical solution for , and obtained by using a 4th order Runge-Kutta method.

To use this solution as a code test, the shock is set to be aligned with the grid in the
direction. We use outflow boundary conditions in this direction. In
multi-dimensional tests, periodic boundary conditions are used in other directions.
The shock solution should be stationary (i.e., a standing shock), thus we evolve
the solution for sufficiently long time () and compare to the initial conditions.
In Figure 1, we further show the absolute error of the shock
profile compared with the semi-analytic solution. Since the shock is
grid-aligned, 1D, 2D and 3D tests essentially produce the same result. In our tests,
the grid resolution is chosen to be 2 and 4 cells per .^{3}^{3}3In comparison,
Mac Low et al. (1995) achieved comparable accuracy as ours using 5 and 10 cells
per , Choi et al. (2009) achieved similar or better accuracy at 6.4 cells per .
We see that our code very accurately resolves the structure of the C-type shock using
only a few cells per . The main source of the error lie in the region where density and
velocity profile vary quickly. In our comparison, the position of the shock is fixed at
the initial place, while in reality, the shock position can shift slightly during numerical
relaxation.

#### 2.3.2 Damping of MHD Waves

Linear MHD waves are damped due to AD. Because the exact eigenvectors in the AD regime can be very complicated, we initialize the problem with ideal MHD wave eigenvectors and measure the damping rate. This means the initial conditions are a linear superposition of more than one eigen-mode, but the averaged damping rate should approach the analytical value for a single mode as long as the AD coefficient is sufficiently small.

The analytical damping rate for various MHD waves due to AD can be found in and/or derived from Balsara (1996), as we summarize below. The damping rate of the Alfvén wave is given by the solution of

(14) |

where , is the Alfvén velocity, is the angle between the magnetic field and the wave vector. The damping rate of fast and slow waves can be obtained by solving the quadratic equation

(15) |

where and are the fast and slow magnetosonic speeds.

When , the damping rate is small and can be found by expanding the ideal MHD dispersion relation to powers of , and to the first order, we find for the damping rate of the Alfvv́en wave

(16) |

The damping rate for fast and slow magnetosonic waves are

(17) |

We perform the linear wave damping test in 1D, 2D and 3D. In 1D, the wave is grid-aligned, whose wave length equals in code unit. In 2D and 3D test problems, the wave vectors are not grid-aligned, with box sizes chosen such that the wave length is also [in 2D, the box size is () and in 3D, the box size is (3, 1.5, 1.5)]. We use isothermal equation of state with . The background gas density is . As before, we choose . In 1D, the wave vector is along the direction, and the adopted magnetic field is , and . In 2D and 3D the background magnetic field vector is rotated with the wave vector accordingly while keeping the same as in 1D (and a vector potential is used to initialize the wave in order to preserve the divergence free condition). From the above we get , and , therefore we find , and .

In practice, we adopt and run our simulation to . By default, the grid size is in 1D, in 2D, and in 3D. Accounting for the box size in each dimension, the effective resolution, characterized by number of cells per wavelength, is , and in 1D, 2D and 3D respectively. The results are shown in Figure 2. From left to right, we show the damping curve from 1D, 2D and 3D simulations in the solid lines, where black, blue and red lines label Alfvén, fast and slow MHD waves. Dashed lines show the theoretical damping curve. We see that the numerical damping rate matches very well with the theoretical damping rate. In the 3D runs, the damping rate for the Alfvén wave is slightly faster than expected, but this may be because the effective resolution is less.

We have also run the simulations with double and half the resolution. With double resolution, the numerical damping curves in 1D, 2D and 3D cases almost match exactly the analytical damping curves (besides some small oscillations due to the initial conditions). At half resolution, however, the numerical damping rate deviates substantially (about to at ). These results indicate that at least cells per wavelength is need to accurately capture the effect of AD.

## 3. Simulations and Results

In this section we describe three groups of simulations and study the effect of AD on the non-linear evolution of the MRI. All our simulations are vertically unstratified by ignoring the disk vertical gravity with fixed box height to be one disk scale height . We initialize our simulations with Keplerian velocity and seed density perturbations of of the background density . We consider three different magnetic field geometries (net vertical flux, net toroidal flux, both vertical and toroidal flux) as described in the following three subsections. Since all our simulations contain net magnetic flux, they are not subject to the issue of convergence found by (Fromang & Papaloizou, 2007) in zero net-flux simulations, and numerical convergence is confirmed in our test simulations (and see Section 3.2.2 for the case of net toroidal flux). For relatively small AD coefficient (large ), MRI grows quickly from the seed perturbations and saturates into turbulence; when the effect of AD is strong, however, MRI does not grow from our seed perturbations. In such cases, we initialize the simulations from a turbulent state which is obtained from simulations with relatively large (see individual subsections for details).

The most important diagnostics are the volume averaged (normalized) Reynolds stress, defined as

(18) |

where the over bar indicates volume averaging, is the azimuthal velocity with the Keplerian velocity subtracted, and the volume averaged (normalized) Maxwell stress, defined as

(19) |

The total stress, namely, the parameter (Shakura & Sunyaev, 1973) is . We also monitor the kinetic and magnetic energy density, which characterize the strength of the MRI turbulence.

The main purpose of this study is to identify the criterion under which sustained turbulence generated by the MRI can be maintained. However, the term “sustained turbulence” is a somewhat ambiguous concept. In the context of 3D shearing box simulations, small-amplitude oscillations left from decayed hydrodynamical turbulence is present (Shen et al., 2006). Such oscillations produce Reynolds stress on the order of with kinetic energy density on the order of , both in normalized unit . Such oscillations are likely to associate with linear modes in the shearing sheet with the origin of acoustic and/or nearly incompressible inertia waves (Balbus, 2003). Being a conservative Godunov MHD code, the Athena code preserves the amplitude of these waves without much damping. Therefore, throughout this paper, the level of turbulence we are interested in are those whose time and volume averaged kinetic energy density is on the order of or higher, and/or whose total stress is no less than . Meanwhile, analysis of all our simulations show that the threshold where the MRI turbulence can be marginally self-sustained is roughly at the same level. (see Section 3.2.4 for further discussion).

Our simulations are run for at least orbits (). A period of
orbits is sufficiently long for the MRI to saturate from initial growth, which typically occurs
in orbits, or for the restart runs to reach a steady state, which typically occurs in
orbits. Our time averaged quantities are mostly taken from after about orbits
(since ) unless otherwise noted. Although a time average over
orbits () is relatively short, it is sufficient for our purpose to judge whether
MRI turbulence can be self-sustained^{4}^{4}4Winters et al. (2003) found that more than
a few hundred orbits as are required to accurately measure the properties of the MRI
turbulence in ideal MHD. This conclusion is based simulations with radial box size being
, while the radial box size in about half of our simulations is , which reduces the
time fluctuations. Also, our Figures 3 and 8 show that the
fluctuations in the Maxwell stress are less severe in the presence of AD than in the ideal
MHD case.. Many of our simulations are run for orbits or longer where better statistics
on the turbulence properties can be obtained.

### 3.1. Net Vertical Flux Simulations

In the first group of simulations, we choose the initial field configuration to be uniform along the vertical axis , characterized by the plasma , where is the background pressure and is the initial field strength. The vertical flux is conserved numerically by remapping the toroidal component of the magnetic field in the ghost zones of the radial boundaries (see Section 4 of Stone & Gardiner (2010) for details). For all simulations, we fix the box size to be in the radial, azimuthal and vertical dimensions, with fixed grid resolution at cells per . We have chosen a relatively large radial box size (), as suggested by Pessah & Goodman (2009), which is needed to properly capture the parasitic modes to break the channel mode into turbulence. It also help substantially reduce the intermittence of the MRI turbulence (Bodo et al., 2008). We note that for local unstratified net vertical flux MRI simulations without explicit dissipation, turbulence properties converge at about cells per (Hawley et al., 1995). The grid resolution in our simulations is two times higher, thus we expect numerical convergence.

All our net vertical flux simulations are listed in Table 1. We first perform a fiducial set of simulations with fixed . We choose a series of values, ranging from down to , and study the critical value of below which MRI turbulence is no longer self-sustained (Section 3.1.1). In the next, we vary the net vertical flux by setting and and run a number of simulations around to study how the critical value of is affected by the vertical flux (Section 3.1.2). Moreover, in Section 3.1.3 we briefly investigate the effect of by varying from the fiducial value to (run Z5a) and (run Z5b) (see equation (4)). Finally, we discuss the properties of the MRI turbulence in the presence of AD (Section 3.1.4).

Our choices of the net vertical flux derive from the linear dispersion relation of the MRI as well as physical considerations. In the case of ideal MHD, the wavelength for the fastest growing linear MRI mode is given by (Hawley et al., 1995). For and , our vertical box size of fits and most unstable wavelengths respectively in ideal MHD. The ideal MHD dispersion relation is considerably modified when . Unstable modes exist for wavelength longer than the critical wavelength (Wardle, 1999)

(20) |

The wavelength for the most unstable mode is about twice larger. An approximate fitting formula that is accurate within for all values of is

(21) |

where . Note that for pure vertical magnetic field and vertical wavenumber, the linear dispersion relation for Ohmic resistivity is exactly the same as that for AD (Wardle, 1999), with replaced by the Elsasser number . For , the most unstable wavelengths at and are , and respectively. Clearly, the most unstable mode does not fit into our simulation box when , and no unstable modes fit into the box for . Since , these modes do fit into our simulation box as we increase to and respectively. In the mean time, since AD tends to be important in the more strongly magnetized upper layers of the protoplanetary disks (Wardle, 2007; Bai, 2011), it is also interesting to study whether the MRI turbulence can be sustained when is relatively small, even if the most (or all) unstable modes do not fit into our simulation box. We have not run simulations with a taller box since we do not include vertical stratification.

Run | Orbits | Restart | Turbulence | |||
---|---|---|---|---|---|---|

Z1 | 1000 | 400 | 0.5 | 48 | No | Yes |

Z2 | 100 | 400 | 0.5 | 24 | No | Yes |

Z3 | 10 | 400 | 0.5 | 24 | No | Yes |

Z4 | 3.33 | 400 | 0.5 | 24 | No | Yes |

Z5 | 1 | 400 | 0.5 | 43 | No | Yes |

Z6 | 0.33 | 400 | 0.5 | 24 | Z5 | No |

Z7 | 0.1 | 400 | 0.5 | 24 | Z5 | No |

Z3s | 10 | 100 | 0.5 | 24 | No | Yes |

Z5s | 1 | 100 | 0.5 | 24 | No | Yes |

Z6s | 0.33 | 100 | 0.5 | 24 | Z5s | No |

Z7s | 0.1 | 100 | 0.5 | 24 | Z5s | No |

Z3w | 10 | 1600 | 0.5 | 24 | No | Yes |

Z5w | 1 | 1600 | 0.5 | 24 | No | Yes |

Z6w | 0.33 | 1600 | 0.5 | 48 | No | Yes |

Z7w | 0.1 | 1600 | 0.5 | 24 | Z5w | No |

Z6e | 0.33 | 0.5 | 48 | No | Yes | |

Z7e | 0.1 | 0.5 | 48 | Z6e | Yes | |

Z5a | 1 | 400 | 0.0 | 24 | No | Yes |

Z5b | 1 | 400 | 1.0 | 24 | No | Yes |

Box size is fixed at , grid resolution is cells per .

whether simulation is initiated by restarting from a turbulent run.

whether turbulence can be self-sustained.

#### 3.1.1 A Fiducial Set of Runs

As the fiducial set of runs, we fix , and run simulations with different values (see Table 1), labeled from with , which essentially corresponds to the ideal MHD case, to with , where the evolution of magnetic field is dominated by AD. Our scan of is more narrowly sampled near where the transition is expected to occur. In Figure 3, we show the time evolution of the Maxwell stress from the fiducial set of runs. We find that for , the growth of the MRI from linear perturbations leads to vigorous MRI turbulence, while for the two runs with , MRI either does not grow from the initial vertical field (Z7), or grows too slowly (Z6), since the most unstable modes do not fit into our simulation box. Therefore, for these two models, we start the simulations from the end of run Z5 (), which is turbulent, and reset to be and respectively. Nevertheless, turbulence continues to decay throughout the span of our simulation in run Z7. Run Z6 is a marginal case where turbulence is neither fully sustained nor decayed continuously (see discussions below).

We first look at run Z1 to Z5 with . The initial growth of the MRI is due to the axisymmetric channel mode (Goodman & Xu, 1994; Pessah & Chan, 2008). The mode becomes non-linear (producing an overshoot in the Maxwell stress up to in the ideal MHD case) until broken down by secondary parasitic modes to produce turbulence. In the turbulent state, it is evident that the Maxwell stress monotonically decreases as decreases, analogous to the Ohmic case (Sano et al., 1998; Fleming et al., 2000). In Table 2 we list the general properties of the turbulence from all our vertical net flux simulations. The quantities are averaged over space and time after saturation (after ). The total stress in run Z1, which agrees with the ideal MHD case (Hawley et al., 1995). It drops slowly with decreasing when , but very rapidly when is around 1. Moreover, as decreases, the ratio of kinetic to the fluctuating part (with background field subtracted) of the magnetic energy increases (see also Figure 5). Similarly, the ratio of Reynolds stress to Maxwell stress increases.

As we have discussed before, the most unstable mode does not fit into our simulation box for runs Z6 and Z7, and our simulations initiated from a turbulent state also show no sign of sustained MRI turbulence. This is not surprising for run Z7, where no unstable MRI mode even exists in the simulation box. Our run Z6 is a marginal case, where a wavelength of is only slightly larger than the critical wavelength for instability () but far from the most unstable wavelength (). This explains the long-term variations in the Figure 3 since the growth rate is only slightly larger than zero. Our analysis in Section 3.2.4 indicates that although some non-zero stress close to is maintained in the simulation, it is unlikely to be due to the MRI turbulence. In real disks, one may expect sustained turbulence to be supported at if it were at the disk midplane, where the density variation over one above and below the midplane is not significant. In the upper layers, may fall off substantially over one , which strongly stabilizes the flow. In sum, we see from our fiducial set of net vertical flux simulations that presence of turbulence mainly depends on whether the most unstable mode of the MRI fits into the simulation box. This aspect will be further explored in the next subsection 3.1.2.

The results reported above are qualitatively different from those observed in HS using two-fluid simulations. One may compare our results with Table 3 of HS, where the value for their four runs Z24, Z17, Z25, Z28 are , , and respectively. The total stress in these simulations does not scale monotonically with , and in particular is on the order of when is as small as , while in our simulations MRI is suppressed. This reflects the difference in the physical assumptions about the two approaches. In the two-fluid limit, ions are coupled to the neutrals only via collisions. When the ion-neutral collisions are infrequent (), the ions and neutrals behave as independent fluid: vigorous MRI is generated in the ion fluid while the neutrals remain quiescent, and the overall is proportional to the ion fraction . In the strong-coupling limit, the ions and the neutrals are coupled not only by collisions, but also via the ionization-recombination reactions. Since the strong-coupling limit requires the recombination time to be much smaller than the orbital time, this means that ions are continuously created and destroyed on a time scale that is much shorter than the time scale for MRI to grow. It is this additional chemical coupling that quenches the MRI in the ion fluid, and suppresses angular momentum transport for .

#### 3.1.2 The Effect of Vertical Field Strength

We select a number of models from the fiducial series and rerun the simulations with three additional initial values: and (e.g., with magnetic field strength two times and half of that in the fiducial models, as well as one case with an very weak field). These simulations are labeled with an additional letter “s” (for strong), “w” (for weak) and “e” (for extremely weak) in Table 1. For strong field simulations with , the wavelength for the fastest growing mode exceeds the vertical box size when . When , there are essentially no unstable mode in the simulation box. Runs Z3s and Z5s are initiated from seed perturbations, while runs Z6s and Z7s are initiated from the turbulent state at the end of run Z5s. For weak field runs with , on the other hand, the most unstable mode can be fitted into the simulation box for all runs with . No unstable mode is fitted into the simulation box when , and as before, Run Z7w is initiated from turbulent state from Z5w to test whether turbulence can be sustained. Our simulations allows the most unstable wavelength to be fitted in the simulation box at small . We conduct two runs in this case with (Z6e) and (Z7e). Run Z7e is initialized from the turbulent state in Z6e to avoid the extremely long time in the linear growth stage. Time averaging in runs Z6w, Z6e and Z7e are taken since (time averaging in other runs are taken since by default).

We find from our simulations that sustained MRI turbulence is present in all models except Z6s, Z7s and Z7w. In particular, the MRI turbulence can be self-sustained even the value is as small as , provided that the net vertical field is sufficiently weak. These results confirm our speculation in the fiducial set of simulations that the MRI turbulence is self-sustained as long as unstable MRI modes fit into the simulation box.

The diagnostic quantities from time and volume averaged quantities in the turbulent state from the weak and strong field series of runs are also listed in Table 2. We see that for , the averaged kinetic energy, Reynolds and Maxwell stress monotonically decreases with increasing . Although not all our simulations are not run long enough for these quantities to be measured accurately, the trend is significant enough and indicates that the MRI saturate at a higher level with higher net vertical flux (small ), in agreement with the ideal MHD case (Hawley et al., 1995). For , the monotonicity trend is still present by comparing our fiducial run Z5 and the weak field run Z5w. The saturation level of the MRI turbulence in the strong field run Z5s is weaker than that for run Z5. This is most likely because the most unstable mode does not fit into our simulation box (but some less unstable modes fit) in run Z5s. The monotonicity trend further preserves at , where the kinetic energy density and total stress from run Z6w is larger than those from run Z6e by about a factor of 2.

We summarize the main results from the net vertical flux simulations in Figure 4. Shown are the total stress from all the simulations that the most unstable mode is properly resolved so that the MRI turbulence is self-sustained and a reliable value of can be obtained. As discussed before, at a fixed , there exists a critical value of below which the most unstable wavelength would exceed and the mode tend to be suppressed due to the vertical stratification. This is effect is illustrated by the colored arrows in the Figure. Equivalently, for turbulence to be sustained at small , must be sufficiently large , as can be obtained from equation (21). At a given , since the stress monotonically increases with the net vertical flux, there exist a maximum stress, corresponding to the largest allowed net flux (smallest allowed value of ). This maximum value of as a function of is illustrated in the dashed line by connecting results from runs Z7e, Z6w, Z5 and Z3s. We see that the maximum drops by a factor of about from the ideal MHD case to , and another factor of about as decreases to . By extrapolating this trend, we expect the MRI turbulence can be self-sustained for arbitrarily small value of , as long as the background magnetic field is sufficiently weak. Nevertheless, the turbulence would seem to be too weak () to produce significant amount of angular momentum transport as required by most astrophysical disks.

#### 3.1.3 The Effect of

The parameter reflects the sensitivity of how the AD coefficient depends on gas density (see Equation (4)). Most of our simulations are run with fixed value of , while can in principal span a range from to . The significance about the effect of largely depends on the level of density fluctuation in the MRI turbulence. In Table 2, we list the rms density fluctuation relative to the background gas density from all vertical net flux runs (see column ). The rms density fluctuation in the ideal MHD case (run Z1) is relatively large, up to , and the largest and smallest densities reach about and times the background density. Since AD reduces the saturation level of the MRI turbulence, the density fluctuations become smaller as decreases. This fact undermines the importance of : when the effect of may be important (large density fluctuations), AD only plays an insignificant role in the MRI turbulence (large ); when AD strongly affect the MRI turbulence (small ), the density fluctuation becomes much smaller and is much less likely to be important. This above implies that variations in the value of should not have a major impact, and in particular, the critical value of below which MRI is suppressed is unlikely to be altered by different choices of .

To confirm our expectations, we perform two additional runs with the same initial conditions as run Z5 (, ), but set to be and respectively. These two runs are named Z5a and Z5b. We see from Table 2 that the turbulence properties from these two runs are essentially identical to those in run Z5. Even though our time averages are taken over relatively short periods, the deviations are generally within . This is understandable since the density fluctuations in these runs are as small as . It appears certain that the value of only plays a very minor role in the MRI turbulence in the strong coupling limit.

Run | |||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|

Z1 | 0.22 | ||||||||||

Z2 |
0.11 | ||||||||||

Z3 |
|||||||||||

Z4 |
|||||||||||

Z5 |
|||||||||||

Z3s |
|||||||||||

Z5s |
|||||||||||

Z3w |
|||||||||||

Z5w |
|||||||||||

Z6w |
|||||||||||

Z6e |
|||||||||||

Z7e |
|||||||||||

Z5a |
|||||||||||

Z5b |

#### 3.1.4 Properties of the MRI turbulence with AD

Besides the general properties of the MRI turbulence listed in Table 2, we study two other aspects of the MRI turbulence with AD.

First, we study the power spectrum density (PSD) of magnetic and kinetic energies by Fourier analysis. The Fourier analysis in the shearing periodic system is performed by the remapping technique before and after Fourier transformation, as described in Section 2.4 of Hawley et al. (1995). Although the PSD is anisotropic in space, it would be beneficial to plot the PSD in one dimensional form by some averaging procedure. Following Davis et al. (2010), we compute shell-integrated power spectrum of the magnetic field , where denotes the average of over shells of constant , and is the Fourier transform of . Here is the volume of the simulation box. The Fourier transformation is of course discrete, but for notational convenience we write the formulas in continuous form. According to Parseval’s theorem, we have

(22) |

Dividing by a factor of we obtain the PSD for the magnetic energy density . Similarly, one can obtain the PSD for the kinetic energy .

In Figure 5, we show the PSDs computed from our runs Z1 and Z5. These two simulations are representative for the MRI turbulence in the ideal MHD and AD dominated regimes respectively, and are run for two times longer than many other simulations (thus giving better statistics). We see that the shape of the PSD obtained from our simulations are very similar. The PSD roughly follows a power-law form at small , with the power law index approximately equals to , which is the index for incompressible Kolmogorov turbulence spectrum. There appears to be a spectral break at , corresponding to a wavelength of about , and the PSD falls off rapidly toward smaller scales. The turbulent power in the case is about times smaller than that in the ideal MHD case. Magnetic energy fluctuations dominate kinetic energy fluctuations in the ideal MHD case, while in the AD dominated regime, more turbulence power resides in the kinetic energy. Moreover, we have also checked the contour plot of vertically integrated PSD (not shown) and found that the turbulence becomes more anisotropic in the AD dominated regime: the turbulent power is more elongated in than in .

Second, we study the effect of AD on the distribution of current in the MRI turbulence. It has been shown that in one and two dimensions, sharp current structure can be developed around magnetic nulls in the presence of AD (Brandenburg & Zweibel, 1994). To examine whether the same effect is present in the MRI turbulence, we show in Figure 6 the cumulative probability distribution of the current density in our simulation runs Z1, Z3 and Z5. If sharp current structure were to form, one would expect to see extended tails in the probability distribution. However, we see that as decreases, the probability distribution shifts leftward since turbulence becomes weaker, but its shape remains largely unchanged. We also note that the current sharpening phenomenon is not observed in simulations by (Brandenburg et al., 1995) either. It is likely the sharpening of current by AD is overwhelmed in 3D MHD turbulence.

AD has also been shown to tend to reduce the current component that is
perpendicular to the magnetic field (Brandenburg et al., 1995) and make the magnetic
configuration more force-free. To examine this effect, we show the cumulative
probability distribution of from our simulation runs Z1, Z2 and Z3 in
Figure 7, where is the angle between the current and the magnetic
field. The cumulative distribution functions from our runs Z4 and Z5 are almost identical
to that from run Z3, where . We confirm that AD makes the distribution more
concentrated toward (i.e., ).^{5}^{5}5We also
report a difference in our results from Brandenburg et al. (1995). In their simulations,
the distribution function peaks at in ideal the MHD case, and AD
concentrate the current toward . In our simulations, we find that the
distribution function is already concentrated at under ideal MHD, and
AD simply makes it more concentrated. Nevertheless, since the distribution of
is already peaked at 1 in the essentially ideal MHD run (Z1), the effect
of AD does not modify the distribution of current orientation substantially.

### 3.2. Net toroidal flux simulations

In the second group of simulations, the initial field configuration is chosen to be
uniform along the azimuthal direction , with strength characterized by
, where is the initial field strength. Following
Simon & Hawley (2009), we fix the box size to be in the radial,
azimuthal and vertical dimensions for all simulation runs in this group^{6}^{6}6We have
also performed our run series Y1 to Y6 using a larger box and
found that the turbulence properties are very similar.. We choose the fiducial resolution
to be cells per in the radial and vertical direction, and cells per in
the azimuthal direction. All our net toroidal flux simulations are listed in Table
3, including one set of fiducial simulations with , one set of
higher resolution simulations, and one set of weak field simulations with .
Unlike net vertical flux, the net toroidal flux is not precisely conserved in our shearing
box simulations. As discussed in (Simon & Hawley, 2009), ensuring strict conservation of
toroidal flux numerically is more complex, and is also less important than conserving net
vertical flux because the saturation level of the MRI turbulence is not very sensitive to the
toroidal flux. Throughout all our simulations in this group, we find that the deviation of net
toroidal flux from the initial value is generally less than .

Run | Resolution | Restart | Turbulence | ||
---|---|---|---|---|---|

Y1 | 1000 | 100 | No | Yes | |

Y2 | 100 | 100 | No | Yes | |

Y3 | 10 | 100 | Y1 | Yes | |

Y4 | 3.33 | 100 | Y1 | Yes | |

Y5 | 1 | 100 | Y1 | No | |

Y6 | 0.33 | 100 | Y1 | No | |

Y7 | 0.1 | 400 | Y1 | No | |

Y1w | 1000 | 100 | No | Yes | |

Y2w | 100 | 100 | No | Yes | |

Y3w | 10 | 400 | Y1w | Yes | |

Y4w | 3.33 | 400 | Y1w | Yes | |

Y5w | 1 | 400 | Y1w | No | |

Y6w | 0.33 | 400 | Y1w | No | |

Y7w | 0.1 | 400 | Y1w | No | |

Y1h | 1000 | 100 | No | Yes | |

Y3h | 10 | 100 | Y1h | Yes | |

Y5h | 1 | 100 | Y1h | No | |

Y6h | 0.33 | 100 | Y1h | No |

Box size is fixed at , is fixed at . All simulations are run for orbits (150).

The linear stability of Keplerian disks in the presence of pure toroidal field is more complex than that for the vertical field case. It requires consideration of non-axisymmetric perturbations (Balbus & Hawley, 1992), and involves the time-dependent amplification of wave modes as the radial wave number swings from leading to trailing. In ideal MHD, pure toroidal MRI favors high wave numbers, and requires relatively large numerical resolution. In the case of Ohmic resistivity, swing amplification of modes is suppressed when the diffusion time of the mode is comparable to the orbital frequency (Papaloizou & Terquem, 1997). A linear stability with non-axisymmetric perturbations in AD dominated regime has yet to be performed. Nevertheless, one might expect that a similar argument holds for AD, with as the boundary for stability.

#### 3.2.1 A Fiducial Set of Runs

We fix and run 7 fiducial simulations with different values, named from Y1 with to Y7 with (see Table 3) similar to the case of net vertical flux runs. Initial growth of the MRI from pure toroidal field is more difficult and is only achieved when is greater than . We initialize the rest of the simulations from the turbulent state at the end of run Y1.

Figure 8 illustrates the time evolution of the Maxwell stress from this fiducial set of runs. The time and volume averaged quantities from these runs are listed in Table 4. We find that at a given value of , the saturation level of the MRI with net toroidal flux is much lower than the net vertical flux case. The turbulent energy density and total stress in run Y1 (essencially ideal MHD) is a few times , about an order of magnitude less than run Z1. As drops below , the saturation level of the MRI turbulence falls off very rapidly. At (run Y4), the total stress falls below . At , although a total stress is maintained at a level of , we do not observe any signature of the MRI turbulence by examining the structure of the velocity field, which is essentially laminar (see further discussion in Section 3.2.4). Since the simulation is initialized from a turbulent state, the low level of stress and kinetic energy are mostly due to the eigen-modes of the shearing box excited from the initial turbulence, and that are not damped due to the low dissipation in the Athena code. Unlike the case for net vertical flux, there appears to exist a critical value of below which MRI turbulence with net toroidal flux is not self-sustained. This critical value of is about . This fact will be further discussed shortly in Section 3.2.3.

HS also performed a number of net toroidal flux two-fluid simulations with and , and in both cases turbulence is self-sustained with total stress on the order of to . Again, these results are no longer valid in the strong-coupling regime and are not directly comparable to our results (see discussion in Section 3.1.1).

We see from Table 4 that the density fluctuations in the net toroidal flux simulations are generally smaller than those in the net vertical flux case. Therefore, following the discussion in Section 3.1.3, we expect the effect of has essentially no impact on the conclusions we have drawn above.

Run | |||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|

Y1 | |||||||||||

Y2 |
|||||||||||

Y3 |
|||||||||||

Y4 |
|||||||||||

Y5 |
|||||||||||

Y1w |
|||||||||||

Y2w |
|||||||||||

Y3w |
|||||||||||

Y4w |
|||||||||||

Y5w |
|||||||||||

Y1h |
|||||||||||

Y3h |
|||||||||||

Y5h |

: These runs are not turbulent.

#### 3.2.2 The Effect of Resolution

Relatively high resolution is needed for net toroidal flux MRI simulations in order to properly capture the amplification of wave modes as they swing from leading to trailing (Simon & Hawley, 2009). In order to justify our results in the previous subsection, we perform a few of the simulations with doubled resolution. These runs are labeled with an additional letter “h” (i.e., high resolution) in Table 3.

The time and volume averaged quantities from these high resolution runs are shown in Table 4. We see that the kinetic and magnetic energy densities from in the high resolution simulations are generally larger than those in the low resolution runs, but only by a small factor. In particular, the difference between low and high resolutions with relatively large AD coefficient is only about (e.g., comparing runs Y3 and Y3h), which strongly indicates numerical convergence. This is not surprising since small scale structures can be largely damped by AD thus higher resolution becomes unnecessary. Run Y5h is also very similar to run Y5, where the initial turbulence is damped with remnant small velocity and magnetic fluctuations unlikely to be associated with the MRI turbulence.

The inferences above are further justified by looking at the power spectrum of magnetic and kinetic energies. Following the same procedure described in Section 3.1.4, we show in Figure 9 the shell integrated PSDs for runs Y1, Y1h and Y3, Y3h. The spectral shapes from the low resolution simulations appear to have a spectral peak at relatively large scales of about , while at higher resolution, a power law spectrum at intermediate scales from approximately to analogous to an inertia range appears to have developed. This observation may suggest high numerical resolution of at least cells per is needed for the toroidal field MRI simulations in order to resolve the inertia range in the turbulent spectrum, although smaller resolution of cells per appears to be sufficient for turbulence properties to converge. The shape of the PSDs at small look very different from the PSDs in the net vertical flux simulations, indicating different energy injection mechanism, while at large , the PSDs from both cases fall off in a similar manner, indicating similar dissipation mechanism.

#### 3.2.3 The Effect of toroidal Field Strength

We have also performed the same set of net toroidal flux simulations with the toroidal magnetic field strength lowered by one half, , labeled with an additional letter “w” (i.e., weak field) in Table 3. General properties from the saturated state of these runs are also shown in Table 4. We see that at the same value of , the kinetic and magnetic energy density from the weak field simulations are smaller than those in our fiducial runs by a factor of -. By inspection of the velocity field as well as the distribution of current density, we find that sustained turbulence is supported in run Y4w but not in run Y5w. Note that the time and volume averaged kinetic energy density from run Y4w is , which is sli