EFFECT OF AMBIPOLAR DIFFUSION ON THE NONLINEAR EVOLUTION OF MAGNETOROTATIONAL INSTABILITY IN WEAKLY IONIZED DISKS

The effect of AD 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 (∂/∂t = 0) shock and work in the shock frame, with upstream gas density ρ₀ moving at velocity v_s. The upstream gas is threaded by a uniform magnetic field ${\boldsymbol{B}}_0$ that lies at an angle θ to the velocity. Let ${\boldsymbol{v}}_s$ be in the $\hat{x}$-direction, and ${\boldsymbol{B}}_0$ be in the $\hat{x}{-}\hat{y}$-plane. For a continuous shock, the jump conditions reduce to ∂/∂x = 0. Assuming the ion density ρ_i is constant (ν = 0), the equations that describe the C-type shock read

$$\begin{equation} \rho v_x=\rho _0v_s, \end{equation} \tag{ 9 }$$

$$\begin{equation} \rho \big(v_x^2+c_s^2\big)+\frac{B_y^2}{8\pi }=\rho _0\left(v_s^2+c_s^2\right) +\frac{B_{0y}^2}{8\pi }, \end{equation} \tag{ 10 }$$

$$\begin{equation} \rho v_xv_y-\frac{B_xB_y}{4\pi }=-\frac{B_xB_{0y}}{4\pi }, \end{equation} \tag{ 11 }$$

$$\begin{equation} v_xB_y-v_yB_x-\frac{B^2}{4\pi \rho \gamma \rho _i} \frac{dB_y}{dx}=v_sB_{0y}. \end{equation} \tag{ 12 }$$

Note that B_x = B_0x is constant.

The shock is characterized by three dimensionless parameters: the sonic Mach number M = v_s/c_s, the Alfvén Mach number A = v_s/v_A (where v²_A = B₀²/4πρ₀), and the angle θ of the magnetic field with the upstream flow. The characteristic length scale of the problem is given by L = v_A/γρ_i. We further define D ≡ ρ/ρ₀ and b ≡ B_y/B₀. After some algebra, we arrive at a dimensionless first-order differential equation for D (Mac Low et al. 1995)

$$\begin{eqnarray} && \bigg (\frac{1}{D^2}-\frac{1}{M^2}\bigg)\frac{dD}{d(x/L)} =(b^2+\cos ^2\theta)^{-1}\nonumber\\ &&\quad\quad\times \frac{b}{A}\bigg [b-D\bigg (\frac{b-\sin \theta }{A^2}\cos ^2\theta +\sin \theta \bigg)\bigg ]. \end{eqnarray} \tag{ 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 M = 50, A = 10, and θ = π/4 obtained by using a fourth-order Runge–Kutta method.

Zoom In Zoom Out Reset image size

To use this solution as a code test, the shock is set to be aligned with the grid in the $\hat{\boldsymbol{x}}$-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 (∼5L/c_s) 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, one-dimensional, two-dimensional, and three-dimensional (1D, 2D, and 3D, respectively) tests essentially produce the same result. In our tests, the grid resolution is chosen to be two and four cells per L.³ We see that our code very accurately resolves the structure of the C-type shock using only a few cells per L. The main source of the error lies in the region where the density and velocity profiles 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.

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 that 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

$$\begin{equation} \omega ^2=k^2v_A^2\cos ^2\theta \bigg (1-{\rm i}\frac{\,\Omega }{\,\Omega _a}\bigg). \end{equation} \tag{ 14 }$$

where  Ω_a ≡ γρ_i, v_A is the Alfvén velocity, and θ 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

$$\begin{equation} \big(\omega ^2-v_f^2\big)\big(\omega ^2-v_s^2\big)+{\rm i} \big(\omega ^2-k^2c_s^2\big)k^2v_A^2\frac{\omega }{\omega _a}=0, \end{equation} \tag{ 15 }$$

where v_f and v_s are the fast and slow magnetosonic speeds

$$\begin{equation} v_{f,s}^{2}=\frac{1}{2}\big(c_s^2+v_A^2\big)\pm \frac{1}{2}\sqrt{\big(c_s^2+v_A^2\big)^2-4c_s^2v_A^2\cos \theta }, \end{equation} \tag{ 16 }$$

with plus (minus) sign corresponding to v_f (v_s).

When ω ≪ Ω_a, the damping rate is small and can be found by expanding the ideal MHD dispersion relation to powers of ω/ω_a, and to the first order, we find for the damping rate of the Alfvén wave

$$\begin{equation} \Gamma _A=\frac{1}{2}\frac{k^2v_A^2\cos ^2\theta }{\omega _a}. \end{equation} \tag{ 17 }$$

The damping rates for fast and slow magnetosonic waves are

$$\begin{equation} \Gamma _f=\frac{1}{2}\bigg (\frac{v_f^2-c_s^2}{v_f^2-v_s^2}\bigg) \frac{k^2v_A^2}{\omega _a}\;,\qquad \Gamma _s=\frac{1}{2}\bigg (\frac{c_s^2-v_s^2}{v_f^2-v_s^2}\bigg) \frac{k^2v_A^2}{\omega _a}. \end{equation} \tag{ 18 }$$

We perform the linear wave damping test in 1D, 2D, and 3D. In 1D, the wave is grid-aligned and its wavelength equals 1 in code unit. In 2D and 3D test problems, the wave vectors are not grid-aligned, and box sizes chosen such that the wave length is also 1 (in 2D, the box size is ($\sqrt{5}, \sqrt{5}/2$) and in 3D, the box size is (3, 1.5, 1.5)). We use isothermal equation of state with c_s = 1. The background gas density is ρ₀ = 1. As before, we choose ν = 0. In 1D, the wave vector is along the x-direction, and the adopted magnetic field is B_0x = 1.0, $B_{0y}=\sqrt{2}$ and B_0z = 0.5. 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 $v_A=\sqrt{13/4}$, v_f = 2, and v_s = 1/2, therefore we find Γ_A = 2π²Q_A, Γ_f = 5.2π²Q_A, and Γ_s = 1.3π²Q_A.

In practice, we adopt ω_a = 100 and run our simulation to t = 5. By default, the grid size is 32 in 1D, 64 × 32 in 2D, and 64 × 32 × 32 in 3D. Accounting for the box size in each dimension, the effective resolution, characterized by number of cells per wavelength, is 32, 28.6, and 21.3 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, respectively. 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.

Zoom In Zoom Out Reset image size

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 15%–30% at t = 5). These results indicate that at least 20 cells per wavelength is needed to accurately capture the effect of AD.

As the fiducial set of runs, we fix β₀ = 400, and run 7 simulations with different Am values (see Table 1), labeled from Z1 with Am = 1000, which essentially corresponds to the ideal MHD case, to Z7 with Am = 0.1, where the evolution of magnetic field is dominated by AD. Our scan of Am is more narrowly sampled near Am = 1 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 Am ⩾ 1, the growth of the MRI from linear perturbations leads to vigorous MRI turbulence, while for the two runs with Am < 1, 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 (Am = 1), which is turbulent, and reset Am to be 0.33 and 0.1, 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 runs Z1 to Z5 with Am ⩾ 1. The initial growth of the MRI is due to the axisymmetric channel mode (Goodman & Xu 1994; Pessah & Chan 2008). The mode becomes nonlinear (producing an overshoot in the Maxwell stress up to 1 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 Am 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 100 Ω⁻¹). The total stress α ≈ 0.2 in run Z1, which agrees with the ideal MHD case (Hawley et al. 1995). It drops slowly with decreasing Am when Am ≫ 1, but very rapidly when Am is around 1. Moreover, as Am decreases, the ratio of kinetic to the fluctuating part (with background field B_z0 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 H is only slightly larger than the critical wavelength for instability (λ_c = 0.89H) but far from the most unstable wavelength (λ_m = 1.72H). This explains the long-term variations in 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 α = 10⁻⁴ 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 Am ≈ 0.3 if it were at the disk midplane, where the density variation over one H above and below the midplane is not significant. In the upper layers, β₀ may fall off substantially over one H, which strongly stabilizes the flow. In sum, we see from our fiducial set of net vertical flux simulations that the 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 Section 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 Am value for their four runs Z24, Z17, Z25, and Z28 are 0.11, 1.1, 11.1, and 111, respectively. The total stress α in these simulations does not scale monotonically with Am, and in particular α is on the order of 10⁻² when Am is as small as 0.11, while in our simulations MRI is suppressed. This reflects the difference in the physical assumptions about the two approaches. In these two-fluid simulations, the ion to neutral mass ratios (f = ρ_i/ρ) were kept to be relatively large (0.001–0.1) to avoid numerical stiffness, thus the ion inertia always plays a role and the ion drift velocity deviates from that in the strong-coupling limit as in Equation (1). The departure is more and more significant as Am decreases, and in particular, for Am ≲ 0.1, ions and neutrals behave as independent fluids in HS: vigorous MRI is generated in the ion fluid while the neutrals remain quiescent, and the overall α is simply proportional to f. In addition, the strong-coupling limit requires the recombination time to be smaller than the orbital time, which means that ions are continuously created and destroyed on a time scale that is shorter than the time scale for MRI to grow. This additional chemical coupling introduces ionization modes (Blaes & Balbus 1994) which were not captured in HS.

We select a number of models from the fiducial series and rerun the simulations with three additional initial β₀ values: β₀ = 100, 1600, and 10⁴ (e.g., with magnetic field strength two times and half of that in the fiducial models, as well as one case with a 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 β₀ = 100, the wavelength for the fastest growing mode exceeds the vertical box size when Am < 10. When Am < 1, 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 β₀ = 1600, on the other hand, the most unstable mode can be fitted into the simulation box for all runs with Am ≳ 0.3. No unstable mode is fitted into the simulation box when Am = 0.1, and as before, Run Z7w is initiated from turbulent state from Z5w to test whether turbulence can be sustained. Our β₀ = 10⁴ simulations allows the most unstable wavelength to be fitted in the simulation box at small Am. We conduct two runs in this case with Am = 0.33 (Z6e) and Am = 0.1 (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 starting at t = 200 Ω⁻¹ (time averaging in other runs are taken starting at t = 100 Ω⁻¹ 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 Am value is as small as 0.1, 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 Am = 10, the averaged kinetic energy, Reynolds and Maxwell stress monotonically decreases with increasing β₀. Although not all of our simulations are run long enough for these quantities to be measured accurately, the trend is significant enough and indicates that the MRI saturates at a higher level with higher net vertical flux (small β₀), in agreement with the ideal MHD case (Hawley et al. 1995). For Am = 1, 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 Am = 0.33, where the kinetic energy density and total stress from run Z6w are larger than those from run Z6e by about a factor of two.

We summarize the main results from the net vertical flux simulations in Figure 4. Shown are the total stress α values from all the simulations for which the most unstable mode is properly resolved so that the MRI turbulence is self-sustained and a reliable value of α can be obtained. As previously discussed, at a fixed β₀, there exists a critical value of Am below which the most unstable wavelength would exceed H and the mode tends to be suppressed due to the vertical stratification. This effect is illustrated by the colored arrows in the figure. Equivalently, for turbulence to be sustained at small Am, β₀ must be sufficiently large β₀ ≳ 100/Am², as can be obtained from Equation (22). At a given Am, since the stress α monotonically increases with the net vertical flux, there exists a maximum stress, corresponding to the largest allowed net flux (smallest allowed value of β₀). This maximum value of α as a function of Am 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 40 from the ideal MHD case to Am = 1, and another factor of about 60 as Am decreases to 0.1. By extrapolating this trend, we expect that the MRI turbulence can be self-sustained for an arbitrarily small value of Am as long as the background magnetic field is sufficiently weak. Nevertheless, the turbulence would seem to be too weak (α < 10⁻⁴) to produce a significant amount of angular momentum transport as required by most astrophysical disks.

Zoom In Zoom Out Reset image size

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 ν = 0.5, while ν can in principal span a range from 0 to 1. 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 0.3, and the largest and smallest densities reach about 4 and 0.2 times the background density, respectively. Since AD reduces the saturation level of the MRI turbulence, the density fluctuations become smaller as Am 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 Am); when AD strongly affects the MRI turbulence (small Am), 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 Am 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 (Am = 1, β₀ = 400), but set ν to be 0 and 1, 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 10%. This is understandable since the density fluctuations in these runs are as small as 0.07. It appears certain that the value of ν only plays a very minor role in the MRI turbulence in the strong coupling limit. This result further implies that the requirement of electron recombination time being shorter than the dynamical time is secondary and may be moderately relaxed.

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 k-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 $B_k^2\equiv 4\pi k^2|\widetilde{\boldsymbol{B}}(k)|^2$, where $|\widetilde{\boldsymbol{B}}(k)|^2$ denotes the average of $|\widetilde{\boldsymbol{B}}({\boldsymbol{k}})|^2$ over shells of constant $k=|{\boldsymbol{k}}|$, and $\widetilde{\boldsymbol{B}}({\boldsymbol{k}})=\int {\boldsymbol{B}}({\boldsymbol{x}})e^{-i{\boldsymbol{k}}\cdot {\boldsymbol{x}}}d^3{\boldsymbol{x}}/V$ is the Fourier transform of ${\boldsymbol{B}}({\boldsymbol{x}})$. Here V is the volume of the simulation box. The Fourier transformation is of course discrete, but for notational convenience we write the formulae in continuous form. According to Parseval's theorem, we have

$$\begin{eqnarray} \frac{1}{V}\int _V |{\boldsymbol{B}}({\boldsymbol{x}})|^2d^3{\boldsymbol{x}}&=&\int |\widetilde{\boldsymbol{B}}({\boldsymbol{k}})|^2 \frac{d^3{\boldsymbol{k}}}{(2\pi)^3}\nonumber\\ &=&\frac{1}{(2\pi)^3}\int 4\pi k^2|\widetilde{\boldsymbol{B}}(k)|^2dk. \end{eqnarray} \tag{ 23 }$$

Dividing by a factor of 8π we obtain the PSD for the magnetic energy density M_k = B²_k/8π. Similarly, one can obtain the PSD for the kinetic energy K_k.

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 is very similar. The PSD roughly follows a power-law form at small k, with the power-law index approximately equals −11/3, which is the index for incompressible Kolmogorov turbulence spectrum. There appears to be a spectral break at kH ≈ 70 in both simulations, corresponding to a wavelength of about 0.1H, and the PSD falls off rapidly toward smaller scales. The turbulent power in the Am = 1 case is about 20 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. We note that although AD provides explicit magnetic dissipation, the scale of the spectral break is similar to the ideal MHD case, which is mainly because no explicit viscosity is included and the viscous damping scale remains unchanged. 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 k_x than in k_y.

Zoom In Zoom Out Reset image size

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 $J=|{\boldsymbol{J}}|$ 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 Am decreases, the probability distribution shifts leftward since turbulence becomes weaker, but its shape remains largely unchanged. We note that the cutoff of current density in the ideal MHD case is likely to be due to numerical dissipation,⁵ while the cutoffs at Am = 1 and 10 are physical. 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.

Zoom In Zoom Out Reset image size

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 |cos α| 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 Am = 10. We confirm that AD makes the distribution more concentrated toward |cos α| = 1 (i.e., ${\boldsymbol{J}}\parallel {\boldsymbol{B}}$).⁶ Nevertheless, since the distribution of |cos α| is already peaked at 1 in the essentially ideal MHD run (Z1), the effect of AD does not substantially modify the distribution of current orientation.

Zoom In Zoom Out Reset image size

We fix β₀ = 100 and run 7 fiducial simulations with different Am values, named from Y1 with Am = 1000 to Y7 with Am = 0.1 (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 Am is greater than 10. 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 Am, 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 (essentially ideal MHD) are a few times 10⁻², about an order of magnitude less than run Z1. As Am drops below 10, the saturation level of the MRI turbulence falls off very rapidly. At Am = 3 (run Y4), the total stress falls below 10⁻³. At Am = 1, although a total stress is maintained at a level of 10⁻⁵, 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 Am below which MRI turbulence with net toroidal flux is not self-sustained. This critical value of Am is about 3. This fact will be further discussed shortly in Section 3.2.3.

Zoom In Zoom Out Reset image size

HS also performed a number of net toroidal flux two-fluid simulations with Am ≈ 1 and Am ≈ 100, and in both cases turbulence is self-sustained with total stress α on the order of 10⁻⁴ to 10⁻³. 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.

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 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 10% (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 0.5H, while at higher resolution, a power-law spectrum at intermediate scales from approximately 0.1H to 0.5H analogous to an inertia range appears to have developed. This observation may suggest high numerical resolution of at least 128 cells per H is needed for the toroidal field MRI simulations in order to resolve the inertia range in the turbulent spectrum, although smaller resolution of 64 cells per H appears to be sufficient for turbulence properties to converge. The shape of the PSDs at small k looks very different from the PSDs in the net vertical flux simulations, indicating different energy injection mechanism, while at large k, the PSDs from both cases fall off in a similar manner, indicating similar dissipation mechanism.

Zoom In Zoom Out Reset image size

We have also performed the same set of net toroidal flux simulations with the toroidal magnetic field strength lowered by one half, β₀ = 400, 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 Am, the kinetic and magnetic energy density from the weak field simulations are smaller than those in our fiducial runs by a factor of 2–3. 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 4.9 × 10⁻⁴, which is slightly below our limit of 10⁻³, but the total stress α = 3.2 × 10⁻⁴ is reasonably large and is unlikely to be caused by inertia waves in the simulation box. Further discussion will be given in Section 3.2.4.

Together with the results from Section 3.2.1, we summarize the results from net toroidal flux simulations in Figure 10. We see that for both values of β₀ (100 and 400), we do not observe any sustained MRI turbulence for Am < 3. Although we have not explored weaker net toroidal flux, based on the trend we see from the figure, even if MRI can be self-sustained at Am ≲ 1 with weaker net toroidal flux, the resulting total stress α is unlikely to be above 10⁻⁴. This is in stark contrast with the net vertical flux simulations and indicates that pure toroidal field geometry is more stable in the AD dominated regime.

Zoom In Zoom Out Reset image size

A key objective of this paper is to study when MRI turbulence can be self-sustained in the presence of AD in the strong coupling limit. For the net vertical flux simulations, the criterion is relatively clear: turbulence can be sustained as long as the most MRI unstable mode fits into the simulation box. The situation is less clear for net toroidal flux simulations, and it remains to be seen whether the latter can be characterized as self-sustained turbulence.

We take run Y5 as an illustrative example in this subsection, but the analysis also applies to other marginal runs including Z6, Y5w, and Y5h. In run Y5, after the initial turbulence has damped, the flow is largely laminar, except in a few very localized regions where some narrow but azimuthally elongated current structure is present and evolves very slowly with time. These features can be better demonstrated by computing the vertically integrated Fourier power spectrum of magnetic and kinetic energies. Following the same procedure as described in Section 3.1.4 but integrating the full 3D PSD over k_z, we show the contour plot of the kinetic energy PSD in the k_x–k_y-plane from our runs Y1, Y4, and Y5 in Figure 11. We see that in ideal MHD (Y1), the vertically integrated PSD has elliptic contours elongated and tilted toward the k_x-axis. The contours are distorted and more elongated when AD is added (Y4). However, in run Y5, we see that the overall shape of the PSD contours are extremely elongated in the k_x-direction. The original elliptic contours are almost destroyed, with irregular fragments distributed around the center. These irregular features in the vertically integrated PSD strongly indicate that the system is not in a turbulent state. We have also found that such irregular features are also present in runs Z6, Y5w, and Y5h, but not present in runs Z5, Y4w, and Y3h.

Zoom In Zoom Out Reset image size

In sum, we conclude that non-zero kinetic energy density and total stress do not prove the existence of the MRI turbulence in shearing box simulations. Transition from self-sustained turbulence to a non-turbulent state can be judged by looking at the vertically integrated PSD of kinetic and magnetic energies.