IONIZATION AND DUST CHARGING IN PROTOPLANETARY DISKS

An accurate calculation of ionization–recombination balance in dense protoplanetary conditions is essential for understanding various fundamental problems, such as the coupling of gas with the magnetic field (Li et al. 2014, pp. 101–123), accretion processes (Turner et al. 2014, pp. 411–432), chemistry (Semenov et al. 2004; Larsson et al. 2012), and dust evolution (Okuzumi et al. 2011b; Akimkin 2015). Both the ionization and recombination processes can arise from several sources. While the treatment of ionization, despite the variety of ionization sources, could be reduced to a single (total) ionization rate, the description of recombination is less straightforward. At sufficiently high densities, the dominant sink of free electrons and ions is dust grains, and the recombination rate non-trivially depends on the properties of the grains. Furthermore, the collection of electrons and ions leads to non-zero grain charges, which effectively changes the grain–grain (Okuzumi 2009) and ion–grain (Weingartner & Draine 1999) interactions as well as the grain dynamics.

Depletion of electrons caused by the presence of dust grains significantly reduces the degree of ionization in dense interstellar conditions (Umebayashi 1983; Umebayashi & Nakano 1990; Nishi et al. 1991): in comparison with dust-free gas, the electron-to-ion ratio may drop by as much as a square root of the effective ion-to-electron mass ratio (which is a factor of 74 for a plasma with dominant ${{\rm{H}}}_{3}^{+}$ ions, or 231 for N₂H⁺/HCO⁺ ions). As ionization controls the coupling of the gas to the magnetic field, and hence the development of magnetorotational instability (MRI; e.g., Velikhov 1959; Balbus & Hawley 1991), dust is an essential ingredient for any MRI model. It has been shown that the grain size critically affects the size of a disk's "dead zone" (Sano et al. 2000; Salmeron & Wardle 2008; Bai 2011a, 2011b; Dudorov & Khaibrakhmanov 2014). Nevertheless, analysis of MRI has usually been carried out assuming that properties of dust are fixed.

In dense protoplanetary environments, the coagulation of submicron interstellar dust particles becomes an important process. The planet formation in protoplanetary disks requires the dust to form larger and larger aggregates, until gravitational forces become dominant (e.g., Testi et al. 2014, pp. 339–361, and references therein). There is clear evidence of grain growth to millimeter and centimeter sizes within protoplanetary disks (e.g., Pérez et al. 2015; van der Marel et al. 2015). However, significant difficulties are found during this coagulation process, such as bouncing barriers (e.g., Zsom et al. 2010) and particle fragmentation (Birnstiel et al. 2012) after initial grain compaction and growth. Many theoretical and laboratory studies have greatly advanced our understanding of grain growth and planetesimal formation in recent years (Dominik et al. 2007, pp. 783–800; Johansen et al. 2014, pp. 547–570, and references therein), with particular attention dedicated to dust traps, now detected with ALMA toward protoplanetary disks (van der Marel et al. 2013, 2015; Flock et al. 2015; Pinilla et al. 2015; Ruge et al. 2016; Zhang et al. 2016). In dust traps, particles are expected to grow more easily due to the locally enhanced dust-to-gas mass ratio (Booth & Clarke 2016; Surville et al. 2016), although the details of this coagulation process are far from being understood, considering the largely unknown dust properties.

As has already been pointed out, ionization not only determines dynamical and chemical processes occurring in protoplanetary disks, but it also leads to dust charging and thus affects coagulation. The collection of electrons and ions results in (on average) negative grain charges due to higher electron velocities. Recently, it has been shown that the coagulation of larger aggregates in protoplanetary disks can be inhibited due to increasing Coulomb repulsion between them—the resulting electrostatic potential barrier is roughly proportional to the aggregate size (Okuzumi 2009; Okuzumi et al. 2011a, 2011b). Along with the plasma charging, other charging mechanisms can operate in protoplanetary disks. In Akimkin (2015) the photoelectric emission from grains, induced by stellar radiation and leading to their positive charging, was considered as a mechanism to overcome the electrostatic barrier in upper disk regions. A similar mechanism—photoelectric charging due to H₂ fluorescence induced by cosmic rays (CRs)—operates in much deeper regions at the disk periphery (Ivlev et al. 2015). However, both mechanisms become negligible in dense regions of the disk. We notice that the (still poorly investigated) effect of charging on the dust evolution has recently received increased interest (Carballido et al. 2016).

It is worth mentioning that the coagulation in protoplanetary disks is accompanied by the formation of porous aggregates characterized by an open, fluffy structure (e.g., Dominik et al. 2007, pp. 783–800; Okuzumi et al. 2009). The porosity has been pointed out to have a strong impact on the ionization in protoplanetary disks (Okuzumi 2009; Dzyurkevich et al. 2013; Mori & Okuzumi 2016). However, while the growth of uncharged aggregates is well studied, an accurate description of their charging as well as of the charging feedback on their further growth poses a serious problem. One of the fundamental difficulties is that, unlike compact spherical grains (whose charging is described using the Orbital Motion Limited (OML) approximation, e.g., Whipple 1981; Fortov et al. 2005), no accurate approximation is known for electron and ion collection by irregular fluffy aggregates. Given these difficulties, here we completely leave aside an in-depth discussion of the porosity effects.

In this paper, we present an analytical model that becomes exact in sufficiently dense astrophysical environments and allows us to self-consistently calculate densities of the charged species, in particular, to obtain the dust charges for arbitrary grain-size distributions. Unlike other known approaches (Ilgner & Nelson 2006; Okuzumi 2009; Fujii et al. 2011; Dzyurkevich et al. 2013; Mori & Okuzumi 2016), our model does not make assumptions on the form of the charge distribution, and yields closed analytical expressions for important limiting cases. The latter enables convenient analysis of results in a general form, in terms of a few dimensionless numbers.

The presented model has only one free parameter (the effective mass of the ions, which we show to have a small effect on results), and can be easily included in numerical simulations following the dust evolution in dense molecular clouds and protoplanetary disks. We employ the model to verify whether the broadly used approximation of low grain charges is valid for typical protoplanetary disks. Furthermore, we identify a "dust–dust" plasma regime, where the grain charge distribution becomes quasi-symmetric with respect to the uncharged state. This leads to removal of the repulsive electrostatic barrier and, for large aggregates, opens a "coagulation window" operating in the  dense inner region of protoplanetary disks. Also, we discuss the importance of self-consistent analysis of the ionization and grain evolution, as these processes are mutually coupled via several mechanisms operating in the disks.

The paper is organized as follows. In Section 2 we consider the overall ionization-recombination balance and introduce a recombination threshold—the gas density above which the electron–ion recombination is dominated by the processes on the dust surface. In Section 3 we present the grain charge distribution determined by the collection of electrons and ions and point out limiting cases of "big" and "small" grains. In Section 4 we derive the governing equations for the dust-phase recombination regime, complemented by the grain charge distribution, which allow us to calculate densities of the charged species in a general form. To reveal generic properties of the solution, we consider "monodisperse" dust (grains of the same size) and investigate separately the big- and small-grain limits. The effect of the grain-size distribution is studied in Section 5. We discuss implications of the proposed model for protoplanetary disks in Section 6 and summarize the results in Section 7.

The ionization–recombination balance for electrons is generally governed by the following equation:

$$\begin{eqnarray}&&\zeta {n}_{{\rm{g}}}={R}_{{\rm{g}}}+{R}_{{\rm{d}}},\end{eqnarray} \tag{ 1 }$$

where ζ is the total ionization rate³ of gas with the number density ${n}_{{\rm{g}}}$. The recombination is represented by the two terms. The first one, ${R}_{{\rm{g}}}={\sum }_{k}{\beta }_{{\rm{g}}}^{(k)}{n}_{{\rm{e}}}{n}_{{\rm{i}}}^{(k)}$, describes the gas-phase recombination of electrons (with the density ${n}_{{\rm{e}}}$) and ions, where ${\beta }_{{\rm{g}}}^{(k)}$ is the rate of recombination for the kth ion species (with the density ${n}_{{\rm{i}}}^{(k)}$; the summation is over all ion species). The second recombination term, ${R}_{{\rm{d}}}={\beta }_{{\rm{d}}}{n}_{{\rm{e}}}{n}_{{\rm{d}}}$, describes the electron collection on dust (with the density ${n}_{{\rm{d}}};$ in equilibrium, the electron collection is equal to the collection of all ion species). The dust-phase recombination rate ${\beta }_{{\rm{d}}}=2\sqrt{2\pi }{a}^{2}{v}_{{\rm{e}}}{e}^{-{\rm{\Psi }}}$ (where ${v}_{{\rm{e}}}=\sqrt{{k}_{{\rm{B}}}T/{m}_{{\rm{e}}}}$ is the thermal velocity scale for electrons) is determined by the normalized potential ${\rm{\Psi }}=| \langle Z\rangle | {e}^{2}/{{ak}}_{{\rm{B}}}T$ of a grain of radius a (Fortov et al. 2005). It is assumed that dust comprises a small mass fraction of gas ${f}_{{\rm{d}}}$ (typically, of the order of 10⁻² in the interstellar medium), i.e., the dust density is proportional to the gas density:

$$\begin{eqnarray}&&{m}_{{\rm{d}}}{n}_{{\rm{d}}}={f}_{{\rm{d}}}{m}_{{\rm{g}}}{n}_{{\rm{g}}},\end{eqnarray} \tag{ 2 }$$

where ${m}_{{\rm{d}}}$ is the mass of a grain and ${m}_{{\rm{g}}}\simeq 2.3{m}_{{\rm{p}}}$ is the mean mass of a gas particle (expressed in units of the proton mass ${m}_{{\rm{p}}}$). We note that generation of local dust traps in the disk (van der Marel et al. 2013; Flock et al. 2015; Ruge et al. 2016) as well as dust settling to the midplane (Zsom et al. 2011) may substantially increase the value of ${f}_{{\rm{d}}}$.

To evaluate the relative contribution of the two recombination terms, one can set ${n}_{{\rm{e}}}={\sum }_{k}{n}_{{\rm{i}}}^{(k)}$ (this is verified in Section 4.1). We substitute Equation (2) in Equation (1) and obtain a quadratic equation for ${n}_{{\rm{e}}}$. The solution suggests that for sufficiently low ${n}_{{\rm{g}}}$ the gas-phase recombination is the dominant process (i.e., ${R}_{{\rm{d}}}$ is negligible), and ${n}_{{\rm{e}}}$ varies with the gas density as $\propto \sqrt{\zeta {n}_{{\rm{g}}}}$.⁴ At higher ${n}_{{\rm{g}}}$ the situation is reversed: when the gas density exceeds a dust-phase recombination threshold,⁵

$$\begin{eqnarray}&&{n}_{{\rm{g}}}\gtrsim {n}_{{\rm{g}}}^{\mathrm{rec}}=\displaystyle \frac{\zeta {\beta }_{{\rm{g}}}{({m}_{{\rm{d}}}/{m}_{{\rm{g}}})}^{2}}{16\pi {f}_{{\rm{d}}}^{2}{v}_{{\rm{e}}}^{2}{a}^{4}}{e}^{2{\rm{\Psi }}},\end{eqnarray} \tag{ 3 }$$

(where ${\beta }_{{\rm{g}}}$ is the characteristic gas-phase recombination rate), the electron/ion density is primarily determined by the recombination on grains. By substituting typical values $\zeta \sim {10}^{-17}$ s⁻¹ and T ∼ 10² K, and taking into account that ${\beta }_{{\rm{g}}}\lesssim {10}^{-7}$ cm³ s⁻¹ (e.g., Okuzumi 2009), for micron-size grains we obtain that the gas-phase recombination in a heavy-ion HCO⁺/N₂H⁺ plasma (${\rm{\Psi }}=3.86$) is negligible for ${n}_{{\rm{g}}}\gg {n}_{{\rm{g}}}^{\mathrm{rec}}\sim {10}^{9}$ cm⁻³; for a ${{\rm{H}}}_{3}^{+}$ plasma (${\rm{\Psi }}=2.94$) the recombination threshold is an order of magnitude lower. The characteristic recombination rate ${\beta }_{{\rm{g}}}$ depends on details of the gas-phase reactions (Oppenheimer & Dalgarno 1974), and therefore generally the value of ${n}_{{\rm{g}}}^{\mathrm{rec}}$ is known only approximately.

Thus, for ${n}_{{\rm{g}}}\gg {n}_{{\rm{g}}}^{\mathrm{rec}}$—below we call this the dust-phase recombination regime—the ionization degree is (practically) not affected by a variety of reactions occurring in the gas phase. The emergence of this regime reflects the growing importance of dust in the global ionization–recombination balance. When the grain-size polydispersity is taken into account, then (depending on the particular shape of the size distribution) the value of ${n}_{{\rm{g}}}^{\mathrm{rec}}$ can decrease significantly (see Section 5).

A stationary discrete charge distribution $N(Z,a)\equiv {N}_{Z}$ is obtained from the detailed equilibrium of the charging master equation (Draine & Sutin 1987; Draine 2011), as presented in Appendix B. The charge distribution, derived for dust grains of a given size a and normalized to the total differential dust density at that size,

$$\begin{eqnarray}&&\displaystyle \sum _{Z}\,{N}_{Z}={{dn}}_{{\rm{d}}}(a)/{da},\end{eqnarray} \tag{ 4 }$$

depends on two dimensionless numbers.

The first number is the effective ion-to-electron mass ratio $\widetilde{m}$, determined by the partial contributions of all ions. It is defined as (Draine & Sutin 1987; Ivlev et al. 2015)

$$\begin{eqnarray}&&\displaystyle \frac{\widetilde{m}}{1836}={\left(\displaystyle \sum _{k}\displaystyle \frac{1}{\sqrt{{A}_{k}}}\displaystyle \frac{{n}_{{\rm{i}}}^{(k)}}{{n}_{{\rm{e}}}}\right)}^{-2}\equiv A{\left(\displaystyle \frac{{n}_{{\rm{e}}}}{{n}_{{\rm{i}}}}\right)}^{2},\end{eqnarray} \tag{ 5 }$$

where ${n}_{{\rm{i}}}\equiv {\sum }_{k}{n}_{{\rm{i}}}^{(k)}$ is the total ion density and A_k is the atomic mass number (in amu) of the kth ion species; the identity in Equation (5) introduces the effective atomic mass of ions A. The value of $\widetilde{m}$ has well-defined upper and lower bounds: when dust does not noticeably affect the overall charge neutrality, we have ${n}_{{\rm{e}}}\simeq {n}_{{\rm{i}}}$ and therefore $\widetilde{m}\simeq 1836A\,( \ggg 1);$ when the contribution of the negatively charged dust is important, electrons become depleted and one can show (see Section 4) that $\widetilde{m}$ asymptotically tends to unity (i.e., ${n}_{{\rm{e}}}/{n}_{{\rm{i}}}$ tends to a small constant of $1/\sqrt{1836A}$), in order to satisfy the charge neutrality.

The second number is the grain potential energy of the unit charge, ${e}^{2}/a$, normalized by ${k}_{{\rm{B}}}T$ (Ivlev et al. 2015),

$$\begin{eqnarray}&&\widetilde{\varphi }=\displaystyle \frac{{e}^{2}}{{{ak}}_{{\rm{B}}}T}=\displaystyle \frac{1.67}{(a/0.1\,\mu {\rm{m}})(T/100\,{\rm{K}})}\end{eqnarray} \tag{ 6 }$$

(equivalently, $\widetilde{\varphi }$ is the inverse reduced temperature; Draine & Sutin 1987). The value of $\widetilde{\varphi }$ can, in principle, be arbitrary small or large. The case $\widetilde{\varphi }\ll 1$ corresponds to situations where grains grow beyond several microns or/and the ambient gas temperature exceeds a few hundred Kelvin—such conditions are at best matched in the inner midplane region of the disk (Armitage 2007). The opposite limit, corresponding to small grains with $a\lesssim 0.1\,\mu {\rm{m}}$ or/and low temperatures of $\lesssim 30$ K, primarily represents the initial coagulation stage or a stage where the particle fragmentation barrier is reached after initial grain growth; this may also represent the outer (colder) disk regions, where the gas density is nevertheless high enough to satisfy condition (3).

From Equations (38) and (39) we obtain the charge distribution for two limiting cases: irrespective of the magnitude of $\widetilde{m}$, for $\widetilde{\varphi }\ll 1$ the distribution tends to a Gaussian form (Draine & Sutin 1987),

$$\begin{eqnarray}&&\widetilde{\varphi }\ll 1:\quad {N}_{Z}\propto {e}^{-{(Z-\langle Z\rangle )}^{2}/2{\sigma }_{Z}^{2}},\end{eqnarray} \tag{ 7 }$$

the average charge $\langle Z\rangle =-{\widetilde{\varphi }}^{-1}{\rm{\Psi }}$ and the charge variance ${\sigma }_{Z}^{2}={\widetilde{\varphi }}^{-1}(1+{\rm{\Psi }})/(2+{\rm{\Psi }})$ are determined by the normalized potential Ψ, which is the solution of the charging equation,

$$\begin{eqnarray}&&(1+{\rm{\Psi }}){e}^{{\rm{\Psi }}}=\sqrt{\widetilde{m}}.\end{eqnarray} \tag{ 8 }$$

For $\widetilde{\varphi }\gg 1$ the distribution is essentially discrete—the singly charged states are

$$\begin{eqnarray}&&\widetilde{\varphi }\gg 1:\quad \displaystyle \frac{{N}_{\pm 1}}{{N}_{0}}\simeq \displaystyle \frac{{\widetilde{m}}^{\mp 1/2}}{\widetilde{\varphi }},\end{eqnarray} \tag{ 9 }$$

i.e., grains with $Z=-1$ are the most abundant; the multiply charged states ($| Z| \geqslant 2$) are usually exponentially small and practically negligible. We notice that the Gaussian charge distribution, usually assumed in ionization models (Okuzumi 2009; Fujii et al. 2011; Dzyurkevich et al. 2013; Mori & Okuzumi 2016), may only be employed for $\widetilde{\varphi }\ll 1$.

In Section 4.2 we identify a narrow range of $\widetilde{\varphi }\,(\sim 1)$ where a gradual transition between the charge triplet (9) and the Gaussian distribution (7) occurs. Also, we discuss the role of the polarization interactions (Draine & Sutin 1987), neglected in our consideration, and show that they have practically no effect on the obtained results.

In the dust-phase recombination regime, the self-consistent ionization degree and the grain charge distribution are determined by a set of two equations: the ionization–recombination balance equation and the charge-neutrality equation.

We start with the derivation of the ionization-recombination equation. The equilibrium charge distribution discussed above is determined by the detailed balance of the electron and ion fluxes on a grain surface, and therefore it does not matter which of these fluxes is used to calculate the recombination. For convenience, we consider the ion collection term, which can be presented in the following general form:

$$\begin{eqnarray}&&{R}_{{\rm{d}}}=2\sqrt{2\pi }\ {v}_{{\rm{i}}}{n}_{{\rm{i}}}\int {da}\ {a}^{2}{ \mathcal N }(a),\end{eqnarray} \tag{ 10 }$$

where ${v}_{{\rm{i}}}=\sqrt{{k}_{{\rm{B}}}T/{{Am}}_{{\rm{p}}}}$ is the effective thermal velocity scale of ions and ${ \mathcal N }(a)$ is the effective number density of grains of radius a, obtained for the ion collection in Appendix C. By substituting ${ \mathcal N }(a)$ from Equation (40), we derive the equation

$$\begin{eqnarray}\zeta {n}_{{\rm{g}}} & = & 2\sqrt{2\pi }\ {v}_{{\rm{i}}}{n}_{{\rm{i}}}\\ & & \times \displaystyle \sum _{Z\geqslant 0}(\sqrt{\widetilde{m}}+{\widetilde{m}}^{-Z})\displaystyle \int {da}\ {a}^{2}{e}^{-Z\widetilde{\varphi }}{N}_{-Z}.\end{eqnarray} \tag{ 11 }$$

Note that the sign of the index in N_Z is inverted, so that the summation is in fact performed over non-positive charge states ($Z\leqslant 0$).

Next, we obtain the charge-neutrality equation. The charge density of dust grains (of a given size) is ${\sum }_{Z}{{ZN}}_{Z}$, where the summation over positive charges can be eliminated by using Equation (38). Then the integration over the size distribution yields

$$\begin{eqnarray}&&{n}_{{\rm{i}}}-{n}_{{\rm{e}}}=\displaystyle \sum _{Z\gt 0}(1-{\widetilde{m}}^{-Z})\int {da}\ {{ZN}}_{-Z},\end{eqnarray} \tag{ 12 }$$

where the summation is over negative charge states.

Equations (11) and (12) are complemented by the normalization condition for the dust density, Equation (4). Again, by eliminating the summation over positive charges we obtain

$$\begin{eqnarray}&&{N}_{0}+\displaystyle \sum _{Z\gt 0}(1+{\widetilde{m}}^{-Z}){N}_{-Z}={{dn}}_{{\rm{d}}}(a)/{da}.\end{eqnarray} \tag{ 13 }$$

The derived set of governing Equations (11)–(13), along with the relation

$$\begin{eqnarray*}&&\widetilde{m}=1836A{({n}_{{\rm{e}}}/{n}_{{\rm{i}}})}^{2}\end{eqnarray*}$$

and Equation (39), allows us to calculate self-consistent charge distributions for arbitrary grain-size distributions and obtain the corresponding electron and (total) ion densities. The grain-size distribution ${{dn}}_{{\rm{d}}}(a)/{da}$, the gas density ${n}_{{\rm{g}}}$, and the ionization rate ζ (which is generally a decreasing function of ${n}_{{\rm{g}}}$) are the input parameters for the derived model. The effective atomic mass A of the ions is the only free parameter (entering through the effective mass ratio $\widetilde{m}$ and the velocity scale ${v}_{{\rm{i}}}\propto 1/\sqrt{A}$). The model is exact as long as the gas-phase recombination plays no role, i.e., when the strong condition (3) is satisfied.

The numerical solution of the derived equations, obtained for conditions of a typical protoplanetary disk and for several characteristic grain-size distributions, is presented in Section 6. However, before discussing these results, in the following sections we study important limiting cases. This allows us to reveal and better understand the critical role of different physical mechanisms controlling ionization and dust charging. First, in Sections 4.1 and 4.2 we study analytically the two limiting cases of $\widetilde{\varphi }\ll 1$ and $\widetilde{\varphi }\gg 1$, assuming that all grains have the same size. Then, in Section 5 we analyze generic effects introduced by the size polydispersity and show how the results derived for monodisperse particles can be generalized for an arbitrary size distribution.

For convenience, in Table 1 we summarize the main notations used throughout the paper.

Table 1.  Notations Used in the Article

Symbol	Meaning
A	effective atomic mass of ions [amu]
a	dust grain radius [cm]
${f}_{{\rm{d}}}$	dust-to-gas mass ratio
${m}_{{\rm{g}}},{m}_{{\rm{d}}},{m}_{{\rm{e}}},{m}_{{\rm{p}}}$	mass of a gas particle, dust grain (of radius a), electron, and proton [g]
$\widetilde{m}$	effective ion-to-electron mass ratio
${N}_{Z}\equiv N(Z,a)$	discrete charge distribution of grains of radius a [cm−4]
${{dn}}_{{\rm{d}}}(a)/{da}$	differential size distribution of grains [cm−4]
${n}_{{\rm{g}}},{n}_{{\rm{e}}},{n}_{{\rm{i}}}$	number density of gas particles, electrons, and the total number density of ions [cm−3]
${n}_{\mathrm{EI}}$	number density of an electron–ion (EI) plasma [cm−3]
${n}_{{\rm{g}}}^{\mathrm{rec}}$	recombination threshold, separating the gas-phase and dust-phase recombination regimes [cm−3]
${n}_{{\rm{g}}}^{\mathrm{dep}}$	electron-depletion threshold, at the transition to a dust–ion (DI) plasma [cm−3]
${n}_{{\rm{g}}}^{\mathrm{asy}}$	asymptotic threshold, at the transition to a dust–dust (DD) plasma [cm−3]
${R}_{{\rm{g}}},{R}_{{\rm{d}}}$	gas-phase and dust-phase recombination rates [cm−3 s−1]
T	gas/dust temperature (same for all species) [K]
${v}_{{\rm{e}}},{v}_{{\rm{i}}},{v}_{{\rm{d}}}$	thermal velocity scale of electrons, ions (of the atomic mass A), and grains (of radius a) [cm s−1]
$Z,\langle Z\rangle$	grain charge state and average charge number
${\beta }_{{\rm{g}}}$	characteristic rate coefficient for the gas-phase recombination [cm3 s−1]
${\beta }_{{\rm{d}}}$	rate coefficient for the dust-phase recombination [cm3 s−1]
ζ	total ionization rate [s−1]
$\widetilde{\varphi }$	normalized grain potential of the unit charge
${\rm{\Psi }}=\langle Z\rangle \widetilde{\varphi }$	normalized grain potential for charge number $\langle Z\rangle$ (limiting case $\widetilde{\varphi }\ll 1$)

Download table as:  ASCIITypeset image

4.1. Case $\widetilde{\varphi }\ll 1$ (Big Grains or/and High Temperature)

The charge distribution for $\widetilde{\varphi }\ll 1$, Equation (7), is a broad Gaussian function (${\sigma }_{Z}\gg 1$), so the summation over Z can be replaced with the integration. Moreover, as we show later, N_Z in this case can be well approximated by a shifted delta function, ${N}_{Z}\propto {n}_{{\rm{d}}}\delta (Z-\langle Z\rangle )$ (where ${n}_{{\rm{d}}}$ is the density of monodisperse particles of radius a; the integration over a in the governing equations is removed). By substituting the delta function in Equations (11) and (12), employing Equation (8), and neglecting small terms ${\widetilde{m}}^{-| \langle Z\rangle | }$, we reduce the governing equations to the following form:

$$\begin{eqnarray}&&\zeta {n}_{{\rm{g}}}=2\sqrt{2\pi }\ {v}_{{\rm{i}}}{n}_{{\rm{i}}}(1+{\rm{\Psi }}){a}^{2}{n}_{{\rm{d}}},\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{n}_{{\rm{i}}}-{n}_{{\rm{e}}}={\rm{\Psi }}{\widetilde{\varphi }}^{-1}{n}_{{\rm{d}}},\end{eqnarray} \tag{ 15 }$$

where the dust and gas densities are related by Equation (2). Together with the charging Equation (8), these equations are solved for ${n}_{{\rm{e}}}$, ${n}_{{\rm{i}}}$, and Ψ; they can be reduced to a single nonlinear equation for Ψ, which depends on the input parameters via the ratio ${n}_{{\rm{g}}}/\zeta ({n}_{{\rm{g}}})$. The presented approach is equivalent to that used by Okuzumi (2009).

The generic behavior of the solution can be easily understood from a simple scaling analysis: as ${n}_{{\rm{d}}}\propto {n}_{{\rm{g}}}$, from Equation (14) it follows that ${n}_{{\rm{i}}}\propto \zeta$. Since ζ is generally a decreasing function of ${n}_{{\rm{g}}}$, we conclude that the lhs of Equation (15) does not increase with ${n}_{{\rm{g}}}$, while the rhs scales as $\propto {\rm{\Psi }}({n}_{{\rm{g}}}){n}_{{\rm{g}}}$. The latter describes a contribution of the negatively charged dust to the overall charge neutrality. As long as ${n}_{{\rm{g}}}$ is sufficiently small, the dust contribution is negligible and the charge neutrality is reduced to ${n}_{{\rm{i}}}-{n}_{{\rm{e}}}\simeq 0$. This corresponds to "regular" electron–ion (EI) plasmas, where ${\rm{\Psi }}={{\rm{\Psi }}}_{\mathrm{EI}}$ does not depend on ${n}_{{\rm{g}}}$ and is determined from Equation (8) with $\widetilde{m}=1836A\equiv {\widetilde{m}}_{\mathrm{EI}}$. The dust contribution becomes crucial at larger ${n}_{{\rm{g}}}$, where Equation (15) can only be satisfied if Ψ and, therefore $\widetilde{m}\propto {({n}_{{\rm{e}}}/{n}_{{\rm{i}}})}^{2}$, decreases with ${n}_{{\rm{g}}}$, i.e., if the electron density is depleted. Hence, there exists a certain electron-depletion threshold ${n}_{{\rm{g}}}^{\mathrm{dep}}$ for the gas density—it identifies a crossover from EI to dust–ion (DI) plasmas, where the charge neutrality is regulated by positive ions and negatively charged dust grains.

We define the electron-depletion threshold ${n}_{{\rm{g}}}^{\mathrm{dep}}$ as ${n}_{{\rm{g}}}$ at which the dust charge density, given by the rhs of Equation (15), becomes equal to the electron density.⁶ This yields

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{n}_{{\rm{g}}}}{\zeta }\right)}_{\mathrm{dep}}=\displaystyle \frac{{({m}_{{\rm{d}}}/{m}_{{\rm{g}}})}^{2}}{4\sqrt{2\pi }{f}_{{\rm{d}}}^{2}{v}_{{\rm{i}}}{a}^{2}}\displaystyle \frac{\widetilde{\varphi }}{{{\rm{\Psi }}}_{\mathrm{EI}}(1+{{\rm{\Psi }}}_{\mathrm{EI}})},\end{eqnarray} \tag{ 16 }$$

where ${({n}_{{\rm{g}}}/\zeta )}_{\mathrm{dep}}\equiv {n}_{{\rm{g}}}^{\mathrm{dep}}/\zeta ({n}_{{\rm{g}}}^{\mathrm{dep}})$ is a function of ${n}_{{\rm{g}}}^{\mathrm{dep}}$. To obtain the magnitude of ${n}_{{\rm{g}}}^{\mathrm{dep}}$, we set for simplicity $\zeta =$ const and assume typical values used to estimate ${n}_{{\rm{g}}}^{\mathrm{rec}}$ from Equation (3); for example, for $\widetilde{\varphi }\sim 0.03$ this yields ${n}_{{\rm{g}}}^{\mathrm{dep}}\sim {10}^{3}{n}_{{\rm{g}}}^{\mathrm{rec}}$. To distinguish between the gas-phase and dust-phase recombination regimes in EI plasmas, below we adopt notations EI[g] and EI[d] for their respective density ranges (of ${n}_{{\rm{g}}}\lesssim {n}_{{\rm{g}}}^{\mathrm{rec}}$ and ${n}_{{\rm{g}}}^{\mathrm{rec}}\lesssim {n}_{{\rm{g}}}\lesssim {n}_{{\rm{g}}}^{\mathrm{dep}}$).

In EI plasmas, the density of electrons (ions) ${n}_{\mathrm{EI}}$ is directly obtained from Equation (14) with ${\rm{\Psi }}={{\rm{\Psi }}}_{\mathrm{EI}}$. This can be expressed in terms of ${({n}_{{\rm{g}}}/\zeta )}_{\mathrm{dep}}$ as

$$\begin{eqnarray}&&\displaystyle \frac{{n}_{\mathrm{EI}}}{\zeta }=2{f}_{{\rm{d}}}\displaystyle \frac{{m}_{{\rm{g}}}}{{m}_{{\rm{d}}}}\displaystyle \frac{{{\rm{\Psi }}}_{\mathrm{EI}}}{\widetilde{\varphi }}{\left(\displaystyle \frac{{n}_{{\rm{g}}}}{\zeta }\right)}_{\mathrm{dep}},\end{eqnarray} \tag{ 17 }$$

i.e., the ratio ${n}_{\mathrm{EI}}/\zeta$ does not depend on ${n}_{{\rm{g}}}$ (note that it is also independent of $\widetilde{\varphi }$). Simultaneously, this equation provides a convenient normalization for the electron and ion densities in DI plasmas. In this case, ${n}_{{\rm{i}}}$ and Ψ can be approximately calculated by neglecting ${n}_{{\rm{e}}}$ in Equation (15). Together with Equations (14), (16), and (17), this leads to the following simple relations for ${n}_{{\rm{g}}}\gtrsim {n}_{{\rm{g}}}^{\mathrm{dep}}$:

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Psi }}(1+{\rm{\Psi }})}{{{\rm{\Psi }}}_{\mathrm{EI}}(1+{{\rm{\Psi }}}_{\mathrm{EI}})}\simeq \displaystyle \frac{{({n}_{{\rm{g}}}/\zeta )}_{\mathrm{dep}}}{({n}_{{\rm{g}}}/\zeta )},\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{n}_{{\rm{i}}}}{{n}_{\mathrm{EI}}}\simeq \displaystyle \frac{1+{{\rm{\Psi }}}_{\mathrm{EI}}}{1+{\rm{\Psi }}}.\end{eqnarray} \tag{ 19 }$$

To derive ${n}_{{\rm{i}}}({n}_{{\rm{g}}})$, one has to substitute ${\rm{\Psi }}({n}_{{\rm{g}}})$ from Equation (18) and ${n}_{\mathrm{EI}}({n}_{{\rm{g}}})$ from Equation (17) into Equation (19). The electron density ${n}_{{\rm{e}}}$ is obtained by substituting ${n}_{{\rm{i}}}$ in Equation (8) and noting that $\sqrt{\widetilde{m}}=({n}_{{\rm{e}}}/{n}_{{\rm{i}}})\sqrt{{\widetilde{m}}_{\mathrm{EI}}}$. This yields

$$\begin{eqnarray}&&\displaystyle \frac{{n}_{{\rm{e}}}}{{n}_{\mathrm{EI}}}\simeq {e}^{-({{\rm{\Psi }}}_{\mathrm{EI}}-{\rm{\Psi }})}.\end{eqnarray} \tag{ 20 }$$

Equations (18)–(20) show that, in a DI plasma, Ψ, and hence the average charge $| \langle Z\rangle | ={\widetilde{\varphi }}^{-1}{\rm{\Psi }}$, monotonically decreases with ${n}_{{\rm{g}}}$ and, when ${\rm{\Psi }}\lesssim 1$, tends to zero as ${\rm{\Psi }}({n}_{{\rm{g}}})\propto \zeta /{n}_{{\rm{g}}}$. Correspondingly, ${n}_{{\rm{i}}}({n}_{{\rm{g}}})$ tends to $(1+{{\rm{\Psi }}}_{\mathrm{EI}})$ times ${n}_{\mathrm{EI}}({n}_{{\rm{g}}})$. The density ratio ${n}_{{\rm{e}}}/{n}_{{\rm{i}}}$ approaches a small (but finite) value of $1/\sqrt{{\widetilde{m}}_{\mathrm{EI}}}$.

We remind that the results presented in this section are obtained by substituting a shifted delta function for the Gaussian charge distribution (Equation (7), entering the governing Equations (11) and (12)). This approximation is formally justified as long as the average charge $| \langle Z\rangle | ={\widetilde{\varphi }}^{-1}{\rm{\Psi }}$ exceeds the charge variance ${\sigma }_{Z}\sim {\widetilde{\varphi }}^{-1/2}$, i.e., for ${\rm{\Psi }}\gtrsim \sqrt{\widetilde{\varphi }}$. Remarkably, it turns out that the derived results also remain valid for ${\rm{\Psi }}\lesssim \sqrt{\widetilde{\varphi }}$: Equations (11) and (12), solved in this case with the Gaussian distribution, directly yield relations (18)–(20).

The condition $| \langle Z\rangle | \sim {\sigma }_{Z}$ (or, equivalently, ${\rm{\Psi }}\sim \sqrt{\widetilde{\varphi }}$) corresponds to another characteristic gas density, the asymptotic threshold ${n}_{{\rm{g}}}^{\mathrm{asy}}$. Using Equation (18), we obtain

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{n}_{{\rm{g}}}}{\zeta }\right)}_{\mathrm{asy}}=\displaystyle \frac{{{\rm{\Psi }}}_{\mathrm{EI}}(1+{{\rm{\Psi }}}_{\mathrm{EI}})}{\sqrt{\widetilde{\varphi }}}{\left(\displaystyle \frac{{n}_{{\rm{g}}}}{\zeta }\right)}_{\mathrm{dep}}.\end{eqnarray} \tag{ 21 }$$

This threshold identifies the next important crossover, now from DI to dust–dust (DD) plasmas, where the charge neutrality is increasingly dominated by a balance between negatively and positively charged grains. For ${n}_{{\rm{g}}}/\zeta \gtrsim {\widetilde{\varphi }}^{-1/2}{(n/\zeta )}_{\mathrm{asy}}$ the average charge becomes less than unity. Since ${\sigma }_{Z}\sim {\widetilde{\varphi }}^{-1/2}\gg 1$ for any ${n}_{{\rm{g}}}$, the Gaussian charge distribution becomes asymptotically symmetric with respect to Z = 0.

We see that the obtained results are completely characterized by the normalized potential ${{\rm{\Psi }}}_{\mathrm{EI}}$ as well as by the characteristic densities ${n}_{\mathrm{EI}}$ and ${n}_{{\rm{g}}}^{\mathrm{dep}}$ (we note that the asymptotic and depletion thresholds are related through Equation (21)). All of these parameters are functions of the effective atomic mass of ions A, which is the only free parameter of our model. The left panel in Figure 1 shows ${{\rm{\Psi }}}_{\mathrm{EI}}$ as well as the asymptotic density ratio ${({n}_{{\rm{i}}}/{n}_{{\rm{e}}})}_{\mathrm{DD}}$, both plotted versus A. In the right panel, functions representing the dependencies ${n}_{\mathrm{EI}}(A)$ and ${n}_{{\rm{g}}}^{\mathrm{dep}}(A)$ are depicted (recall that ${v}_{{\rm{i}}}\propto 1/\sqrt{A}$). In principle, A may vary between unity (hydrogen ions) and some large numbers (e.g., 56 for iron ions), but astrochemical models (Semenov et al. 2004) suggest that molecular and metal ions, such as HCO⁺/N₂H⁺ (A = 29) and Mg⁺ (A = 24), usually dominate in very dense molecular clouds and protoplanetary disks. From Figure 1 we see that the uncertainty in ${n}_{\mathrm{EI}}$ and ${n}_{{\rm{g}}}^{\mathrm{dep}}$, associated with the plasma composition in this case, does not exceed a few dozen percent. This effect is practically negligible next to uncertainties introduced by poorly known grain-size distribution (e.g., Kim et al. 1994; Weingartner & Draine 2001) and grain morphology (the latter determines dependence of the grain mass on its effective size).

Zoom In Zoom Out Reset image size

Figure 2 summarizes the behavior of the average dust charge and of the electron and ion densities in the entire dust-phase recombination regime, and also illustrates modification of the dust charge distribution with increasing ${n}_{{\rm{g}}}$. The densities are normalized by ${n}_{\mathrm{EI}}({n}_{{\rm{g}}})$, so for ${n}_{{\rm{g}}}\ll {n}_{{\rm{g}}}^{\mathrm{dep}}$ we have ${n}_{{\rm{e}}}/{n}_{\mathrm{EI}}={n}_{{\rm{i}}}/{n}_{\mathrm{EI}}=1$. The value of ${{\rm{\Psi }}}_{\mathrm{EI}}$ weakly depends on the average atomic mass of ions A (the plotted curves are for HCO⁺/N₂H⁺ ions). The solid red lines show the exact solution of the governing Equations (14) and (15) together with the charging Equation (8), which are equivalent to Equations (32)–(34) of Okuzumi (2009). With the normalization used, the curves have a universal form applicable for arbitrary set of parameters in the dust-phase recombination regime. The behavior at ${n}_{{\rm{g}}}\gtrsim {n}_{{\rm{g}}}^{\mathrm{dep}}$ is well reproduced by approximate relations (18)–(20), shown by the red dashed lines, with the maximum deviation at ${n}_{{\rm{g}}}={n}_{{\rm{g}}}^{\mathrm{dep}}$ (where the ion density is underestimated by $\simeq 20 \%$, while the grain potential and the electron density are overestimated by $\simeq 30 \%$ and a factor of $\simeq 2.5$, respectively). Note that $\widetilde{\varphi }$ is an arbitrary small parameter in the considered case, and therefore the point ${\rm{\Psi }}\sim \sqrt{\widetilde{\varphi }}$ (which identifies the crossover to a DD plasma) could in principle be located below the asymptotic value of ${n}_{{\rm{e}}}/{n}_{\mathrm{EI}}\to {e}^{-{{\rm{\Psi }}}_{\mathrm{EI}}}$.

Zoom In Zoom Out Reset image size

So far we have assumed that the gas temperature is constant. In protoplanetary disks, the relative temperature increase toward the center is usually not as strong as the increase of the gas density. Nevertheless, this is a noticeable effect that can in fact be straightforwardly incorporated in our model—the temperature simply becomes an additional input parameter. Moreover, it turns out that in many cases, for instance in the disk midplane, T and ${n}_{{\rm{g}}}$ are related by a simple power-law dependence, $T\propto {n}_{{\rm{g}}}^{\epsilon }$ (e.g., with $\epsilon =2/11$; see Armitage 2007). We notice that the gas temperature enters ${v}_{{\rm{i}}}\propto \sqrt{T}$ and $\widetilde{\varphi }\propto {T}^{-1}$, thus affecting the electron-depletion threshold ${n}_{{\rm{g}}}^{\mathrm{dep}}$. From this we immediately infer that the ratio ${({n}_{{\rm{g}}}/\zeta )}_{\mathrm{dep}}$ on the lhs of Equations (16) should be replaced with ${({n}_{{\rm{g}}}^{1+3\epsilon /2}/\zeta )}_{\mathrm{dep}}$ (then the rhs is properly normalized with the respective density scale entering the $T\propto {n}_{{\rm{g}}}^{\epsilon }$ dependence); both ratios on the rhs of Equation (18) are replaced in the same way. With regard to ${n}_{\mathrm{EI}}$, the lhs of Equation (17) should be multiplied with ${({n}_{{\rm{g}}}{/n}_{{\rm{g}}}^{\mathrm{dep}})}^{\epsilon /2}$, i.e., ${n}_{\mathrm{EI}}$ falls off with gas density as $\propto \zeta {/n}_{{\rm{g}}}^{\epsilon /2}$.

When considering recombination on grains, we have so far also implicitly assumed that only free electrons and ions contribute to this process. Such an approach is natural since the thermal velocities of plasma species are much larger than typical thermal velocities of grains, and therefore the terms describing recombination due to mutual dust collisions (Umebayashi 1983; Umebayashi & Nakano 1990; Marchand et al. 2016) are omitted in Equation (11). On the other hand, in a DD plasma grains become the most abundant charged species and, thus, collisions between them may provide an important contribution to the net recombination rate. This occurs when the product ${n}_{{\rm{d}}}{v}_{{\rm{d}}}$ becomes comparable to the thermal ion flux ${n}_{{\rm{i}}}{v}_{{\rm{i}}}$, where ${v}_{{\rm{d}}}$ is the relevant (thermal or non-thermal) scale for the relative velocity of grains. By substituting the asymptotic expression for ${n}_{{\rm{i}}}({n}_{{\rm{g}}})$ we conclude that the recombination mechanism due to mutual dust collisions should be important at ${n}_{{\rm{g}}}/\zeta \gtrsim {({v}_{{\rm{i}}}/{v}_{{\rm{d}}})({n}_{{\rm{g}}}/\zeta )}_{\mathrm{dep}}$. If the thermal motion dominates dust dynamics, the corresponding gas density exceeds ${n}_{{\rm{g}}}^{\mathrm{dep}}$ by many orders of magnitudes: e.g., for micron-size grains in a HCO⁺ plasma, ${n}_{{\rm{g}}}$ should be of the order of $3\times {10}^{5}{n}_{{\rm{g}}}^{\mathrm{dep}}\,(\sim {10}^{4}{n}_{{\rm{g}}}^{\mathrm{asy}}$). However, the condition can be significantly relaxed for grains exhibiting strong non-thermal motion, e.g., due to differential drift or sedimentation (Testi et al. 2014, pp. 339–361). Furthermore, since ${({n}_{{\rm{g}}}/\zeta )}_{\mathrm{dep}}\propto {f}_{{\rm{d}}}^{-2}$, an increase of the dust fraction due to various processes operating in the disk midplane (see Section 6.1) also promotes this mechanism of recombination.

4.2. Case $\widetilde{\varphi }\gg 1$ (Small Grains or/and Low Temperature)

In this case the charge distribution is very different from the Gaussian form—practically, it is limited by the singly charged and neutral states, as it follows from Equation (9). By substituting the charge distribution in Equation (13) we derive

$$\begin{eqnarray}&&{N}_{0}\simeq {(1+\sqrt{\widetilde{m}}/\widetilde{\varphi })}^{-1}{n}_{{\rm{d}}},\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&{N}_{-1}\simeq {(1+\widetilde{\varphi }/\sqrt{\widetilde{m}})}^{-1}{n}_{{\rm{d}}},\end{eqnarray} \tag{ 23 }$$

where small terms $1/(\sqrt{\widetilde{m}}\widetilde{\varphi })$ were neglected. The governing Equations (11) and (12) are reduced to

$$\begin{eqnarray}&&\zeta {n}_{{\rm{g}}}=2\sqrt{2\pi }\ {v}_{{\rm{i}}}{n}_{{\rm{i}}}(1+\sqrt{\widetilde{m}}\ ){a}^{2}{N}_{0},\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&{n}_{{\rm{i}}}-{n}_{{\rm{e}}}=(1-{\widetilde{m}}^{-1}){N}_{-1},\end{eqnarray} \tag{ 25 }$$

where $\widetilde{m}={({n}_{{\rm{e}}}/{n}_{{\rm{i}}})}^{2}{\widetilde{m}}_{\mathrm{EI}}$. By substituting Equation (23) in Equation (25) we obtain the average charge

$$\begin{eqnarray}&&| \langle Z\rangle | =\displaystyle \frac{1-{\widetilde{m}}^{-1}}{1+\widetilde{\varphi }/\sqrt{\widetilde{m}}},\end{eqnarray} \tag{ 26 }$$

which is always smaller than unity and, as expected, tends to zero when $\widetilde{m}\to 1$.

For an EI plasma, one should set $\widetilde{m}={\widetilde{m}}_{\mathrm{EI}}$ and neglect the rhs of the charge-neutrality Equation (25), exactly as in the case $\widetilde{\varphi }\ll 1$. The corresponding electron-depletion threshold is determined by

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{n}_{{\rm{g}}}}{\zeta }\right)}_{\mathrm{dep}}=\displaystyle \frac{{({m}_{{\rm{d}}}/{m}_{{\rm{g}}})}^{2}}{4\sqrt{2\pi }{f}_{{\rm{d}}}^{2}{v}_{{\rm{i}}}{a}^{2}}\displaystyle \frac{{(1+\widetilde{\varphi }/\sqrt{{\widetilde{m}}_{\mathrm{EI}}})}^{2}}{\widetilde{\varphi }},\end{eqnarray} \tag{ 27 }$$

and the plasma density for ${n}_{{\rm{g}}}\lesssim {n}_{{\rm{g}}}^{\mathrm{dep}}$ is

$$\begin{eqnarray}&&\displaystyle \frac{{n}_{\mathrm{EI}}}{\zeta }=2{f}_{{\rm{d}}}\displaystyle \frac{{m}_{{\rm{g}}}}{{m}_{{\rm{d}}}}{(1+\widetilde{\varphi }/\sqrt{{\widetilde{m}}_{\mathrm{EI}}})}^{-1}{\left(\displaystyle \frac{{n}_{{\rm{g}}}}{\zeta }\right)}_{\mathrm{dep}}.\end{eqnarray} \tag{ 28 }$$

We take into account that ${\widetilde{m}}_{\mathrm{EI}}$ is very large, so the terms $\sim {\widetilde{m}}_{\mathrm{EI}}^{-1}$ and $\sim 1/\sqrt{{\widetilde{m}}_{\mathrm{EI}}}$ are omitted in both expressions (whereas the retained terms $\widetilde{\varphi }/\sqrt{{\widetilde{m}}_{\mathrm{EI}}}$ may be arbitrarily large in the considered case).

For a DI plasma, ${n}_{{\rm{g}}}\gtrsim {n}_{{\rm{g}}}^{\mathrm{dep}}$, we derive the following relations:

$$\begin{eqnarray}&&(1-{\widetilde{m}}^{-1})(1+1/\sqrt{\widetilde{m}}){\left(\displaystyle \frac{1+\widetilde{\varphi }/\sqrt{{\widetilde{m}}_{\mathrm{EI}}}}{1+\widetilde{\varphi }/\sqrt{\widetilde{m}}}\right)}^{2}\\ &&\quad \simeq \displaystyle \frac{{({n}_{{\rm{g}}}/\zeta )}_{\mathrm{dep}}}{({n}_{{\rm{g}}}/\zeta )},\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{n}_{{\rm{i}}}}{{n}_{\mathrm{EI}}}\simeq {(1+1/\sqrt{\widetilde{m}})}^{-1}\left(\displaystyle \frac{1+\widetilde{\varphi }/\sqrt{\widetilde{m}}}{1+\widetilde{\varphi }/\sqrt{{\widetilde{m}}_{\mathrm{EI}}}}\right).\end{eqnarray} \tag{ 30 }$$

Equation (29) yields the solution for $\widetilde{m}({n}_{{\rm{g}}})$, and hence for ${n}_{{\rm{e}}}/{n}_{{\rm{i}}}=\sqrt{\widetilde{m}/{\widetilde{m}}_{\mathrm{EI}}};$ the dependence ${n}_{{\rm{i}}}({n}_{{\rm{g}}})$ is obtained by substituting $\widetilde{m}({n}_{{\rm{g}}})$ into Equation (30). Asymptotically we obtain $\widetilde{m}-1\propto | \langle Z\rangle | \propto \zeta /{n}_{{\rm{g}}}$ and ${n}_{{\rm{i}}}\propto \zeta$.

Thus, the qualitative behavior of ${n}_{{\rm{e}}}$, ${n}_{{\rm{i}}}$, and $| \langle Z\rangle |$ remains exactly the same as in the case $\widetilde{\varphi }\ll 1$. On the other hand, the values of ${n}_{{\rm{g}}}^{\mathrm{dep}}$ (Equations (16) or (27)) and ${n}_{\mathrm{EI}}$ (Equations (17) or (28)) can be quite different in the two cases—their relative magnitudes depend on ${{\rm{\Psi }}}_{\mathrm{EI}}$ and $\widetilde{\varphi }$ (for the corresponding case). Nevertheless, the curves for the normalized electron and ion densities, plotted versus the normalized gas density in the case $\widetilde{\varphi }\gg 1$, look similar to those in Figure 2. Interestingly, the absolute value of the asymptotic ion density ${n}_{{\rm{i}}}({n}_{{\rm{g}}})$ for $\widetilde{\varphi }\gg 1$ is exactly half of that derived for $\widetilde{\varphi }\ll 1$. The dust charge distribution gradually changes at ${n}_{{\rm{g}}}\gtrsim {n}_{{\rm{g}}}^{\mathrm{dep}}$, from the asymmetric triplet with the maximum at $Z=-1$, as given by Equation (9) with $\widetilde{m}={\widetilde{m}}_{\mathrm{EI}}$, to a quasi-symmetric triplet ${N}_{\pm 1}/{N}_{0}\simeq {\widetilde{\varphi }}^{-1}$ with the peak at Z = 0. This latter asymptotic form represents a DD plasma, which is discussed in Section 4.1.

Above, we have completely neglected the polarization interactions of electrons and ions with dust grains. These interactions noticeably increase the electron/ion collection cross sections by uncharged (or weakly charged) grains in the case $\widetilde{\varphi }\gg 1$ (Draine & Sutin 1987). As a result, the relative abundances ${N}_{\pm 1}/{N}_{0}$ in Equation (9) are increased by a factor of $\simeq \sqrt{\pi \widetilde{\varphi }/8}$. Correspondingly, $\sqrt{\widetilde{m}}$ in Equations (22) and (23) should be multiplied with this factor, while $\sqrt{\widetilde{m}}$ in Equation (24) should be multiplied with $\sqrt{\pi \widetilde{\varphi }/2}$, and Equation (25) is left unchanged. In practice, such modification does not affect the characteristic densities given by Equations (27) and (28), as this becomes important only for extremely large values of $\widetilde{\varphi }\gtrsim \widetilde{m}$ (∼5 × 10⁴ for a HCO⁺/N₂H⁺ plasma).

To conclude this section, let us make a note on a crossover between the limiting cases $\widetilde{\varphi }\ll 1$ and $\widetilde{\varphi }\gg 1$. In an EI plasma, the magnitude of the average charge given by Equation (26) is always significantly smaller than $| \langle {Z}_{\mathrm{EI}}\rangle | \,(={{\rm{\Psi }}}_{\mathrm{EI}}/\widetilde{\varphi }$) for $\widetilde{\varphi }\ll 1$. This fact does not allow us to smoothly match these two cases: from Equations (16) and (27) we conclude that the value of ${({n}_{{\rm{g}}}/\zeta )}_{\mathrm{dep}}$ for $\widetilde{\varphi }\gg 1$ is a factor of ${{\rm{\Psi }}}_{\mathrm{EI}}(1+{{\rm{\Psi }}}_{\mathrm{EI}})$ larger than that for $\widetilde{\varphi }\ll 1$, when compared at the formal "matching point" $\widetilde{\varphi }=1$. Similarly, from Equations (17) and (28), we obtain that ${n}_{\mathrm{EI}}/\zeta$ calculated at the matching point for $\widetilde{\varphi }\gg 1$ is $1+{{\rm{\Psi }}}_{\mathrm{EI}}$ times the value for $\widetilde{\varphi }\ll 1$. Setting for simplicity $\zeta =$ const and taking a HCO⁺/N₂H⁺ plasma as an example, we obtain the relative mismatch of $\simeq 19$ for ${n}_{{\rm{g}}}^{\mathrm{dep}}$ and $\simeq 4.8$ for ${n}_{\mathrm{EI}}$. The reason for the discrepancy is obvious: Equations (37) and (38) imply that the negatively charged states with $| Z| \geqslant 2$, neglected in Equations (24) and (25), rapidly become dominant when $\widetilde{\varphi }$ decreases below a value of $\simeq \mathrm{ln}\sqrt{{\widetilde{m}}_{\mathrm{EI}}}$ ($\simeq 5.4$ for a HCO⁺/N₂H⁺ plasma). Hence, there is a relatively narrow range (of $1\lesssim \widetilde{\varphi }\lesssim \mathrm{ln}\sqrt{{\widetilde{m}}_{\mathrm{EI}}}$) where a gradual transition between the triplet (9) and the Gaussian distribution (7) takes place.

The grain polydispersity, i.e., the presence of grains of different sizes, plays an exceptionally important role in the discussed processes. The reason for that is twofold, as one can directly see from the governing equations. (i) The integral in the ionization–recombination Equation (11) determines the magnitude of the plasma recombination rate on the grain surface and, hence, the equilibrium plasma (ion) density. The integral can be dominated by different parts of the size distribution, depending on its particular form. (ii) In the same way, the size distribution affects the integral determining the contribution of charged grains in the charge-neutrality Equation (12).

Let us first consider the case $\widetilde{\varphi }\ll 1$ which, remarkably, allows us to obtain rigorous results for an arbitrary size distribution from the solution derived in Section 4.1 for single-size grains. Indeed, keeping the integrals in Equations (11) and (12) yields in this case Equations (14) and (15) with a modified rhs: the "recombination factor" ${a}^{2}{n}_{{\rm{d}}}$ in the former equation is replaced by $\int {{dn}}_{{\rm{d}}}\ {a}^{2}$, and, similarly, the "charge-neutrality factor" ${\widetilde{\varphi }}^{-1}{n}_{{\rm{d}}}$ in the latter equation is replaced by $\int {{dn}}_{{\rm{d}}}\ {\widetilde{\varphi }}^{-1}$ (we remind readers that ${\widetilde{\varphi }}^{-1}\propto a$). Thus, the effects of polydispersity are factorized, leading to a simple renormalization of the characteristic densities.

To illustrate, we employ a widely used power-law dependence for the differential dust density,

$$\begin{eqnarray}&&{{dn}}_{{\rm{d}}}(a)/{da}={{Ca}}^{-p},\end{eqnarray} \tag{ 31 }$$

defined for the size range ${a}_{\min }\leqslant a\leqslant {a}_{\max }$ (here, we naturally assume that the condition $\widetilde{\varphi }\ll 1$ is satisfied for all sizes down to ${a}_{\min }$). The constant C is determined from the dust-to-gas density ratio, Equation (2), where ${m}_{{\rm{d}}}{n}_{{\rm{d}}}$ should be replaced by the integral $\int {{dn}}_{{\rm{d}}}\ {m}_{{\rm{d}}}$. For certainty, we assume that the power-law index p does not exceed the MRN value of 3.5 (Kim et al. 1994; Weingartner & Draine 2001), so the mass integral is always dominated by the upper-size cutoff ${a}_{\max }$. For this size distribution, the integral $\int {{dn}}_{{\rm{d}}}\ {a}^{2}$ in the modified Equation (14) is equal to ${a}_{\max }^{2}{n}_{{\rm{d}}}$ multiplied by the following renormalization factor:

$$\begin{eqnarray}\displaystyle \frac{4-p}{| 3-p| }\left\{\begin{array}{cc}1, & p\lt 3;\\ {\tilde{a}}^{p-3}, & p\gt 3.\end{array}\right.\end{eqnarray} \tag{ 32 }$$

Here, $\tilde{a}\equiv {a}_{\max }/{a}_{\min }\gg 1$ and ${n}_{{\rm{d}}}$ is the effective dust density determined by the condition $\int {{dn}}_{{\rm{d}}}\ {m}_{{\rm{d}}}={m}_{{\rm{d}},\max }{n}_{{\rm{d}}}$, with ${m}_{{\rm{d}},\max }\equiv {m}_{{\rm{d}}}({a}_{\max })$. For the sake of clarity we do not consider situations with $p\simeq 3$ (where the contributions of the upper and lower cutoffs are comparable), allowing us to neglect small terms $\propto {\tilde{a}}^{-| p-3| }$. Similarly, the integral ${\widetilde{\varphi }}^{-1}{n}_{{\rm{d}}}$ in the modified Equation (15) is equal to ${\widetilde{\varphi }}_{\max }^{-1}{n}_{{\rm{d}}}$ multiplied by the renormalization factor

$$\begin{eqnarray}\displaystyle \frac{4-p}{| 2-p| }\left\{\begin{array}{cc}1, & p\lt 2;\\ {\tilde{a}}^{p-2}, & p\gt 2,\end{array}\right.\end{eqnarray} \tag{ 33 }$$

where ${\widetilde{\varphi }}_{\max }\equiv \widetilde{\varphi }({a}_{\max })$. We see that the contribution of smaller grains can become dominant when ${{dn}}_{{\rm{d}}}(a)/{da}$ decreases sufficiently steeply with size. This lowers the plasma density (due to enhanced recombination) and increases the average grain charge, as described by Equations (32) and (33), respectively.

From Equations (32) and (33) we infer that the characteristic densities can be dramatically decreased for polydisperse grains: from Equation (3) we conclude that the resulting recombination threshold ${({n}_{{\rm{g}}}/\zeta )}_{\mathrm{rec}}$ is inversely proportional to the squared dust-phase recombination rate and hence is reduced by the squared factor (32). For the MRN distribution this implies a decrease by a factor of about ${a}_{\max }/{a}_{\min }\simeq 50$ (Mathis et al. 1977). As for the electron-depletion threshold ${({n}_{{\rm{g}}}/\zeta )}_{\mathrm{dep}}$, its value is inversely proportional to the product of the two factors, i.e., the rhs of Equation (16) should be divided by

$$\begin{eqnarray}\displaystyle \frac{{(4-p)}^{2}}{| (2-p)(3-p)| }\left\{\begin{array}{cc}1, & p\lt 2;\\ {\tilde{a}}^{p-2}, & 2\lt p\lt 3;\\ {\tilde{a}}^{2p-5}, & p\gt 3.\end{array}\right.\end{eqnarray} \tag{ 34 }$$

Again, taking the MRN distribution as an extreme example, we obtain a reduction by almost three orders of magnitude! Finally, the density of an EI plasma, Equation (17), is directly derived from Equation (14). Therefore, ${n}_{\mathrm{EI}}/\zeta$ for polydisperse grains should be divided by the factor (32).

In the case $\widetilde{\varphi }\gg 1$, the effects of polydispersity are generally no longer factorized—the factors ${a}^{2}{N}_{0}$ and ${N}_{-1}$ on the rhs of Equations (24) and (25), respectively, should be replaced with integrals. Using Equations (22) and (23) we obtain the respective integrals $\int {{dn}}_{{\rm{d}}}\ {a}^{2}{(1+\sqrt{\widetilde{m}}/\widetilde{\varphi })}^{-1}$ and $\int {{dn}}_{{\rm{d}}}\ {(1+\widetilde{\varphi }/\sqrt{\widetilde{m}})}^{-1}$. One can see that a renormalization of the characteristic densities is only possible when the range of $\widetilde{\varphi }$ (corresponding to a given range of grain sizes) does not overlap with the value of $\sqrt{\widetilde{m}}$. As the latter decreases monotonically with ${n}_{{\rm{g}}}$, from $\sqrt{{\widetilde{m}}_{\mathrm{EI}}}\,(=231$ for a HCO⁺/N₂H⁺ plasma) to unity, such a situation appears unlikely, and thus the governing equations should be solved numerically. Nevertheless, in the following section (where the numerical results are presented), we demonstrate that the qualitative effects of polydispersity remain similar to those discussed above for the case $\widetilde{\varphi }\ll 1$.

For convenience, in Table 2 we summarize characteristic values of the threshold gas densities. Using these values and taking into account the following scaling dependencies (for compact grains): ${n}_{{\rm{g}}}^{\mathrm{rec}}\propto \zeta {a}^{2}/T$, ${n}_{{\rm{g}}}^{\mathrm{dep}}\propto \zeta {a}^{3}/{T}^{3/2}$, and ${n}_{{\rm{g}}}^{\mathrm{asy}}\propto \zeta {a}^{7/2}/T$, one can deduce the thresholds for arbitrary parameters.

Table 2.  Characteristic Values of the Threshold Gas Densities for $\zeta ={10}^{-17}$ s⁻¹, T = 100 K, and Two "Standard" Models of the Grain-size Distribution

Grain-size Distribution	${n}_{{\rm{g}}}^{\mathrm{rec}}$, cm−3	${n}_{{\rm{g}}}^{\mathrm{dep}}$, cm−3	${n}_{{\rm{g}}}^{\mathrm{asy}}$, cm−3
MRN	2 × 105	3 × 106	1 × 109
$a=0.1\,\mu {\rm{m}}$	8 × 106	2 × 109	2 × 1011

Download table as:  ASCIITypeset image

In this section we employ a typical protoplanetary disk model to numerically calculate the gas ionization and dust charging. In particular, the results are aimed to answer the following questions:

• 1.  
  Can the "coagulation window" be opened, to overcome the electrostatic barrier?
• 2.  
  How accurate is the (conventional) assumption of low dust charges?

We restrict ourselves by considering the disk midplane region, where the gas and dust mass is concentrated. The gas surface density of the disk is assumed to obey a power-law dependence on the disk radius, ${\rm{\Sigma }}{(R)={{\rm{\Sigma }}}_{\mathrm{AU}}(R/\mathrm{au})}^{-q}$ with ${{\rm{\Sigma }}}_{\mathrm{AU}}=200$ g cm⁻² and q = 1, which leads to the disk mass of $0.02\,{M}_{\odot }$ ($0.04\,{M}_{\odot })$ for the outer radius of 150 au (300 au). We consider a sharp cutoff, as a tapered power-law profile does not qualitatively change the results. Furthermore, we assume that the disk thermal structure is determined by re-radiation of the central star emission. The gas/dust temperature is calculated from ${T}^{4}{(R)={T}_{* }^{4}({R}_{* }/R)}^{2}\sin \phi$ (Brauer et al. 2008) with the effective star temperature ${T}_{* }=4000$ K, stellar radius ${R}_{* }=2.6\,{R}_{\odot }$, and $\sin \phi =0.05$ for the sine of the grazing angle.⁷ The disk is assumed to be vertically isothermal, with the vertical density structure derived from the hydrostatic equilibrium for the central star mass of ${M}_{* }=0.7\,{M}_{\odot }$. We consider CRs, stellar X-rays, and radionuclides as primary sources of ionization in the midplane (Turner & Sano 2008),

$$\begin{eqnarray}&&\zeta (R)={\zeta }_{\mathrm{CR}}{e}^{-\tfrac{{\rm{\Sigma }}(R)}{2{{\rm{\Sigma }}}_{\mathrm{CR}}}}+\displaystyle \frac{{\zeta }_{\mathrm{XR}}}{{(R/\mathrm{au})}^{2}}{e}^{-\tfrac{{\rm{\Sigma }}(R)}{2{{\rm{\Sigma }}}_{\mathrm{XR}}}}+{\zeta }_{\mathrm{RA}},\end{eqnarray} \tag{ 35 }$$

where ${\zeta }_{\mathrm{CR}}={10}^{-17}$ s⁻¹ is the (unattenuated) CR-ionization rate with the attenuation surface density of ${{\rm{\Sigma }}}_{\mathrm{CR}}=96$ g cm⁻², ${\zeta }_{\mathrm{XR}}=5.2\times {10}^{-15}$ s⁻¹ is the X-ray ionization rate at 1 au (which corresponds to the X-ray luminosity of the central star of $\simeq 2\times {10}^{30}$ erg s⁻¹) with the attenuation surface density of ${{\rm{\Sigma }}}_{\mathrm{XR}}=8$ g cm⁻², and ${\zeta }_{\mathrm{RA}}={10}^{-21}$ s⁻¹ is the ionization rate due to long-lived radioactive nuclei.

An important ingredient of any disk model are dust properties, i.e., the dust density, size distribution, and grain morphology. Below we consider four characteristic models for the dust-size distribution—two "monodisperse" (single-size) populations of grains, and two MRN-like distributions (Mathis et al. 1977), given by Equation (31) with p = 3.5:

• 1.  
  monodisperse, $a=0.1\,\mu {\rm{m}}$;
• 2.  
  monodisperse, $a=10\,\mu {\rm{m}}$;
• 3.  
  MRN: ${a}_{\min }=0.005\,\mu {\rm{m}}$, ${a}_{\max }=0.25\,\mu {\rm{m}}$;
• 4.  
  evolved MRN: ${a}_{\min }=0.5\,\mu {\rm{m}}$, ${a}_{\max }=25\,\mu {\rm{m}}$.

For all models, grains are supposed to be compact spheres, and the dust-to-gas ratio is the same and equal to ${f}_{{\rm{d}}}=0.01$ at any location in the disk. We adopt 3.5 g cm⁻³ for the solid mass density of dust grains, and assume that N₂H⁺ or HCO⁺ are the dominant ions (A = 29).

6.1. Can the "Coagulation Window" be Opened?

The electrostatic barrier against the dust growth, caused by the mutual repulsion of the negatively charged grains (Okuzumi 2009), presents a fundamental but still poorly investigated issue of the modern dust evolution models. The issue stems from estimates (Okuzumi 2009) showing that, in most regions of a protoplanetary disk, the average electrostatic energy of interaction between grains at their contact can exceed their relative kinetic energy. Indeed, for "big" grains ($\widetilde{\varphi }\ll 1$) in EI plasmas, the ratio of the electrostatic energy to the thermal energy, $\simeq \tfrac{1}{2}| \langle Z\rangle | {{\rm{\Psi }}}_{\mathrm{EI}}$, is always very large.

The two obvious ways to overcome the electrostatic barrier are to decrease the magnitude of the grain charges or to increase the relative velocity of the grains. Recently it was shown that the photoelectric emission caused by stellar UV radiation may drive grain charges to positive values and thus allow dust coagulation at intermediate heights of the protoplanetary disks (Akimkin 2015); similarly, the photoelectric grain charging due to CR-induced H₂ fluorescence can operate in the much deeper, outer midplane regions of the disk (Ivlev et al. 2015). None of these mechanisms, however, is able to affect dust charging in dense disk regions. Okuzumi et al. (2011a) suggested that the presence of a large number of small grains may remove free electrons from the gas in these regions and thus make larger grains less charged. Let us elaborate on the latter mechanism.

Under low-density/high-ionization conditions with ${n}_{{\rm{g}}}/\zeta \ \lesssim {({n}_{{\rm{g}}}/\zeta )}_{\mathrm{dep}}$ (EI plasmas), grains are highly charged due to the abundance of free electrons. The average grain charge in this case tends to the maximum possible value of ${\widetilde{\varphi }}^{-1}{{\rm{\Psi }}}_{\mathrm{EI}}$, determined by the ion mass and temperature, and therefore the coagulation of micron-size (or larger) grains is usually hampered. In denser regions of the disk, where ${({n}_{{\rm{g}}}/\zeta )}_{\mathrm{dep}}\lt {n}_{{\rm{g}}}/\zeta \lt {({n}_{{\rm{g}}}/\zeta )}_{\mathrm{asy}}$ (DI plasmas), the grain charges are lower due to depletion of free electrons. A specific feature of DI plasmas is that charging of grains of a given size is determined by the entire dust ensemble. By moving into the densest disk regions, where ${n}_{{\rm{g}}}/\zeta \gt {({n}_{{\rm{g}}}/\zeta )}_{\mathrm{asy}}$ (DD plasmas), the depletion of electrons eventually becomes so strong that their accretion onto a neutral grain is practically equal to the ion accretion. Thus, the grain charge distribution becomes practically symmetric with respect to zero, opening up the opportunity for barrier-free coagulation.

The threshold parameters ${({n}_{{\rm{g}}}/\zeta )}_{\mathrm{rec}}$, ${({n}_{{\rm{g}}}/\zeta )}_{\mathrm{dep}}$, and ${({n}_{{\rm{g}}}/\zeta )}_{\mathrm{asy}}$, identifying boundaries between different plasma regions, vary across the disk. As introduced above, ${({n}_{{\rm{g}}}/\zeta )}_{\mathrm{rec}}$ is determined by equal gas and dust contributions to the recombination rates, ${({n}_{{\rm{g}}}/\zeta )}_{\mathrm{dep}}$ corresponds to equal number densities of free electrons and electrons carried by dust grains, while at ${({n}_{{\rm{g}}}/\zeta )}_{\mathrm{asy}}$ the total positive charge is equally distributed between free ions and grains. We express these parameters in terms of the gas density, with the ionization rate according to Equation (35), and depict the resulting plasma regions EI[g], EI[d], DI, and DD in Figure 3 for the four dust models.

Zoom In Zoom Out Reset image size

Figure 3 shows that for the monodisperse 0.1 μm dust model (upper left panel), the DD plasma region is limited to $R\lesssim 4\,\mathrm{au}$ and the EI region is at $R\gtrsim 20\,\mathrm{au}$ (with the DI region filling the gap between them). Introducing the grain polydispersity leads to a strong shift of the plasma boundaries: for the MRN model (lower left panel), the DD plasma region extends up to $R\sim 30\,\mathrm{au}$, while the DI region continues beyond $R\sim 300\,\mathrm{au}$. These results confirm conclusions of Section 5, demonstrating that an excess of small particles in a broad size distribution may dramatically reduce the values of ${({n}_{{\rm{g}}}/\zeta )}_{\mathrm{dep}}$ and ${({n}_{{\rm{g}}}/\zeta )}_{\mathrm{asy}}$. The increase of the overall grain size has the opposite effect: according to the results of Sections 4.1 and 4.2, the boundaries between different plasma regions in this case are shifted to much higher ${n}_{{\rm{g}}}$ (i.e., to smaller R), with the strongest effect being on ${({n}_{{\rm{g}}}/\zeta )}_{\mathrm{dep}}\propto {a}^{4}$ (for monodisperse grains). Indeed, the right panels of Figure 3, depicting the results for monodisperse 10 μm grains and for the evolved MRN model, show no presence of a DD plasma.

Thus, the (initial) MRN size distribution ensures that the coagulation window is opened (DD plasma state) in the entire inner disk, promoting rapid dust growth. However, the latter process rapidly modifies the overall charge balance toward DI and EI plasmas: for the evolved dust distributions (only by two orders of magnitude in size, as in the right panels of Figure 3) the electrostatic barrier is completely restored (EI plasma state) for $R\gtrsim 1\,\mathrm{au}$. The excess of small grains makes the coagulation conditions more favorable, but this factor alone is insufficient for further efficient dust growth.

We conclude that some feedback is necessary to facilitate the coagulation. This could be dust fragmentation leading to, e.g., a bimodal small/big size distribution (where small grains provide conditions for a reduced potential barrier between big grains; Okuzumi et al. 2011a), or this could be an increase in the dust-to-gas ratio ${f}_{{\rm{d}}}$ (Booth & Clarke 2016; Surville et al. 2016), since both the EI–DI and DI–DD plasma boundaries are determined by the value of ${({n}_{{\rm{g}}}/\zeta )}_{\mathrm{dep}}\propto {f}_{{\rm{d}}}^{-2}$. Such feedback may be achieved due to the turbulent motion outside of the dead zone (see Section 6.3), which simultaneously increases the relative grain velocities (Testi et al. 2014, pp. 339–361) and hence further reduces the effect of the electrostatic barrier. Turbulence generates local dust traps (van der Marel et al. 2013; Flock et al. 2015; Ruge et al. 2016), where ${f}_{{\rm{d}}}$ can become as high as a few (Johansen & Youdin 2007; Surville et al. 2016). With such values of ${f}_{{\rm{d}}}$, the coagulation window opens up at $R\sim 1\,\mathrm{au}$ even for the monodisperse 10 μm dust (while for the evolved MRN, ${f}_{{\rm{d}}}\simeq 0.2$ is sufficient).

6.2. Low Charges for Dust Grains: Is this Always Justified?

In the contemporary MRI and astrochemical models of protoplanetary disks, it is routinely assumed that the charge states of dust grains are around zero ($0,\pm 1,\pm 2$). This is generally true for small dust grains ($a\lesssim 0.1\,\mu {\rm{m}}$) in low-temperature gas ($T\lesssim 100$ K) (Tielens 2005, see Equation (5.75)). However, such small grains cannot be representative of the dense protoplanetary environments. As the average charge of larger grains scales linearly with the size, the coagulation inevitably leads to the breakdown of the low-charge assumption (see also Perez-Becker & Chiang 2011; Ilgner 2012).

In Figure 4 we demonstrate the grain charge distributions for the MRN and evolved MRN models, plotted at the radial distance of 1 au (left) and 100 au (right). Note that different dust models correspond to independent simulations and are only shown on the same plot to facilitate the comparison.

Zoom In Zoom Out Reset image size

In the MRN case shown in Figure 4, most grains indeed carry low charges. For R = 1 au, the charge distribution is practically symmetric with respect to the zero state. This symmetry is a characteristic feature of DD plasmas, where the singly charged positive and negative grains are the dominant charge carriers. For R = 100 au, the charge distribution is slightly shifted to the negative values, representing DI plasmas. Low grain charges in the MRN case are due to a large number of small ($a\lesssim 0.01\,\mu {\rm{m}}$) grains—they effectively reduce the abundance of free electrons and prevent larger (sub-micron) grains from being multiply charged. However, as small grains are expected to disappear rapidly during the initial stages of coagulation, larger grains become dominant in the size distribution.

For the evolved MRN distribution in Figure 4, the average grain charge exhibits a linear scaling with the size; the charge may be as high as −50 for 1 μm grains. We notice that the results for R = 1 au and 100 au represent, respectively, DI plasmas and EI plasmas (as follows from the lower right panel of Figure 3), so their direct comparison may appear counterintuitive: for otherwise the same parameters, dust charges in DI plasmas should be lower than in EI plasmas. The observed "discrepancy" is, however, due to the fact that temperatures are higher at smaller R.

The derived charge distributions show that one should be careful assuming low-charge states for protoplanetary disk conditions. A moderate increase of the average grain size (above $\sim 1\,\mu {\rm{m}}$) can easily break this assumption.

6.3. Discussion

The "dead zone" is a region of the protoplanetary disk where the development of MRI is suppressed by non-ideal MHD effects (e.g., Armitage 2007). It is well known that the size of the dead zone should strongly depend on dust properties (Sano et al. 2000; Bai 2011a, 2011b; Dudorov & Khaibrakhmanov 2014). The inner boundary of the dead zone is associated with the thermal ionization near the central star and is usually located at distances well below 1 au (Chatterjee & Tan 2014). To roughly estimate the position of the outer boundary, one can use the condition that the magnetic Reynolds number,

$$\begin{eqnarray}&&{\mathrm{Re}}_{\mathrm{mag}}=\displaystyle \frac{{\alpha }^{1/2}{c}_{{\rm{s}}}^{2}}{{\eta }_{{\rm{O}}}{{\rm{\Omega }}}_{{\rm{K}}}},\end{eqnarray} \tag{ 36 }$$

exceeds unity outside of the dead zone (see, e.g., Equation (170) of Armitage 2015). Here, $\alpha \sim {10}^{-3}$ is the turbulent alpha-parameter (outside of the dead zone), ${c}_{{\rm{s}}}=\sqrt{{k}_{{\rm{B}}}T/(2.3{m}_{{\rm{p}}})}$ is the sound speed, and ${{\rm{\Omega }}}_{{\rm{K}}}=\sqrt{{{GM}}_{* }/{R}^{3}}$ is the Keplerian angular velocity. For the estimates, we assume that the magnetic diffusivity ${\eta }_{{\rm{O}}}$ at high gas densities is determined by the Ohmic term (Xu & Bai 2016), which can be approximated by ${\eta }_{{\rm{O}}}=234({n}_{{\rm{g}}}/{n}_{{\rm{e}}})\sqrt{T/{\rm{K}}}$ for a weakly ionized gas (Spitzer 1962; Armitage 2011).⁸

To demonstrate the influence of the grain-size distribution on the size of the dead zone, in Figure 5 we plot the abundances of the major charged species, calculated for different dust models. One can see that even a moderate increase in the average grain size leads to significant changes in the ionization degree which, in turn, strongly affects the position of the outer dead-zone boundary, marked by the vertical dashed line in each panel. The resulting size of the dead zone is about 30, 10, 3, and 2 au for the MRN, monodisperse 0.1 μm, evolved MRN, and monodisperse 10 μm dust-size distribution, respectively.

Zoom In Zoom Out Reset image size

As was pointed out in Section 6.1, dust coagulation increases the ionization degree, which, in turn, leads to higher grain charges and prevents further coagulation due to a growing electrostatic barrier. On the other hand, higher ionization favors the development of MRI and thus stimulates the coagulation via turbulent dust motion. The situation becomes even more complicated when the turbulence leads to the fragmentation of dust aggregates—apart from destruction, this process generates new populations of small dust, which may reduce the electrostatic barrier and promote coagulation of larger grains. Altogether, this suggests the existence of positive and negative feedback loops that may unpredictably halt or accelerate the coagulation at different locations in the disk, and also highlights the importance of a self-consistent analysis of the ionization and dust evolution processes.

We have developed an exact analytical model that describes ionization and dust charging in dense protoplanetary disk conditions for an arbitrary grain-size distribution. Unlike previously developed approaches (Ilgner & Nelson 2006; Okuzumi 2009; Fujii et al. 2011; Dzyurkevich et al. 2013; Mori & Okuzumi 2016), our model does not make assumptions about the form of the grain charge distribution and enables convenient analysis of results in a general form, in terms of a few dimensionless numbers, which allows us to identify universality in the behavior of the charged species. The governing equations for the different cases are summarized in Appendix A, Table 3.

Table 3.  Summary of Governing Equations

Dust-size Distribution; Limiting Case	Parameters	Governing Equations	Auxiliary Relations	Solution for
Polydisperse; arbitrary $\widetilde{\varphi }$	${n}_{{\rm{g}}},\zeta ,{f}_{{\rm{d}}},T,A,{{dn}}_{{\rm{d}}}/{da}$	(11)–(13), (38), (39)	(2), (4), (5), (6)	${n}_{{\rm{i}}},{n}_{{\rm{e}}},{N}_{Z}$
Monodisperse; $\widetilde{\varphi }\ll 1$	${n}_{{\rm{g}}},\zeta ,{f}_{{\rm{d}}},T,A,a$	(8), (14), (15)	(2), (5), (6), (7)	${n}_{{\rm{i}}},{n}_{{\rm{e}}},\langle Z\rangle$
Monodisperse; $\widetilde{\varphi }\gg 1$	${n}_{{\rm{g}}},\zeta ,{f}_{{\rm{d}}},T,A,a$	(22)–(24), (25)	(2), (5), (6), (9), (26)	${n}_{{\rm{i}}},{n}_{{\rm{e}}},{N}_{0},{N}_{-1}$

Download table as:  ASCIITypeset image

For given dust properties and conditions of the disk, the presented model has only one free parameter (the effective mass of the ions A), and is developed for the regime where the dust-phase recombination of free electrons and ions dominates over the gas-phase recombination. A transition to this regime occurs in an EI plasma (where charged grains still do not play any role in the overall charge neutrality), and is characterized by the dust-phase recombination threshold ${({n}_{{\rm{g}}}/\zeta )}_{\mathrm{rec}}$ for the gas density.

At higher gas densities, ${n}_{{\rm{g}}}/\zeta \gtrsim {({n}_{{\rm{g}}}/\zeta )}_{\mathrm{rec}}$, charged grains play an increasingly important role in the charge neutrality. We have determined two characteristic parameters, the electron-depletion threshold ${({n}_{{\rm{g}}}/\zeta )}_{\mathrm{dep}}\gg {({n}_{{\rm{g}}}/\zeta )}_{\mathrm{rec}}$ and the asymptotic threshold ${({n}_{{\rm{g}}}/\zeta )}_{\mathrm{asy}}\gg {({n}_{{\rm{g}}}/\zeta )}_{\mathrm{dep}}$, marking, respectively, transitions from the EI to DI plasma state, and then to the dust–dust (DD) state. The thresholds are determined in such a way that at ${n}_{{\rm{g}}}/\zeta ={({n}_{{\rm{g}}}/\zeta )}_{\mathrm{dep}}$ electrons and negative grains equally contribute to the total negative charge, while at ${n}_{{\rm{g}}}/\zeta ={({n}_{{\rm{g}}}/\zeta )}_{\mathrm{asy}}$ ions and positive grains provide equal contribution to the total positive charge.

The immediate important implications of the derived results for protoplanetary disks are as follows.

• 1.  
  Unless the dust-size distribution is dominated by grains much smaller than $\sim 1\,\mu {\rm{m}}$, larger grains are typically multiply (negatively) charged. In this case, irrespective of the location in the disk, the average grain charge scales linearly with the size. As the size distribution in protoplanetary disk conditions moves toward bigger grains, the conventional approximation of low grain charges may only be used for (very) initial stages of the disk evolution. The presented results are obtained assuming compact dust, and therefore the implication for porous aggregates (which are expected to carry higher charges due to bigger effective sizes) is even stronger. In situations where an aggregate is approximated by a sphere, the effects of porosity can be straightforwardly included by adopting a fractal scaling law ${m}_{{\rm{d}}}\propto {a}^{D}$, relating the dust mass and the effective size (with the appropriate fractal dimensionality $D\lt 3$, which is known to vary during the dust evolution; Okuzumi et al. 2009).
• 2.  
  The asymptotic transition to a DD plasma implies that the grain charge distribution becomes quasi-symmetric with respect to Z = 0. This completely removes the repulsive electrostatic barrier and opens a "coagulation window" for large aggregates operating in the dense inner region of protoplanetary disks. On the other hand, the DD plasma state only "delays" the formation of the barrier: the coagulation itself leads to decreasing dust number density and a gradual shift back to DI/EI plasmas. The (re)appearance of the electrostatic barrier in this case, with the maximum achieved in the EI state (where the energy of the barrier normalized to the thermal energy of grains, $\sim | \langle Z\rangle | {{\rm{\Psi }}}_{\mathrm{EI}}$, is always very large), may completely inhibit further dust growth. The effect of the barrier can be reduced by various feedback mechanisms operating in the disk and leading to increased local gas-to-dust ratio, such as the dust trapping or moderate fragmentation.

We point out that the dust evolution, change in the charged species abundances, and development of MRI are strongly interrelated processes whose mutual effect is poorly understood: the dust coagulation increases the ionization degree which, in turn, leads to higher grain charges and prevents further coagulation due to a growing electrostatic barrier. On the other hand, higher ionization favors the development of MRI, making a disk turbulent; moderate turbulence facilitates dust growth by increasing the relative grain velocities, while strong turbulence leads to dust fragmentation. The latter generates small grains, which may decrease the electron fraction (asymptotically, by a factor of $\sqrt{{m}_{{\rm{i}}}/{m}_{{\rm{e}}}}\sim {10}^{2}$) and lead to MRI quenching. A complex interplay of these nonlinear processes suggests the existence of multiple positive and negative feedback loops that may dramatically affect the ultimate dust evolution.

A rigorous treatment of the ionization fraction and dust evolution could be critical during all stages in the process that links molecular clouds to stellar systems—this is the motivation behind our work. The exact analytical model presented here can be easily implemented in non-ideal MHD simulations, to properly follow the ionization fraction and the dust growth during the process of protoplanetary disk formation (e.g., Zhao et al. 2016) and evolution (e.g., Dullemond et al. 2007; Armitage 2011, and references therein).

We provide a FORTRAN source code, applicable for arbitrary dust-size distributions and also for porous grains, to calculate abundances of the charged species: http://www.inasan.ru/~akimkin/codes.html.

This research made use of NASA's Astrophysics Data System. The authors thank an anonymous referee for detailed comments and suggestions, and Dr. Bo Zhao and Dr. Sergey Khaibrakhmanov for useful discussions. V.V.A. acknowledges financial support from the Russian Foundation for Basic Research (16-32-00012 mol_a) and from a grant of the President of the Russian Federation (NSh-9576.2016.2). P.C. acknowledges financial support from the European Research Council (ERC Advanced Grant PALs 320620).

For convenience, in Table 3 we list the equations to be used in a general case, as well as in different limiting cases.

When the photo-emission from grains as well as other emission mechanisms are negligible, N_Z is determined by the collection of electrons and ions from the ambient plasma. The sticking probabilities of electrons and ions are both assumed to be equal to unity. The collection cross sections, determined by the electrostatic interactions with a charged grain, are derived from the OML approximation (Whipple 1981; Fortov et al. 2005). For positive charge states, the detailed equilibrium yields (Draine & Sutin 1987; Draine 2011)

$$\begin{eqnarray}&&Z\geqslant 0:\quad \displaystyle \frac{{N}_{Z+1}}{{N}_{Z}}=\left(\displaystyle \frac{{e}^{-Z\widetilde{\varphi }}}{1+(Z+1)\widetilde{\varphi }}\right)\displaystyle \frac{1}{\sqrt{\widetilde{m}}},\end{eqnarray} \tag{ 37 }$$

negative charge states are related to the respective positive states via

$$\begin{eqnarray}&&{N}_{-Z}={\widetilde{m}}^{Z}{N}_{Z}.\end{eqnarray} \tag{ 38 }$$

From Equation (37) we derive

$$\begin{eqnarray}&&Z\gt 0:\quad \displaystyle \frac{{N}_{Z}}{{N}_{0}}=\displaystyle \frac{{e}^{-\tfrac{1}{2}Z(Z-1)\widetilde{\varphi }}}{{\widetilde{m}}^{Z/2}{\displaystyle \prod }_{Z^{\prime} =1}^{Z}(1+Z^{\prime} \widetilde{\varphi })},\end{eqnarray} \tag{ 39 }$$

${N}_{-Z}$ is readily obtained by using Equation (38). Note that the charge distribution is slightly modified when the polarization interactions are taken into account (e.g., Equations (3.3) and (3.4) in Draine & Sutin 1987). If needed, this effect can easily be included in Equations (37)–(39) (Ivlev et al. 2015). To include the sticking probabilities of electrons, ${s}_{{\rm{e}}}(Z)$, and ions, ${s}_{{\rm{i}}}(Z)$, the rhs of Equation (37) should be multiplied by the ratio ${s}_{{\rm{i}}}(Z)/{s}_{{\rm{e}}}(Z+1)$.

The effective number density ${ \mathcal N }(a)$ of grains of radius a, entering the dust-phase recombination term ${R}_{{\rm{d}}}$ in Equation (10), takes into account the electrostatic interaction between ions by charged grains. Depending on the sign of the grain charge, the (geometrical) cross section of the ion collection by the grain is increased or decreased; the corresponding factors are derived from the OML approximation. Summing up partial contributions of all charged states yields

$$\begin{eqnarray*}&&{ \mathcal N }(a)=\displaystyle \sum _{Z\lt 0}(1-Z\widetilde{\varphi }){N}_{Z}+\displaystyle \sum _{Z\geqslant 0}{e}^{-Z\widetilde{\varphi }}{N}_{Z}.\end{eqnarray*}$$

We use Equation (38) to eliminate the summation over positive charges, and employ the recurrent relation (37) to rewrite ${ \mathcal N }$ in an identical form,

$$\begin{eqnarray}&&{ \mathcal N }(a)=\displaystyle \sum _{Z\geqslant 0}(\sqrt{\widetilde{m}}+{\widetilde{m}}^{-Z}){e}^{-Z\widetilde{\varphi }}{N}_{-Z},\end{eqnarray} \tag{ 40 }$$

which is more convenient for the analysis (see Sections 4.1 and 4.2).