INTERDEPENDENCE OF ELECTRIC DISCHARGE AND MAGNETOROTATIONAL INSTABILITY IN PROTOPLANETARY DISKS

There are three diffusion terms in magnetohydrodynamic (MHD); they are ohmic diffusion, Hall diffusion, and ambipolar diffusion. In protoplanetary disks, any one of the three modes can be the dominant mode depending on dust and gas density (e.g., Wardle 2007), and interaction between the different modes may alter the MRI (Wardle & Salmeron 2012). In this paper, we only focus on the ohmic diffusion because it is the most studied one in the context of the MRI. We leave the treatment of other diffusion modes for future studies.

The electric discharge taken into account, we construct a simple model of the discharge as follows, in terms of an appropriate η₀ and J_crit:

$$\begin{eqnarray} {\mathbf {E^{\prime }}} &=& \frac{4\pi }{c^2} \eta \left(J\right){\mathbf {J}}, \end{eqnarray} \tag{ 1 }$$

and

$$\begin{eqnarray} \eta \left(J\right) &=& \eta _0\ \ \ \ \mathrm{if} \ J< J_{\mathrm{crit}}, \nonumber \\ &=& \frac{J_{\mathrm{crit}}}{J}\eta _0\ \ \ \ \mathrm{if} \ J> J_{\mathrm{crit}}. \end{eqnarray} \tag{ 2 }$$

This nonlinear diffusivity model states that the electric field on the fluid comoving frame never exceeds a critical value, E'_crit = 4πc⁻²η₀J_crit, and thus the magnetic diffusivity η varies depending on electric current J, and Ohm's law become nonlinear.

Note that the smallest space scale dealt with in this paper is on the order of 10⁻² AU. The actual scale of the discharge structures can be much smaller than this. The estimate for the lower limit of the size of such structures is their thickness, which is on the order of 5000 times the electron mean free path (Pilipp et al. 1992).

However, we can derive the macroscopic discharge model from the microscopic discharge relation |E'| < E'_crit that holds everywhere in the plasma. Since, e.g., the x-component of the discretized electric field 〈E'_x〉 is obtained from the line integral of the real field over the discretization length ΔV,

$$\begin{eqnarray} |\langle {E_x^{\prime }} \rangle | &=& \left| \frac{\int {E_x^{\prime }} \mathit {dV}}{\Delta V}\right| \nonumber \\ &\le & \frac{\int |{E_x^{\prime }}| \mathit {dV}}{\Delta V} \nonumber \\ &\le & \frac{\int {E_{\mathrm{crit}}^{\prime }}\mathit {dV}}{\Delta V} \nonumber \\ &=&{E_{\mathrm{crit}}^{\prime }}. \end{eqnarray} \tag{ 3 }$$

Therefore, we can use Equations (1) and (2) as a "coarse-grained model" where we can interpret E' and J as spatial averages. If electrical breakdown occurs in a scale smaller than the grid size, the spatially averaged electric field is smaller than E'_crit. Thus, in general, the electrical breakdown may occur even in the region where 〈E'_crit〉 is smaller than E'_crit. Therefore, if we adopt Equations (1) and (2), we may underestimate, but not overestimate, the occurrence of electric discharges.

We adopted a local, Kepler-rotation shearing box that has radial (x), azimuthal (y), and vertical (z) axes and solved the following resistive MHD equations numerically:

$$\begin{eqnarray} \frac{\partial \rho }{\partial t} + \nabla \cdot \left(\rho {\mathbf {v}}\right) &=& 0 , \end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray} \frac{\partial \mathbf {v}}{\partial t} + {\mathbf {v}} \cdot \nabla {\mathbf {v}} &=& - \frac{1}{\rho }\nabla \left(P + \frac{B^2}{8\pi }\right) + \frac{1}{4\pi \rho }\left({\mathbf {B}}\cdot \nabla {\mathbf {B}}\right)\nonumber \\ & & - 2{\mathbf {\Omega }}\times {\mathbf {v}}+3\Omega ^2x{\hat{\mathbf {x}}} , \end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray} \frac{\partial {\mathbf {B}}}{\partial t} &=& - c \nabla \times {\mathbf {E}} , \end{eqnarray} \tag{ 6 }$$

with isothermal equation of state (EOS)

$$\begin{eqnarray} &&P = c_s^2 \rho . \end{eqnarray} \tag{ 7 }$$

Had we used the adiabatic EOS, the internal energy would have kept growing as the linear function of time (HGB95). Therefore, we use the isothermal EOS to approximate the steady state attained by the cooling processes present in the protoplanetary disks.

The nonlinear Ohm's law reads:

$$\begin{eqnarray} {\mathbf {E}} &=& -\frac{1}{c}{\mathbf {v}}\times {\mathbf {B}}+\frac{4\pi }{c^2}\eta \left(J\right){\mathbf {J}}, \end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray} {\mathbf {J}} &=& \frac{c}{4\pi }\nabla \times {\mathbf {B}}. \end{eqnarray} \tag{ 9 }$$

Here, we have studied three different models for nonlinear diffusivity η(J). In addition to our fiducial model (fid), Equations (1) and (2), we have studied the following two models:

$$\begin{eqnarray} \left(\mathtt {p2}\right):\ \eta \left(J\right) &=& \left(1 + \left(\frac{J_{\mathrm{crit}}}{J}\right)^{-2}\right)^{-{1}/{2}} \end{eqnarray} \tag{ 10 }$$

and

$$\begin{eqnarray} \left(\mathtt {p4}\right):\ \eta \left(J\right) &=& \left(1 + \left(\frac{J_{\mathrm{crit}}}{J}\right)^{-4}\right)^{-{1}/{4}}. \end{eqnarray} \tag{ 11 }$$

The electric field as a function of current density in these three models is shown in Figure 1.

Zoom In Zoom Out Reset image size

Following HGB95, we set up our numerical initial conditions as follows. We use the disk scale height H as the unit length. The box size is (L_x, L_y, L_z) = (1, 2π, 1). First,we set the average values ρ₀ = 1 and P₀ = 10⁻⁶ to every mesh and let fluid velocity be at rest in the shearing box frame; (v_x, v_y, v_z = 0, −(3/2)Ωx, 0). Here, c_s = 10⁻³, and also Ω = 10⁻³.

Next, we introduce random perturbations in density, pressure, and velocity. The density and pressure perturbations are in proportion so that the isothermal EOS is met, and the amplitude is δρ/ρ₀ = δP/P₀ = 2.5 × 10⁻². We perturb the velocity componentwise, with the amplitude δv_i = 5 × 10⁻³c_s for each.

We use the following units for the magnetic field, electric current, and scale height, respectively:

$$\begin{eqnarray} B_{\mathrm{eqp}}&\equiv & \sqrt{8 \pi P_0} , \end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray} {J_{\mathrm{eqp}}}&\equiv & \frac{{c B_{\mathrm{eqp}}} }{ 4\pi H} , \end{eqnarray} \tag{ 13 }$$

and

$$\begin{eqnarray} H&=& \frac{c_s}{\Omega } . \end{eqnarray} \tag{ 14 }$$

We set the uniform magnetic field in the z-direction and express the initial field strength by the plasma beta,

$$\begin{eqnarray} &&\beta \equiv {B_{\mathrm{eqp}}}^2/{B_{z0}}^2 = 8\pi P_0/ {B_{z0}}^2 . \end{eqnarray} \tag{ 15 }$$

The plasma beta satisfy the following relation between the sound speed c_s and the Alfvén velocity along the magnetic field $v_{\mathrm{{}_{AZ}}}$:

$$\begin{eqnarray} &&\beta = \frac{2 c_s^2}{{v_{\mathrm{{}_{AZ}}}}^2} . \end{eqnarray} \tag{ 16 }$$

We define the magnetic Reynolds number as

$$\begin{eqnarray} &&R_{\mathrm{M}}\equiv {v_{\mathrm{{}_{AZ}}}}^2 / \eta _0 \Omega , \end{eqnarray} \tag{ 17 }$$

using the Alfvén velocity $v_{\mathrm{{}_{AZ}}}= B_{z0}/\sqrt{4\pi \rho }$ set by the initial vertical magnetic field. This is in accordance with SMUN00 and IS05,  though some literature adopts a different definition (e.g., R_M ≡ c_s²/η₀Ω in Fleming et al. 2000).

We have used Athena (Gardiner & Stone 2005, 2008), an open-source MHD code, for our simulations.

We vary the initial magnetic field strength and the diffusivity models and we classify each set of parameters as either ◯active zone, ×dead zone, or ▵sustained zone. The experiment method and the definition of the three classes are given in this section.

The parameters we have investigated are the initial vertical field strength (represented by plasma β), the linear diffusivity (represented by magnetic Reynolds number R_M), and the critical current J_crit. The range of the survey was 400 ⩽ β ⩽ 25600, 0.002 ⩽ R_M ⩽ 2, and 0.01 ⩽ J_crit/J_eqp ⩽ 100. In addition, the limiting cases of R_M = ∞ and J_crit/J_eqp = ∞, which correspond to ideal MHD models and linear Ohm's law models, respectively, are studied for comparison with the literature.

For each value of β, we prepared the initial condition as described in Section 2.1 and continued the simulation for 10 orbits (t = 20π/Ω), at first with magnetic diffusivity turned off (η(J) = 0). While running the simulations, we created restart data for each periodic point (t = 2nπ/Ω, where n is an integer). The condition allowed the MRI to grow and saturate in about 5 orbits (t = 10π/Ω).

Then, for each pair of (R_M, J_crit/J_eqp), we turned on the diffusivity and restarted the simulation either from the initial laminar flow (t = 0) or the saturated MRI states at 8, 9, and 10 orbits (t = 16π/Ω, 18π/Ω, 20π/Ω). We numerically evolved them until they reach 20 orbits (t = 40π/Ω). The reason why we have adopted three different MRI saturated initial conditions (t = 16π/Ω, 18π/Ω, 20π/Ω) for each set of the parameters is that a "turbulent initial condition" is not unique, and therefore we need to test if our results depend on the choice of the initial condition or not.

During each simulation run, we recorded the space averages of physical quantities as functions of time, such as magnetic energy density B², the Reynolds and Maxwell stress ρv_xδv_y, −B_xB_y/4π, and the squared current J². After the simulations, we studied the time average of the quantities. For a physical quantity A, we denote its space and time average by 〈A〉 and $\overline{A}$, respectively. Their definitions are as follows:

$$\begin{eqnarray} \langle A \rangle &\equiv & \frac{\int \mathit {dx}\int \mathit {dy}\int \mathit {dz}A}{\int \mathit {dx}\int \mathit {dy}\int \mathit {dz}} \end{eqnarray} \tag{ 18 }$$

and

$$\begin{eqnarray} \overline{A} &\equiv & \frac{\int \mathit {dt} A}{\int \mathit {dt}}. \end{eqnarray} \tag{ 19 }$$

The space average is taken for the entire computational domain (0 < x < L_x, 0 < y < L_y, 0 < z < L_z), the numerical resolution is (N_x, N_y, N_z) = (64, 128, 64), and the time average is taken for the last five orbits (30π < t < 40π) unless otherwise mentioned. Statistics for the important sets of parameters are presented at the end of the paper.

Using the average values, we classify each set of the parameter (β, R_M, J_crit/J_eqp) as follows (cf. Table 2): First, a parameter is in ◯active zone if the MRI is observed both in the simulation started from the laminar flow as well as in all of the three simulations started from the saturated MRI states. Second, a parameter is in ×dead zone if the instability is observed neither in the simulation started from the laminar flow nor in any of the three simulations started from the saturated MRI states. Finally, a parameter is in ▵sustained zone if the MRI is observed in all the three simulations started from the saturated MRI states, but not at the end of the simulation started from the laminar flow.

The typical behavior of the current for the active zone, dead zone, and sustained zone is in Figure 2. From the simulations, we have observed that the three classes (active, dead, and sustained) are exhaustive: the runs started from 8, 9, and 10 orbits always agree in terms of the classification and if the MRI dies when starting from a saturated initial condition, it also does not activate starting from a laminar initial condition.

Zoom In Zoom Out Reset image size

To classify the active, dead, and sustained zones, we need to assess the magnetorotational instability of the system, so we introduce the criteria below for quantitative assessment. We say that the system is magnetorotationally unstable if the averaged current 〈J²〉^1/2 is greater than 0.1J_eqp, and stable if otherwise. Here, the time average is taken for the last five orbits (30π < t < 40π). We have learned from the simulations that the quantity 〈J²〉^1/2 is a good indicator of stability, since it either fluctuates around the mean value 〈J²〉^1/2 ≃ 10J_eqp (unstable) or goes under 〈J²〉^1/2 < 0.1J_eqp almost monotonically with little vibration (stable). There is no ambiguity between the two (cf. Figure 2).

However, as we make R_M smaller, the diffusion timescale becomes shorter and the wall clock time for the simulations until t = 40π becomes impractically larger. Therefore, for the parameter range R_M < 0.01, we terminate the simulations at 1.5 × 10⁶ cycles and determine the class by extrapolations.

Figures 3 and 4 show the result of our parameter survey. We found that sustained zones exist—the MRI does exhibit hysteresis behavior for a certain set of parameters.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To study possible influences of the numerical resolutions, we have performed simulations using the different numerical resolutions N_x = 32, 48, 64, (fiducial) and 80 while keeping the aspect ratio N_x: N_y: N_z = 1: 2: 1, for three different sets of parameters (β, R_M, J_crit/J_eqp) = (400,0.6,1), (400,0.2,1), and (400,0.2,1). Figure 5 summarizes the convergence tests. Within these sets of parameters, we observe that our classification does not depend on the numerical resolution.

Zoom In Zoom Out Reset image size

We also studied if our result depends on the model of the nonlinear Ohm's law. In addition to our fiducial (fid) model (Equation (2)), we have studied two smoothly transiting models (Equations (10) and (11)). To distinguish the three models, see Figure 1.

Figure 6 summarizes the time evolution of the current density in the simulations with the three typical sets of parameters (β, R_M, J_crit/J_eqp) = (400,0.6,1), (400,0.2,1), and (400,0.2,1). Within the regime we have tested, the hysteresis behavior does not depend on the detail of the nonlinear resistivity models.

Zoom In Zoom Out Reset image size

From Figures 3 and 4, we can see the following condition for the ▵sustained zone:

$$\begin{eqnarray} \frac{J_{\mathrm{crit}}}{{J_{\mathrm{eqp}}}} \frac{1}{R_{\mathrm{M}}} &\le & f_{\mathrm{whb}}, \end{eqnarray} \tag{ 20 }$$

where f_whb is a proportionality constant that satisfies f_whb ≃ 5–15 for 400 ⩽ β ⩽ 1600, and f_whb ≃ 15–50 for β ⩾ 3200. Hereafter, we interpret this in terms of the work-heat balance in the resistive MRI.

We also remark that within parameter regions that is in active zone, R_M < 1, and with larger J_crit, we observed the large-amplitude time-variability of physical quantities, such as magnetic fields, due to repeated growth and reconnection of the channel solutions. This phenomenon is reported by Fleming et al. (2000).

In this section, we show that Equation (20) can be understood as a condition of the balance between the magnetic energy dissipated by Joule heating per unit volume (W_J) and the work done by shearing motion per unit volume (W_sh).

Let us define Q_whb as the left-hand side of Equation (20):

$$\begin{eqnarray} &&{Q_{\mathrm{whb}}}\equiv \frac{J_{\mathrm{crit}}}{{J_{\mathrm{eqp}}}} \frac{1}{R_{\mathrm{M}}}. \end{eqnarray} \tag{ 21 }$$

The condition for self-sustained MRI is Q_whb < f_whb, which will be explained in this section.

First, substitute $R_{\mathrm{M}}\equiv {v_{\mathrm{{}_{AZ}}}}^2 / \eta _0 \Omega$:

$$\begin{eqnarray} &&{Q_{\mathrm{whb}}}= \frac{ J_{\mathrm{crit}}}{{J_{\mathrm{eqp}}}} \frac{\eta _0 \Omega }{ {v_{\mathrm{{}_{AZ}}}}^2}. \end{eqnarray} \tag{ 22 }$$

Next, when the MRI is active, 〈J²〉^1/2 ≃ f_satJ_eqp where f_sat is of the order of 10. Since this average current strength lies in the super-critical regime of Ohm's law (J > J_crit), the electric field is E'_crit = 4πc⁻²η₀J_crit as modeled in Equations (1) and (2). Therefore, the Joule heating per unit volume is estimated to be

$$\begin{eqnarray} W_{\mathrm{J}}&=& {E_{\mathrm{crit}}^{\prime }}\cdot \langle J^2 \rangle ^{1/2} \nonumber \\ &=& {E_{\mathrm{crit}}^{\prime }}\cdot f_{\mathrm{sat}}{J_{\mathrm{eqp}}}\nonumber \\ & = & 4 \pi f_{\mathrm{sat}}c^{-2}\eta _0J_{\mathrm{crit}}{J_{\mathrm{eqp}}}. \end{eqnarray} \tag{ 23 }$$

To explain the work-heat balance qualitatively, this estimate needs correction due to the high space variability of the current field under the discharge conditions. We introduce f_fill, or the filling factor, which is the ratio of the volume that contributes to the Joule heating to the total volume. Formally, f_fill is defined as the ratio between volume averages of the actual Joule heat generated and the Joule heat estimated by this method:

$$\begin{eqnarray} f_{\mathrm{fill}}&\equiv & \frac{\langle E^{\prime }J \rangle }{ {E_{\mathrm{crit}}^{\prime }}\cdot \langle J^2 \rangle ^{1/2} } \nonumber \\ &=&\frac{\langle f\left(J\right) J\rangle }{\langle J^2 \rangle ^{1/2}}, \end{eqnarray} \tag{ 24 }$$

where

$$\begin{eqnarray} f\left(J\right) = \left\lbrace \begin{array}{@{} ll} 1 & \mathrm{if} \; J> J_{\mathrm{crit}}\\ \frac{J}{J_{\mathrm{crit}}} & \mathrm{if} \; J< J_{\mathrm{crit}} \end{array}. \right. \end{eqnarray} \tag{ 25 }$$

Using this f_fill, Equation (23) is rewritten as

$$\begin{eqnarray} W_{\mathrm{J}}& \equiv & 4 \pi f_{\mathrm{fill}}f_{\mathrm{sat}}c^{-2}\eta _0J_{\mathrm{crit}}{J_{\mathrm{eqp}}}. \end{eqnarray} \tag{ 26 }$$

Substituting η₀ in Equation (22) with Equation (26) gives

$$\begin{eqnarray} {Q_{\mathrm{whb}}}&=&\frac{c^2 W_{\mathrm{J}}\Omega }{ 4 \pi f_{\mathrm{fill}}f_{\mathrm{sat}}{J_{\mathrm{eqp}}}^2 {v_{\mathrm{{}_{AZ}}}}^2} \nonumber , \\ W_{\mathrm{J}}&=& \frac{ 4 \pi f_{\mathrm{fill}}f_{\mathrm{sat}}{J_{\mathrm{eqp}}}^2 {v_{\mathrm{{}_{AZ}}}}^2 {Q_{\mathrm{whb}}}}{c^2 \Omega }. \end{eqnarray} \tag{ 27 }$$

Substituting $v_{\mathrm{{}_{AZ}}}$ with Equation (16) and then c_s with Equation (14), J_eqp with Equation (13) and then B_eqp with Equation (12) gives

$$\begin{eqnarray} W_{\mathrm{J}}&=& \frac{4 f_{\mathrm{fill}}f_{\mathrm{sat}}P_0 \Omega {Q_{\mathrm{whb}}}}{ \beta }. \end{eqnarray} \tag{ 28 }$$

On the other hand,

$$\begin{eqnarray} W_{\mathrm{sh}}&\equiv & \frac{3}{2} \Omega \langle w_{xy} \rangle \nonumber \\ &=& \frac{3}{2} \alpha \Omega P, \end{eqnarray} \tag{ 29 }$$

where α is Shakura & Sunyaev's (1973) α parameter. Substituting Equation (29) into Equation (28), one obtains

$$\begin{eqnarray} W_{\mathrm{J}}&=& \frac{8 f_{\mathrm{fill}}f_{\mathrm{sat}}{Q_{\mathrm{whb}}}}{3 \alpha \beta } W_{\mathrm{sh}}. \end{eqnarray} \tag{ 30 }$$

For the MRI to sustain itself by the discharge process, the Joule heating W_J needs to be equal to or smaller than W_sh:

$$\begin{eqnarray} W_{\mathrm{J}}&\lesssim & W_{\mathrm{sh}}. \end{eqnarray} \tag{ 31 }$$

Therefore, we have the following constraint on the left-hand side of Equation (22):

$$\begin{eqnarray} {Q_{\mathrm{whb}}}&=& \frac{3 \alpha \beta }{8 f_{\mathrm{fill}}f_{\mathrm{sat}}} \frac{W_{\mathrm{J}}}{W_{\mathrm{sh}}} \nonumber , \\ \frac{J_{\mathrm{crit}}}{{J_{\mathrm{eqp}}}} \frac{1}{R_{\mathrm{M}}} &\lesssim & \frac{3 \alpha \beta }{8 f_{\mathrm{fill}}f_{\mathrm{sat}}} \equiv f_{\mathrm{whb}}\left(\beta \right). \end{eqnarray} \tag{ 32 }$$

Thus, the work-heat balance imposes an upper limit on the product of J_crit/J_eqp and 1/R_M, provided that f_fill, f_sat, and α are constants that do not depend on J_crit and R_M, but only on β. This explains the inverse-proportional relations observed in Figures 3 and 4. The f_whb(β) values calculated with this interpretation using the experimental data are in Table 3. We list the statistics for selected set of parameters from Figures 3 and 4 in Table 4.

Table 3. The f_whb(β) Calculated from the Experimental Data

β	$\overline{\alpha }$	$\displaystyle \overline{\langle { B^2}/{8 \pi P_0} \rangle }$	$\overline{f_{\mathrm{sat}}}$	$\overline{f_{\mathrm{fill}}}$	fwhb(β)
400	0.176 ± 0.036	0.222 ± 0.040	14.5 ± 1.7	0.252 ± 0.001	7.19 ± 0.81
800	0.143 ± 0.035	0.206 ± 0.057	17.4 ± 2.7	0.259 ± 0.001	9.33 ± 1.03
1600	0.104 ± 0.027	0.166 ± 0.047	17.4 ± 1.8	0.265 ± 0.002	13.3 ± 2.3
3200	0.0523 ± 0.0206	0.0916 ± 0.0411	12.2 ± 3.9	0.262 ± 0.008	19.0 ± 3.9
6400	0.0236 ± 0.0041	0.0386 ± 0.0072	9.58 ± 1.04	0.263 ± 0.003	22.2 ± 1.6
12800	0.0182 ± 0.0032	0.0300 ± 0.0060	9.84 ± 0.88	0.270 ± 0.000	32.6 ± 2.8
25600	0.0185 ± 0.0063	0.0315 ± 0.0128	10.7 ± 1.5	0.275 ± 0.001	58.5 ± 12.4

Notes. We first calculated the time and space averaged quantities $\overline{\alpha }$ and $\overline{f_{\mathrm{sat}}}$ for each run. Then for the ensemble of runs, we calculated the means and standard deviations of the quantities. The ensemble constitutes runs that (1) belong to the ▵sustained zone, (2) are restarted runs (t = 16π/Ω, 18π/Ω, 20π/Ω) so that they are magnetorotationally unstable, and (3) have the largest product J_crit/J_eqp · 1/R_M so that they face the ▵sustained zone—×dead zone boundary.

Download table as:  ASCIITypeset image

Table 4. Statistics of the Local Simulations Abridged

runID	t	β	$\displaystyle {J_{\mathrm{crit}}}/{{J_{\mathrm{eqp}}}}$	RM	$\displaystyle \langle { B^2}/{8 \pi P_0} \rangle$	$\displaystyle \langle { -B_x B_y}/{4 \pi P_0} \rangle$	$\displaystyle \langle {\rho v_x \delta v_y}/{P_0} \rangle$	$\displaystyle \langle {J^2}/{{{J_{\mathrm{eqp}}}}^2} \rangle ^{0.5}$
3160	0	400	∞	∞	(1.99 ± 0.60) × 10−1	(1.05 ± 0.37) × 10−1	(3.01 ± 1.49) × 10−2	(2.04 ± 0.24) × 101
10	0	800	∞	∞	(1.64 ± 0.69) × 10−1	(8.27 ± 3.65) × 10−2	(2.26 ± 0.96) × 10−2	(1.94 ± 0.30) × 101
220	0	1600	∞	∞	(1.94 ± 0.61) × 10−1	(9.47 ± 3.19) × 10−2	(2.34 ± 0.84) × 10−2	(2.07 ± 0.26) × 101
3180	0	3200	∞	∞	(1.30 ± 0.27) × 10−1	(6.08 ± 1.47) × 10−2	(1.68 ± 0.52) × 10−2	(1.85 ± 0.18) × 101
3612	0	6400	∞	∞	(7.16 ± 2.57) × 10−2	(3.14 ± 0.85) × 10−2	(8.28 ± 2.51) × 10−3	(1.51 ± 0.15) × 101
3616	0	12800	∞	∞	(3.61 ± 0.83) × 10−2	(1.70 ± 0.35) × 10−2	(6.15 ± 2.40) × 10−3	(1.24 ± 0.10) × 101
3620	0	25600	∞	∞	(3.61 ± 0.84) × 10−2	(1.63 ± 0.36) × 10−2	(5.40 ± 1.84) × 10−3	(1.22 ± 0.09) × 101
3292	0	400	1.0	0.6	(2.92 ± 1.56) × 10−1	(1.67 ± 1.04) × 10−1	(4.81 ± 3.63) × 10−2	(2.10 ± 0.43) × 101	⌉
3293	16π/Ω	400	1.0	0.6	(2.74 ± 1.43) × 10−1	(1.62 ± 0.99) × 10−1	(4.63 ± 3.05) × 10−2	(2.06 ± 0.48) × 101	(a)
3294	18π/Ω	400	1.0	0.6	(2.07 ± 0.71) × 10−1	(1.21 ± 0.44) × 10−1	(3.68 ± 1.62) × 10−2	(1.88 ± 0.27) × 101	◯
3295	20π/Ω	400	1.0	0.6	(2.34 ± 0.82) × 10−1	(1.36 ± 0.54) × 10−1	(3.86 ± 1.80) × 10−2	(2.00 ± 0.30) × 101	⌋
3352	0	400	1.0	0.2	(2.50 ± 0.00) × 10−3	(8.24 ± 6.55) × 10−13	(4.83 ± 2.58) × 10−7	(4.11 ± 0.94) × 10−3	⌉
3353	16π/Ω	400	1.0	0.2	(2.03 ± 1.27) × 10−1	(1.27 ± 0.86) × 10−1	(4.07 ± 2.99) × 10−2	(1.50 ± 0.60) × 101	(b)
3354	18π/Ω	400	1.0	0.2	(2.21 ± 0.98) × 10−1	(1.33 ± 0.50) × 10−1	(4.68 ± 2.68) × 10−2	(1.66 ± 0.41) × 101	▵
3355	20π/Ω	400	1.0	0.2	(2.52 ± 1.65) × 10−1	(1.56 ± 1.02) × 10−1	(4.57 ± 2.72) × 10−2	(1.67 ± 0.59) × 101	⌋
3348	0	400	10.0	0.2	(2.50 ± 0.00) × 10−3	(8.24 ± 6.55) × 10−13	(4.83 ± 2.58) × 10−7	(4.11 ± 0.94) × 10−3	⌉
3349	16π/Ω	400	10.0	0.2	(2.50 ± 0.00) × 10−3	(1.77 ± 1.01) × 10−7	(9.90 ± 7.50) × 10−5	(3.60 ± 0.54) × 10−2	(c)
3350	18π/Ω	400	10.0	0.2	(2.50 ± 0.00) × 10−3	(2.31 ± 1.64) × 10−6	(6.06 ± 9.16) × 10−5	(3.49 ± 0.54) × 10−2	×
3351	20π/Ω	400	10.0	0.2	(2.50 ± 0.00) × 10−3	(4.85 ± 3.58) × 10−7	(6.62 ± 8.64) × 10−5	(3.43 ± 0.36) × 10−2	⌋
3352	0	400	1.0	0.2	(2.50 ± 0.00) × 10−3	(8.24 ± 6.55) × 10−13	(4.83 ± 2.58) × 10−7	(4.11 ± 0.94) × 10−3	}▵
3353	16π/Ω	400	1.0	0.2	(2.03 ± 1.27) × 10−1	(1.27 ± 0.86) × 10−1	(4.07 ± 2.99) × 10−2	(1.50 ± 0.60) × 101
3452	0	400	0.9	0.06	(2.50 ± 0.00) × 10−3	(3.05 ± 2.55) × 10−13	(4.29 ± 2.47) × 10−7	(2.22 ± 0.52) × 10−3	} ×
3453	16π/Ω	400	0.9	0.06	(2.50 ± 0.00) × 10−3	(1.22 ± 0.75) × 10−8	(2.29 ± 1.58) × 10−4	(3.02 ± 0.46) × 10−2
3448	0	400	0.3	0.06	(2.50 ± 0.00) × 10−3	(3.05 ± 2.55) × 10−13	(4.29 ± 2.47) × 10−7	(2.22 ± 0.52) × 10−3	}▵
3449	16π/Ω	400	0.3	0.06	(2.29 ± 2.53) × 10−1	(1.43 ± 1.71) × 10−1	(4.61 ± 4.95) × 10−2	(1.39 ± 0.90) × 101
3436	0	400	0.3	0.02	(2.50 ± 0.00) × 10−3	(9.13 ± 12.07) × 10−13	(1.20 ± 1.54) × 10−7	(3.43 ± 2.16) × 10−4	} ×
3437	16π/Ω	400	0.3	0.02	(2.50 ± 0.00) × 10−3	(2.39 ± 1.08) × 10−9	(8.02 ± 6.85) × 10−5	(1.17 ± 0.23) × 10−2
3432	0	400	0.1	0.02	(2.50 ± 0.00) × 10−3	(9.55 ± 12.38) × 10−13	(8.32 ± 8.82) × 10−8	(2.91 ± 1.54) × 10−4	}▵
3433	16π/Ω	400	0.1	0.02	(2.00 ± 2.42) × 10−1	(1.14 ± 1.50) × 10−1	(3.54 ± 4.61) × 10−2	(1.23 ± 0.98) × 101
3372	0	3200	1.0	0.2	(1.11 ± 0.62) × 10−1	(5.22 ± 2.79) × 10−2	(1.39 ± 0.68) × 10−2	(1.58 ± 0.31) × 101	}◯
3373	16π/Ω	3200	1.0	0.2	(1.00 ± 0.41) × 10−1	(4.29 ± 1.61) × 10−2	(9.90 ± 3.04) × 10−3	(1.52 ± 0.24) × 101
3460	0	3200	0.9	0.06	(3.12 ± 0.00) × 10−4	(3.02 ± 1.86) × 10−12	(4.68 ± 2.21) × 10−7	(1.89 ± 0.39) × 10−3	} ×
3461	16π/Ω	3200	0.9	0.06	(3.35 ± 0.24) × 10−4	(9.82 ± 10.64) × 10−6	(3.36 ± 3.24) × 10−5	(5.07 ± 2.78) × 10−2
3456	0	3200	0.3	0.06	(3.12 ± 0.00) × 10−4	(3.02 ± 1.86) × 10−12	(4.68 ± 2.21) × 10−7	(1.89 ± 0.39) × 10−3	}▵
3457	16π/Ω	3200	0.3	0.06	(1.19 ± 0.48) × 10−1	(5.44 ± 1.88) × 10−2	(1.42 ± 0.51) × 10−2	(1.65 ± 0.26) × 101
3444	0	3200	0.3	0.02	(3.12 ± 0.00) × 10−4	(8.60 ± 6.12) × 10−13	(5.18 ± 2.54) × 10−7	(1.55 ± 0.35) × 10−3	}▵
3445	16π/Ω	3200	0.3	0.02	(9.38 ± 4.06) × 10−2	(4.50 ± 1.88) × 10−2	(1.15 ± 0.46) × 10−2	(1.29 ± 0.27) × 101
3440	0	3200	0.1	0.02	(3.12 ± 0.00) × 10−4	(8.60 ± 6.12) × 10−13	(5.18 ± 2.54) × 10−7	(1.55 ± 0.35) × 10−3	}▵
3441	16π/Ω	3200	0.1	0.02	(8.15 ± 3.91) × 10−2	(3.79 ± 1.70) × 10−2	(1.05 ± 0.50) × 10−2	(1.41 ± 0.23) × 101

Notes. Each run is labeled by an integer. The restart time is in the second column. The next three columns indicate the initial magnetic field strength, the critical current, and the magnetic Reynolds number. The physical quantities are represented in terms of the time average and standard deviation of the space average, i.e., A is in the format $\overline{\langle A \rangle } \pm (\overline{{\langle A \rangle }^2} - {\overline{\langle A \rangle }}^2)^{0.5}$. In this table, there are runs for ideal MHD; runs in Figure 2 that represent the behavior in (a) active, (b) sustained, and (c) dead zones; and runs that constitute sustained-dead zone boundaries for β = 400, 3200.

Download table as:  ASCIITypeset image

We can further simplify Equation (32) by using the saturation predictor proposed by HGB95. The proposed predictors read

$$\begin{eqnarray} \alpha P &=& 0.61 \pm 0.06 \overline{\left\langle \frac{B^2}{8 \pi } \right\rangle } \end{eqnarray} \tag{ 33 }$$

and

$$\begin{eqnarray} \overline{\left\langle \frac{B^2}{8 \pi } \right\rangle } &=& \left(1.21 \pm 0.29\right) \cdot 2 \pi \sqrt{\frac{16}{15} \cdot \frac{2}{\beta }} P_0. \end{eqnarray} \tag{ 34 }$$

Using this, Equation (32) is rewritten as follows:

$$\begin{eqnarray} \frac{J_{\mathrm{crit}}}{{J_{\mathrm{eqp}}}} \frac{1}{R_{\mathrm{M}}} &\lesssim & f_{\mathrm{whb}}\left(\beta \right) \nonumber \\ &\simeq &\left(2.54 \pm 0.66\right) \frac{\beta ^{1/2}}{f_{\mathrm{fill}}f_{\mathrm{sat}}}. \end{eqnarray} \tag{ 35 }$$

By ignoring the dependence of f_fill and f_sat on β, we assume that $\overline{f_{\mathrm{fill}}}= 0.264 \pm 0.007$ and $\overline{f_{\mathrm{sat}}} = 13.1 \pm 3.1$. This further simplifies Equation (32) as

$$\begin{eqnarray} f_{\mathrm{whb}}\left(\beta \right) & \simeq &\left(0.74 \pm 0.26\right) {\beta ^{1/2}}. \end{eqnarray} \tag{ 36 }$$

This is in agreement with our experimental data (Table 3) within a factor of two.

In the previous section, we performed the shearing-box simulations of MRI with nonlinear Ohm's law and found the three classes of MRI behavior; we named them active, dead and sustained zones. We also found the condition for the MRI to be self-sustained. In this section, we apply the findings to the global model of protoplanetary disks and analyze how they are divided into the three zones.

We use the minimum-mass solar nebula model introduced by Hayashi et al. (1985) as the fiducial disk model:

$$\begin{eqnarray} &&\Sigma = {f_\Sigma }{\Sigma _0}{\left(\frac{r}{{\mathrm{AU}}}\right)}^ {-q} , \end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray} &&T= T_0{\left(\frac{r}{{\mathrm{AU}}}\right)}^ {-{{1}/{2}}}. \end{eqnarray} \tag{ 38 }$$

Here, Σ₀ = 1.7 × 10³ g cm⁻³ and T₀ = 2.8 × 10² K are the surface density and the temperature at 1AU, respectively. ${f_\Sigma }$ is the nondimensional surface density parameter. The fiducial value for the surface density power index is q = 3/2. Since we assume the isotropic EOS, the ratio of specific heats γ = 1 in our model and the thermal velocities for gas molecules and electrons are

$$\begin{eqnarray} c_s\left(r\right) &=& \sqrt{\frac{ {k_B}T}{\mu {m_H}}} \end{eqnarray} \tag{ 39 }$$

$$\begin{eqnarray} v_e\left(r\right) &=& \sqrt{\frac{ {k_B}T}{{m_e}}}, \end{eqnarray} \tag{ 40 }$$

respectively, where μ is the mean molecular weight of the gas.

Since we neglect the self gravity of the disk, the disk is in Kepler rotation and its orbital angular velocity is

$$\begin{eqnarray} \Omega \left(r\right) &=& \sqrt{\frac{{G}M_{*}}{ r^3 }}. \end{eqnarray} \tag{ 41 }$$

From the equilibrium between the vertical pressure gradient and vertical component of the stellar gravity, the disk density and pressure distributions are

$$\begin{eqnarray} \rho \left(r,z\right) &=& \frac{\Sigma }{ \sqrt{2 \pi } H} \exp \left({\frac{- z^2 }{2 H^2}}\right) \end{eqnarray} \tag{ 42 }$$

and

$$\begin{eqnarray} P\left(r,z\right) &=& \frac{\rho {k_B}T}{\mu {m_H}}, \end{eqnarray} \tag{ 43 }$$

where H is the definition of the disk scale height in this paper (cf. Equation (14)).

For simplicity, we assume every dust particle to be a solid sphere of equal radius a_d and density ρ_s. The mass m_d and geometrical cross section σ_d of the dust particle are

$$\begin{eqnarray} {m_d}&=& \frac{4 \pi }{3} \rho _sa_d^3 , \end{eqnarray} \tag{ 44 }$$

$$\begin{eqnarray} \sigma _d&=& \pi a_d^2, \end{eqnarray} \tag{ 45 }$$

respectively. The fiducial values are a_d = 0.1 μm and ρ_s = 3 g cm⁻³.

Using this m_d and dust-to-gas density ratio f_d = 0.01, the number densities of the dust and gas components are

$$\begin{eqnarray} n_g\left(r,z\right) &=& \frac{\rho }{\mu {m_H}} \end{eqnarray} \tag{ 46 }$$

and

$$\begin{eqnarray} n_d\left(r,z\right) &=& \frac{f_d\hspace{2.84526pt} \rho }{{m_d}} . \end{eqnarray} \tag{ 47 }$$

Methods for calculating the charge equilibrium of the dust-plasma in protoplanetary disks has been studied (Umebayashi & Nakano 1980, 2009; Fujii et al. 2011). Among those we use O09's method because of its numerical efficiency and generality. First, the gas column density above ($\chi _{{}_\top }$) and below ($\chi _{{}_\bot }$) the coordinate (r, z) in the disk are

$$\begin{eqnarray} \chi _{{}_\top }\left(r,z\right) &=& \int _z^{\infty } \rho \mathit {dz} \nonumber \\ &=& \frac{\Sigma }{2} \left[1 - {\mathrm{erf}}\left(\frac{z}{\sqrt{2} H}\right) \right] \end{eqnarray} \tag{ 48 }$$

and

$$\begin{eqnarray} \chi _{{}_\bot }\left(r,z\right) &=& \int _{-\infty }^z \rho \mathit {dz} \nonumber \\ &=& \frac{\Sigma }{2} \left[1 + {\mathrm{erf}}\left(\frac{z}{\sqrt{2} H}\right) \right], \end{eqnarray} \tag{ 49 }$$

respectively, where

$$\begin{eqnarray} &&{\mathrm{erf}}\left(x\right) \equiv \frac{2}{\sqrt{\pi }} \int _x^{\infty } e^{-t^2} \mathit {dt} \end{eqnarray} \tag{ 50 }$$

is the error function.

According to SMUN00, the effective ionization rate in the disk is

$$\begin{eqnarray} \zeta \left(r,z\right) &\approx & \frac{\zeta _{\mathrm{CR}}}{2} \left\lbrace \exp \left(-\frac{\chi _{{}_\top }}{ \chi _{\mathrm{CR}}}\right) \right. \nonumber \\ && + \left. \exp \left(-\frac{\chi _{{}_\bot }}{\chi _{\mathrm{CR}}}\right) \right\rbrace + \zeta _{\mathrm{RA}}, \end{eqnarray} \tag{ 51 }$$

where χ_CR = 96 g cm⁻², ζ_CR = 1.0 × 10⁻¹⁷ s⁻¹, and ζ_RA = 6.9 × 10⁻²³ s⁻¹.

Using this, O09's nondimensional parameter Θ is calculated as

$$\begin{eqnarray} &&\Theta = \frac{\zeta n_ge^2}{s_ic_s\sigma _da_dn_d^2 {k_B}T} \end{eqnarray} \tag{ 52 }$$

and Γ is defined as the solution of the equation

$$\begin{eqnarray} &&\frac{1}{1+\Gamma } - \left(\frac{s_i}{s_e}\sqrt{\frac{{m_e}}{\mu {m_H}}} \exp \Gamma + \frac{\Gamma }{\Theta }\right) = 0. \end{eqnarray} \tag{ 53 }$$

Once Equation (53) is numerically solved for Γ, we can calculate the number density of ions and electrons, n_i and n_e, as well as the rms of the charge per dust particle, $\sqrt{{\langle Z^2 \rangle }} e$, as

$$\begin{eqnarray} n_i\left(r,z\right) &=& \frac{\zeta n_g}{ s_ic_s\sigma _dn_d\left(1+\Gamma \right)}, \end{eqnarray} \tag{ 54 }$$

$$\begin{eqnarray} n_e\left(r,z\right) &=& \frac{\zeta n_g\exp \Gamma }{ s_ev_e\sigma _dn_d}, \end{eqnarray} \tag{ 55 }$$

and

$$\begin{eqnarray} {\langle Z^2 \rangle }&=& \left(\frac{\Gamma a_d}{\lambda }\right)^2 + \frac{1+\Gamma }{2+\Gamma }\frac{a_d}{\lambda }\\ \mathrm{where} \lambda &=&\frac{e^2}{{k_B}T}.\nonumber \end{eqnarray} \tag{ 56 }$$

We assume sticking probabilities s_i = 1 and s_e = 0.3 as assumed by O09.

Next, we estimate the plasma conductivity using the method of SMUN00. The rate coefficient for the collision between the neutrals and the ions is

$$\begin{eqnarray} &&{\langle \sigma v\rangle _i}= 2.41 \pi \left(\frac{\alpha e^2 }{{\mu m_H}}\right)^{{1}/{2}}. \end{eqnarray} \tag{ 57 }$$

We use α = 7.66 × 10⁻²⁵ cm³ as an averaged polarizability. 〈σv〉_e is the rate coefficient for collision between neutrals and electrons, whose form is given in SMUN00. The rate coefficient for dust particles is

$$\begin{eqnarray} &&{\langle \sigma v\rangle _d}= \frac{4\pi }{3}a_d^2 c_s. \end{eqnarray} \tag{ 58 }$$

This expression is valid as long as a_d is much smaller than the mean free path of the gas molecules.

With these, the magnetic diffusivity is calculated componentwise:

$$\begin{eqnarray} {\eta _e}&=& \frac{c^2 {m_e}n_g{\langle \sigma v\rangle _e}}{4 \pi e^2 n_e} , \end{eqnarray} \tag{ 59 }$$

$$\begin{eqnarray} {\eta _i}&=& \frac{c^2 {\mu m_H}n_g{\langle \sigma v\rangle _i}}{4 \pi e^2 n_i} , \end{eqnarray} \tag{ 60 }$$

$$\begin{eqnarray} {\eta _d}&=& \frac{c^2 {\mu m_H}n_g{\langle \sigma v\rangle _d}}{4 \pi {\langle Z^2 \rangle }e^2 n_d} , \end{eqnarray} \tag{ 61 }$$

and

$$\begin{eqnarray} \eta _0&=& \big ({{\eta _e}}^{-1} + {{\eta _i}}^{-1} + {{\eta _d}}^{-1}\big)^{-1}. \end{eqnarray} \tag{ 62 }$$

The value of E'_crit is set by the condition that the kinetic energy of the electrons accelerated by the electric field is large enough to initiate an electron avalanche:

$$\begin{eqnarray} \Delta W& = &f_\mathrm{DP}\hspace{2.84526pt} e\hspace{2.84526pt} {E_{\mathrm{crit}}^{\prime }}\hspace{2.84526pt} l_{\mathrm{mfp}} ; \nonumber \\ {E_{\mathrm{crit}}^{\prime }}&=& \frac{\Delta W}{f_\mathrm{DP}\hspace{2.84526pt} e\hspace{2.84526pt} l_{\mathrm{mfp}}} \nonumber \\ &=& \frac{\Delta Wn_g{\langle \sigma v\rangle _e}}{f_\mathrm{DP}e v_e}. \end{eqnarray} \tag{ 63 }$$

Here, $f_\mathrm{DP}= 0.43 \sqrt{{\mu m_H}/ {m_e}}$ is the coefficient for the average energy of electrons in weakly ionized plasma (IS05), and l_mfp = v_e/(n_g〈σv〉_e) is the mean free path of electrons. For ionization energy, we use the value for a hydrogen molecule ΔW = 15.4 eV.

With this, the critical current is

$$\begin{eqnarray} &&J_{\mathrm{crit}}= \frac{c^2}{4 \pi \eta _0} {E_{\mathrm{crit}}^{\prime }}. \end{eqnarray} \tag{ 64 }$$

Note that the discharge current, Equation (63), is calculated using the strong electric field limit of the electron distribution function, while we used the weak field limit formulae for the charge distributions, Equations (52)–(58). In this paper, we adopt this treatment for simplicity. A more consistent treatment will be the topic of another paper in preparation.

Now we study the distribution of active, sustained, and dead zones in the protoplanetary disk models. SMUN00 give the condition for the MRI unstable region as follows:

$$\begin{eqnarray} &&\frac{2\pi v_{\mathrm{{}_{AZ}}}}{\Omega } \le \sqrt{2} H\bigwedge \frac{2\pi \eta }{v_{\mathrm{{}_{AZ}}}} \le \sqrt{2} H . \end{eqnarray} \tag{ 65 }$$

Combining this with the work-heat balance model of Equation (20) gives the following conditions for active, sustained, and dead zones, respectively:

$$\begin{eqnarray} & &\hspace{-19.91684pt}\frac{2\pi v_{\mathrm{{}_{AZ}}}}{\Omega } \le \sqrt{2} H \bigwedge \frac{2\pi \eta }{v_{\mathrm{{}_{AZ}}}} \le \sqrt{2} H, \end{eqnarray} \tag{ 66 }$$

$$\begin{eqnarray} & &\hspace{-19.91684pt} \frac{2\pi v_{\mathrm{{}_{AZ}}}}{\Omega } \le \sqrt{2} H \bigwedge \frac{2\pi \eta }{v_{\mathrm{{}_{AZ}}}} > \sqrt{2} H \bigwedge \frac{J_{\mathrm{crit}}}{{J_{\mathrm{eqp}}}} \frac{1}{R_{\mathrm{M}}} \le f_{\mathrm{whb}}, \nonumber \\ \end{eqnarray} \tag{ 67 }$$

and

$$\begin{eqnarray} & &\hspace{-19.91684pt} \frac{2\pi v_{\mathrm{{}_{AZ}}}}{\Omega } > \sqrt{2} H \bigvee \left(\frac{2\pi \eta }{v_{\mathrm{{}_{AZ}}}} > \sqrt{2} H \bigwedge \frac{J_{\mathrm{crit}}}{{J_{\mathrm{eqp}}}} \frac{1}{R_{\mathrm{M}}} > f_{\mathrm{whb}}\right). \nonumber \\ \end{eqnarray} \tag{ 68 }$$

Using Equations (66)–(68), we plot the active, sustained, and dead zones for various global disk model Equations (37) and (38).

First, Figure 7 shows the unstable zones for varying disk surface density, ${f_\Sigma }= 0.3$, ${f_\Sigma }= 1$ (the fiducial model), ${f_\Sigma }= 3$, and ${f_\Sigma }= 10$. In this figure and following figures, the thick curves are the boundary of the work-heat balance model Equation (20), while the thin curves are the boundary of the instability condition in the resistive limit, i.e., the second condition in Equation (65). The solid and dashed curves correspond to the plasma beta at the midplane β = 100 and β = 1000, respectively. The active zones are marked by meshes and the sustained zones are marked by horizontal stripes.

Zoom In Zoom Out Reset image size

Figure 8 shows the active and sustained zones for dust-to-gas ratio f_d = 1, f_d = 0.1, f_d = 0.01 (the fiducial model), and f_d = 10⁻⁴. The work-heat balance condition Equation (20) is not affected by changing the dust properties such as the dust-to-gas ratio f_d or dust size a_d. One can understand this by rewriting the condition Equation (22) in the following form:

$$\begin{eqnarray} &&\frac{ {E_{\mathrm{crit}}^{\prime }}\Omega }{4 \pi c^{-2} {J_{\mathrm{eqp}}}{v_{\mathrm{{}_{AZ}}}}^2} \le f_{\mathrm{whb}}. \end{eqnarray} \tag{ 69 }$$

This form does not include a term affected by the dust properties, such as the magnetic diffusivity. On the other hand, E'_crit is inversely proportional to the electron mean free path l_mfp, and is proportional to the gas number density.

Zoom In Zoom Out Reset image size

In Figure 9, we study the evolution of the active and sustained zones as the gas density becomes lower while the dust density is kept constant. We change the set $({f_\Sigma }, f_d)$ from (1, 0.01) (the fiducial model) to (0.1, 0.1), (0.01, 1), and (10⁻⁴, 100). In Figure 9, zones are marked for β = 1000. The midplane of the disk between the radii 2 AU–20 AU becomes the sustained zone as the gas density becomes 10⁻⁴ times the fiducial model.

Zoom In Zoom Out Reset image size