THE NONLINEAR OHM'S LAW: PLASMA HEATING BY STRONG ELECTRIC FIELDS AND ITS EFFECTS ON THE IONIZATION BALANCE IN PROTOPLANETARY DISKS

How protoplanetary disks form and evolve is a key question of planet formation studies. It is generally accepted that magnetic fields play many important roles in these processes. Coupling between a disk and a magnetic field induces magnetorotational instability (MRI; Balbus & Hawley 1991), and turbulence driven by MRI provides a high effective viscosity that allows disk accretion (e.g., Hawley et al. 1995; Fromang & Nelson 2006; Flock et al. 2011). A large-scale magnetic field threading a disk also drives outflows from disk surfaces via the recurrent breakup of MRI modes (Suzuki & Inutsuka 2009, 2014) and/or the magnetocentrifugal mechanism (e.g., Blandford & Payne 1982; Spruit 1996; Lesur et al. 2013; Fromang et al. 2013; Bai & Stone 2013a). These outflows can significantly affect disk structure and planet formation in inner disk regions (Suzuki et al. 2010).

In protoplanetary disks, these magnetic activities are strongly subject to non-ideal magnetohydrodynamics (MHD) effects simply because the ionization fraction of the disk gas is low. Thermal ionization is effective only in innermost disk regions where the gas temperature exceeds 1000 K (Umebayashi 1983). Further out, the disk gas is only weakly ionized by external ionizing sources including cosmic rays (Umebayashi & Nakano 1981) and stellar X-rays (Glassgold et al. 1997). The resulting high Ohmic conductivity yields an MRI stable "dead zone" near the midplane (Gammie 1996; Sano et al. 2000). Ambipolar diffusion also suppresses MRI near the disk surface (Perez-Becker & Chiang 2011a; Bai 2011a; Simon et al. 2013a, 2013b), and the combined effect of Ohmic and ambipolar diffusion can render MRI inactive at all altitudes in inner disk regions (Bai & Stone 2013b; Bai 2013; Lesur et al. 2014). The high diffusivities also cause the loss of large-scale magnetic fields that are required for magnetocentrifugal outflow (e.g., Lubow et al. 1994; Okuzumi et al. 2014; Guilet & Ogilvie 2014). Hall drift can either enhance or suppress the magnetic activity of the disks depending on the polarity of the large-scale magnetic field relative to disk's rotation axis (Wardle & Salmeron 2012; Lesur et al. 2014; Bai 2014).

At the same time, these magnetic activities can influence the ionization state of the disk gas. MRI turbulence transports ionized gas to less ionized regions, and this process revives MRI in dead zones under favorable conditions (Inutsuka & Sano 2005; Turner et al. 2007; Ilgner & Nelson 2008). Joule heating by MRI turbulence can change the temperature profile of the disks and even the location of the dead-zone inner edge (Latter & Balbus 2012; Faure et al. 2014). On smaller scales, strong current sheets produced by MRI can locally heat up the disk gas, which can even affect the ionization state if the background temperature is near the thermal ionization threshold (Hubbard et al. 2012; McNally et al. 2013, 2014).

This study focuses on the role of strong electric fields in the ionization balance. Since the electric conductivity of the disk gas is finite, the coupling between the moving gas and magnetic fields inevitably produces a nonzero electric field in the comoving frame of the gas. The field induces systematic drift motions of ionized gas particles, whose effect is expressed by the conventional Ohm's law in which the electric current is linearly proportional to the electric field strength. However, what is largely unappreciated is that in an weakly ionized plasma, an electric field also induces random motion of plasma particles if the field is sufficiently strong (Druyvesteyn & Penning 1940; Golant et al. 1980; Lifshitz & Pitaevskii 1981). In principle, the heating of plasmas affects the chemical reactions of the plasmas, thereby affecting the ionization balance of the gas. This effect has been ignored by all previous models of disk ionization.

Inutsuka & Sano (2005) first pointed out that this electric plasma heating can occur in protoplanetary disks when MRI drives disk turbulence. They noted that MRI-driven turbulence produces a strong electric field in the neutral comoving frame when the ionization degree is low (but not too low for MRI to be active). They estimated the random energy of electrically heated free electrons, and concluded that the energy can be high enough to cause electric discharge (or electron avalanche) in the disk gas. They suggested that MRI in protoplanetary disks can be self-sustained: MRI turbulence can provide sufficient ionization to keep MRI active even in the conventional dead zones. This scenario has been recently tested by Muranushi et al. (2012) using local MHD simulations with a toy resistivity model that mimics electron discharge at high electric field strengths. They found that self-sustained MRI is realized when the work done by the turbulence exceeds the energy consumed by Joule heating.

In this paper, we study in detail the effects of plasma heating by electric fields on the ionization balance of a weakly ionized gas. While the discharge is only produced by very high electric fields, weaker fields still can heat up plasmas and can change the reaction balance. To reveal the consequences of electric plasma heating over a wide range of electric field strengths, we construct a charge reaction model properly taking into account the kinetics of weakly ionized plasmas under an electric field (Druyvesteyn & Penning 1940; Golant et al. 1980; Lifshitz & Pitaevskii 1981). Our model also includes plasma capture by small dust grains, which is essential to study the ionization balance in dense protoplanetary disks (e.g., Umebayashi 1983; Sano et al. 2000; Ilgner & Nelson 2006; Wardle 2007; Bai 2011a). As a first step, we neglect the effects of magnetic fields on the kinetics of plasmas, which means that we only treat Ohmic conductivity and neglect ambipolar diffusion and Hall drift. An extension of our model to the non-Ohmic resistivities will be done in future work.

The paper is organized as follows. In Section 2, we present an order-of-magnitude estimate of the electric field strength in MRI-driven turbulence to highlight the potential importance of electric plasma heating in weakly ionized protoplanetary disks. Our charge reaction model is described in Section 3, and results are presented in Sections 4 and 5. Section 6 discusses important implications for MHD in protoplanetary disks. A summary is given in Section 7.

Before presenting our charge reaction model, we briefly describe the basic physics of electric plasma heating in a weakly ionized gas. We will then demonstrate by simple estimations that the plasma heating can occur in protoplanetary disks under realistic conditions.

Let us consider a weakly ionized gas in which neutral gas particles are much more abundant than plasma particles. In such a gas, plasma particles collide with neutral particles much more frequently than with themselves. Therefore, if there is no externally applied field, the plasma particles tend to be thermally equilibrated with the neutrals, and their mean kinetic energy approaches that of the neutrals, 3k_BT/2, where T is the neutral temperature and k_B is the Boltzmann constant. However, if there is an applied electric field, the field accelerates the plasma particles, and some part of the gained energy is converted to their random energy after collisions with neutrals. This electric heating is particularly significant for electrons because of their high mobility and low-energy transfer efficiency in collisions with neutrals. In equilibrium, the random energy of electrons greatly exceeds that of neutrals (3k_BT/2) when E is well above the threshold

$$\begin{equation} E_{\rm crit}\equiv \sqrt{\frac{6m_e}{m_n}}\frac{k_{\rm B}T}{e\ell _e}, \end{equation} \tag{ 1 }$$

where e is the elementary charge and m_e and m_n are the masses of an electron and a neutral, respectively (Druyvesteyn & Penning 1940; Golant et al. 1980; Lifshitz & Pitaevskii 1981). The electron mean free path ℓ_e is determined by the collisions with neutrals and is given by ℓ_e = 1/(n_nσ_en), where n_n is the neutral number density and σ_en is the momentum-transfer cross section for electron–neutral collisions. Equation (1) neglects energy losses due to inelastic electron–neutral collisions (i.e., collisions that involve electronic/vibrational/rotational excitation of the neutrals) and radiative energy losses upon collisions with positive ions. The former is negligible at least at the onset of electron heating (for details, see Section 3.1) and the latter is generally negligible in weakly ionized plasmas where electron–ion collision are rare. The small factor $\sqrt{m_e/m_n}$ in Equation (1) comes from the fact that electrons lose only a small fraction (∼m_e/m_n) of their kinetic energy in a single elastic collision with a neutral (for details, see Appendix A). For an H₂ gas, σ_en ≈ 10⁻¹⁵ cm² at electron energies <10 eV (Frost & Phelps 1962; Yoon et al. 2008), so we have

$$\begin{equation} E_{\rm crit} \approx 1\times 10^{-9} \left ({\displaystyle \frac{T}{100 \,{\rm K}}} \right) \left ({\displaystyle \frac{n_n}{10^{12}\, {\rm cm^{-3}}}} \right) \,{\rm esu\, cm^{-2}}. \end{equation} \tag{ 2 }$$

Let us see whether MRI turbulence in protoplanetary disks can provide such a strong electric field. We denote the mean amplitude of the electric field in MRI turbulence, as measured in the comoving frame of the neutral gas, by E_MRI.³ To evaluate E_MRI, we use the finding by Muranushi et al. (2012) that the mean current density in MRI turbulence is insensitive to the strength of the Ohmic resistivity. They performed local unstratified resistive MHD simulation and found that the current density in MRI turbulence has a mean amplitude

$$\begin{equation} J_{\rm MRI} \approx f_{\rm sat}J_{\rm eqp}, \end{equation} \tag{ 3 }$$

where f_sat ≈ 10 is a numerical constant and

$$\begin{equation} J_{\rm eqp} \equiv \sqrt{\frac{\rho _g}{2\pi }} c\Omega \end{equation} \tag{ 4 }$$

depends only on the gas mass density ρ_g = m_nn_n and orbital frequency Ω.⁴ The value of f_sat is independent of the Ohmic resistivity as long as sustained MRI turbulence is realized.⁵ The criterion for sustained turbulence can be given in terms of the Elsasser number

$$\begin{equation} \Lambda \equiv \frac{v_{Az}^2}{\eta \Omega }, \end{equation} \tag{ 5 }$$

where η is the Ohmic resistivity and v_Az is the Alfvén speed in the direction perpendicular to the disk's midplane. MRI grows when Λ > Λ_crit, where Λ_crit ≈ 0.1–1 (e.g., Sano et al. 1998; Turner et al. 2007; Muranushi et al. 2012).

Given J_MRI, one can estimate E_MRI by using Ohm's law E = (4πη/c²)J. Here it is useful to rewrite the Ohmic diffusivity as $\eta = v_{Az}^2/\Lambda \Omega = 2c_s^2/\beta _z\Lambda \Omega$, where $\beta _z = 2c_s^2/v_{Az}^2$ is the plasma beta of the vertical magnetic field and $c_s = \sqrt{k_{\rm B}T/m_n}$ is the sound speed. If we use this expression, Ohm's law can be rewritten as

$$\begin{equation} E = \frac{8\pi }{\beta _z\Lambda \Omega } \left ({\displaystyle \frac{c_s}{c}} \right) ^2 J. \end{equation} \tag{ 6 }$$

This expression is useful because β_z ∼ 100–1000 for fully saturated MRI turbulence (e.g., Hawley et al. 1995; Miller & Stone 2000; Fromang & Papaloizou 2006; Suzuki & Inutsuka 2009). From Equation (6), we obtain

$$\begin{eqnarray} E_{\rm MRI} &=& \frac{8\pi }{\beta _z\Lambda \Omega } \left ({\displaystyle \frac{c_s}{c}} \right) ^2 J_{\rm MRI} \nonumber \\ &\approx & {\displaystyle \frac{2\times 10^{-7}}{\Lambda }} \left ({\displaystyle \frac{100}{\beta _z}} \right) \left ({\displaystyle \frac{f_{\rm sat}}{10}} \right) \left ({\displaystyle \frac{T}{100 \,{\rm K}}} \right)\nonumber\\ &&\times\, \left ({\displaystyle \frac{n_n}{10^{12}\, {\rm cm^{-3}}}} \right) ^{1/2} \,{\rm esu\, cm^{-2}}. \end{eqnarray} \tag{ 7 }$$

Thus, the ratio of E_MRI to E_crit is

$$\begin{equation} \frac{E_{\rm MRI}}{E_{\rm crit}} \approx \frac{200}{\Lambda } \left ({\displaystyle \frac{100}{\beta _z}} \right) \left ({\displaystyle \frac{f_{\rm sat}}{10}} \right) \left ({\displaystyle \frac{n_n}{10^{12}\, {\rm cm^{-3}}}} \right) ^{-1/2}. \end{equation} \tag{ 8 }$$

Note that the ratio is independent of the gas temperature T. We find that E_MRI exceeds E_crit if

$$\begin{equation} \Lambda \lesssim 200 \left ({\displaystyle \frac{100}{\beta _z}} \right) \left ({\displaystyle \frac{f_{\rm sat}}{10}} \right) \left ({\displaystyle \frac{n_n}{10^{12}\, {\rm cm^{-3}}}} \right) ^{-1/2}. \end{equation} \tag{ 9 }$$

At the same time, Λ ≳ Λ_crit is required for MRI turbulence to be sustained. Combining these two criteria, we arrive at the condition for plasma heating in MRI turbulence,

$$\begin{equation} \Lambda _{\rm crit} \lesssim \Lambda \lesssim 200 \left ({\displaystyle \frac{100}{\beta _z}} \right) \left ({\displaystyle \frac{f_{\rm sat}}{10}} \right) \left ({\displaystyle \frac{n_n}{10^{12}\, {\rm cm^{-3}}}} \right) ^{-1/2}. \end{equation} \tag{ 10 }$$

Since Λ_crit ≈ 0.1–1, both the criteria are satisfied when n_n ≲ 10¹⁴–10¹⁸ cm⁻³ for f_sat ≈ 10 and β_z ∼ 100–1000.

For protoplanetary disks, the neutral gas density at the midplane is generally supposed to be 10⁹–10¹⁵ cm⁻³. Therefore, in the presence of MRI-driven turbulence, significant plasma heating can occur in some parts of protoplanetary disks.

In this section, we introduce a charge reaction model that takes into account heating of plasmas by strong electric fields. We consider a partially ionized dusty gas consisting of neutrals, singly charged positive ions, electrons, and dust grains. The gas is assumed to be so weakly ionized that the neutral number density n_n can be approximated as a constant. The neutral gas temperature T is also assumed to be a constant by neglecting the Joule heating of the neutral gas (see McNally et al. 2014 for the possibility of significant local heating and thermal ionization of the neutral gas in MRI current sheets). We simplify the reaction network by representing the positive ions by a single dominant species (denoted by i). This one-component approach allows us to compute the ionization fraction of the gas without going into the details of the chemical composition of the ions, and also with a good accuracy in particular when dust dominates the recombination process (Okuzumi 2009; Ilgner 2012). We will assume i = HCO⁺ in this study. Dust grains are spheres of single radius a and are allowed to charge up by capturing plasma particles. For simplicity, we do not consider size distribution of the grains in this study, but it is straightforward to do so because the charge reaction only depends on a handful of moments of the size distribution (see Okuzumi 2009).

Charge reactions in the gas–dust mixture depend on the kinetic states of ions and electrons, which are described by the velocity distribution functions f_i and f_e, respectively. Unlike previous resistivity models, we give the distribution functions as a function of the electric field strength E in the neutral rest frame. We assume a steady state where acceleration by the electric field balances with energy/momentum losses upon collision with neutrals. This assumption is valid for protoplanetary disks because collisions with neutral gas particles are much more frequent than charge reactions (which are collisions between plasma particles themselves or with dust grains) and also than the evolution of the electric field (which occurs on disk's dynamical timescale). The velocity distribution functions will be presented in Section 3.1.

The charge reactions we consider are ionization by external high-energy particles (e.g., cosmic rays and X-rays), impact ionization by heated electrons, recombination of plasma particles in the gas phase, and plasma capture by dust grains. The latter three reactions depend on the velocity distributions of the plasmas and hence on the electric field strength E. Inclusion of impact ionization is essential to study electric discharge at high-field strengths. The rate equations and rate coefficients for the charge reactions will be given in Section 3.2.

One goal of this study is to reveal how the conventional Ohm's law is modified by strong electric fields. This can be done by calculating the current density as a function of the electric field strength E, and we do this in the following way (see Figure 1 and also Section 3.2). First, we calculate the mean drift velocities of plasma particles (denoted by $\langle {\boldsymbol v}_{i||} \rangle$ and $\langle {\boldsymbol v}_{e||} \rangle$) and the charge reaction rates as a function of E. We then calculate the ionization balance and obtain the number densities n_i and n_e of ions and electrons and the charge Z of dust grains in equilibrium. Finally, we obtain the current density as

$$\begin{equation} {\boldsymbol J} ={\boldsymbol J}_i + {\boldsymbol J}_e, \end{equation} \tag{ 11 }$$

$$\begin{equation} {\boldsymbol J}_\alpha \equiv q_\alpha n_\alpha \langle {\boldsymbol v}_{\alpha ||} \rangle \quad (\alpha =i,e), \end{equation} \tag{ 12 }$$

where q_i = e and q_e = −e are the charges of ions or electrons, respectively. We have neglected the contribution of charged grains to the Ohmic conductivity because it is usually small (see, e.g., Bai 2011b). Note that the resulting Ohm's law is nonlinear in E because both n_α and $\langle {\boldsymbol v}_{\alpha ||} \rangle$ depend on E,

Zoom In Zoom Out Reset image size

3.1. Velocity Distribution Functions and Their Moments

Here we describe the velocity distribution functions of plasma particles and some averaged quantities (or "velocity moments") that will be used in later steps. We denote the velocity distribution functions for ions and electrons as $f_{i}({\boldsymbol E},{\boldsymbol v}_i)$ and $f_{e}({\boldsymbol E},{\boldsymbol v}_e)$, respectively, where ${\boldsymbol v}_\alpha$ (α = i, e) is the velocity of each ionized particle. The first- and second-order moments of the distribution functions give the mean drift velocity parallel to the electric field, $\langle {\boldsymbol v}_{\alpha ||} \rangle$, and the mean kinetic energy, 〈_α〉, as

$$\begin{equation} \langle {\boldsymbol v}_{\alpha ||} \rangle = \hat{\boldsymbol E} \int (\hat{\boldsymbol E}\cdot {\boldsymbol v}_\alpha) f_\alpha ({\boldsymbol E},{\boldsymbol v}_\alpha)\, d^3v_\alpha \end{equation} \tag{ 13 }$$

and

$$\begin{equation} \langle \epsilon _\alpha \rangle = \int \epsilon _\alpha f_\alpha ({\boldsymbol E},{\boldsymbol v}_\alpha)\, d^3v_\alpha, \end{equation} \tag{ 14 }$$

respectively, where $\hat{\boldsymbol E} = {\boldsymbol E}/E$ and $\epsilon _\alpha = m_\alpha v_\alpha ^2/2$ ($v_\alpha = |{\boldsymbol v}_\alpha |$). Note that the drift velocity perpendicular to ${\boldsymbol E}$ is zero in the absence of magnetic fields (see, e.g., Nakano & Umebayashi 1986; Wardle 1999).

As mentioned earlier, we assume that acceleration by electric fields is balanced with energy/momentum losses upon collisions with neutrals. In principle, a collision with a neutral is either "elastic" or "inelastic," depending on whether the kinetic energy in the center-of-mass frame of the colliding particles is conserved or not (note, however, that both types of collisions can lead to energy loss of the charged particle in the neutral rest frame). Inelastic energy losses are due to impact excitation (rotational/vibrational/electronic) and ionization of the neutrals. However, these inelastic losses only enhance the efficiency of energy transfer from plasmas to neutrals by a factor of ≲ 10 (which is equivalent to increasing the electron mass by the same factor; see Appendix A) as long as the collision energy is ≲ 1 eV (e.g., see Figure 15 of Engelhardt & Phelps 1963). This effect is particularly negligible at the onset of plasma heating (i.e., E ∼ E_crit) where the collision energy is ∼k_BT  ∼  10⁻² eV. For this reason, we neglect all inelastic losses and only consider elastic collisions in determining the velocity distributions of plasmas. This assumption allows us to use analytic expressions for the velocity distribution functions, which we will introduce below.

3.1.1. Electrons

Having neglected inelastic energy losses, one can analytically obtain the velocity distribution function for electrons in a weakly ionized gas using the Fokker–Planck (diffusion) approximation (see Golant et al. 1980; Lifshitz & Pitaevskii 1981). In the steady state, the distribution function is given by (Davydov 1935)

$$\begin{equation} f_e({\boldsymbol E},{\boldsymbol v}_e) = \left(1-\frac{eE\ell _e}{k_{\rm B}T}\frac{\epsilon _e \hat{\boldsymbol E}\cdot \hat{\boldsymbol v}_e}{\epsilon _e+\chi k_{\rm B}T}\right)f_{e0}(E,v_e), \end{equation} \tag{ 15 }$$

$$\begin{equation} \chi \equiv \left ({\displaystyle \frac{E}{E_{\rm crit}}} \right) ^2 \end{equation} \tag{ 16 }$$

where $\hat{\boldsymbol v}_e = {\boldsymbol v}_e/v_e$ and f_e0 is the "symmetric" part of f_e that depends on the magnitudes of ${\boldsymbol E}$ and ${\boldsymbol v}_e$ but not on the angle between them ($\cos ^{-1}(\hat{\boldsymbol E}\cdot \hat{\boldsymbol v}_e)$). The exact expression of f_e0 is

$$\begin{equation} f_{e0} = \left ({\displaystyle \frac{m_e}{2\pi k_{\rm B}T}} \right) ^{3/2} \frac{({\epsilon _e}/{k_{\rm B}T}+\chi)^\chi }{W(\chi)}\exp \left(-\frac{\epsilon _e}{k_{\rm B}T}\right), \end{equation} \tag{ 17 }$$

$$\begin{equation} W(\chi) \equiv \chi ^{3/2+\chi }U\left(\frac{3}{2},\frac{5}{2}+\chi,\chi\right), \end{equation} \tag{ 18 }$$

where $U(x, y, z) \equiv \Gamma (x)^{-1}\int _0^\infty t^{x-1}(1+t)^{y-x-1}\exp (-zt)dt$ is the confluent hypergeometric function of the second kind and $\Gamma (x) \equiv \int _0^\infty t^{t-1}\exp (-t)dt$ is the Gamma function. The electron mean free path is given by ℓ_e = 1/(n_nσ_en), where we will take the momentum transfer cross section σ_en to be σ_en = 10⁻¹⁵ cm² by assuming that H₂ dominates the gas (Frost & Phelps 1962; Yoon et al. 2008). Equation (15) assumes that ℓ_e (or σ_en) is independent of v_e, which is a good assumption when the electron energy is less than 10 eV (see, e.g., Figure 2 of Frost & Phelps 1962). In the limit of weak electric fields (E ≪ E_crit), f_e0 reduces to the familiar Maxwell distribution

$$\begin{equation} f_{e0}^{\rm (M)} = \left ({\displaystyle \frac{m_e}{2\pi k_{\rm B}T}} \right) ^{3/2}\exp \left(-\frac{\epsilon _e}{k_{\rm B}T}\right). \end{equation} \tag{ 19 }$$

In the opposite limit (E ≫ E_crit), f_e0 reduces to the Druyvesteyn distribution (Druyvesteyn & Penning 1940)

$$\begin{equation} f_{e0}^{\rm (D)} = \frac{1}{\pi \Gamma \left(\frac{3}{4}\right)} \left ({\displaystyle \frac{3m_e^3}{4m_n(eE\ell _e)^2}} \right) ^{3/4} \exp \left(-\frac{3m_e \epsilon _e^2}{m_n(eE\ell _e)^2}\right). \end{equation} \tag{ 20 }$$

In Figure 2, we plot f_e0 for E = 0 and 100E_crit as a function of _e/k_BT.

Zoom In Zoom Out Reset image size

Substituting Equations (15) and (17) into Equations (13) and (14), the mean velocity and energy of electrons are analytically obtained as

$$\begin{equation} \langle {\boldsymbol v}_{e||} \rangle = -\frac{\Gamma (1+\chi,\chi)\exp \chi }{W(\chi)} \frac{e{\boldsymbol E}\ell _e}{3k_{\rm B}T} \sqrt{ {\displaystyle \frac{8 k_{\rm B}T}{\pi m_e}} }, \end{equation} \tag{ 21 }$$

and

$$\begin{equation} \langle \epsilon _e \rangle = \frac{\chi U \left(\frac{5}{2},\frac{7}{2}+\chi,\chi\right)}{U\left(\frac{3}{2},\frac{5}{2}+\chi,\chi\right)} \frac{3k_{\rm B}T}{2}, \end{equation} \tag{ 22 }$$

where $\Gamma (x,z) \equiv \int _z^\infty t^{x-1}\exp (-t)dt$ is the incomplete Gamma function. Figure 3 plots Equation (22) as a function of E/E_crit for T = 100 K.

Zoom In Zoom Out Reset image size

For later convenience, let us see how $\langle {\boldsymbol v}_{e||} \rangle$ and 〈_e〉 behave in the limits of weak and strong electric fields. Substituting Equations (19) and (20) into Equations (13) and (14), we have

$$\begin{equation} \langle {\boldsymbol v}_{e||} \rangle \approx \left\lbrace \begin{array}{@{}ll@{}} - {\displaystyle \frac{e{\boldsymbol E}\ell _e}{3k_{\rm B}T}} \sqrt{ {\displaystyle \frac{8 k_{\rm B}T}{\pi m_e}} }, & E \ll E_{\rm crit}, \\ [4mm] - {\displaystyle \frac{\sqrt{2\pi }}{3^{3/4}\Gamma \left(\frac{3}{4}\right)}} {\displaystyle \frac{\sqrt{e E \ell _e} }{(m_em_n)^{1/4}}} \hat{{\boldsymbol E}}, & E \gg E_{\rm crit}, \end{array}\right. \end{equation} \tag{ 23 }$$

$$\begin{equation} \langle \epsilon _e \rangle \approx \left\lbrace \begin{array}{@{}ll@{}} {\displaystyle \frac{3k_{\rm B}T}{2}}, & E \ll E_{\rm crit}, \\ [3mm] {\displaystyle \frac{\Gamma \left(\frac{5}{4}\right)}{\Gamma \left(\frac{3}{4}\right)}} \sqrt{ {\displaystyle \frac{m_n}{3m_e}} } eE\ell _e, & E \gg E_{\rm crit}. \end{array}\right. \end{equation} \tag{ 24 }$$

We can see three important properties of the electron velocity distribution in the strong field limit. First, the drift speed $|\langle {\boldsymbol v_{e||}} \rangle |$ is proportional to $\sqrt{E}$, not to E. The reason is that the mean free time $\Delta t_e \sim \ell _e/\sqrt{\langle v_e^2 \rangle }$ is inversely proportional to $\sqrt{E}$ in the strong field limit (this can be clearly seen by looking at the momentum conservation law of electrons; see Appendix A). We will see in the following section that the nonlinearity of Ohm's law partly comes from the nonlinearity of $|\langle {\boldsymbol v_{e||}} \rangle |$. Second, the mean electron energy is approximately given by 〈_e〉 ≈ 1.04(E/E_crit)k_BT ≈ (E/E_crit)k_BT. Thus, if T ∼ 100 K, a field of E ≈ 100E_crit gives a mean electron energy of 〈_e〉 ∼ 1 eV (see also Figure 3). Third, the kinetic energy associated with the drift motion, $m_e\langle {\boldsymbol v_{e||}} \rangle ^2/2$, is smaller than the total kinetic energy 〈_e〉 by factor (m_e/m_n)^1/2 ∼ 0.01. This means that electrons's random motion dominates over systematic motion even in the strong field limit. Thus, in weakly ionized plasmas, electric fields "heat" rather than "accelerate" electrons.

3.1.2. Ions

Unlike for electrons, there is no closed expression for the velocity distribution function of ions at high electric fields. The difference arises from the fact that ion's momentum transfer cross section depends on the ion–neutral collision velocity (instead, the mean collision time is approximately constant) owing to the polarization force between ions and neutrals (Wannier 1953). For this reason, we approximate f_i by the offset Maxwell distribution (Hershey 1939)

$$\begin{equation} f_i({\boldsymbol E},{\boldsymbol v}_i) = \left ({\displaystyle \frac{m_i}{2\pi k_{\rm B}T_i}} \right) ^{3/2} \exp \left (-\frac{m_i({\boldsymbol v}_i-\langle {\boldsymbol v}_{i||} \rangle)^2}{2k_{\rm B}T_i}\right), \end{equation} \tag{ 25 }$$

where the ion drift velocity $\langle {\boldsymbol v}_{i||} \rangle$ and ion temperature T_i are given by

$$\begin{equation} \langle {\boldsymbol v}_{i||} \rangle = \frac{m_i+m_n}{m_im_n}e{\boldsymbol E}\Delta t_i, \end{equation} \tag{ 26 }$$

$$\begin{equation} \frac{3}{2}k_{\rm B}T_i = \frac{3}{2}k_{\rm B}T + \frac{1}{2}m_n\langle {\boldsymbol v}_{i||} \rangle ^2, \end{equation} \tag{ 27 }$$

respectively. The mean free time Δt_i is the inverse of the frequency of collisions with neutrals, and is given by Δt_i = 1/K_inn_n, where K_in is the momentum transfer rate coefficient for ion–neutral collisions (assumed to be a constant). We take K_in = 1.6 × 10⁻⁹ cm³s⁻¹ following Nakano & Umebayashi (1986). From Equation (14), the mean kinetic energy is

$$\begin{eqnarray} \langle \epsilon _i \rangle &=& \frac{3}{2}k_{\rm B}T_i + \frac{1}{2}m_i \langle {\boldsymbol v}_{i||} \rangle ^2 \nonumber \\ &=& \frac{3}{2}k_{\rm B}T + \frac{1}{2}(m_n+m_i)\langle {\boldsymbol v_{i||}} \rangle ^2 \nonumber \\ &=& \frac{3}{2}k_{\rm B}T + \frac{(m_n+m_i)^3(eE\Delta t_i)^2}{2(m_im_n)^2}. \end{eqnarray} \tag{ 28 }$$

Although the distribution function given by Equation (25) is approximate, Equations (26) and (28) are the exact expressions for $\langle {\boldsymbol v}_{i||} \rangle$ and 〈_i〉 of ions having a constant mean free time (Wannier 1953, see also Appendix A).

As we will show below, electric heating is much less efficient for ions than for electrons. In the second (or third) line of Equation (28), the first and second terms account for heating by neutrals and electric fields, respectively. The second term dominates when E > E_{crit, i}, where the threshold field strength E_{crit, i} is given by

$$\begin{equation} E_{{\rm crit},i} \equiv \frac{m_i m_n \sqrt{3k_{\rm B}T}}{(m_i+m_n)^{3/2}e\Delta t_i}. \end{equation} \tag{ 29 }$$

However, E_{crit, i} is much larger than E_crit because

$$\begin{eqnarray} \frac{E_{{\rm crit},i}}{E_{\rm crit}}&=& \sqrt{\frac{m_n}{2m_ek_{\rm B}T}} \frac{K_{in} m_i m_n}{\sigma _{en} (m_i+m_n)^{3/2}} \nonumber \\ &\approx & 300 \left ({\displaystyle \frac{m_i}{30\ {\rm amu}}} \right) ^{-1/2} \left ({\displaystyle \frac{T}{100 \,{\rm K}}} \right) ^{-1/2}, \end{eqnarray} \tag{ 30 }$$

where we have assumed m_i ≫ m_n, as is the case for dominant ions in protoplanetary disks like HCO⁺. Therefore, for T ∼ 10–1000 K, ion heating becomes significant only at E ≳ 100–1000E_crit. As an example, in Figure 3, we compare 〈_i〉 with 〈_e〉 at T = 100 K. We see that ions start to be heated up only after electrons are heated to ∼1 eV.

3.2. Charge Reactions

We consider two ionizing mechanisms. One is the conventional "external" ionization by high-energy particles. The sources may include galactic cosmic rays (Umebayashi & Nakano 1981), stellar X-rays (Glassgold et al. 1997) and FUV (Perez-Becker & Chiang 2011b), and/or γ rays from radionuclides (Umebayashi & Nakano 2009). This process is characterized by a constant ionization rate ζ (the rate at which a single neutral gas particle is ionized). The second mechanism is impact ionization by electrically heated electrons. This is an "internal" ionization process in the sense that its rate is proportional to the electron number density n_e. Its rate also depends on the energy distribution of the electrons, and consequently on the strength E of the applied electric field (see Section 3.2.3). We neglect impact ionization by ions since electrons are always hotter than ions (see Section 3.1.2). We also neglect thermal ionization by assuming that the temperature T of the neutral gas is much lower than 1000 K (Umebayashi 1983). Secondary electron emission from dust grains is also neglected since it becomes important only when the electron energy is above 100 eV (Chow et al. 1993; Walch et al. 1995). Photoelectric emission from grains can become important when strong UV irradiation is present (Spitzer 1941; Weingartner & Draine 2001), but we do not consider this in this study.

Ionized particles are removed from the gas through gas-phase recombination and sticking to dust grains. By the latter process, dust grains on average obtain a negative charge because electrons have a higher random velocity than ions. In this study, we express the mean charge of the grains by eZ, where Z < 0. We will also express the mean charge in terms of the grain surface potential

$$\begin{equation} \phi = \frac{eZ}{a}, \end{equation} \tag{ 31 }$$

where a is the grain radius. Because of Coulomb interaction, the plasma accretion rates of the grains depend not only on E but also on ϕ (see Section 3.2.2).

The charge reactions mentioned above determine how n_i, n_e, and Z evolve with time t. This is described by the rate equations

$$\begin{equation} \frac{dn_i}{dt} = \zeta n_n - K_{di}(E,\phi)n_dn_i - K_{\rm rec}(E) n_i n_e + K_*(E)n_nn_e, \end{equation} \tag{ 32 }$$

$$\begin{equation} \frac{dn_e}{dt} = \zeta n_n - K_{de}(E,\phi)n_dn_e - K_{\rm rec}(E) n_i n_e + K_*(E)n_nn_e, \end{equation} \tag{ 33 }$$

$$\begin{equation} \frac{dZ}{dt} = K_{di}(E,\phi)n_i - K_{de}(E,\phi)n_e, \end{equation} \tag{ 34 }$$

where K_dα(α = i, e), K_rec, and K_* are the rate coefficients for plasma accretion by grains, gas-phase recombination, and impact ionization by electrons, respectively. Plasma heating affects the solution of these equations through the E dependences of the rate coefficients.

The set of Equations (32)–(34) has a constant of integration, ρ_c ≡ e(n_i − n_e + Zn_d), which is the net charge of the gas–dust mixture. In this study, we assume ρ_c = 0 and obtain the charge neutrality condition

$$\begin{equation} n_i - n_e + Zn_d = 0. \end{equation} \tag{ 35 }$$

It is important to note here that the plasma gas is generally nonneutral, i.e., $n_i \not= n_e$, because the grains contribute to the overall charge neutrality of the gas–dust mixture. This is particularly true when the ionization rate is low and/or small dust grains are abundant.

3.2.1. Gas-phase Recombination

Gas-phase recombination is dissociative for molecular ions like HCO⁺. For HCO⁺, Ganguli et al. (1988) provide an empirical fit to the experimental data of the recombination rate coefficient

$$\begin{equation} K_{\rm rec} = 2.4 \times 10^{-7} \left ({\displaystyle \frac{T_e}{300 \,{\rm K}}} \right) ^{-0.69}\, {\rm cm^3\, s^{-1}}, \end{equation} \tag{ 36 }$$

where T_e is the electron temperature. In this study, we use Equation (36) but we replace T_e by 2〈_e〉/3k_B. If metal ions like Mg⁺ are dominant, gas-phase recombination is radiative, and therefore K_rec becomes much lower than that given by Equation (36). However, such a difference is unimportant when plasma accretion by dust grains dominates over gas-phase recombination.

3.2.2. Plasma Sticking to Dust Grains

The rate coefficient for plasma accretion by grains is given by

$$\begin{equation} K_{d\alpha }(E,\phi) = \int f_\alpha ({\boldsymbol E}, {\boldsymbol v}_\alpha) \sigma _{d\alpha }(\epsilon _\alpha,\phi)v_\alpha d^3v_\alpha, \end{equation} \tag{ 37 }$$

where σ_dα is the effective collision cross section. In this study, we adopt σ_dα of the form (Spitzer 1941; Shukla & Mamun 2002)

$$\begin{equation} \sigma _{d\alpha }(\epsilon _\alpha,\phi) = \left\lbrace \begin{array}{@{}ll@{}} \pi a^2 \left(1- {\displaystyle \frac{q_\alpha \phi }{\epsilon _\alpha }} \right), & \epsilon _\alpha > q_\alpha \phi, \\ [1mm] 0, & \epsilon _\alpha < q_\alpha \phi. \end{array}\right. \end{equation} \tag{ 38 }$$

It should be noted that the above expression assumes that plasma particles perfectly stick to grains upon a collision. Some ionization models in the astronomical literature (e.g., Umebayashi 1983; Nishi et al. 1991; Ilgner & Nelson 2006; Bai 2011a) assume that the electron sticking probability rapidly decreases as the electron energy increases beyond ∼100 K ∼ 10⁻² eV. However, such an assumption is inconsistent with the results of laboratory experiments. There is ample evidence that dust grains in plasmas are highly negatively charged even if the electron temperature is as high as 0.1–10 eV (e.g., Melzer et al. 1994; Barkan et al. 1994; Walch et al. 1995; Ratynskaia et al. 2004), which is well reproduced by models assuming perfect sticking (Khrapak et al. 2005). Therefore, perfect sticking is a more natural assumption as long as the secondary electron emission from dust grains is negligible (i.e., _e ≲ 100 eV). Equation (38) also neglects the polarization force between grains and charged particles, which is valid for aT ≳ 10 μm K (Draine & Sutin 1987).

For electrons, we use Equation (15) and obtain

$$\begin{eqnarray} K_{de} &=& \pi a^2 \sqrt{ {\displaystyle \frac{8 k_{\rm B}T}{\pi m_e}} } \frac{1}{W(\chi)} [(\psi +\chi)^{1+\chi }\exp (-\psi) \nonumber \\ && +(1-\psi)\Gamma (1+\chi,\psi +\chi) \exp \chi ], \end{eqnarray} \tag{ 39 }$$

where ψ ≡ −eϕ/k_BT is the negative surface potential of the grains normalized by k_BT (note that we assume ϕ < 0 and hence ψ > 0). In the absence of photoelectric and secondary electron emissions, dust grains tend to be negatively charged because electrons move much faster than ions. For neutral grains (ϕ → 0), K_de is simply given by the product of grain's geometric cross section and electron's mean speed,

$$\begin{equation} K_{de}(E,\phi =0) = \pi a^2 \langle v_e \rangle, \end{equation} \tag{ 40 }$$

with

$$\begin{eqnarray} \langle v_e \rangle &\equiv & \int |{\boldsymbol v}_e| f_{e0} d^3v_e \nonumber \\ &=& \sqrt{ {\displaystyle \frac{8 k_{\rm B}T}{\pi m_e}} } \frac{1}{W(\chi)} [ \chi ^{1+\chi } +\Gamma (1+\chi,\chi) \exp \chi ]. \end{eqnarray} \tag{ 41 }$$

Therefore, the dimensionless quantity

$$\begin{equation} {\cal C}(E,\phi) \equiv \frac{K_{de}}{\pi a^2 \langle v_e \rangle } \end{equation} \tag{ 42 }$$

measures how much the electron–grain collision rate is reduced by the Coulomb repulsion. We will call ${\cal C}$ the Coulomb reduction factor for electron–grain collisions. By using $f_{e0}^{\rm (M)}$ and $f_{e0}^{\rm (D)}$ instead of Equations (15), one can obtain the asymptotic expressions of 〈v_e〉 and ${\cal C}$ in the weak and strong field limits,

$$\begin{equation} \langle v_e \rangle \approx \left\lbrace \begin{array}{@{}ll@{}} \sqrt{ {\displaystyle \frac{8 k_{\rm B}T}{\pi m_e}} }, & E \ll E_{\rm crit}, \\ [4mm] {\displaystyle \frac{\sqrt{2}}{3^{1/4}\Gamma \left(\frac{3}{4}\right)}} \left ({\displaystyle \frac{m_n}{m_e}} \right) ^{1/4} \sqrt{ {\displaystyle \frac{eE\ell _e}{m_e}} }, & E \gg E_{\rm crit}, \end{array}\right. \end{equation} \tag{ 43 }$$

and

$$\begin{equation} {\cal C} \approx \left\lbrace \begin{array}{@{}ll@{}} \exp \left({\displaystyle \frac{-e|\phi |}{k_{\rm B}T}} \right), & E \ll E_{\rm crit}, \\ [4mm] \ \exp (-X^2) - \sqrt{\pi }X \, {\rm erfc}(X), & E \gg E_{\rm crit}, \end{array}\right. \end{equation} \tag{ 44 }$$

where erfc(x) is the complementary error function and

$$\begin{equation} X \equiv \sqrt{ {\displaystyle \frac{3m_e}{m_n}} } {\displaystyle \frac{|\phi |}{E\ell _e}}. \end{equation} \tag{ 45 }$$

The form of ${\cal C}$ in the weak field limit is well known (Spitzer 1941; Shukla & Mamun 2002). Because e|ϕ|/k_BT ≈ 1.5e|ϕ|/〈_e〉 for E ≪ E_crit and X ≈ 0.74e|ϕ|/〈_e〉 for E ≫ E_crit, Equation (44) indicates that ${\cal C}$ is determined by the ratio between the electric and kinetic energies e|ϕ|/〈_e〉. This can also be seen in Figure 4, where we plot the exact form of ${\cal C}$ as a function of E/E_crit for three cases eϕ = 0, −3k_BT, and −2〈_e〉. As we see, ${\cal C}$ is nearly constant for ϕ∝〈_e〉, while ${\cal C}$ increases toward unity for constant ϕ.

Zoom In Zoom Out Reset image size

For ions, we use Equation (25) and obtain

$$\begin{eqnarray} K_{di} &=& \pi a^2 \left [ \sqrt{ {\displaystyle \frac{2 k_{\rm B}T_i}{\pi m_i}} } \exp \left(-\frac{m_i\langle {\boldsymbol v}_{i||} \rangle ^2}{2 k_{\rm B}T_i} \right)\right. \nonumber \\ &&\left. + |\langle {\boldsymbol v}_{i||} \rangle | \left(1 + \frac{k_{\rm B}T_i + 2e|\phi | }{m_i \langle {\boldsymbol v}_{i||} \rangle ^2} \right) {\rm erf} \left(\sqrt{\frac{m_i}{2k_{\rm B}T_i}}|\langle {\boldsymbol v}_{i||} \rangle | \right) \right ],\nonumber\\ \end{eqnarray} \tag{ 46 }$$

where erf(x) is the error function. In the limits of E ≪ E_{crit, i} and E ≫ E_{crit, i}, Equation (46) reduces to

$$\begin{equation} K_{di} \approx \left\lbrace \begin{array}{@{}ll@{}} \pi a^2 \sqrt{ {\displaystyle \frac{8 k_{\rm B}T}{\pi m_i}} }\left(1 + {\displaystyle \frac{e|\phi |}{k_{\rm B}T}} \right), & E \ll E_{{\rm crit},i}, \\ [4mm] \pi a^2 |\langle {\boldsymbol v}_{i||} \rangle | \left(1 + {\displaystyle \frac{2e|\phi |}{m_i\langle {\boldsymbol v}_{i||} \rangle ^2}} \right), & E \gg E_{{\rm crit},i}, \end{array}\right. \end{equation} \tag{ 47 }$$

where we have used m_i ≫ m_n for the high-field expression.

3.2.3. Impact Ionization

The impact ionization rate coefficient K_* is given by

$$\begin{equation} K_{*}(E) = \int f_e({\boldsymbol E},{\boldsymbol v}_e) \sigma _*(\epsilon _e) v_e d^3 v_e, \end{equation} \tag{ 48 }$$

where σ_* is the impact ionization cross section of neutrals. For σ_*, we adopt Thomson's expression (Thomson 1912)

$$\begin{equation} \sigma _*(\epsilon _e) = \left\lbrace \begin{array}{@{}ll@{}} {\displaystyle \frac{\pi N_* e^4}{\epsilon _e^2}} \left({\displaystyle \frac{\epsilon _e}{{\rm IP}}} -1\right), & \epsilon _e > {\rm IP}, \\ 0, & \epsilon _e < {\rm IP}, \end{array}\right. \end{equation} \tag{ 49 }$$

where N_* is the number of bound electrons in the outermost shell of the neutrals and IP is the ionization potential of the outermost bound electrons. Equation (49) well approximates experimentally obtained ionization cross sections unless _e is much larger than IP (Lotz 1967). We only consider the impact ionization of H₂ molecules (IP = 15.4 eV) because they dominate the gas of protoplanetary disks. However, if there are a considerable number of metal atoms having a low ionization energy (e.g., K and Ca) in the gas phase, they would effectively lower the value of IP in Equation (49).

Since impact ionization is only important for E ≫ E_crit, it is sufficient to evaluate K_* using the Druyvesteyn distribution $f_{e0}^{\rm (D)}$ (Equation (20)). This allows us to analytically perform the integration in Equation (48), yielding

$$\begin{equation} K_*(E) = \frac{\sqrt{2}\pi^{3/2} N_*e^4}{\Gamma \left(\frac{3}{4}\right)\sqrt{m_e}\,{\rm IP}^{3/2}} \sqrt{Y} \left({\rm erfc}(Y) + \frac{Y}{\sqrt{\pi }}{\rm Ei}(-Y^2) \right) \end{equation} \tag{ 50 }$$

with

$$\begin{equation} Y \equiv \sqrt{\frac{3m_e}{m_n}}\frac{{\rm IP}}{eE\ell _e}, \end{equation} \tag{ 51 }$$

where ${\rm Ei}(x) \equiv -\int _{-x}^\infty t^{-1}\exp (-t) dt$ is the exponential integral. Note that K_* is determined by the ratio 〈_e〉/IP because Y ≈ 0.74IP/〈_e〉. Since the Druyvesteyn distribution neglects inelastic losses, Equation (50) must be taken as a very crude estimate for K_*.

Figure 5 shows K_* for hydrogen molecules as a function of 〈_e〉. As we can see, K_* abruptly increases before 〈_e〉 reaches the ionization threshold IP = 15.4 eV. This means that electrons at the high-energy tail of the energy distribution significantly contribute to the impact ionization. This would remain true, at least qualitatively, even if inelastic ionization losses are included.

Zoom In Zoom Out Reset image size

3.2.4. Charge Equilibrium Solution

In this study, we follow Okuzumi (2009) and calculate the equilibrium solutions to the rate equations in an analytic way. First we solve Equations (32) and (33) with respect to n_i and n_e under the equilibrium condition dn_i/dt = dn_e/dt = 0. The solution $n_\alpha \equiv n_\alpha ^{\rm (eq)}$ (α = i, e) is then a function of E and ϕ. It is easy to show that the solution is given by

$$\begin{equation} n_\alpha ^{(\rm eq)}(E, \phi) = \frac{\zeta n_n}{K_{d\alpha }(E,\phi)n_d} \left(\sqrt{\frac{1}{\cal S} + \left(\frac{1-{\cal I}}{2}\right)^2} + \frac{1-{\cal I}}{2} \right)^{-1}, \end{equation} \tag{ 52 }$$

where the dimensionless quantities ${\cal S}$ and ${\cal I}$ are defined by

$$\begin{equation} {\cal S}(E,\phi) = \frac{K_{di}(E,\phi)K_{de}(E,\phi)n_d^2}{K_{\rm rec}(E) \zeta n_n}, \end{equation} \tag{ 53 }$$

$$\begin{equation} {\cal I}(E,\phi) = \frac{K_*(E)n_n}{K_{de}(E,\phi)n_d}. \end{equation} \tag{ 54 }$$

The parameter ${\cal S}$ indicates which of gas-phase recombination and plasma sticking onto grains dominates (the latter dominates if ${\cal S} > 1$), while ${\cal I}$ indicates which of electron sticking onto grains and impact ionization dominates (the latter dominates if ${\cal I} > 1$). In this paper, we will call ${\cal S}$ the grain recombination parameter.

In order to determine ϕ as a function of E, we substitute Equation (52) into the charge neutrality condition (Equation (35)) to obtain

$$\begin{equation} n_i^{(\rm eq)}(E, \phi) - n_e^{(\rm eq)}(E, \phi) + \frac{a n_d}{e}\phi = 0, \end{equation} \tag{ 55 }$$

where we have used Equation (31) to rewrite Z as aϕ/e. Equations (52)–(55) reduce to Equations (27)–(30) of Okuzumi (2009) in the limit E → 0. As shown by Okuzumi (2009), one can apply our Equations (52)–(55) to arbitrary grain size distribution dn_d/da (number of grains per unit grain radius) if one replaces an_d in Equation (55) with ∫a(dn_d/da)da, and πa²n_d in K_dαn_d with ∫πa²(dn_d/da)da.

We solve Equation (55) with respect to ϕ using the Newton–Raphson method. If impact ionization is neglected (${\cal I}=0$), the left-hand side of Equation (55) monotonically increases with ϕ, so Equation (55) has only one root for each value of E. In this case, the Newton–Raphson procedure converges to the single root with an arbitrary initial guess. With impact ionization, Equation (55) can possess three roots for a certain range of E (see Section 5). In this case, we search for all the roots by varying the initial guess gradually from ϕ = −0.1V to −10V.

3.3. Model Parameters

Impact ionization is important only when the electric field strength E is so high that the mean electron energy exceeds a few eV. However, heating of electrons still occurs at lower E and affects the rates of gas-phase recombination and electron sticking to dust grains. We here ignore impact ionization (${\cal I} = 0$; models A–D) and focus on how plasma heating changes the ionization state at such low E. The effect of impact ionization will be studied in Section 5.

4.1. Ionization State

Figure 6 shows the ionization state of the four models as a function of E. Here we plot the abundances of plasma particles, x_e(= n_e/n_n) and x_i(= n_i/n_n), and the negative grain charge −Z in the abundance form −Zx_d (=   − Zn_d/n_n). From the charge neutrality, x_i is equal to the sum of x_e and −Zx_d. Since impact ionization is not treated here, the equilibrium ionization state is determined by the balance among external ionization, gas-phase recombination, and charging of dust grains. This balance can be characterized by two dimensionless quantities. The first one is ${\cal S}$ already introduced in Section 3.2.4 (Equation (53)). This quantity is a diagnostic of the dominant recombination process: recombination mainly takes place in the gas phase for ${\cal S} <1$, and in the "solid phase" (i.e., on dust grains) for ${\cal S} >1$. The second one is given by

$$\begin{equation} {\cal P} \equiv \frac{|Z|n_d}{n_e}, \end{equation} \tag{ 56 }$$

which is known as the Havnes parameter (Havnes 1984) in the field of dusty plasma physics. This is a diagnostic of the charge neutrality in the gas–dust mixture. If ${\cal P} \ll 1$, gas-phase free electrons dominate over negatively charged grains, and the charge neutrality is approximately established within the gas phase (i.e., n_i ≈ n_e). If ${\cal P} >1$, negatively charged grains dominate, and the number of positive ions in the gas approximately balances with the number of free electrons on grains (n_i ≈ −Zn_d). The two dimensionless quantities measure how strongly dust grains affect the ionization state.

Zoom In Zoom Out Reset image size

Figure 7 plots these two diagnostics for the four models as a function of E. By definition, ${\cal S}$ and ${\cal P}$ increase as the amount of dust f_dg is increased (see models A, B, and C). They also increase as the ionization rate ζ is decreased (see models C and D) because the presence of grains becomes more and more important as the ionized degree decreases. For fixed f_dg and ζ, ${\cal S}$ and ${\cal P}$ increase with E, because more and more electrons are transferred from the gas to grains as the random velocity of electrons (≈collision velocity between electrons and grains) is increased.

Zoom In Zoom Out Reset image size

Model A is characterized by the low dust-to-gas ratio f_dg = 10⁻⁶ and is an example where ${\cal S} < 1$ and ${\cal P} < 1$ over the entire range of E under consideration. In this model, the gas phase dominates both recombination and charge neutrality, and dust grains have essentially no effect on the ionization state of the gas. The balance between external ionization rate ζn_n and gas-phase recombination rate K_recn_in_e gives an approximate expression for the plasma density n_e (≈ n_i) in the case of ${\cal S} \ll 1$,

$$\begin{equation} n_e \approx \sqrt{\frac{\zeta n_n}{K_{\rm rec}(E)}}. \end{equation} \tag{ 57 }$$

Note that n_e increases with E because the gas-phase recombination rate coefficient K_rec is a decreasing function of 〈_e〉. Since n_e ≈ n_i, the electron flux n_e〈v_e〉 is much higher than the ion flux n_i〈v_i〉 (note that in general 〈v_e〉 ≫ 〈v_i〉), so individual dust grains tend to be charged up so that Coulomb repulsion between the grains and electrons becomes effective. This can be seen in Figure 8, where we plot |ϕ| as a function of E. In model A, |ϕ| increases with E in the way that the relation e|ϕ| ∼ 2〈_e〉 is approximately satisfied, i.e., in the way that the Coulomb repulsion energy between a grain and an electron upon collision is comparable to their collision energy. As a result, the Coulomb reduction factor ${\cal C}$ is much less than unity (${\cal C} <0.1$) for all E as shown in Figure 9.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Model B is a more dusty case where f_dg is 100 times larger than in model A. As a result, ${\cal S}$ now exceeds unity (i.e., solid-phase recombination becomes the dominant recombination process) at E ≳ 10E_crit ≈ 10⁻⁸ esu cm⁻². In the case of ${\cal S} \gg 1$, n_e is determined by the balance between external ionization rate ζn_n and electron capture rate K_recn_dn_e, i.e.,

$$\begin{equation} n_e \approx \frac{\zeta n_n}{K_{de}n_d} = \frac{\zeta n_n}{\pi a^2 n_d \langle v_e \rangle {\cal C}(E,\phi)}. \end{equation} \tag{ 58 }$$

Note that n_e decreases with increasing E because both 〈v_e〉 and ${\cal C}$ increase with E (for ${\cal C}$, see Figure 9). If ${\cal P} < 1$, as is the case in model B, n_i is also given by Equation (58).

Model C is an even more dusty case where the dust-to-gas ratio is interstellar. In this model, ${\cal S} > 1$ over the entire range of E. In addition, at high E, ${\cal P}$ exceeds unity, i.e., dust grains become the dominant negative charge carriers. We can see that n_e rapidly decreases when ${\cal P}$ crosses unity. This is a positive feedback effect of grain's negative charging on electron depletion. As ${\cal P}$ exceeds 1, n_e becomes smaller than n_i, and the electron-to-ion flux ratio n_e〈v_e〉/n_i〈v_i〉 becomes closer to unity. For this reason, individual dust grains tend to be less negatively charged than in the case of ${\cal P} \ll 1$. This can be seen in Figure 8, where we see that e|ϕ| falls below 〈_e〉 after ${\cal P}$ exceeds unity. As a result, the Coulomb repulsion between the grains and electrons become ineffective (${\cal C} \approx 1$), leading to a further decrease in the electron number density according to Equation (58). Note that n_i is constant at ${\cal P} > 1$ as long as ion heating is insignificant (i.e., E ≪ E_{crit, i}). This constant value is given by n_i ≈ ζn_n/πa²v_{i, T}, where $v_{i,T} = \sqrt{8 k_{\rm B}T/\pi m_i}$ is the mean thermal speed of ions at T_i = T.

In model D, ζ is decreased by factor 100 from model C, and we see that ${\cal S} > 1$ and ${\cal P} > 1$ over the entire range of E. The electron abundance decreases at E > E_crit, but more slowly than in model C since the Coulomb repulsion factor ${\cal C}$ is already close to unity from the beginning (see Figure 9).

4.2. Current Density

Figure 10 show the magnitude of the current density $J = |{\boldsymbol J}|$ as a function of E for models A–D. This is the sum of the ion and electron currents, $J_i = |{\boldsymbol J}_i| = e n_i |\langle \boldsymbol v_{i||} \rangle |$ and $J_e = |{\boldsymbol J}_e| = e n_e |\langle \boldsymbol v_{e ||} \rangle |$, which are also plotted in Figure 10.

Zoom In Zoom Out Reset image size

As mentioned earlier, the J–E relations are nonlinear in E because the plasma heating changes n_α and $|\langle \boldsymbol v_{\alpha ||} \rangle |$ (α = i, e) at high E. The nonlinearity is, however, weak in model A. In this model, the electric current is dominated by J_e, which is proportional to the product of n_e and $|\langle {\boldsymbol v}_{e||} \rangle |$. At E ≫ E_crit, $|\langle {\boldsymbol v}_{e||} \rangle |$ increases more slowly than at E ≪ E_crit due to the enhanced frequency of electron–neutral collisions (see Equation (23)). Meanwhile, n_e increases with E because the recombination rate coefficient K_rec decreases with ${\cal \epsilon }_e$ (see Equations (36) and (57)). These two opposing effects partially cancel out in the product $n_e|\langle {\boldsymbol v}_{e||} \rangle |$.

In model B, J behaves in the same way as in model A as long as ${\cal S}< 1$. The behavior changes when ${\cal S}$ crosses unity because n_e becomes a decreasing function of E (see Section 4.1). We see that the current is approximately constant (precisely speaking, decreases very slowly with E) at ${\cal S} > 1$. This trend can be explained as follows. When ${\cal S} \gg 1$, n_e is inversely proportional to the mean electron speed 〈v_e〉 (see Equation (58)), which is, at E ≫ E_crit, proportional to the electron drift speed $\langle {\boldsymbol v}_{e||} \rangle$. Equations (23) and (43) imply that the ratio of the two velocities is

$$\begin{equation} |\langle {\boldsymbol v}_{e||} \rangle | \approx \sqrt{\frac{\pi m_e}{3m_n}} \langle v_{e} \rangle. \end{equation} \tag{ 59 }$$

Therefore, in $J_e = {{en}_e} |\langle {\boldsymbol v}_{e||} \rangle |$, the dependences on 〈v_e〉 cancel out, resulting in

$$\begin{equation} J_e \approx \frac{J_{e,\infty }}{{\cal C}}, \end{equation} \tag{ 60 }$$

where J_{e, ∞} is a constant defined by

$$\begin{eqnarray} J_{e,\infty } &\equiv & \sqrt{\frac{\pi m_e}{3m_n}} \frac{\zeta e n_n}{\pi a^2 n_d} \nonumber \\ &\approx & 5 \times 10^{-5} \left ({\displaystyle \frac{10^{-4}}{f_{dg}}} \right) \left ({\displaystyle \frac{a}{1\, {\mu {\rm m}}}} \right) \left ({\displaystyle \frac{\zeta }{10^{-17} {\rm s^{-1}}}} \right) \,{\rm esu\, cm^{-2} s^{-1}}.\nonumber\\ && \end{eqnarray} \tag{ 61 }$$

Hence, if ${\cal C}$ is independent of E, so is J_e. As we have already seen in Section 4.1, ${\cal C}$ varies only slowly with E unless ${\cal P}$ crosses unity. This explains why in model B the current is approximately constant at ${\cal S}>1$.

The result for model C is more complex, but can also be explained in a similar way. We see that J_e starts to decrease with E at the point where ${\cal P}$ crosses unity. This is because the Coulomb reduction factor ${\cal C}$ increases from 0.025 to unity as ${\cal P}$ goes from ≪1 to ≫1. This leads to a 40-fold decrease in J_e across ${\cal P} = 1$ as predicted by Equation (60). By contrast, the ion current J_i continues increasing with E because n_i is approximately constant for ${\cal S}>1$. As a consequence, the net current J = J_e + J_i forms an N-shaped curve in the J–E diagram.

In model D, J_e immediately approaches J_{e, ∞} at E > E_crit since ${\cal C} \approx 1$ from the beginning. The net current monotonically increases with E because of the presence of J_i. At E ≳ E_{crit, i}, J_i also gets saturated at a constant value.

We now include impact ionization (models B* and C*) and see how it changes the ionization balance at very high electric field strengths.

Figure 11 shows the plasma abundances x_i and x_e and current density J as a function of E for models B* and C*. As expected by previous studies (Inutsuka & Sano 2005; Muranushi et al. 2012), impact ionization dramatically changes the ionization state at large E. In both models, we observe an "electric discharge," an abrupt increase in the plasma abundance, when the mean electron energy 〈_e〉 reaches ≈3 eV. This is due to the rapid increase in the impact ionization rate coefficient K_* around that electron energy (see Figure 5). It is interesting to see that the discharge current appears much earlier than 〈_e〉 reaches the ionization potential IP = 15.4 eV. This means that the high-energy tail of the energy distribution function is responsible for this ionization. At lower E, impact ionization has no effect on the ionization state, so the curves' left ends in models B* and C* are identical to those in models B and C, respectively.

Zoom In Zoom Out Reset image size

However, the nature of the discharge current is much more complex than assumed in the previous studies. In model C*, we find that Equation (55) has three equilibrium solutions for a single value of E when the mean electron energy falls within the narrow range 3.0 eV ≲ 〈_e〉 ≲ 3.3 eV. The triple solution forms an S-curve in the J–E space as shown in the small panel of Figure 11. By contrast, in model B*, the discharge current is single-valued for all E. In the following subsections, we analyze the structure of the equilibrium solutions in more detail.

5.1. Classification of Equilibria

First let us see how the reaction balance is changed by impact ionization. In Figure 12, we plot the rates (per unit volume) of external ionization (ζn_n), impact ionization (K_*n_nn_e), gas-phase recombination (K_recn_in_e), and electron sticking onto grains (K_den_dn_e) as a function of E. In the limit of low E, external ionization balances with plasma sticking onto grains (i.e., ζn_n ≈ K_den_dn_e), and impact ionization is negligible as well as gas-phase recombination. In the opposite limit, impact ionization balances with gas-phase recombination (i.e., K_*n_nn_e ≈ K_recn_in_e), and external ionization and plasma sticking are subdominant (in general, gas-phase recombination becomes more and more important as the ionization degree is increased). In the following, we will refer to the former ionization state as the low state, and to the latter as the high state. In model B*, the low state is smoothly connected to the high state at E ≈ 300E_crit.

Zoom In Zoom Out Reset image size

In model C*, we can see the third type of ionization state. Of the three equilibrium solutions appearing at 340 ≲ E/E_crit ≲ 360, the top and bottom solutions are merely an extension of the low and high states, respectively, but the middle solution is characterized by the balance between impact ionization and plasma capture by dust grains (i.e., K_*n_nn_e ≈ K_den_dn_e). We will call this the middle state.

5.2. Emergence of Multiple Equilibria

It is easy to explain why no multiple equilibrium solution appears in model B*. We first note that the grain surface potential ϕ in equilibrium must satisfy the relation

$$\begin{equation} K_{di}(E,\phi)n_i = K_{de}(E,\phi)n_e \end{equation} \tag{ 62 }$$

(see Equation (34)). In model B*, gas-phase electrons are so abundant that n_i ≈ n_e (or equivalently, ${\cal P} \ll 1$) even before the onset of impact ionization. For fixed E, K_di is a decreasing function of ϕ, while K_de is an increasing function of ϕ. Therefore, in the case of n_i ≈ n_e, Equation (62) uniquely specifies ϕ, and in turn n_e (≈n_i), for each value of E. In fact, if we look at ϕ of model B*, there is no appreciable change in ϕ before and after the onset of impact ionization (see the left panel of Figure 13). This is the reason why no multiple solution emerges in model B*. By contrast, in model C*, the condition n_i ≈ n_e is violated (or equivalently, ${\cal P} \gg 1$) before the onset of impact ionization. In this case, Equation (62) can have multiple solutions, because the ratio n_e/n_i can vary with ϕ and hence the monotonicity of Equation (62) as a function of ϕ is not ensured. This argument suggests that multiple solutions emerge only when the condition ${\cal P} > 1$ is satisfied before the onset of electric discharge.

Zoom In Zoom Out Reset image size

5.2.1. Instability of the Intermediate (Middle) State

Emergence of a triple equilibrium or an S-shaped equilibrium curve can be seen in many physical systems. In most cases, two extreme equilibria are stable, while the middle equilibrium is unstable against perturbation. We find that this is also the case for our multiple solutions. We numerically solved the time-dependent rate equations (Equations (32)–(34)) for various initial conditions and looked at which of the equilibria is reached at late times. An example of such tests is shown in Figure 14. Here we plot the time evolution of the electron abundance x_e for different initial conditions in the case of model C* with E = 350E_crit. The low, middle, and high equilibrium states correspond to x_e = 1 × 10⁻¹⁵, 2 × 10⁻¹², and 2 × 10⁻⁷, respectively. As we can see, all the time-dependent solutions converge toward either the high or low state, while the middle state is never reached even if the initial state is very close to it. For completeness, in Appendix B, we perform a linear stability analysis of the middle state by using simplified rate equations, and show that the middle solution is indeed unstable. From these facts, we conclude that the middle equilibrium of the discharge current is never realized in real systems.

Zoom In Zoom Out Reset image size

In this section, we discuss important implications of our model calculations for the MHD of protoplanetary disks.

6.1. Negative Differential Resistance and Its Instability

Negative differential resistance refers to the property of some electric circuits (e.g., Gunn diodes) that an increase in the applied voltage causes a decrease in the electric current. Interestingly, some of our model calculations yield J–E relations that have negative differential resistance (i.e., ${dJ}/{dE}<0$) in some range of E. For example, in model C, we see that J decreases by an order of magnitude when going from E ≈ 3E_crit to E ≈ 30E_crit. As we will discuss below, negative differential resistance has many important implications for the evolution of electric fields and for the MHD of protoplanetary disks.

The most important consequence of negative differential resistance is that the displacement current neglected in Ampère's law ceases to be negligible. To see this, let us consider the Maxwell–Ampère equation

$$\begin{equation} \frac{\partial {\boldsymbol E}}{\partial t} = c\nabla \times {\boldsymbol B} - 4\pi {\boldsymbol J} \end{equation} \tag{ 63 }$$

with ${\boldsymbol J} = J(E) \hat{\boldsymbol E}$, where J is a nonlinear function of the electric field strength $E = |{\boldsymbol E}|$. Note again that all the quantities are defined in the comoving frame of neutrals. We assume that the background electric field ${\boldsymbol E}_0$ is imposed by external sources and is approximately steady over a long timescale. Then, the background magnetic field ${\boldsymbol B}_0$ is related to ${\boldsymbol E}_0$ by the classical Ampère's law

$$\begin{equation} c\nabla \times {\boldsymbol B}_0 = 4\pi {\boldsymbol J}({\boldsymbol E}_0). \end{equation} \tag{ 64 }$$

We examine the stability of this relation by considering a perturbation ${\boldsymbol E}(t) = {\boldsymbol E}_0 + {\boldsymbol E}_1(t)$, where $|{\boldsymbol E}_1| \ll |{\boldsymbol E}_0|$. For simplicity, we drop the perturbation of the $c \nabla \times {\boldsymbol B}$ term by assuming that the wavelength of the perturbed field is sufficiently long.⁶ Substituting this into Equation (63) we obtain

$$\begin{equation} \frac{\partial {\boldsymbol E}_1}{\partial t} = - 4\pi \sigma _{\rm diff}{\boldsymbol E}_1, \end{equation} \tag{ 65 }$$

where

$$\begin{equation} \sigma _{\rm diff} \equiv \frac{{dJ}}{{dE}} (|{\boldsymbol E}_0|) \end{equation} \tag{ 66 }$$

is the differential conductivity evaluated at $E = |{\boldsymbol E}_0|$. To order of magnitude, |σ_diff|⁻¹ ∼ E/J, which is ∼10⁻⁴–10⁻²s for models B* and C* in the region of negative σ_diff. If σ_diff is positive, as is the case for the conventional linear Ohm's law, the perturbation in E decays on a timescale of (4πσ_diff)⁻¹ (which is known as the Faraday time). However, if σ_diff is negative, the perturbation grows exponentially with time, meaning that the Ampère's law (Equation (64)) is unstable.

Figure 15 summarizes the instabilities of the nonlinear Ohm's law identified in this study. As already seen in Section 6.1, the middle branch of the discharge current has, whenever present, an unstable ionization balance. In Figure 15, the unstable middle branch is indicated by the dotted line, and the instability by the vertically diverging arrows. In addition to this, we have found here that a current with a negative ${dJ}/{dE}$ is unstable to perturbations in E. A negative ${dJ}/{dE}$ can appear at E ≳ E_crit and at the discharge current, which are indicated by the dashed and dotted line segments in Figure 15, respectively. The associated instability is represented by the horizontally diverging arrows.

Zoom In Zoom Out Reset image size

The instability of negative differential resistance has a significant impact on the MHD of the system. The fundamental assumption of the standard non-relativistic MHD is that the displacement current is negligible (i.e., Ampère's law approximately holds) on a dynamical timescale. As shown above, this assumption is always valid for linear Ohm's laws, but not for nonlinear Ohm's laws that exhibit negative differential resistance over some range of E. In the latter case, one must in principle treat the dynamics of the system by using the fully time-independent Maxwell–Ampére's law (Equation (63)) instead of the quasi-steady Ampére's law. Of course, such a task is computationally challenging since one then needs to treat unwanted electromagnetic modes that appear at the same time.

It should be noted, however, that the linear analysis presented above does not directly apply to real protoplanetary disks. In a typical disk environment, the predicted timescale of the instability (≪1s) is much shorter than the relaxation timescales of the velocity distributions and charge reactions of the plasmas. This invalidates the use of the relation J = J(E), which assumes that both n_α and $\langle {\boldsymbol v}_{||\alpha } \rangle$ (α = i, e) instantaneously reach an equilibrium for given E. To go further, we need to consider non-equilibrium evolution of J instead of using the quasi-equilibrium relation J = J(E). Such a task is beyond the scope of this paper, but will be addressed in our future work.

Nevertheless, we may argue that inclusion of the displacement current is essential to treat the dynamics of the system with negative differential resistance. When we neglect the displacement current, we lose the ability to treat the electric field ${\boldsymbol E}$ as an independent dynamical quantity since the displacement current is responsible for the time evolution of ${\boldsymbol E}$. This is not an issue when ${\boldsymbol J}$ is monotonic in ${\boldsymbol E}$, because ${\boldsymbol E}$ is then uniquely determined as a function of ${\boldsymbol J}$. However, in the presence of negative differential resistance, ${\boldsymbol E}$ becomes multivalued as a function of ${\boldsymbol J}$, and therefore cannot be determined by instantaneous relations. In this case, the state of the system depends on the history of ${\boldsymbol E}$, and this can only be determined by Maxwell–Ampére's equation with the $\partial {\boldsymbol E}/\partial t$ term. In a forthcoming paper, we will show that the displacement current naturally solves this issue by allowing hysteresis for the relation between ${\boldsymbol J}$ and ${\boldsymbol E}$.

6.2. Implications for MRI Turbulence

One important finding of this study is that plasma heating reduces the electric conductivity J/E before the onset of impact ionization. This was not considered in our previous studies (Inutsuka & Sano 2005; Muranushi et al. 2012), which only assumed that plasma heating enhances the conductivity via impact ionization. The reduction of the conductivity might lead to self-regulation of magnetohydrodynamic motion in protoplanetary disks: coupling between a moving gas and a magnetic field generates an electric field in the comoving frame, but this causes a reduction of the conductivity and hence of the coupling between the gas and magnetic field.

Such an effect is of particular importance to MRI turbulence in protoplanetary disks as it could limit or even determine the saturation level of the turbulence. However, in order to prove this, we would have to perform a resistive MHD simulation including the effect of plasma heating on the resistivity, which is clearly beyond the scope of this paper. Below, we shall only speculate, by using two illustrative examples, how MRI turbulence will develop in a protoplanetary disk under the effect of plasma heating.

As mentioned in Section 2, fully saturated MRI turbulence is characterized by the universal average current density J_MRI (Equation (3)). Therefore, MRI tends to grow until the current density J reaches this saturation value. However, MRI does not grow but decays when the Elsasser number Λ (Equation (5)) is less than the critical value Λ_crit ∼ 0.1–1. Taken together, we may assume that MRI grows until either J reaches J_MRI or Λ falls below Λ_crit.

Let us map these criteria onto the J–E diagrams presented in this study. From Equations (3), J_MRI can be evaluated as

$$\begin{eqnarray} J_{\rm MRI} &\approx& 2\times 10^{-3} \left ({\displaystyle \frac{f_{\rm sat}}{10}} \right) \left ({\displaystyle \frac{n_n}{10^{12}\, {\rm cm^{-3}}}} \right) ^{1/2}\nonumber\\ &&\times \left ({\displaystyle \frac{30 \,{\rm yr}}{t_K}} \right)\, {\rm esu\, cm^{-2} s^{-1}}, \end{eqnarray} \tag{ 67 }$$

where t_K = 2π/Ω is the local orbital period. The MRI stability criterion Λ < Λ_crit can be written as the condition for the conductivity J/E (see Equation (6)),

$$\begin{eqnarray} \frac{J}{E} & < & \frac{\beta _z\Lambda _{\rm crit}\Omega }{8\pi } \left ({\displaystyle \frac{c}{c_s}} \right) ^2 \nonumber \\ &\approx & 2\times 10^{4}\Lambda _{\rm crit} \left ({\displaystyle \frac{\beta _z}{100}} \right) \left ({\displaystyle \frac{100 \,{\rm K}}{T}} \right) \left ({\displaystyle \frac{30 \,{\rm yr}}{t_K}} \right) {\rm s^{-1}}. \end{eqnarray} \tag{ 68 }$$

The orbital period of t_K = 30 yr corresponds to the distance of ≈10 AU from a solar-mass star. Below we assume f_sat = 10 and Λ_crit = 1.

In Figure 16, we plot the equilibrium J–E relations for models B* and C* together with the saturation criterion J = J_MRI and stability criterion Λ < 1. Here we assume t_K = 30 yr, so that the parameter set (n_n, T, t_K) approximately corresponds to the orbital distance of 10 AU in the minimum-mass solar nebula. The value of β_z is taken to be 100, which corresponds to strong turbulence induced by a large net poloidal flux (see, e.g., Hawley et al. 1995). For comparison, the conventional linear Ohm's law is also plotted. In both models, the instability criterion Λ > 1 is satisfied at E < E_crit(≈ 10⁻⁹ esu cm⁻²), so MRI is active at least in its early growth stages. However, at E > E_crit, J starts to decrease and cross the Λ = 1 line before it reaches J_MRI. Since MRI grows at Λ > 1 and decays at Λ < 1, we expect that the turbulence will saturate on the Λ = 1 line. In particular, the saturation level in model B* is expected to be much lower than that of fully developed MRI turbulence because the value of J at Λ = 1 is two orders of magnitude smaller than J_MRI (i.e., the value of f_sat introduced in Section 2 is as small as 0.1). Of course, a quantitative estimate of the saturation level requires MHD simulations.

Zoom In Zoom Out Reset image size

MRI generates strong electric fields, which can significantly heat up plasmas in weakly ionized protoplanetary disks. To study how this affects the ionization state and MHD of the disks, we have formulated a charge reaction model that takes into account plasma heating and impact ionization by heated electrons as well as plasma accretion onto dust grains. The output of our model is the electric current density J as a function of the electric field strength E in the comoving frame of the neutral gas. Because the plasma heating changes the ionization degree of the gas, the resulting Ohm's law is nonlinear in E.

We have presented some model calculations to illustrate the effects of plasma heating on the ionization balance of a dusty gas. The key findings are summarized as follows.

• 1.  
  When impact ionization is negligible, the ionization states are characterized by (1) which of gas-phase recombination and plasma sticking to grains ("solid-phase recombination") dominates the reaction balance, and by (2) which of gas-phase electrons and charged grains are the dominant negative charge carriers. These two conditions are quantified by the dimensionless grain recombination parameter ${\cal S}$ (Equation (53)) and Havnes parameter ${\cal P}$ (Equation (56)), respectively. For both conditions, the presence of dust becomes more and more important (both ${\cal S}$ and ${\cal P}$ increase) as the number of small dust grains increases, the external ionization rate ζ decreases, and/or the electric field strength E increases (Figure 7). The field strength is relevant because electrons hit and stick to dust grains more and more frequently as they are heated up.
• 2.  
  When plasma accretion by dust grains dominates over gas-phase recombination (${\cal S}>1$), the electron abundance decreases with increasing E (Section 4.1, Figure 6) because of the electron–grain collisions facilitated by the electron heating. The current density J also decreases until the electron current is taken over by the ion current (Section 4.2 and Figure 10). In particular, J rapidly decreases when charged grains replace free electrons as the dominant negative charge carriers of the system (i.e., when ${\cal P}$ crosses unity). These results have very important implications for the MHD of the system. First, the decrease of the electron abundance implies that MRI turbulence can be self-regulating: as MRI grows, the magnetic resistivity increases, which prevents further growth of MRI (Section 6.2). Furthermore, the N-shaped J–E curve violates the fundamental assumption of the standard non-relativistic MHD that a single value of J corresponds to a single value of the comoving field strength E. In fact, our simple linear analysis suggests that the negative differential resistance $({dJ}/{dE} < 0)$ should destabilize the electric field via the displacement current, which implies that the dynamical evolution of the system should depend on the history of the electric field (Section 6.1).
• 3.  
  Impact ionization by hot electrons sets in when the mean electron energy exceeds a few eV. This results in an abrupt increase in the electric current as previously investigated by Inutsuka & Sano (2005) and Muranushi et al. (2012) (Section 5). We find that this discharge current is triple-valued as a function of E (i.e., the J–E curve is S-shaped) when charged dust grains dominate the charge neutrality (${\cal P} > 1$) at low E. Furthermore, the middle branch of the S-shape current is found to be unstable to perturbations to the ionization balance. Therefore, the MHD near the discharge current could be more complex than self-sustained turbulence as envisaged by Inutsuka & Sano (2005) and Muranushi et al. (2012).

Plasma heating could also have a significant influence on the collisional growth of dust grains. Since grains in a plasma have a nonzero (and negative) mean charge, their collisional cross section is on average smaller than their geometric cross section. This "charge barrier" can slow down the growth of small dust grains even in weakly ionized protoplanetary disks (Okuzumi 2009; Okuzumi et al. 2011; Matthews et al. 2012). The maximum value of the grain negative surface potential −ϕ is given by the energy balance −eϕ ∼ 〈_e〉 (see also our Figure 8). The maximum negative potential is ∼−10 mV in a cool gas of T ∼ 100 K, but can exceed ∼−1 V when the plasma is heated by a strong electric field. Because the Coulomb repulsion energy between two grains scales as ϕ², plasma heating can lead to a ∼10⁴ fold enhancement of the repulsion energy. Therefore, in a disk region where plasma heating is effective, the growth of dust grains could be more strongly suppressed than previously thought.

A major limitation of this study is that the velocity distribution functions adopted here neglect the effects of magnetic fields on the kinetics of the plasmas. In terms of non-ideal MHD, we have only considered Ohmic diffusivity and neglected ambipolar diffusion and Hall drift. Such a treatment is only valid in dense gases where the frequency of plasma–neutral collisions is much higher than the gyration frequency of the plasmas (e.g., Nakano & Umebayashi 1986; Wardle 1999). If the neutral drag is weak, plasma particles undergo gyromotion, which prevents plasma heating when the electric field is perpendicular to the magnetic field (Golant et al. 1980). This effect is non-negligible over a wide region of protoplanetary disks where Hall drift or ambipolar diffusion dominates over Ohmic diffusion (Balbus & Terquem 2001; Wardle 2007; Bai 2011a). Our future modeling will take into account this effect.

We thank Takayuki Muranushi, Xuening Bai, Takeru Suzuki, Shigenobu Hirose, Hidekazu Tanaka, Naoki Watanabe, Shota Nunomura, Takayuki Muto, and Shoji Mori for useful comments and inspiring discussions. We also thank the referee, Neal Turner, for his prompt and detailed report that greatly improved our presentation. This work is supported by Grant-in-Aid for Research Activity Start-up (#25887023) from JSPS, and by Grants-in-Aid for Scientific Research (#23103005 and 23244027) from MEXT.

In this Appendix, we briefly review the kinetics of weakly ionized plasmas under an applied electric field. A more comprehensive review can be found in Wannier (1953) and in Golant et al. (1980).

We consider the motion of charged particles in a neutral gas in the presence of an applied E-field. We assume that the number density n_α of the charged particles are so low that collisions between the charged particles are rare. In this case, the charged particles gain and lose their momentum $m_\alpha {\boldsymbol v}_\alpha$ and kinetic energy _α through collision with neutrals and acceleration by the electric fields. The equations of momentum and energy balance are given by (e.g., Hershey 1939)

$$\begin{equation} m_\alpha \left\langle \left\langle \frac{\Delta {\boldsymbol v}_\alpha }{\Delta t_\alpha } \right\rangle \right\rangle + q_\alpha {\boldsymbol E} = 0, \end{equation} \tag{ A1 }$$

$$\begin{equation} \left\langle \left\langle \frac{\Delta \epsilon _\alpha }{\Delta t_\alpha } \right\rangle \right\rangle + q_\alpha {\boldsymbol E}\cdot \langle {\boldsymbol v}_{\alpha ||} \rangle = 0. \end{equation} \tag{ A2 }$$

Here, $\Delta {\boldsymbol v}_\alpha$ and Δ_α are the change in the changes in ${\boldsymbol v}_\alpha$ and _α upon individual collisions averaged over the scattering angle, respectively, Δt_α is the mean free time of the charged particles, and the double brackets 〈〈⋅⋅⋅〉〉 stand for the average over the charged and neutral particle velocities. Note that Δt_α is generally a function of the relative speed $|{\boldsymbol v}_\alpha - {\boldsymbol v}_n|$. If the collision with neutrals is elastic and isotropic, $\Delta {\boldsymbol v}_\alpha$ and Δ_α can be written as (e.g., Golant et al. 1980)

$$\begin{equation} \Delta {\boldsymbol v}_\alpha = - \lambda _{\alpha n}({\boldsymbol v}_\alpha - {\boldsymbol v}_n), \end{equation} \tag{ A3 }$$

$$\begin{equation} \Delta \epsilon _\alpha = -\kappa _{\alpha n}\left[\epsilon _\alpha -\epsilon _n -\frac{1}{2}(m_\alpha -m_n){\boldsymbol v}_\alpha \cdot {\boldsymbol v}_n \right], \end{equation} \tag{ A4 }$$

where

$$\begin{equation} \lambda _{\alpha n} \equiv \frac{m_n}{m_\alpha +m_n}, \end{equation} \tag{ A5 }$$

$$\begin{equation} \kappa _{\alpha n} \equiv \frac{2m_\alpha m_n}{(m_\alpha +m_n)^2} \end{equation} \tag{ A6 }$$

are the momentum and energy transfer efficiencies for the elastic collisions, respectively. Note that λ_αn ≪ 1 for heavy charged particles (m_α ≫ m_n), while λ_αn ≈ 1 for intermediate-mass and light charged particles (m_α ≲ m_n). By contrast, κ_αn ≪ 1 for both light and heavy charged particles, and κ_αn ∼ 1 only for m_α ∼ m_n. For electrons (n_e ≪ m_n), a single collision perfectly isotropizes the velocity distribution of the electrons but hardly affects their energy distribution. For heavy ions (m_i ≫ m_n), a single collision with a neutral hardly changes the momentum and energy of the ions.

It is known that the mean free time of ions is approximately constant (i.e., independent of $|{\boldsymbol v}_\alpha - {\boldsymbol v}_n|$) because of the polarization force acting between ions and neutrals (see, e.g., Wannier 1953). In this case, Equations (A1) and (A2) are closed with respect to $\langle {\boldsymbol v}_{\alpha ||} \rangle$ and 〈_α〉, and we obtain

$$\begin{equation} {-}m_\alpha \lambda _{\alpha n} \langle {\boldsymbol v}_{\alpha ||} \rangle + q_\alpha {\boldsymbol E} \Delta t_\alpha = 0, \end{equation} \tag{ A7 }$$

$$\begin{equation} \kappa _{\alpha n} \left(\frac{3k_{\rm B}T}{2} - \langle \epsilon _\alpha \rangle \right) + q_\alpha {\boldsymbol E}\cdot \langle {\boldsymbol v}_{\alpha ||} \rangle \Delta t_\alpha = 0, \end{equation} \tag{ A8 }$$

where we have used that $\langle {\boldsymbol v}_n \rangle = 0$ and 〈_n〉 = 3k_BT/2. Solving these equations for $\langle {\boldsymbol v}_{\alpha ||} \rangle$ and 〈_α〉, we obtain

$$\begin{equation} \langle {\boldsymbol v}_{\alpha ||} \rangle = {\displaystyle \frac{q_\alpha {\boldsymbol E}\Delta t_\alpha }{m_\alpha \lambda _{\alpha n} }} = {\displaystyle \frac{m_\alpha +m_n}{m_\alpha m_n}} q_\alpha {\boldsymbol E}\Delta t_\alpha, \end{equation} \tag{ A9 }$$

$$\begin{eqnarray} \langle \epsilon _\alpha \rangle &=& {\displaystyle \frac{3}{2}} k_{\rm B}T + {\displaystyle \frac{(q_\alpha E \Delta t_\alpha)^2}{m_\alpha \kappa _{\alpha n}\lambda _{\alpha n}}} \nonumber \\ &=& {\displaystyle \frac{3}{2}} k_{\rm B}T + {\displaystyle \frac{(m_\alpha +m_n)^3}{2(m_\alpha m_n)^2}} (q_\alpha E \Delta t_\alpha)^2. \end{eqnarray} \tag{ A10 }$$

In Equation (A10), the first term comes from the collisional heating with neutrals, while the second term from the field heating. Equation (A10) implies that the field heating dominates when

$$\begin{equation} E > \frac{m_\alpha m_n \sqrt{3k_{\rm B}T}}{(m_\alpha +m_n)^{3/2}|q_\alpha |\Delta t_\alpha } \equiv E_{{\rm crit},\alpha }. \end{equation} \tag{ A11 }$$

If we eliminate q_αEΔt_α in Equation (A10) by using Equation (A9), we obtain an interesting relation

$$\begin{equation} \langle \epsilon _\alpha \rangle = \frac{3}{2}k_{\rm B}T + \frac{1}{2}m_n \langle {\boldsymbol v}_{\alpha ||} \rangle ^2 + \frac{1}{2}m_\alpha \langle {\boldsymbol v}_{\alpha ||} \rangle ^2. \end{equation} \tag{ A12 }$$

Since the third term in Equation (A12) is the energy associated with the mean drift motion, the sum of the first and second terms represents the energy of random motion. It follows that, in the limit of high E (or low T), the random energy is exactly m_n/m_α times the drift energy, and the ratio of the root-mean-squared speed to the drift speed approaches

$$\begin{equation} \frac{\left\langle v_\alpha ^2 \right\rangle ^{1/2}}{|\langle {\boldsymbol v}_{\alpha ||} \rangle |} \rightarrow \sqrt{1+\frac{m_n}{m_\alpha }} \quad (E \gg E_{\rm crit,\alpha }). \end{equation} \tag{ A13 }$$

If Δt_α depends on $|{\boldsymbol v}_\alpha -{\boldsymbol v}_n|$, Equations (A1) and (A2) are not closed with respect to $\langle {\boldsymbol v}_{\alpha ||} \rangle$ and 〈_α〉. For example, electrons have a mean free time that is nearly inversely proportional to $|{\boldsymbol v}_e-{\boldsymbol v}_n|$, because its momentum-transfer cross section (or equivalently, its mean free path $\ell _e = |{\boldsymbol v}_e-{\boldsymbol v}_n|\Delta t_e$) depends on $|{\boldsymbol v}_e-{\boldsymbol v}_n|$ very weakly (e.g., Frost & Phelps 1962; Yoon et al. 2008). Nevertheless, one can use Equations (A7) and (A8) (and hence Equations (A9) and (A10)) for electrons to a good approximation if one approximates the electron mean free time Δt_e as $\ell _e/\sqrt{\langle v_e^2 \rangle } = \ell _e/\sqrt{2\langle \epsilon _e \rangle /m_e}$ (Wannier 1953). The approximate equations can be solved with respect to $\langle {\boldsymbol v}_{e||} \rangle$ and 〈_e〉; in the limit of high fields, the result is

$$\begin{equation} \langle {\boldsymbol v}_{e||} \rangle \approx - \left ({\displaystyle \frac{\kappa _{en}}{2}} \right) ^{1/4} \sqrt{\frac{eE\ell _e}{m_e}}\hat{{\boldsymbol E}} \approx - \frac{\sqrt{eE\ell _e}}{(m_em_n)^{1/4}}\hat{{\boldsymbol E}}, \end{equation} \tag{ A14 }$$

$$\begin{equation} \langle \epsilon _e \rangle \approx \frac{eE\ell _e}{\sqrt{2\kappa _{en}}} \approx \frac{1}{2}\sqrt{\frac{m_n}{m_e}}eE\ell _e, \end{equation} \tag{ A15 }$$

where we have used that λ_en ≈ 1 and κ_en ≈ 2m_e/m_n. Comparing Equation (A15) with the weak-field expression 〈_e〉 ≈ 3k_BT/2, we find that significant electron heating occurs when

$$\begin{equation} E \gg \frac{\sqrt{\kappa _{en}}k_{\rm B}T}{e \ell _e} \sim E_{\rm crit}. \end{equation} \tag{ A16 }$$

Equations (A14) and (A15) agree with the asymptotic expressions of $\langle {\boldsymbol v}_{e||} \rangle$ and 〈_e〉 in the limit of E ≫ E_crit (Equations (23) and (24)) within an relative error of only 11% and 17%, respectively. Equation (A13) also holds within an error of only 3% (see Equation (59)).

Here, we prove the instability of the middle solution using a standard linear perturbation theory. The middle solution is characterized by the balance between the ionization by high-energy electrons and the sticking of ions and electrons onto grains. With this fact in mind, we may simplify Equations (32) and (33) as

$$\begin{equation} \frac{dn_i}{dt} = -K_{di}(Z)n_dn_i + K_* n_n n_e, \end{equation} \tag{ B1 }$$

$$\begin{equation} \frac{dn_e}{dt} = -K_{de}(Z)n_dn_e + K_* n_n n_e, \end{equation} \tag{ B2 }$$

respectively. The equilibrium conditions dn_i/dt = 0 and dn_e/dt = 0 give

$$\begin{equation} K_{di}n_dn_i = K_*n_nn_e, \end{equation} \tag{ B3 }$$

$$\begin{equation} K_{de}n_d = K_*n_n, \end{equation} \tag{ B4 }$$

respectively.

To analyze the stability of the equilibrium solution, we consider perturbations δn_i, δn_e, and δZ around the equilibrium values n_i, n_e, and Z. Up to the first order in the perturbations, Equations (B1), (B2), and (35) are written as

$$\begin{equation} \frac{d}{dt}\delta n_i = -K_{di}n_d\delta n_i - K^{\prime }_{di}n_dn_i\delta Z + K_* n_n \delta n_e, \end{equation} \tag{ B5 }$$

$$\begin{equation} \frac{d}{dt}\delta n_e = - K^{\prime }_{de}n_dn_e\delta Z, \end{equation} \tag{ B6 }$$

$$\begin{equation} \delta n_i - \delta n_e + n_d \delta Z = 0. \end{equation} \tag{ B7 }$$

where $K^{\prime }_{d\alpha } \equiv dK_{d\alpha }(Z)/dZ$. In deriving Equation (B6), we have used Equation (B4) to eliminate (− K_den_d + K_*n_n)δn_e. Substituting δn_i, δn_e, δZ∝exp (st) into Equations (B5)–(B7), we obtain the equation for s,

$$\begin{equation} s^2 + s(K_{di}n_d + K^{\prime }_{de}n_e - K^{\prime }_{di}n_i) + (K_i - K_*)K^{\prime }_{de} n_en_d=0. \end{equation} \tag{ B8 }$$

The equilibrium solution is stable against the charge perturbations if the roots of Equation (B8) are all negative. However, Equation (B3) and (B4) imply that

$$\begin{equation} \frac{K_i}{K_*} = \frac{n_e}{n_i}. \end{equation} \tag{ B9 }$$

Since n_e < n_i for Z < 0, Equation (B9) suggests K_i − K_* < 0. Furthermore, $K^{\prime }_{de}$ is generally positive because the collisional cross section between grains and electrons increases with increasing Z (or decreasing −Z). Since K_i − K_* < 0 and $K^{\prime }_{de}>0$, the third term in Equation (B8) is negative, implying that the equation has one positive and one negative root. The existence of a positive s means that the middle solution is unstable.