EFFECTS OF DUST GROWTH AND SETTLING ON THE IONIZATION BY RADIONUCLIDES. I. FORMULATION AND RESULTS IN A QUIESCENT STATE OF PROTOPLANETARY DISKS

In the core accretion paradigm of planet formation, planetesimal formation in protoplanetary disks is the first step for the formation of rocky planets and the rocky cores of gas giant planets. A standard scenario begins with the collisional growth of submicron-sized dust particles to dust aggregates. They grow further by mutual sticking, gradually settle to the midplane of the disk, and finally form a dense dust layer. The formation of planetesimals in itself remains an open issue whether the gravitational instability of the thin layer or the direct growth of the aggregates occurs.

Observationally, some evidence of dust size growth has been proposed based on dust continuum emission from protoplanetary disks, such as power-law indices in spectral energy distribution and a feature of 10 μm silicate emission (e.g., Beckwith & Sargent 1991; Miyake & Nakagawa 1993, 1995; Beckwith et al. 2000; Kitamura et al. 2002; van Boekel et al. 2003; Lommen et al. 2010; Ricci et al. 2010; Sicilia-Aguilar et al. 2011; Furlan et al. 2011). Dust evolution in which micron-sized interstellar dust particles grow into centimeter-sized particles with settling toward the midplane of the disk has been studied theoretically both in a quiescent state (Weidenschilling 1980; Nakagawa et al. 1981, 1986) and in a local turbulence induced by the shear between dust and gas layers near the midplane (Cuzzi et al. 1993; Weidenschilling 1997, 2004).

In addition to the shear-induced turbulence, there exists global turbulent motion of the gas in protoplanetary disks, which is caused mainly by magnetorotational instability (MRI; e.g., Balbus & Hawley 1991). The global motion of turbulent gas can affect the dust evolution processes extensively: relative motion induced by turbulence increases the mutual collision rate, i.e., the growth rate of dust particles (e.g., Völk et al 1980); turbulent mixing motion lets particles move diffusively (e.g., Cuzzi et al. 1993); and turbulent eddies trap dust particles (e.g., Klahr & Henning 1997; Cuzzi et al. 2001; Johansen et al. 2004). On the other hand, the disk instability induced by MRI depends on the spatial and size distribution of dust particles (e.g., Mizuno et al. 1988; Sano et al. 2000; Nomura 2004). Therefore, unified investigation of the evolution of dust particles and the turbulent motion (disk instability) is needed to understand the very beginning of the planet formation in the disk where dust particles settle toward the midplane and form a dust layer.

Magnetic field is expected to play a major role in the self-consistent evolution of dust particles and disk turbulence at this stage. Coupling of the gas with magnetic field that depends mainly on sorts of charged particles and their relative densities is the most important factor in determining which parts of the disk become magnetorotationally unstable. Many authors have examined whether protoplanetary disks can be ionized enough to sustain MHD turbulence (e.g., Gammie 1996; Sano et al. 2000; Ilgner & Nelson 2006a, 2006b, 2006c; Wardle 2007). One of the important findings is that MRI in protoplanetary disks is strongly controlled by the concentration of dust material since surface areas of dust particles provide important places for sticking and reaction (especially recombination) of various chemical (ionized) species in gas phase. Thus, the settling and size growth of dust particles affect extensively the disk instability, and vice versa.

Charge state of dust particles directly affects the collisional growth of dust particles in protoplanetary disks. Okuzumi (2009) pointed out that electrostatic repulsion between dust aggregates becomes strong enough to prevent the mutual collisions, and consequently the negative charging of dust particles in a weakly ionized disk could significantly slow down the coalescence process. Furthermore, Okuzumi et al. (2011a, 2011b) investigated how the charging affects the evolution of the size distribution and found out a location where the charging stops dust growth at the fractal growth stage where the settling of dust particles remains of little importance.

In addition, the chemistry of protoplanetary disks is strongly related to the ionization process and the charge states of gas and dust particles. The chemical evolution in the disk has been investigated extensively by using highly sophisticated models of the chemical reaction network (e.g., Finocchi & Gail 1997; Ilgner & Nelson 2006a; Bergin et al. 2007; Walsh et al. 2012). Although the model of chemical evolution was constrained by the chemical record gathered from meteorites, planetary atmospheres, and cometary comae, large questions such as the origin of cometary ices and the organic contents of gas/solid reservoirs that produced life remain regarding the chemical evolution of gas and dust particles in protoplanetary disks.

In the part of protoplanetary disks where the temperature is below several hundred kelvin and thermal ionization is inefficient, the gas is mainly ionized by non-thermal particles such as cosmic rays, X-rays, and UV radiation from central stellar objects (e.g., Matsumura & Pudritz 2003; Bergin et al. 2007). As the depth χ (in column density) from the surface of the disk increases, the flux of these particles decreases extensively. A very thin layer with the characteristic depth smaller than 10⁻¹ g cm⁻² is adequate to attenuate UV photons (Perez-Becker & Chiang 2011). X-ray observations of star-forming regions revealed a surprisingly strong X-ray activity of young stellar objects (YSOs). Based on the observations (ROSAT, Einstein, and ASCA) of soft X-ray emission from YSOs, Igea & Glassgold (1999) investigated the ionization rate by low-energy X-rays (<20 keV) in the inner region (≲ 5 AU) of the disk and found that X-rays penetrate and ionize the gas at the depth from the surface up to ∼10 g cm⁻². Recently, the Chandra and XMM-Newton observations show that hard X-ray emission up to ∼8–10 keV is common in disk-bearing YSOs (e.g., Preibisch et al. 2005). Furthermore, superhot flares from pre-main-sequence stars observed in the Chandra Orion Ultradeep Project are definitely hotter than 100 MK, with derived peak temperature above 200 MK (Getman et al. 2008). Although we need detailed calculations of hard X-ray penetration in the disk, the ionization by X-rays can be efficient in the upper zone with χ ≲ 100 g cm⁻² (e.g., Ilgner & Nelson 2006c). The ionization rate by highly energetic (cosmic-ray) particles including the effects of secondary particles such as electrons and γ-rays decreases exponentially with the characteristic (attenuation) length λ_CR ≈ 96 g cm⁻² (Umebayashi & Nakano 1981). At depth χ larger than some critical value, matter is ionized primarily by the activity of radionuclides contained in the disk. As shown by Umebayashi & Nakano (1981), both the long-lived radionuclides, which have half-lives t_1/2 much longer than 10⁸ yr and survive in many objects of the present solar system, and the short-lived ones, which have t_1/2 ≪ 10⁸ yr and are already extinct in the objects of the present solar system but alive to some extent in the primitive solar nebula and in the other dense objects, contribute the ionization of matter. Actually, there is evidence in meteorites that short-lived nuclides existed in the formation stage of the Sun and the solar system about 4.6 billion years ago (e.g., Wadhwa et al. 2007).

Many γ-ray spectroscopy instruments have also detected γ-ray emission in a narrow line at 1809 keV from the equatorial plane of the Galaxy (Mahoney et al. 1984; Harris et al. 1990, 1997; Plüschke et al. 2001; Smith 2004; Diehl et al. 2006). This line is interpreted as resulting from the decay of ²⁶Al in the present interstellar medium (ISM), whose mass ratio ²⁶Al/²⁷Al is about six times smaller than the ratio estimated by the analysis of meteoritic material (Wadhwa et al. 2007). Detection of the ⁶⁰Fe decay γ-ray line from the inner Galaxy was also reported (Wang et al. 2007). These short-lived radionuclides must be generated continuously by nucleosynthesis in the Galaxy and exist universally in the ISM from which stars and exoplanetary systems are formed.

The ionization rate by the decay of radionuclides was also investigated by many authors (Cameron 1962; Nakano & Tademaru 1972; Umebayashi & Nakano 1981; Stepinski 1992; Finocchi & Gail 1997). Umebayashi & Nakano (2009, hereafter UN09) have reinvestigated comprehensively the ionization rates both by long- and short-lived nuclides whose relative abundance has been measured in some meteorites and/or by γ-ray astronomy. Because most of these radionuclides can be naturally condensed in dust material, the ionization rate by radionuclides is greatly affected by the change in size and spatial distribution of dust particles. Thus, UN09 pointed out that the growth and settling of dust particles, which are believed to be important processes leading to planet formation in the disk, affect the ionization rate by radionuclides. Inversely, the ionization rate is one of the essential factors for the ionization state in the disk, and it can affect the coupling of the gas with magnetic field and the charge state of dust particles. A unified approach to the physical and chemical evolution of the disk is needed to understand the early stage of the planet formation that could lead to planetesimal formation.

We investigate in this series of papers the fundamental process that the ionization rate by radionuclides is affected by the change in size and spatial distribution of dust particles. The size of dust particles must increase significantly by coalescence through mutual collisions, and the spatial distribution of dust particles must change considerably by settling toward the midplane of the disk. Because it is difficult to solve the dynamics of the disk simultaneously, we calculate the change of the ionization rate by radionuclides in various models of protoplanetary disks.

In this paper, we investigate the elementary processes about the settling and size growth of dust particles and the energy loss of energetic particles emitted by radionuclides and apply them to the disk models similar to the solar nebula in a quiescent state. Because we want to clarify the mechanism for changing the ionization rate, we simply assume in numerical calculations represented in this paper that a coalescence model deals with perfect sticking of the compact sphere. In Section 2, we summarize the model of protoplanetary disks and the motion of dust particles and formulate the settling and size growth (coalescence) equation of dust particles generally. In Section 3, we investigate the elementary loss processes of available energy emitted by the decay of radionuclides, i.e., energy loss in the solid material of dust particles and absorption (shielding) of emitted energies by other dust particles located nearby. Using the coalescence equation for settling dust particles obtained in Section 2, we numerically simulate the evolution of the dust size distribution, calculate the change of the ionization rates by radionuclides at representative orbits of various disk models, and find the maximum ionization rates around the midplane in Section 4. Discussion of some related problems, such as adequacy for the coalescence model of dust particles, contribution of minor radionuclides, and effects of internal structure of dust aggregates, occurs in Section 5. In Appendices A and B, we summarize the mean free path in a gas mixture and calculate the traveling path length of the emitted particle in the dust particle, the results of which are essential to the elementary processes. The ionization rate by radionuclides under a turbulent motion of the disk and/or various models of non-compact dust aggregates will be investigated in the forthcoming paper (T. Umebayashi et al., in preparation, hereafter referred to as Paper II).

2.1. Disk Model

There seems to be no really acceptable model of protoplanetary disks in the evolutionary stage of dust settling and size growth. In view of such a situation, we adopt a simplified model used previously in order to investigate the effects of the dust evolution on the ionization rate by radionuclides. As the gas surface density Σ and the gas and dust temperature T of the disk we assume a power-law distribution given by

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

$$\begin{eqnarray} T(R) & = T_0 \left(\frac{R}{1\,\mathrm{AU}} \right)^{-p}, \end{eqnarray} \tag{ 2 }$$

where R is the radial distance in the conventional cylindrical coordinate system (R, ϕ, z) anchored at the central star and p and q are arbitrary positive parameters (e.g., Cassen 1994). We adopt this temperature distribution even for the opaque region near the midplane because it may be regarded as being in a thermal bath of this temperature. Hayashi (1981) proposed a model with Σ₀ = Σ_H = 1700 g cm⁻², T₀ = 280(L_⋆/1 L_☉)^1/4 K, p = 1/2, and q = 3/2 for the minimum mass solar nebula (MMSN; e.g., Safronov 1969; Hayashi et al. 1985), where L_⋆ is the luminosity of the central star. In the following, we primarily investigate the evolution of the ionization rate in the models with q = 3/2 and $\Sigma _0 = f_\Sigma \Sigma _{\rm H}$, where f_Σ (≳ 1) is an arbitrary parameter. Indeed, Zhu et al. (2010) suggested that at least the inner part (dead zone) of the solar nebula can be more massive than the MMSN.

The surface density of dust particles is

$$\begin{equation} \Sigma ^\mathrm{dust}(R) = \Sigma _0^\mathrm{dust} \left(\frac{R}{1\,\mathrm{AU}} \right)^{-q}, \end{equation} \tag{ 3 }$$

where Σ^dust₀ depends on the temperature of the dust particles at distance R. For the MMSN we have Σ^dust₀ ≈ 7.1 and 30 g cm⁻² for R ≲ 2.7 and ≳ 2.7 AU, respectively, taking into account the effect of water ice sublimation (Hayashi 1981). Note that Σ^dust₀ is also proportional to f_Σ and the border radius between different Σ^dust₀ depends on T₀ and p. The sound velocity is given by

$$\begin{eqnarray} c_\mathrm{s}(R) &=& \left(\frac{kT(R)}{\mu m_\mathrm{H}} \right)^{1/2} = 9.93 \times 10^4 \left(\frac{T_0}{280\,\mathrm{K}} \right)^{1/2}\nonumber\\ &&\times \left(\frac{\mu }{2.34} \right)^{-1/2} \left(\frac{R}{1\,\mathrm{AU}} \right)^{-p/2}\,\mathrm{cm\,s}^{-1}, \end{eqnarray} \tag{ 4 }$$

where μ is the mean molecular weight of the gas and k and m_H are the Boltzmann constant and the mass of atomic hydrogen, respectively.

In a gravitationally stable disk, the gravity perpendicular to the disk is mainly contributed by the central star. The balance of the stellar gravity and the pressure force in the vertical direction determines the density at height z from the midplane as

$$\begin{equation} \rho (R, z) = \rho _\mathrm{c}(R) \exp \left(-\frac{z^2}{H^2(R)} \right), \end{equation} \tag{ 5 }$$

where

$$\begin{eqnarray} H(R) &=& \frac{\sqrt{2}c_\mathrm{s}(R)}{\Omega _\mathrm{K}(R)} = 4.72 \times 10^{-2} \left(\frac{T_0}{280\,\mathrm{K}} \right)^{1/2} \left(\frac{\mu }{2.34} \right)^{-1/2}\nonumber\\ &&\times \left(\frac{M_\star }{1\,M_\odot } \right)^{-1/2} \left(\frac{R}{1\,\mathrm{AU}} \right)^{(3-p)/2}\,\mathrm{AU} \end{eqnarray} \tag{ 6 }$$

is the scale height of the disk and

$$\begin{eqnarray} \Omega _\mathrm{K}(R) &=& \left(\frac{GM_\star }{R^3} \right)^{1/2}\nonumber\\ &=& 1.99 \times 10^{-7}\left(\frac{M_\star }{1\,M_\odot } \right)^{1/2} \left(\frac{R}{1\,\mathrm{AU}} \right)^{-3/2}\,\mathrm{s}^{-1}\qquad \end{eqnarray} \tag{ 7 }$$

is the Keplerian angular velocity, with M_⋆ and G being the mass of the central star and the gravitational constant, respectively. The density at the midplane (z = 0) is given by

$$\begin{eqnarray} \rho _\mathrm{c}(R) &=& \frac{\Sigma (R)}{\sqrt{\pi }H(R)} = 1.36 \times 10^{-9}f_\Sigma \left(\frac{T_0}{280\,\mathrm{K}} \right)^{-1/2} \left(\frac{\mu }{2.34} \right)^{1/2}\nonumber\\ &&\times \left(\frac{M_\star }{1\,M_\odot } \right)^{1/2} \left(\frac{R}{1\,\mathrm{AU}} \right)^{(p-2q-3)/2}\,\mathrm{g\,cm^{-3}}. \end{eqnarray} \tag{ 8 }$$

The disk is gravitationally stable as long as Toomre's (1964) Q = c_sΩ/πGΣ > 1. In our case stable disks must fulfill a condition

$$\begin{equation} f_\Sigma < 55.5\! \left(\frac{T_0}{280\,\mathrm{K}} \right)^{1/2}\! \left(\!\frac{\mu }{2.34}\! \right)^{-1/2}\! \left(\!\frac{M_\star }{1\,M_\odot }\! \right)^{1/2}\! \left(\!\frac{R}{1\,\mathrm{AU}}\! \right)^{q-(p+3)/2}. \end{equation} \tag{ 9 }$$

For the disk with q = 3/2 and f_Σ ≲ 10 this condition is satisfied in the part R ≲ 100 AU, which is far beyond the outermost region of planet formation.

2.2. Equations of Motion for Gas and Dust Particles

The equations of motion for gas and dust particles in the disk are given by

$$\begin{eqnarray} \frac{d\boldsymbol{u}}{dt} & = A\rho ^\mathrm{dust}(\boldsymbol{U}-\boldsymbol{u})-\frac{G M_\star }{r^3}\boldsymbol{r}-\frac{\nabla P^\mathrm{gas}}{\rho ^\mathrm{gas}}, \end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray} \frac{d\boldsymbol{U}}{dt} & = -A\rho ^\mathrm{gas}(\boldsymbol{U}-\boldsymbol{u})-\frac{G M_\star }{r^3}\boldsymbol{r}, \end{eqnarray} \tag{ 11 }$$

where u, U, ρ^gas, and ρ^dust are the velocities (in the inertial frame of reference) and the mass densities of gas and dust particles, respectively, r = (R, ϕ, z) is the position vector in the cylindrical coordinate system anchored at the central star, and P^gas is the gas pressure. Following Epstein's law or Stokes' law, the coefficient of the drag force A is given by

$$\begin{equation} A = \left\lbrace \begin{array}{@{}ll}c_\mathrm{th}/\rho _\mathrm{s}a & \quad \mbox{for } a \lesssim l^\mathrm{gas} \\ 3c_\mathrm{th}l^\mathrm{gas}/2\rho _\mathrm{s}a^2 & \quad \mbox{for } a \gtrsim l^\mathrm{gas}, \end{array} \right. \end{equation} \tag{ 12 }$$

where a and ρ_s are the radius and the inner solid density of dust particles, respectively, and $c_\mathrm{th}=\sqrt{8kT(R)/\pi \mu m_\mathrm{H}}=1.60\,c_\mathrm{s}(R)$ and l^gas are the thermal speed and the mean free path of gas particles, respectively. Around the disk midplane, we have

$$\begin{eqnarray} l^\mathrm{gas} &=& \frac{\mu m_\mathrm{H}}{\overline{\sigma }^\mathrm{\,gas}\rho ^\mathrm{gas}} = 2.87 f_\Sigma ^{-1} \left(\frac{\overline{\sigma }^\mathrm{\,gas}}{10^{-15}\,\mathrm{cm^2}} \right)^{-1} \left(\frac{T_0}{280\,\mathrm{K}} \right)^{1/2}\nonumber \\ &&\times \left(\frac{\mu }{2.34} \right)^{1/2} \left(\frac{M_\star }{1\,M_\odot } \right)^{-1/2} \left(\frac{R}{1\,\mathrm{AU}} \right)^{(3-p+2q)/2}\,\mathrm{cm},\nonumber\\ \end{eqnarray} \tag{ 13 }$$

where $\overline{\sigma }^\mathrm{\,gas}$ is the average cross section for mutual collisions of gas particles. For the neutral gas composed of hydrogen molecules (H₂) and helium (He) atoms with solar abundance, we have $\overline{\sigma }^\mathrm{\,gas} = 2.24 \times 10^{-15}\,\mathrm{cm}^2$ on average (for further details, see Appendix A) and obtain l^gas = 1.3 cm in Earth's region of the MMSN. In this paper, we simply assume that the shape of dust particles is a compact sphere. We note that the area-to-mass ratio, which strongly affects the drag force coefficient and the sticking rate of dust aggregates, depends on their shape (the fractal structure of dust aggregates; see also Section 5).

Furthermore, we assume that the disk is axisymmetric and rotates around the central star nearly with circular Keplerian velocity v_K = RΩ_K(R)e_ϕ, where e_ϕ is a unit vector in the azimuthal (ϕ) direction. We use the gas and dust velocities relative to the Keplerian velocity $\boldsymbol{v} = \boldsymbol{u}-\boldsymbol{v}_\mathrm{K}$ and $\boldsymbol{V} = \boldsymbol{U}-\boldsymbol{v}_\mathrm{K}$. Because the time for dust particles to approach the equilibrium velocity owing to the gas drag force is much smaller than the Keplerian time and the mutual collision time of dust particles as long as the dust radii are small enough (e.g., Nakagawa et al. 1981), the mean motion of dust particles becomes steady soon. By setting ∂/∂t = 0 in Equations (10) and (11), we obtain (for details, see Nakagawa et al. 1986) the terminal velocities in the radial (R) and the vertical (z) directions as

$$\begin{eqnarray} V_R-v_R &=& -\frac{2D\Omega _\mathrm{K}(R)}{D^2+\Omega _\mathrm{K}(R)^2}\eta v_\mathrm{K}, \nonumber\\ V_R &=& -\frac{\rho ^\mathrm{gas}}{\rho }\frac{2D\Omega _\mathrm{K}(R)}{D^2+\Omega _\mathrm{K}(R)^2}\eta v_\mathrm{K}, \end{eqnarray} \tag{ 14 }$$

and

$$\begin{equation} V_z-v_z = -\frac{\Omega _\mathrm{K}(R)^2}{D^\mathrm{gas}} z, \end{equation} \tag{ 15 }$$

where ρ = ρ^gas + ρ^dust, D^gas = Aρ^gas, D = Aρ, and

$$\begin{eqnarray} \eta &=&-\frac{1}{2R\Omega _\mathrm{K}(R)^2}\frac{1}{\rho ^\mathrm{gas}}\frac{\partial P^\mathrm{gas}}{\partial R} \nonumber \\ &=& 2.78 \times 10^{-4}(p+2q+3) \left(\frac{T_0}{280\,\mathrm{K}} \right) \left(\frac{\mu }{2.34} \right)^{-1} \nonumber\\ &&\times\left(\frac{M_\star }{1\,M_\odot } \right)^{-1} \left(\frac{R}{1\,\mathrm{AU}} \right)^{-p+1} \end{eqnarray} \tag{ 16 }$$

is half the ratio of the pressure force to the stellar gravity in the R-direction around the disk midplane (z ≈ 0).

2.3. Coalescence Equation for Settling Dust Particles

We solve the settling and size growth (coalescence) equation of dust particles to investigate the effects of dust evolution on the ionization rate by radionuclides. Letting n(m, z)dm be the number density of dust particles with mass between m and m+dm at the height z from the disk midplane, the size growth of settling dust particles is represented by the equation

$$\begin{eqnarray} &&\frac{\partial n(m, z)}{\partial t}+\frac{\partial [n(m, z)V_z(m, z)]}{\partial z} \nonumber\\ &&\quad= -n(m, z)\int _0^\infty \langle \sigma v\rangle _{m, m^{\prime }}n(m^{\prime }, z)\,dm^{\prime }\nonumber\\ &&\quad\phantom{=} +\frac{1}{2}\int _0^m \langle \sigma v\rangle _{m-m^{\prime }, m^{\prime }}n(m-m^{\prime }, z)n(m^{\prime }, z)\,dm^{\prime }, \end{eqnarray} \tag{ 17 }$$

where the rate coefficient of the dust particle of mass m and radius a colliding to attach with one of m' and a' is given by

$$\begin{equation} \langle \sigma v\rangle _{m, m^{\prime }} = p_\mathrm{s}\pi (a+a^{\prime })^2\delta V, \end{equation} \tag{ 18 }$$

p_s and δV being the probability for sticking together and the relative velocity for mutual collision of dust particles (Nakagawa et al. 1981). We note that p_s depends on size, relative velocity of collision, chemical composition, and/or the shape of dust particles (e.g., Weidenschilling & Cuzzi 1993; Weidenschilling 2004; see also references therein). Owing to the lower values of δV, we simply neglect in Equation (17) the fragmentation of dust particles in the course of mutual collisions, which could occur if the dust particles become large and their velocities relative to small particles become high enough (e.g., Dullemond & Dominik 2005).

The second term on the left-hand side of Equation (17) shows the transport of dust particles in the vertical direction. Because the mean vertical velocity of the gas is negligibly small (v_z ≈ 0) in a quiescent state of the disk, the dust particles of mass m and radius a settle toward the midplane with the mean velocity

$$\begin{equation} V_z(m, z) = \left\lbrace \begin{array}{@{}ll}-(\rho _\mathrm{s}a/\rho ^\mathrm{gas}c_\mathrm{th})\Omega _\mathrm{K}(R)^2z & \quad \mbox{for } a \lesssim l^\mathrm{gas}, \\ -(2\rho _\mathrm{s}a^2/3\rho ^\mathrm{gas}c_\mathrm{th}l^\mathrm{gas})\Omega _\mathrm{K}(R)^2z & \quad \mbox{for } a \gtrsim l^\mathrm{gas}. \end{array} \right. \end{equation} \tag{ 19 }$$

We completely neglect the transport of dust particles in the radial direction since we are mainly interested in the evolutionary stage that a mass fraction of gas in a certain position is larger than that of dust particles where the dust particles scarcely move in the radial direction (Nakagawa et al. 1986).

Numerically, we solve the following discrete equation of coalescence:

$$\begin{eqnarray} \frac{\partial \varphi (i, z_l)}{\partial t}&+&\frac{\partial [\varphi (i, z_l)V_z(i, z_l)]}{\partial z} = -m_i\varphi (i, z_l)\!\sum _{j=1}^n \beta (i,j)\varphi (j, z_l) \nonumber \\ &+&\frac{1}{2}m_i\sum _{j=1}^{i-1}\sum _{k=1}^{i-1}\gamma _{ijk}\beta (j, k)\varphi (j, z_l)\varphi (k, z_l), \end{eqnarray} \tag{ 20 }$$

where

$$\begin{equation} \varphi (i, z_l) = \int _{m_{i-1/2}}^{m_{i+1/2}}m^{\prime }\,n(m^{\prime }, z_l)\,dm^{\prime } \end{equation} \tag{ 21 }$$

is the integrated mass density of dust particles whose mass ranges from m_{i − 1/2} to m_{i + 1/2} in the position of z = z_l, $\beta (j, k) = \langle \sigma v\rangle _{m_j, m_k}/(m_jm_k)$ is the average rate coefficient of the dust particle of mass m_j = (m_{j − 1/2} + m_{j + 1/2})/2 colliding to attach with one of m_k, and γ_ijk is the factor between 0 and 1 that represents the branching ratio that the dust particle of mass m_i is formed through the collision between particles of m_j and m_k. This factor is induced because the integrated form of mass density φ(i, z_l) is used to calculate Equation (17) numerically (Nakagawa et al. 1981; Nakagawa & Kohno 1999; Nomura & Nakagawa 2006). The total mass density of dust particles at a given position z_l can be obtained from

$$\begin{equation} \rho ^\mathrm{dust}(z_l) = \int _0^\infty m^{\prime }\,n(m^{\prime }, z_l)\,dm^{\prime } = \sum _{i=1}^\infty \varphi (i, z_l). \end{equation} \tag{ 22 }$$

The transport term in Equation (20) is calculated by using the spatial grids in the z-direction indexed by discrete number (integer) l.

For the relative collision velocity between dust particles, δV, we take into account the thermal Brownian motion

$$\begin{equation} \delta V_\mathrm{B} = \left(\frac{8kT}{\pi } \right)^{1/2}\left(\frac{1}{m_j}+\frac{1}{m_k} \right)^{1/2}, \end{equation} \tag{ 23 }$$

the velocity differences in the R- and z-directions δV_R = V_R(j) − V_R(k) and δV_z = V_z(j) − V_z(k), and the turbulence-induced relative velocity δV_t, which can be neglected in a quiescent disk. We neglect the velocity difference in the ϕ-direction δV_ϕ because the azimuthal velocity V_ϕ of dust particles considered in this paper hardly depends on their radius a (e.g., Nakagawa et al. 1986). We adopt the square root of the sum of the squares of these velocity differences as relative collision velocity:

$$\begin{equation} \delta V = \big(\delta V_\mathrm{B}^2+\delta V_R^2+\delta V_z^2+\delta V_\mathrm{t}^2\big)^{1/2}. \end{equation} \tag{ 24 }$$

3.1. General Remarks on the Ionization by Radionuclides

As for the long-lived radionuclides, UN09 shows that the dominant nuclide for the ionization of matter is ⁴⁰K and that ²³²Th (Th series), ²³⁵U (Ac series), and ²³⁸U (U series) contribute about 20% of the total ionization rate by the long-lived nuclides. Among the short-lived radionuclides, UN09 reveals that ²⁶Al is the dominant ionization source overwhelming the ionization by the long-lived nuclides and that the ionization rates by ⁶⁰Fe and ³⁶Cl are more efficient than those by the long-lived nuclides. Cosmochemically, elements Al, Th, and U are refractory, Fe is a siderophile element whose condensation temperature for iron–nickel mixtures is about 1700 K, and K (potassium) is a volatile element that condenses at a temperature between 600 K and 1300 K. The only (major radioactive) element that is so highly volatile that its condensation temperature is below 600 K is Cl. Because the ionization by radionuclides becomes important only in the region of the protoplanetary disk with T ≲ 600 K, where the thermal ionization is inefficient (Umebayashi 1983), we can regard that most of these radionuclides are substances of the solid material in dust particles. In the course of investigating the effects of dust settling and size growth on the ionization by radionuclides, we can safely consider that ⁴⁰K and ²⁶Al are the primary sources of the long- and short-lived nuclides for ionization.

The average energy released in the decay series starting at parent nuclide X and used for ionization per unit time per unit hydrogen nucleus is given by

$$\begin{equation} \mathcal {E}_\mathrm{ion}(\mathrm{X}) = \overline{E}_\mathrm{em}(\mathrm{X})w(\mathrm{X})x(\mathrm{X}), \end{equation} \tag{ 25 }$$

where $\overline{E}_\mathrm{em}(\mathrm{X})$, w(X) = ln 2/t_1/2(X), and x(X) are the energy actually used for ionization, the decay rate constant, and the abundance relative to hydrogen nucleus by number, respectively. The energy deposited into the medium is eventually used for ionization and electric excitation of the surrounding gas that is mainly composed of H₂ molecules and He atoms. The ionization rate of an H₂ molecule by the decay of X is then given by

$$\begin{equation} \zeta _\mathrm{R}^{\mathrm{H}_2}(\mathrm{X}) \approx \frac{\overline{E}_\mathrm{em}(\mathrm{X})}{W_{\mathrm{H}_2}}\frac{w(\mathrm{X})\,x(\mathrm{X})}{x(\mathrm{H}_2)+x(\mathrm{He})}, \end{equation} \tag{ 26 }$$

where $W_{\mathrm{H}_2} \approx 36\,\mathrm{eV}$ is the average energy required for energetic particles to produce an ion–electron pair passing through molecular hydrogen gas and x(H2) = 0.5 and x(He) are the relative abundance of H2 molecules and that of He atoms by number (see UN09 for details). Because the ionization rate of an He atom can be obtained through the equation $\zeta _\mathrm{R}^\mathrm{He}(\mathrm{X}) \approx 0.84 \zeta _\mathrm{R}^{\mathrm{H}_2}(\mathrm{X})$ (see UN09), we deal only with ζH2R(X) in the following (hereafter we simply write it as ζ(X) without confusion).

The size growth of settling dust particles affects the ionization by radionuclides through the following three processes: (1) the change of the relative abundance of radionuclides x(X) due to the movement (mainly settling toward the disk midplane) of dust particles, (2) the energy loss of the emitted particles along the traveling path in the solid material of dust particles before they emerge from the surface of dust particles (hereafter we simply call it energy loss), and (3) the absorption (shielding) of the emitted energetic particles by other dust particles located nearby (hereafter shielding). In the following subsection, we investigate the processes of energy loss and shielding in detail.

3.2. Energy Loss in the Solid Material of Dust Particles

In the ordinary environment of the protoplanetary disks we can safely regard that all the radionuclides considered in UN09 are found in the dust particles, as mentioned above. The decay of the parent radionuclide X and subsequent processes emit α-particles, electrons, positrons, and/or photons, and the (kinetic) energies of these emitted particles are at most several MeV. The ranges λ of emitted α-particle, electron (positron), and photon with energy ∼1 MeV are about 0.001, 0.5, and 20 g cm⁻², respectively, in the solid material of dust particles, and these ranges decrease steeply as their energies decrease (Berger 1993; Hubbell & Seltzer 1996). In case the mass thickness χ^dust of the dust particle becomes close to or over the range λ of emitted energetic particles, these particles must lose a certain amount of their (kinetic) energies before they emerge from the surface of dust particle and ionize the gas. Thus, we investigate the energy loss of emitted particles along the traveling path in the solid material and obtain the residual energy available to ionize the gas $\overline{E}_\mathrm{res}(j)$ for each emitted particle j.

Charged particles and photons emitted with energies less than several MeV lose their energies mainly through the ionization of surrounding matter and the Compton effects. Since the energy-loss rates of these processes are approximately taken to be constant around the energy region considered in this paper, the residual energy $\overline{E}_\mathrm{res}(j)$ can depend mainly on the mass thickness χ^dust along the traveling path of the emitted particle in the solid material of the dust particle. For the sake of simplicity, we calculate the average thickness 〈χ^dust〉 for the dust particle of mass m (radius a), which contains a parent nuclide X based on the following assumptions: (1) All the dust particles have the shape of a compact sphere, and their inner solid density ρ_s is constant. (2) The parent radionuclide X is uniformly distributed throughout the solid material of the dust particle (see also Section 5). (3) The direction of emitted energetic particle j is isotropic in solid angle around the decay place of X as a reference point. (4) The traveling path of the emitted particle j is straight, and the change in its direction due to the collisions with atomic electrons on the way can be neglected. The last assumption is not accurate enough, especially for Compton scattering of photons. However, we rest upon this assumption because emitted photons scatter at most a few times even in the largest dust particles considered in this paper.

We calculate the average traveling path length 〈s〉 of the emitted energetic particle in the solid material under these assumptions: we first average the path length over the emitted direction to obtain $\overline{s}$ at a fixed place of particle emission; then the length $\overline{s}$ is averaged over the distribution of the emitted place, which corresponds to the decay position of nuclide X (for further details, see Appendix B). The resultant length for the dust particle of radius a is 〈s〉 = (3/4)a. Since the average mass thickness is simply related to this length as 〈χ^dust〉 = ρ_s〈s〉, we have 〈χ^dust〉 ≈ (81ρ_s²m/256π)^1/3 for the dust particle of mass m.

The available energies (per particle) $\overline{E}(j)$ and the branching ratios R(j) of various energetic particles j with their types y(j) emitted due to the decay of parent nuclides ²⁶Al and ⁴⁰K are summarized in Table 1. We estimate the average energy loss in the solid material using the thickness 〈χ^dust〉. In the case of electron and positron (y(j) = e⁻ and e⁺) the residual kinetic energy $\overline{E}_\mathrm{res}(j)$ to ionize the gas on the surface of the dust particle can be obtained by integrating the mean rate of energy loss (or stopping power) over 〈χ^dust〉. The stopping power is essentially given by the Bethe–Bloch equation, which for the atomic component (absorber) i at the low-energy region (β ≲ 1) approximately has the form

$$\begin{equation} \left(\!-\frac{dE}{d\chi }\right)_i\! = 4\pi N_\mathrm{A} r_\mathrm{e}^2 m_\mathrm{e}c^2\frac{Z_i}{A_i}\frac{1}{\beta ^2}\!\left[\frac{1}{2}\ln \frac{2m_\mathrm{e}c^2\beta ^2\gamma ^2T_\mathrm{max}}{I_i^2}-\beta ^2\!\right], \end{equation} \tag{ 27 }$$

where N_A, r_e, and m_e are Avogadro's number, classical electron radius, and electron mass, respectively, Z_i, A_i, and I_i are atomic number, atomic mass, and mean excitation energy of component i, respectively, and T_max is the maximum kinetic energy that can be imparted to a free electron in a single collision. Stopping powers and the resultant ranges of electrons obtained from Equation (27) agree well with those calculated with a more sophisticated computer program ESTAR (Berger 1993). Concerning the chemical compositions of the dust particle, we use the fractions by weight f(i) summarized in Table 2. Consequently, the total energy-loss rate in the solid material is given by (− dE/dχ)^dust = ∑_if(i)(− dE/dχ)_i.

Table 1. The Average (Kinetic) Energies (Per Particle) $\overline{E}(j)$ and the Branching Ratios R(j) of the Particles j with Their Types y(j) Emitted by the Decay of Parent Nuclides X

Nuclide	Half-life	Abundance	Decay Mode	Type	Energy	Branching
X	t1/2(X) (yr)	x(X)		y(j)	$\overline{E}(j)$ (keV)	Ratio R(j)
26Al	7.4 × 105	(1.8–2.5) × 10−10a	β+ (82%)	e+	473	0.82
				γ	1808b	0.82
				2γ	512c	0.82
			EC (18%)	γ	2938b	0.0027
				γ	1808b	0.1773
				γ	1130b	0.0247
40K	1.28 × 109	2.2 × 10−10	β− (89.3%)	e−	534	0.893
			EC (10.7%)	γ	1461d	0.107

Notes. ^aThe relative abundance ²⁶Al/²⁷Al is taken to be (5–7) × 10⁻⁵. ^bThe energy of a photon emitted from the excitation state of ²⁶Mg. ^cThe energy of two photons produced by positron annihilation. ^dThe energy of a photon emitted from the excitation state of ⁴⁰Ar.

Download table as:  ASCIITypeset image

Table 2. Fractions by Weight f(i) for Atomic Components i in Dust Particles

Element	Chemical	Abundancea	Fraction	Solar
i	Substances	Ndust(i)	f(i)	Abundanceb
Si	SiO2	1.00 × 106	0.167	1.000 × 106
Mg	MgO	1.02 × 106	0.147	1.020 × 106
Fe	FeO, FeS	8.38 × 105	0.277	8.380 × 105
S	FeS	4.45 × 105	0.0848	4.449 × 105
O	SiO2, MgO, FeOc	3.41 × 106d	0.324	1.413 × 107

Notes. ^aThe relative abundance by number based on the assumption that each element (except O) is completely included in solid material of dust particles. ^bThe relative abundance by number (Lodders 2003). ^cAll Fe atoms other than those included in FeS are in the form of FeO. ^dAn amount of O in solid material is estimated on the assumption that dust particles consist mainly of the chemical substances in this table.

Download table as:  ASCIITypeset image

As for the photons (y(j) = γ) emitted from the excited states of the daughter nuclides, the residual energy $\overline{E}_\mathrm{res}(j)$ to ionize the gas on the surface of the dust particle is simply related to the emitted photon energy $\overline{E}(j)$ as

$$\begin{equation} \overline{E}_\mathrm{res}(j) \approx \overline{E}(j)\exp \left(-\frac{\langle \chi ^\mathrm{dust} \rangle }{\lambda ^\mathrm{dust}(\overline{E}(j))}\right), \end{equation} \tag{ 28 }$$

where λ^dust(E) is the effective energy-absorption length of the photon emitted with energy E. The length λ^dust(E) in the solid material can be described as

$$\begin{equation} \lambda ^\mathrm{dust}(E) = 1/\sum _i f(i)(\mu _\mathrm{en}/\rho)_i, \end{equation} \tag{ 29 }$$

where f(i) and (μ_en/ρ)_i are the fraction by weight in the solid material and the mass energy-absorption coefficients (in units of g⁻¹ cm²) for atomic component i, respectively. In calculating $\overline{E}_\mathrm{res}(j)$ from Equations (28) and (29), we also use f(i) summarized in Table 2 and the National Institute of Standards and Technology (NIST) data of X-ray mass attenuation coefficients as the coefficients (μ_en/ρ)_i (Hubbell & Seltzer 1996). The NIST data of attenuation coefficients take into account not only the secondary photon radiation produced at the initial photon–atomic interaction site but also the quanta of radiation from the annihilation of positrons originating in the initial pair production. These coefficients further involve the emission of radiation produced by the secondary charged particles traveling through the medium.

In the case of two photons emitted by the positron annihilation (y(j) = 2 γ) we must carefully estimate the residual energy (per positron) $\overline{E}_\mathrm{res}(2\,\gamma)$ of two respective photons (here we can denote j = 2 γ without confusion) available to ionize the gas on the surface of the dust particle. Because the rate coefficient of positron annihilation is much smaller than that of atomic electron scattering or ionization (see Equation (27)), it is reasonable to assume that the pair annihilation occurs when a positron comes to rest, i.e., its kinetic energy becomes zero. As long as the residual kinetic energy available to ionize the gas $\overline{E}_\mathrm{res}({e^+})$ is greater than zero, such a positron can leave the dust particle and the energy of photons emitted due to a subsequent pair annihilation in the gas phase is exactly given by $\overline{E}_\mathrm{res}(2\,\gamma)=2 \overline{E}(2\,\gamma) =2 \times 512\,\mathrm{keV}$. In case the kinetic energy $\overline{E}_\mathrm{res}(e^+)$ becomes zero on the way to the surface of the dust particle, a positron annihilation occurs on the inside of the dust particle and the emitted photons lose a certain amount of energy while they go through the solid material of this particle. We simply calculate the traveling path lengths of emitted photons under the following assumptions: (1) A pair annihilation occurs at the position where its distance r_pa from the center of the dust particle is equal to the range 〈χ^dust(e⁺)〉/ρ_s (=〈s(e⁺)〉) of the positron. (2) The directions of emitted photons are isotropic in solid angle around the spot of positron annihilation. (3) Two photons transit each other in the exactly opposite senses. With the average path lengths $\overline{s}$ of two emitted photons in the dust particle (for further details, see Appendix B), the residual energy per pair annihilation available to ionize the gas is given by

$$\begin{eqnarray} \overline{E}_\mathrm{res}(2\,\gamma) &\approx& \overline{E}(2\,\gamma)\left[\exp \left(-\frac{\overline{\chi }_\mathrm{long}}{\lambda ^\mathrm{dust}(\overline{E}(2\,\gamma))} \right)\right.\nonumber\\ &&\left.+\exp \left(-\frac{\overline{\chi }_\mathrm{short}}{\lambda ^\mathrm{dust}(\overline{E}(2\,\gamma))} \right) \right], \end{eqnarray} \tag{ 30 }$$

where $\overline{\chi }_\mathrm{long}$ ($=\rho _\mathrm{s}\overline{s}_\mathrm{long}$) and $\overline{\chi }_\mathrm{short}$ ($=\rho _\mathrm{s}\overline{s}_\mathrm{short}$) represent the mass thicknesses for the long and short traveling path lengths of emitted photons in the solid material, respectively, and $\lambda ^\mathrm{dust}(\overline{E}(2\,\gamma))$ is obtained by Equation (29) with $\overline{E}(2\,\gamma)=512\,\mathrm{keV}$.

We calculate numerically the residual energies available to ionize the gas $\overline{E}_\mathrm{res}(j)$ for various emitted particles j as functions of the dust radius a (or the corresponding mass m). The total residual energy for the decay of X is given by $\overline{E}_\mathrm{res}(\mathrm{X})=\sum _jR(j)\overline{E}_\mathrm{res}(\mathrm{X};\, j)$, where R(j) is the branching ratio of emitted particle j due to the decay of radionuclide X.

The results with the inner solid density ρ_s = 2 g cm⁻³ for the decays of nuclides X = ²⁶Al and ⁴⁰K are shown in Figures 1 and 2, respectively. The energies $\overline{E}_\mathrm{res}(\mathrm{X};\, j)$ of various emitted particles j are also represented by dashed and dot-dashed lines in these figures.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The results for ²⁶Al and ⁴⁰K resemble well each other because the types of emitted energetic particles are positron and photons, and electron and photon, respectively. As long as the radius of dust particle a is smaller than about 1 × 10⁻² cm, the kinetic energy losses of electron and positron due to the ionization of solid material are completely inefficient. Above this radius, which corresponds to the mass thickness 〈χ^dust〉 ≳ a few × 10⁻² g cm⁻², the energy loss of electron (positron) becomes efficient and the emitted electron (positron) stops within the dust particle of a ≳ 10⁻¹cm, or 〈χ^dust〉 ≳ (2–3) × 10⁻¹ g cm⁻². As for ⁴⁰K, the total residual energy $\overline{E}_\mathrm{res}({}^{40}\mathrm{K})$ decreases abruptly above this radius because most of the energy available to ionize the gas is carried by the emitted electron.

The absorption of photons emitted from the excited states of daughter nuclides and that of two photons emitted due to the positron annihilation become efficient only at a ≳ 10 cm, or 〈χ^dust〉 ≳ (2 – 3) × 10 g cm⁻², because the effective energy-absorption length of a photon with E ≳ 1 MeV is greater than 20 g cm⁻² for the chemical compositions of the dust particle. Thus, the total residual energies $\overline{E}_\mathrm{res}(\mathrm{X})$ for both nuclides remain constant values as long as a ≲ 10 cm. The residual energies become abruptly close to zero for dust particles of a ≳ several × 10 cm, or 〈χ^dust〉 ≳ 10² g cm⁻², since the emitted energy of photons is at most 3 MeV. In conclusion, the size growth of the dust particle affects the total residual energy $\overline{E}_\mathrm{res}({}^{26}\mathrm{Al})$ moderately as long as a ≲ 10 cm.

3.3. Shielding by Other Dust Particles

When the settling of dust particles proceeds extensively, the mass density of dust particles ρ^dust around the midplane becomes close to or exceeds that of gas ρ^gas (see Section 4). In such a situation, the energetic particles such as α-particles, electrons, positrons, and photons that have emerged from the surface of a certain dust particle collide with other dust particles located nearby and lose considerable amounts of the (kinetic) energies before they ionize the gas particles sited in the same place.

Because the column density of the disk where this effect becomes efficient is much larger than the ranges of the energetic particles, we can regard that these particles deposit their energies at the same place where they leave from the surfaces of dust particles. The part of the deposit energy used for ionizing the gas depends on the densities and the energy-loss rates of the gas and dust particles at that site.

Charged particles with energies less than a few MeV lose their kinetic energies primarily through the ionization and atomic excitation. We can see from Equation (27) that the mean energy-loss rate is essentially proportional to the ratio Z_i/A_i of absorber i and that the variation of loss rate due to the atomic property of the absorber can be insignificant because the rate depends on the logarithm of its excitation energy I_i. As for photons with energies less than a few MeV the primary energy-loss process is incoherent scattering (Compton scattering off an electron). The atomic photoelectric effect, which is characterized by discontinuities due to thresholds for photoionization of various atomic levels, becomes predominant in the lower energy region of E ≲ 100 keV. Photon attenuation length for E ≳ 100 keV depends slightly on the atomic property of absorber i. Thus, we can safely regard that the mass attenuation coefficient of the emitted photon is proportional to the number of electrons Z_i in the atom (absorber). The effective mass attenuation coefficient is approximately given by the sum of the coefficients that are proportional to the ratio Z_i/A_i weighted with respect to the fraction by weight f(i) for atomic component i in the chemical compound.

In such a situation the portion of the (kinetic) energy used for ionizing the gas particles sited in a certain place is given by

$$\begin{equation} p^\mathrm{gas} = \frac{\alpha ^\mathrm{gas}}{\alpha ^\mathrm{gas}+\alpha ^\mathrm{dust}} \approx \frac{(Z_i/A_i)^\mathrm{gas}\rho ^\mathrm{gas}}{(Z_i/A_i)^\mathrm{gas}\rho ^\mathrm{gas}+(Z_i/A_i)^\mathrm{dust}\rho ^\mathrm{dust}}, \end{equation} \tag{ 31 }$$

where α^gas and α^dust are the energy-absorption probabilities by gas and dust particles in this place, which are approximately proportional to products of the partial mass densities and the average ratios Z_i/A_i of these particles. Absorbers consist mainly of elements H and He for gas particles and of heavy elements summarized in Table 2 for dust particles. We can easily calculate the average ratios and use (Z_i/A_i)^gas = 0.86 and (Z_i/A_i)^dust = 0.49 in the following calculations.

We solve numerically the discrete equation for size growth (coalescence) of settling dust particles (Equation (20)) at a typical orbit (R = constant) of the disk and obtain the dust size distribution and the change of the gas ionization rates by the decay of radionuclides as functions of time t and the height z from the midplane. As for the physical conditions and/or the evolutionary stages of the protoplanetary disk, we consider primarily a quiescent disk model (evolutionary stage) in this paper. The effects of the (global) turbulence on the ionization rates by radionuclides will be investigated in Paper II.

Energetic particles emitted by the decay of both long-lived radionuclide ⁴⁰K and short-lived one ²⁶Al can penetrate the neutral gas (composed of H₂ molecules and He atoms) of thickness χ_R (∼10 g cm⁻²; see also Section 3). Thus, we only consider the inner region where the surface density Σ(R) exceeds about 100 g cm⁻² and the emitted particles deposit their energies into the medium located in the same place as the parent nuclides. Moreover, the ionization by radionuclides becomes significant only in the zone where the ionization rate by cosmic rays decreases to a certain extent (Umebayashi & Nakano 1981), and the X-ray ionization from the central star must be unimportant. The ionization by low-energy X-rays becomes efficient only at the surface layer whose column density is below ∼10 g cm⁻² (Igea & Glassgold 1999). Because hard X-ray emission from YSOs (e.g., Preibisch et al. 2005; Getman et al. 2008) can penetrate and ionize the gas above this column density, the ionization by radionuclides can be efficient only in the inner region with surface density ≫10 g cm⁻². In the MMSN, the region inside Jupiter's orbit (R = 5.2 AU) satisfies this condition.

Furthermore, we restrict ourselves to investigating the ionization rate by radionuclides in the zone from the midplane to the scale height H(R) within which most of the gas and dust particles for a given R exist. Outside this zone the emitted particles can escape freely from the disk to the surrounding interstellar space without ionizing the gas of the disk.

UN09 calculated the ionization rates by the long- and short-lived nuclides based on the following assumptions, which were not mentioned explicitly: (1) The dust particles that contain these nuclides are mixed so well with the gas that relative abundance of dust particles to gas remains unchanged everywhere in the disk. (2) The size of dust particles is small enough for emitted energetic particles to move freely from the solid material of dust particle to the surrounding space, and consequently the loss of the available energies to ionize the gas is negligible. We perform numerical calculations of the dust evolution under the initial conditions similar to these assumptions: dust particles are initially mixed well with the gas and their abundance relative to gas is homogeneous throughout the vertical zones below z = H(R). As for the initial dust radius we use the same initial value a_init = 10⁻⁴ cm, independent of the orbital radius R and the factor $f_\Sigma$. The same initial radius for various regions of the disks is different from those used by Nakagawa et al. (1986) and by Nomura & Nakagawa (2006) simply because we clarify the evolutionary changes of the ionization rates at different orbits under the same conditions about the initial model and abundance of dust particles. These conditions affect little the following results obtained numerically. Moreover, we simply assume the sticking probability of mutual dust collisions p_s = 1, which affects primarily the evolutionary timescale of dust particles.

In the numerical calculations the spatial grid in the vertical direction is taken to be equally spaced into 20 intervals within 0 ⩽ z ⩽ H(R), which is the same as that used by Nomura & Nakagawa (2006). On the other hand, we do not take additional logarithmically spaced subintervals within the lowest interval 0 ⩽ z ⩽ H(R)/20 because we need not resolve the dust evolution near the midplane of the disk in further details. The dust mass is binned into n intervals logarithmically as m_{i + 1/2} = m_{i − 1/2} with n = 160 and = 2^1/2, unlike Nomura & Nakagawa (2006), who adopted n = 320 and = 2^1/4. This change causes no essential difference in the dust evolution as described below.

We can safely regard that the emitted particles deposit most of their energies in such an interval that the decay of radionuclides has occurred. The matter thickness of a spatial interval in the vertical direction is about 4 g cm⁻² at the orbit R = 5 AU of the MMSN. This thickness is a few times smaller than the range of energetic particles χ_R. Because the particles emitted isotropically from the decay place of the parent nuclides, half of the emitted particles run inward. The matter thickness of the spatial intervals in the radial and azimuthal directions is considered to be close to the surface density for a given R. Moreover, some of emitted particles in the innermost interval go through the midplane and deposit their energies in its counterpart on the opposite side. Thus, the matter thickness of an interval penetrated by the energetic particles is generally much larger than that in the vertical direction. A notable exception is the outermost interval where the emitted particles can escape upward from the disk. Although the disk matter does not end abruptly around the position z ≈ H(R), the ionization rate in this part of the disk with Σ(R) ≈ 100 g cm⁻² is treated as an upper limit to the actual values.

The motion and distribution of dust particles are characterized by typical heights Z_k (k = 1, ..., 4) from the midplane (Nakagawa et al. 1986; Nomura & Nakagawa 2006). The spatial grid in the numerical calculations can resolve clearly the zone z ≳ Z₁ = 2ηR where dust particles move mainly in the vertical direction (V_z ≫ V_R). Thus, we can simply assume the velocity difference in the R-direction δV_R = 0 in calculating the relative collision velocity of dust particles δV from Equation (24). Because the position Z₁ is close to the upper boundary of the innermost zone H(R)/20, the radial velocity V_R is not absolutely small. However, this assumption is not worse for calculating the size distribution around the midplane of the disk in a quiescent state, since the velocity difference δV_R is found to be so small that it affects little the relative velocity δV.

In the zone $z \gtrsim Z_2=(\sqrt{\pi }/2)(\Sigma _0^\mathrm{dust}/\Sigma _0)H(R)$, the gas density ρ^gas is always larger than the dust density ρ^dust. The radial velocity of dust particles in this zone reaches its maximum value V_R ≈ ηv_K (see Equation (14)) when the dust radius becomes above a critical value a_cr. Setting the friction time between gas and dust particles τ_f ≈ 1/D^gas equal to the characteristic time of the radial migration τ_mig ≈ 1/Ω_K(R), the critical radius of dust particles is roughly given by

$$\begin{equation} a_\mathrm{cr} = \left\lbrace \begin{array}{@{}ll}c_\mathrm{th}\rho ^\mathrm{gas}/\rho _\mathrm{s}\Omega _\mathrm{K}(R) &\quad \mbox{for } a_\mathrm{cr} \lesssim l^\mathrm{gas} \\ {[3c_\mathrm{th}\rho ^\mathrm{gas}l^\mathrm{gas}/2\rho _\mathrm{s}\Omega _\mathrm{K}(R)]}^{1/2} &\quad \mbox{for } a_\mathrm{cr} \gtrsim l^\mathrm{gas}. \end{array} \right. \end{equation} \tag{ 32 }$$

In the inner region of the disk a_cr is generally given by the latter formula, i.e., dust particles can migrate efficiently after they grow larger than the mean free path of gas particles (a ≳ l^gas). Because dust particles with a ≳ a_cr migrate inward very rapidly with the timescale R/V_R ≈ t_K/(2πη), we artificially remove from the numerical calculations the dust particles whose radius reaches above a_cr.

Although the dust fluid whose velocity is an average over a mass of dust particles does not oscillate around the midplane, Equation (11) shows that each individual dust particle begins to oscillate if D^gas = Aρ^gas becomes smaller than 2Ω_K(R) (Nakagawa et al. 1986). Setting the friction time τ_f ≈ 1/D^gas equal to the characteristic time of the oscillation τ_osc ≈ 1/2Ω_K(R), the critical radius for oscillations a_osc is roughly given by either a_cr/2 for a_osc ≲ l^gas or $a_\mathrm{cr}/\sqrt{2}$ for a_osc ≳ l^gas. In a quiescent disk, however, dust particles always settle toward the midplane in the mean motion (Nomura & Nakagawa 2006). The oscillation of a large particle that satisfies a > a_osc is a fluctuation of motion that is dampened with a timescale of 2τ_f, while a smaller particle (a < a_osc) settles without oscillation. Thus, we neglect the oscillation of large dust particles around the midplane of the disk in numerical calculations.

4.1. A Fiducial Model: The Minimum Mass Solar Nebula

First, we investigate the evolution of dust particles and its effects on the ionization rates by radionuclides in the MMSN at a quiescent stage: $f_\Sigma =1$, q = 3/2, and p = 1/2. Because the innermost orbital radius R of the solar nebula is about 0.35 AU (Hayashi 1981) where the temperature becomes close to 500 K but the thermal ionization remains inefficient (Umebayashi 1983), we investigate the dust evolution and its effects on the ionization rates in the region of 0.35 AU ≲ R ≲ 5 AU.

4.1.1. Dust Size Distribution

Figures 3, 4, and 5 show the numerical results for the size distribution of the mass density of dust particles normalized by the total mass density at its initial stage (t = 0), φ(a, z)/ρ^dust, for R = 1, 0.5, and 5 AU, respectively. In these figures, we represent the size distribution as a function of the dust radius a (or the corresponding dust mass m) instead of the numerical bin numbered by i (see Equation (20)) and the normalized height from the midplane z/H(R). Panels (a)–(d) of Figure 3 show the results at four characteristic stages of the dust evolution: (1) the early stage that the size growth of dust particles proceeds at their initial places, and the settling occurs little in most of the vertical zones (t = 3.2 × 10² yr); (2) the stage that both the settling and the size growth proceed extensively (t = 1 × 10³ yr); (3) the stage that most of the dust particles settle in the zone z ≲ Z₁, where the radial migration of dust particles becomes efficient (t = 3.2 × 10³ yr); (4) the later stage that the dust layer can already be gravitationally unstable to fragment into planetesimals if the disk continues to be in a quiescent state (t = 1 × 10⁴ yr). In Figures 4 and 5, the size distributions at typical times corresponding to stages (2) and (3) are shown in panels (a) and (b), respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

At stage (1) the size growth of dust particles occurs faster in the zones of larger z/H(R) and the settling from their initial location remains inefficient, as can be seen from Figure 3(a). Because the terminal velocity of dust particles in the vertical direction V_z is proportional to the height z from the midplane (see Equation (15)), the relative collision velocity δV of dust particles in the upper part becomes larger than that in the lower part. Larger dust particles that have grown around the surface layer of the disk settle more rapidly toward the midplane and grow further in size. This tendency leads to bimodal size distribution in the zone close to the midplane at the later stages (2)–(4) (see, e.g., Weidenschilling 1997; Nomura & Nakagawa 2006). Above the innermost zone, on the other hand, single-peaked and almost flat distribution is achieved owing to this fast settling of larger dust particles toward the midplane. At stages (3) and (4), the local minimum between the peaks of bimodal distribution around the midplane becomes shallower and amounts of floating dust particles in the residual zones continue to decrease due to the settling of smaller particles as time goes by.

The numerical results for R = 1 and 5 AU agree well with those obtained by Nomura & Nakagawa (2006). Because we use the interpolation formula for the coefficient A of the gas drag force (see Equation (12)) in numerical calculations, no gaps in size distribution appear around a ≈ l^gas where A begins to follow Stokes' law, rather than Epstein's law. We note that the initial radius of dust particles a_init = 10⁻⁴ cm affects only the early evolutionary stages where the size growth of dust particles in their initial positions is efficient. Thus, the times to reach stages (2)–(4) for R = 5 AU become a few thousand years later than those obtained by Nomura & Nakagawa (2006).

We find the critical dust radii a_cr = 21, 32, and 97 cm at the orbits of R = 0.5, 1, and 5 AU, respectively. Because dust particles with a ≳ a_cr migrate inward very rapidly with the timescales R/V_R ≈ 44, 88, and 4.4 × 10² yr for R = 0.5, 1, and 5 AU (e.g., Adachi et al. 1976), we artificially remove the dust particles whose radii reach above a_cr from the numerical calculations. Some dust particles exceed the critical values, but total mass fractions Σ^dust_rm(R)/Σ^dust(R) of the removed particles (a > a_cr) are negligible even at the later stages (see Table 3). Thus, the removal of larger dust particles affects little the numerical results of the size distribution and the ionization rates by radionuclides.

Table 3. The Maximum Ratios ζ_max(X)/ζ_init(X) of the Ionization Rates by the Decay of Radionuclides X around the Midplane at the Orbits R of the Protoplanetary Disks with the Surface Density Factor $f_\Sigma$

R	(AU)	0.5	1	5	10	30
$H(R)/(Z_\mathrm{R}f_\Sigma)$		240	85	7.6	2.7	0.52
$\Sigma ^\mathrm{dust}(R)/f_\Sigma$	(g cm−2)	20	7.1	2.7	0.95	0.18
$a_\mathrm{s}/f_\Sigma$	(cm)	1.3	0.44	0.34	0.12	0.023
(a)   $f_\Sigma =1$	(MMSN)
Σdustrm(R)/Σdust(R)		0.031	2 × 10−8	2 × 10−17
$E_\mathrm{av}({}^{26}\mathrm{Al})/\overline{E}_\mathrm{em}({}^{26}\mathrm{Al})$		0.66	0.73	0.81
$E_\mathrm{av}({}^{40}\mathrm{K})/\overline{E}_\mathrm{em}({}^{40}\mathrm{K})$		0.23	0.26	0.32
pgas		0.64	0.83	0.93
ζmax(26Al)/ζinit(26Al)		100	51	14a
ζmax(40K)/ζinit(40K)		36	19	5.4a
ζlim/ζinit		150	70	19a
(b)   $f_\Sigma =3$
Σdustrm(R)/Σdust(R)		0.34	0.012	0.037	4 × 10−18
$E_\mathrm{av}({}^{26}\mathrm{Al})/\overline{E}_\mathrm{em}({}^{26}\mathrm{Al})$		0.42	0.60	0.58	0.81
$E_\mathrm{av}({}^{40}\mathrm{K})/\overline{E}_\mathrm{em}({}^{40}\mathrm{K})$		0.16	0.21	0.22	0.32
pgas		0.47	0.63	0.82	0.92
ζmax(26Al)/ζinit(26Al)		140	95	11	13a
ζmax(40K)/ζinit(40K)		53	34	4.1	5.3a
ζlim/ζinit		220	160	18	18a
(c)   $f_\Sigma =10$
Σdustrm(R)/Σdust(R)		0.63	0.28	0.32	0.028	8 × 10−24
$E_\mathrm{av}({}^{26}\mathrm{Al})/\overline{E}_\mathrm{em}({}^{26}\mathrm{Al})$		0.23	0.38	0.30	0.53	0.84
$E_\mathrm{av}({}^{40}\mathrm{K})/\overline{E}_\mathrm{em}({}^{40}\mathrm{K})$		0.10	0.15	0.13	0.21	0.35
pgas		0.32	0.41	0.66	0.79	0.95
ζmax(26Al)/ζinit(26Al)		180	130	15	11	14a
ζmax(40K)/ζinit(40K)		79	51	6.5	4.4	5.8a
ζlim/ζinit		280	250	34	20	19a

Note. ^aThe ratio in the innermost zone 0 ⩽ z ⩽ H(R)/20 obtained from the numerical calculation.

Download table as:  ASCIITypeset image

It should also be noted that shear motion between the dust layer and the gas induces a local turbulence near the midplane as the settling of dust particles proceeds extensively (see, e.g., Cuzzi et al. 1993; Weidenschilling & Cuzzi 1993; Weidenschilling 1997). A quiescent state cannot maintain its ideal situation of non-turbulence. Computing the Richardson number, Nomura & Nakagawa (2006) found that the shear instability occurs in the zone z ≲ Z₂ at each orbit. Dust particles do not concentrate in the midplane and migrate inward very rapidly once they grow large enough, a ≳ a_cr.

4.1.2. Ionization Rates

Using the size distribution of dust particles described above, we calculate the evolution of the ionization rates by the decay of either short-lived nuclide ²⁶Al or long-lived one ⁴⁰K. The ionization rates at some typical evolutionary stages normalized by their initial values ζ(X)/ζ_init(X) at the orbits R = 1, 0.5, and 5 AU are shown in Figures 6, 7, and 8, respectively, as functions of the relative height z/H(R) from the midplane of the disk. Panels (a) and (b) of Figure 6 show the rates by the decay of nuclides X = ²⁶Al and ⁴⁰K. We only represent the results of ²⁶Al at R = 0.5 (Figure 7) and 5 AU (Figure 8). As can be seen from these figures, the characteristics of the evolution of the rates at different orbits R resemble each other. Thus, we can easily estimate from these figures how the ratios ζ(⁴⁰K)/ζ_init(⁴⁰K) at R = 0.5 and 5 AU change gradually.

Zoom In Zoom Out Reset image size

At the early stage (1) that the size growth of dust particles proceeds at their initial places, the ionization rate remains unchanged except for the uppermost zone where the settling is most efficient. The size growth of dust particles scarcely affects the ionization rate because most of the dust particles are still so small that the energetic particles emitted due to the decay of radionuclides can freely penetrate the solid material of dust particles and the loss of their (kinetic) energies is negligible (see Figures 1 and 2). As the settling proceeds, the rate in the innermost zone increases considerably, that in the intermediate zone continues to be a constant value close to its initial rate, and at the outer zone near the scale height H(R) the ratio decreases noticeably. At the later stage that most of the dust particles have settled in the innermost zone, the ratio ζ(X)/ζ_init(X) in this zone reaches a large value and remains unchanged afterward. The ratio at the outer zone, however, keeps on decreasing in proportion to the amount of the floating dust particles as time goes by.

The energy available to ionize the gas E_av(X), which is given by the residual energy $\overline{E}_\mathrm{res}(\mathrm{X})$ averaged over the size distribution φ(a, z) at the height z, is essential to calculate the ratio ζ(X)/ζ_init(X). Both the energy loss due to the size growth of dust particles and the amount of particles removed from the calculation decrease the value of E_av(X), although the latter affects little the results in the present situation (MMSN). As for the ionization rate by the decay of ²⁶Al, the size growth of dust particles decreases the energy E_av(²⁶Al) only by a factor of about 0.7 in the innermost zone (see also Table 3). This can be easily understood from Figure 1 that the radii of the most abundant dust particles in Figures 3, 4, and 5 are still smaller than the radius at which the energy loss of emitted energetic particles becomes efficient. Because most of the dust material is settled into the innermost zones of thickness Δz ≈ H(R)/20, the rates ζ(²⁶Al) in these zones are about 14 times greater than their initial values.

The size growth of dust particles can be efficient to decrease the ionization rate by the decay of ⁴⁰K. Thus, the rate ζ(⁴⁰K) in the innermost zone for R = 1 AU at the later stage becomes only several (≈5) times larger than its initial value. We can see from Figure 2 that the energy available to ionize the gas $\overline{E}_\mathrm{res}({}^{40}\mathrm{K})$ decreases to about 1/3 times the total emitted energy for dust particles with a ≳ 0.1 cm. The residual energy that is carried by photons is almost constant up to a ≈ 10 cm. Since this value is larger than the radius of the most abundant particles at the later stage, the ratio ζ(⁴⁰K)/ζ_init(⁴⁰K) in the innermost zone remains unchanged. The ratios in the outer zones, on the other hand, change in a similar way to those for ²⁶Al as the settling proceeds. Because most of the dust particles floating in these zones are smaller than 10⁻² cm as can be seen from panels (c) and (d) of Figure 3, the energy loss in the solid material is inefficient. Thus, the ratios ζ(X)/ζ_init(X) that are proportional to the amount of dust particles take the same value for both nuclides X = ²⁶Al and ⁴⁰K (Figures 6(a) and (b)). Similar results are obtained in the case of R = 0.5 and 5 AU, and the ratios by the decay of ⁴⁰K in the innermost zones are about 4.5 and 5.4, respectively. The difference in the ratios is partly due to the inner solid density ρ_s of dust particles.

The concentration of dust particles in the innermost zones and the resultant absorption (shielding) of emitted energetic particles by other dust particles located nearby are ineffective in decreasing the ionization rates. Because the relative abundance ρ^dust/ρ^gas becomes at most 0.08 for R ≲ 2.7 AU and 0.35 for R ≳ 2.7 AU, the shielding decreases the large ratios shown in Figures 6, 7, and 8 only by 5%, 5%, and 20%, respectively (see Equation (31)).

As for the outer orbit R = 5 AU, the surface density within H(R)/20 is about 4 g cm⁻², which is already smaller than the stopping range of the emitted particles χ_R. Simply the rates in Figure 8 decrease by a factor of about 2.5, i.e., the ratio of the range of energy deposition to the surface thickness. However, this factor must decrease and approach close to one because the energetic particles are emitted isotropically from the decay place of the parent nuclides, as mentioned earlier. This indicates that the ratios ζ(²⁶Al)/ζ_init(²⁶Al) shown in Figure 8 are close to the precise values both in the innermost zone and in the outer zones. The ratio can be overestimated only in the outermost zone z ≈ H(R). Moreover, the results of the ionization rates by ⁴⁰K are much the same.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.2. Disks with Various Surface Densities

We investigate quiescent disk models with various surface densities rotating around a star of M_⋆ = 1 M_☉. First, we consider the model of different surface density factor $f_\Sigma$ and the other parameters are the same as those of the MMSN. Panels (a) and (b) of Figure 9 show the numerical results of the ratios ζ(²⁶Al)/ζ_init(²⁶Al) for $f_\Sigma = 10$ at the orbits R = 1 and 30 AU, respectively. It should be noted that the disk of $f_\Sigma \lesssim 10$ is gravitationally stable in the orbital region R ≲ 100 AU (see Equation (9)) and that the effect of dust evolution on the ionization rates by radionuclides becomes significant in the region of Σ(R) ≳ 100 g cm⁻².

Zoom In Zoom Out Reset image size

Characteristics of the results are broadly similar to those in the MMSN. (1) A significant increase of the ionization rate occurs only in the innermost zone around the midplane. (2) In the intermediate zones and the outer zones, the ionization by radionuclides becomes completely inefficient as the settling of dust particles proceeds.

The ratio ζ(²⁶Al)/ζ_init(²⁶Al) for R = 1 AU in the numerical calculation reaches up to 7.4, which is about half of the ratio in the innermost zone of the MMSN. Because there exists sufficient dust material in the disk of large $f_\Sigma$, many dust particles grow large and some exceed the critical radius a_cr. Note that a_cr is independent of $f_\Sigma$ as long as a_cr ≳ l^gas (see Equation (32)). Consequently, about 30% of dust material is removed from the calculation. The energy loss of emitted particles in the solid material also increases as the fraction of large dust particles increases and the relative energy available to ionize the gas $E_\mathrm{av}({}^{26}\mathrm{Al})/\overline{E}_\mathrm{em}({}^{26}\mathrm{Al})$ decreases to about 4.

Similar results are obtained in the inner orbital region R ≲ 13 AU of the disk with $f_\Sigma = 10$. The condition a_cr ≳ l^gas is satisfied in this region, and a_cr increases moderately as R increases. The ratios ζ(²⁶Al)/ζ_init(²⁶Al) in the innermost zones become about 4.6, 5.3, and 8.9 for R = 0.5, 5, and 10 AU, respectively. These values depend mainly on the energies available to ionize the gas E_av(²⁶Al), which are affected by total mass fraction of the removed dust particles Σ^dust_rm(R)/Σ^dust(R) and/or by the energy loss in the solid material. For R = 0.5 and 5 AU considerable amounts (∼60% and 30%) of dust material are removed from the calculations, and as a result E_av(²⁶Al) decreases significantly. At R ≈ 10 AU, on the other hand, a fraction of removed dust material becomes small (≈3%). Because dust particles can grow close to the critical radius a_cr ≈ 190 cm (>l^gas), the energy loss in the solid material becomes the major process and E_av(²⁶Al) decreases by about 50%. As for the case of long-lived nuclide ⁴⁰K, situations are similar and the ratios in the innermost zones decrease to about 2.0, 2.9, 2.3, and 3.5 for R = 0.5, 1, 5, and 10 AU, respectively. The shielding by other particles decreases the ratios of both nuclides by about 2%, 3%, 12%, and 16% at the same orbits R. The reduction of the factor p^gas in the inner region R ≲ 5 AU is owing to a large amount of removed dust material.

At the orbit R = 30 AU where the surface density Σ(R) is close to 100 g cm⁻², we find numerically the large ratios ζ(X)/ζ_init(X) ≈ 14 and 5.8 for X = ²⁶Al and ⁴⁰K around the midplane (see Figure 9(b)). Because the dust particles scarcely grow above the critical value a_cr = 66 cm (<l^gas), a fraction of removed dust particles Σ^dust_rm(R)/Σ^dust(R) is negligible (∼8 × 10⁻²⁴) and the available energy E_av(X) decreases moderately. This situation resembles that in the inner orbital region (R ≲ 1 AU) of the MMSN. Since the surface density within H(R)/20 is about 2.5 g cm⁻², the real ratios ζ(X)/ζ_init(X) can decrease at most by a factor of four. However, for the same reason as mentioned earlier, both ratios are close to their precise values.

In the outer orbital region R ≳ 13 AU where a_cr becomes smaller than l^gas, the critical radius a_cr, which is proportional to $f_\Sigma$ (or the gas density ρ^gas), decreases as R increases. Because the amount of dust material in this region is insufficient to make massive dust particles (a ≳ a_cr), characteristics of the results closely resemble those in the inner orbital region of the MMSN. As R increases, a fraction of removed dust material becomes insignificant. The energy available to ionize the gas E_av(X) increases as the radius of the most abundant dust particles decreases gradually. Thus, the ratios ζ(X)/ζ_init(X) approach the values of about 14 and 6 for X = ²⁶Al and ⁴⁰K obtained in the outermost region R = 30 AU.

We also investigate the disk model of $f_\Sigma =3$ numerically and find much the same characteristics of results. Quantitatively, the ratios for ²⁶Al in the innermost zones are found to be about 8.2, 11, 9.6, and 13 for R = 0.5, 1, 5, and 10 AU, respectively, and those for ⁴⁰K are 3.1, 4.1, 3.7, and 5.3 at the same orbits. Because the amount of dust material in the disk is proportional to $f_\Sigma$, a small portion of dust particles can grow large enough to exceed the critical radius a_cr in this model. Thus, a fraction Σ^dust_rm(R)/Σ^dust(R) becomes smaller and E_av(X) increases, compared with the disk for $f_\Sigma =10$ at the same orbit. These situations adequately change both ratios to the intermediate values.

The ratios ζ(X)/ζ_init(X) in the intermediate zones and the outer zones at a given R are almost independent of the surface density of the disk model, i.e., $f_\Sigma$ and R. Because the ratios are directly proportional to the amounts of dust particles in these zones at the later stages, this tendency becomes obvious as time goes by.

Furthermore, we investigate the disk models with different power index q of the surface density and find that characteristics of results are similar among all the disk models. The surface density Σ(R) at a given R is essential for an increase in the ratio ζ(X)/ζ_init(X) around the midplane. As long as Σ(R) is close to that in the MMSN, a fraction of removed dust material is small and the settling of dust particles increases the ionization rates by ²⁶Al and ⁴⁰K extensively. The largest ratios are about 14 and 5, almost independent of the radial distance R. In the region that Σ(R) is large enough to make massive dust particles with a ≳ a_cr, the increase in the ratio is moderate because the reduction of dust material due to the radial migration toward the central star becomes significant. The condition that a_cr ≳ l^gas is satisfied or not affects little the results summarized above.

4.3. The Maximum Ionization Rates around the Disk Midplane

The large ionization rates by radionuclides in the innermost zones obtained numerically are restricted to the thickness of interval Δz ≈ H(R)/20 used in the calculations. In order to estimate the maximum ratios ζ_max(X)/ζ_init(X) in the thin gaseous layer around the midplane, we must consider additional settling of dust particles toward the midplane.

The minimum thickness of the gaseous layer into which the energetic particles emitted due to the decay of nuclide X can deposit their (kinetic) energies is given by χ_R. This thickness corresponds to the layer of the height Z_R ≈ [2χ_R/Σ(R)] H(R) from the midplane. As long as $f_\Sigma \lesssim 10$, the boundary Z_R is not far away from the typical height Z₂, and is located between Z₁ and Z₂ for a given R in most of the disk models considered in the previous subsections. It should be noted that both heights Z₁ and Z₂ are independent of $f_\Sigma$. As the settling in the innermost zone proceeds further, most of the dust material is concentrated in the thin layer of thickness Δz (≳ Z_R). The amount of dust material in this layer at the later stage is at most H(R)/Δz times larger than its initial value. Even if dust particles settle further into a very thin region close to the midplane (Δz ≪ Z_R), energetic particles emitted from such a region deposit their energies into the gaseous layer of the matter thickness χ_R.

Nomura & Nakagawa (2006) calculated additional settling and size growth of dust particles at R = 1 and 5.2 AU in a quiescent disk (MMSN) numerically. Their results show that in the inner part where most of the dust particles are concentrated the size distribution of dust particles has essentially a similar form and the radius of the most abundant particles remains unchanged. This part becomes thinner as the settling toward the midplane proceeds further. Features of the size distribution in the innermost zone obtained from our calculations resemble their results. Only a minor difference is the radius of the most abundant dust particles. Because dust particles can grow up to the critical radius a_cr in our calculations, the radius of the most abundant particles is several times larger than that obtained by Nomura & Nakagawa (2006).

We use a fraction of energy available to ionize the gas $E_\mathrm{av}(\mathrm{X})/\overline{E}_\mathrm{em}(\mathrm{X})$ obtained numerically in the innermost zone for estimating the maximum ratios ζ_max(X)/ζ_init(X) in the thin gaseous layer. Both the energy loss of emitted energetic particles and the reduction of dust material due to the radial migration of massive dust particles (a ≳ a_cr) have effects on E_av(X). Since additional settling changes little the size distribution of the dust particles in the inner layer where most of the dust material is concentrated, the fraction of dust particles removed in the course of further settling is small. The radius of the most abundant particles is at most comparable to the length above which the energy loss in the solid material becomes abruptly effective, and the energy loss in smaller dust particles is roughly constant (see Figures 1 and  2). Thus, we can safely consider that further settling affects little the energy loss of emitted particles and that a fraction $E_\mathrm{av}(\mathrm{X})/\overline{E}_\mathrm{em}(\mathrm{X})$ in the thin layer is approximately given by that found in the numerical calculations.

As dust particles settle further toward the midplane, the abundance ratio ρ^dust/ρ^gas increases gradually and the shielding by other particles located nearby becomes efficient. Because most of the dust particles remaining in the numerical calculations can be settled into the thin gaseous layer of thickness Δz, the maximum ratio of ρ^dust/ρ^gas is virtually given by [Σ^dust(R) − Σ^dust_rm(R)]/2χ_R. The portion of the (kinetic) energy to ionize the gas p^gas is obtained from Equation (31) with this ρ^dust/ρ^gas. Consequently, the maximum ratio of the ionization rate by the decay of X around the midplane is given by

$$\begin{equation} \frac{\zeta _\mathrm{max}(\mathrm{X})}{\zeta _\mathrm{init}(\mathrm{X})} = \frac{H(R)}{\Delta z}\frac{E_\mathrm{av}(\mathrm{X})}{\overline{E}_\mathrm{em}(\mathrm{X})}p^\mathrm{gas} \approx \frac{\Sigma (R)}{2 \chi _\mathrm{R}}\frac{E_\mathrm{av}(\mathrm{X})}{\overline{E}_\mathrm{em}(\mathrm{X})}p^\mathrm{gas}, \end{equation} \tag{ 33 }$$

where we take Δz = Z_R in the last expression.

Using $E_\mathrm{av}(\mathrm{X})/\overline{E}_\mathrm{em}(\mathrm{X})$ and Σ^dust_rm(R)/Σ^dust(R) obtained numerically, we calculate the maximum ratios ζ_max(X)/ζ_init(X) by the decay of nuclides X = ²⁶Al and ⁴⁰K. The results with Δz = Z_R for the disk models of $f_\Sigma =1$, 3, and 10 are summarized in Table 3. In the outermost orbital region Σ(R) ≲ 200 g cm⁻² the surface density of interval Δz ≈ H(R)/20 becomes smaller than the stopping range of the emitted particles χ_R. Because the increase of the ionization rates around the midplane due to the additional settling of dust particles scarcely occurs in this region, we represent the results of the numerical calculations in the columns of the outermost orbit R for a given $f_\Sigma$, which are close to the maximum ionization rates, as mentioned earlier.

In the case of (a) $f_\Sigma =1$ (MMSN) the maximum ratios ζ_max(²⁶Al)/ζ_init(²⁶Al) reach up to about 100 and 50 at the orbits R = 0.5 and 1 AU, respectively. The surface thickness χ_R corresponds to the layers of heights Z_R ≈ H(R)/240 and H(R)/85. Because fractions of removed dust particles Σ^dust_rm(R)/Σ^dust(R) are minor (≲ 0.03), dust material in these layers increases extensively. Both the energy loss due to the size growth of dust particles and the shielding by other particles become efficient in this situation. Thus, the maximum ratios increase moderately, compared with large amounts of dust material in the thin layers. The resultant ionization rates ζ_max(²⁶Al) are close to or larger than the ionization rate by energetic particles (cosmic rays) in the ISM (ζ_CR ≈ (3–10) × 10⁻¹⁷ s⁻¹; van der Tak & van Dishoeck 2000; van der Tak et al. 2006), where the damping of the flux is negligible.

As for the ionization by ⁴⁰K, the characteristics are similar, but the maximum ratios ζ_max(⁴⁰K)/ζ_init(⁴⁰K) decrease by a factor of about 2.7 in comparison to ²⁶Al. The resultant rates ζ_max(⁴⁰K) exceed the ionization rates by cosmic rays around the midplane of the orbital region R ≲ 2 AU, where the attenuation of the cosmic-ray flux becomes efficient.

Compared to the results in the MMSN, the maximum ratios ζ_max(X)/ζ_init(X) for the disk models of (b) $f_\Sigma =3$ and (c) 10 increase only by a factor of about 1.4–1.9 and 1.8–2.7, respectively, at the orbit R ≲ 1 AU. Because many dust particles can grow massive in the disks of large $f_\Sigma$, considerable amounts (30%–60%) of dust material are removed due to the rapid migration toward the central star. Thus, fractions of energy available to ionize the gas $E_\mathrm{av}(\mathrm{X})/\overline{E}_\mathrm{em}(\mathrm{X})$ decrease significantly. The residual amounts of dust material in the thin layers, however, are sufficient to decrease the portion of energy to ionize gas, p^gas ≲ 0.6.

In the intermediate orbital regions 1 AU ≲ R ≲ 2.7 AU, the characteristics of maximum ratios ζ_max(X)/ζ_init(X) are similar to those in R ≲ 1 AU for all the disk models. As R increases, both the available energy $E_\mathrm{av}(\mathrm{X})/\overline{E}_\mathrm{em}(\mathrm{X})$ and the factor p^gas increase, and a fraction Σ^dust_rm(R)/Σ^dust(R) becomes minor even for the disk of $f_\Sigma \approx 10$. Thus, the maximum ratios ζ_max(X)/ζ_init(X) decrease and approach gradually the values proportional to the relative amount of dust material in the thin layer Σ(R)/2χ_R.

The relative abundance and the inner solid density ρ_s of dust particles change suddenly at the orbit R ≈ 2.7 AU where the temperature is below ∼180 K and the water vapor begins to condense into ice. Hence, the maximum ratio ζ_max(X)/ζ_init(X) decreases abruptly by a factor of 2–3. In the outer orbital region R ≳ 2.7 AU, ζ_max(X)/ζ_init(X) decreases slightly and approaches constant values 11 and 4 for X = ²⁶Al and ⁴⁰K, respectively. Because the minimum heights of the thin layer Z_R are $H(R)/7.6 f_\Sigma$ and $H(R)/2.7 f_\Sigma$ for R = 5 and 10 AU, the concentration of dust material in the thin layers and fractions Σ^dust_rm(R)/Σ^dust(R) (≲ 0.04) become of little importance. As R increases further, the radius of the most abundant dust particles and the energy loss in the solid material decrease gradually. Thus, the ratios ζ_max(X)/ζ_init(X) increase to about 14 and 6 in the outermost region for a given $f_\Sigma$.

As the additional settling toward the midplane proceeds, the relative abundance of dust material in the outer part z ≲ Δz (≈ Z_R) decreases considerably. The ratio ζ(X)/ζ_init(X) in this part begins to decrease and comes closer to the typical value outside the innermost zone obtained numerically (see Figures 6, 7, and 9(a)). This tendency is independent of the radial distance R and the disk models $f_\Sigma$ and q. Similar to the intermediate zone and the outermost zone near the disk scale height, the ionization by radionuclides in this part becomes completely inefficient as the settling of dust particles proceeds.

In numerical calculations, we implicitly include the radial motion of dust particles by allowing the growth of particles above the critical radius a_cr. If dust particles cannot move in the R-direction completely, the maximum radius of dust particles is given by the total amount of dust material in the z-direction for a given R, a_s = p_sΣ^dust(R)/8ρ_s (Safronov 1969; see Table 3 for definite values). A similar kind of dust radius is obtained from the criteria for the coalescence of dust particles in a possible physical situation of the disk (Ohtsuki et al. 1990; Nakagawa & Kohno 1999): the collision and growth of particular dust particles whose number in the narrow area (≈10³ cm²) for a certain vertical interval becomes less than one are prohibited in calculating Equation (20) numerically. However, there exists the fluctuating motion of dust particles beyond this limited area because the movement of particles in the R-direction is not so small even in the zone of z ≳ Z₁. This allows further growth of such scarce particles, with the result that the dust radius reaches up to a_cr in our calculations. The required distance in the movement of particle growth is still negligibly small, compared with the radial migration of particles during the sedimentation (Nakagawa et al. 1986).

Most of the other short- and long-lived radionuclides such as ³⁶Cl, ⁶⁰Fe, and decay series of Th, Ac, and U (²³²Th, ²³⁵U, and ²³⁸U), which are mainly contained in the solid material of dust particles, also contribute the ionization of the gas. As for the short-lived nuclide ³⁶Cl the available energy ∼270 keV is mainly emitted in the form of the kinetic energy of the electron by β⁻ decay. In the decay of ⁶⁰Fe two photons with energies around 1.3 MeV and electrons with kinetic energies ∼60–400 keV are emitted. The energies of these photons are somewhat smaller than the one emitted in the decay of ⁴⁰K, whose residual energy $\overline{E}_\mathrm{res}({}^{40}\mathrm{K}; \gamma)$ is shown in Figure 2. Because the stopping ranges of electrons with these energies are less than ∼0.2 g cm⁻², they can scarcely penetrate the dust particles of a ≳ 0.1 cm. Thus, the residual energy $\overline{E}_\mathrm{res}({}^{60}\mathrm{Fe})$ changes in a similar way, and the settling of dust particles increases the ratio ζ(⁶⁰Fe)/ζ_init(⁶⁰Fe) by almost the same factor as that of ²⁶Al. Its contribution, however, is minor since the initial ionization rate ζ_init(⁶⁰Fe) is at least 10 times smaller than ζ_init(²⁶Al) (UN09).

The long-lived nuclides Th, Ac, and U series, on the other hand, mainly emit α-particles and electrons in the course of sequential decays toward the stable nuclides (²⁰⁸Pb, ²⁰⁷Pb, and ²⁰⁶Pb). The kinetic energies of α-particles are roughly several MeV, and their stopping ranges are ∼0.01 g cm⁻² or less. The energies of electrons are at most a few MeV, and the electron with the average energy (≲ 1 MeV) has a range around 0.5 g cm⁻². Because these ranges are smaller than the mass thickness of the most abundant dust particles in the thin layer around midplane, an increase in the ratios for these series is insignificant, compared with the ratio ζ(⁴⁰K)/ζ_init(⁴⁰K) obtained in Section 4.

In the outer region of the disk (R ≳ 2.7 AU) most of the dust particles have the icy mantel and their inner solid density decreases to about 1 g cm⁻³ on average. In calculating the ratio ζ(X)/ζ_init(X) in this region, we simply use ρ_s = 1 g cm⁻³ and do not consider the detailed structure of dust particles. Because the radionuclides mentioned above are overwhelmingly contained within the core of the dust particle, these nuclides are not distributed homogeneously in the solid material of the dust particle. However, we can safely regard that the parent nuclides are uniformly distributed throughout the solid material because the energy loss of the emitted particles becomes important after an extensive number of dust particles are aggregated together. Furthermore, the mean energy-loss rates are essentially proportional to the ratio Z_i/A_i of absorber i and the atomic properties of the absorber affect slightly the loss rates (see Section 3). Both icy materials (mantles) and minerals (cores) of dust particles have almost the same value of Z_i/A_i. Thus, we can use with sufficient accuracy the residual energies $\overline{E}_\mathrm{res}(\mathrm{X})$ for ρ_s = 1 g cm⁻³, which resemble the results shown in Figures 1 and 2.

It should also be noted that the decay of radionuclides in the dust particles affects little the sublimation of icy material in dust particles. For a spherical dust particle with radius a, specific heat C, and thermal conductivity κ, the diffusion time of heat and the maximum rise in temperature are approximately given by t_d = ρ_sCa²/κ and ΔT(X) = Q₀(X)a²/6κ, respectively, where Q₀(X) is the energy per unit volume and time emitted by the decay of nuclide X. In the case of icy dust aggregates for a ≈ a_cr around R = 5 AU of the MMSN most of the emitted energies are deposited in the solid material, and we obtain t_d ≈ 6 × 10⁵ s and ΔT(²⁶Al) ≈ 2 × 10⁻⁵ K if the abundance ratio is taken to be the canonical value ²⁶Al/²⁷Al ≈ 5 × 10⁻⁵ (UN09). Thus, we can safely neglect the effect of the emitted energies on an increase in the dust temperature. The decay of radionuclides becomes an essential heat source when the radius exceeds ∼7 × 10⁶ cm, which corresponds to planetesimal rather than dust aggregate.

Following most of the previous studies on the settling and size growth of dust particles in the protoplanetary disks, we neglect the internal structure of dust aggregates and assume compact aggregates with a constant density. At the early stage of dust growth, however, dust particles collide mutually with a significantly low velocity (≲ 1 mm s⁻¹), and the resultant aggregates can simply stick without any restructuring and form clusters of ballistic cluster–cluster aggregation (BCCA). Dominik & Tielens (1997) investigated the property of such dust aggregations in the mutual collisions by performing two-dimensional numerical simulations of aggregate collisions. Wada et al. (2007, 2008) and Suyama et al. (2008) performed more extensively two- and three-dimensional numerical simulations of the collisions of BCCA clusters to clarify the property of the resultant aggregates.

The important results obtained in these investigations are that the dust aggregates can retain their porous structure after the sticking growth and that the fractal dimension of the aggregates is two or less for BCCA clusters. Even if the aggregates collide continuously, the aggregates compressed have the inner density of about 10⁻⁴ g cm⁻³. Because the mass thickness of such dust aggregates decreases at least by a factor of 500, the energetic particles emitted by the decay of radionuclides can escape freely from the aggregates and ionize the surrounding gas. Such situations must change the ionization rate significantly, and both nuclides ²⁶Al and ⁴⁰K have the similar ratio ζ/ζ_init because the energy loss of the emitted particles in the solid material becomes negligible. In the case of extremely low ρ_s, the critical radius of the dust particle for radial migration always satisfies the condition a_cr ≳ l^gas, and its value is given by the latter expression of Equation (32), a_cr ≈ [3c_thρ^gasl^gas/2ρ_sΩ_K(R)]^1/2. Thus, the critical mass and radius of aggregates increase by a factor of ∼10² for ρ_s ≈ 10⁻⁴ g cm⁻³. Because a_s is inversely proportional to ρ_s, dust aggregates can grow larger than a_cr and a fraction of removed dust material Σ^dust_rm(R)/Σ^dust(R) becomes significant. Consequently, the limit of the ratio of the ionization rates in the thin layer is given by

$$\begin{eqnarray} \frac{\zeta _\mathrm{lim}}{\zeta _\mathrm{init}} &=& \frac{H(R)}{\Delta z}\left[1-\frac{\Sigma _\mathrm{rm}^\mathrm{dust}(R)}{\Sigma ^\mathrm{dust}(R)}\right]p^\mathrm{gas}\nonumber\\ &\approx& \frac{\Sigma (R)}{2 \chi _\mathrm{R}}\left[1-\frac{\Sigma _\mathrm{rm}^\mathrm{dust}(R)}{\Sigma ^\mathrm{dust}(R)}\right]p^\mathrm{gas}, \end{eqnarray} \tag{ 34 }$$

where we also take Δz = Z_R in the last expression.

Since this situation is analogous to that investigated in Section 4.3, we use Σ^dust_rm(R)/Σ^dust(R) obtained numerically and calculate the limits of the ratios ζ_lim/ζ_init. The results with Δz = Z_R for the disk models of $f_\Sigma =1$, 3, and 10 are also represented in Table 3.

The limits of the ratios become at most twice as large as ζ_max(²⁶Al)/ζ_init(²⁶Al). They do not take extreme values even in the inner orbital region R ≲ 1 AU for the disk model $f_\Sigma =10$ because the shielding factor p^gas becomes efficient to decrease the ionization rates. This indicates that the limits ζ_lim/ζ_init in a particular condition of Σ^dust_rm(R) = 0 are about 360, 280, and 43 for R = 0.5, 1, and 5 AU of the disk with $f_\Sigma =10$, and 270 for R = 0.5 AU of the disk with $f_\Sigma =3$, respectively. At the outer orbits R ≳ 2.7 AU an increase in the ratios remains large, but their absolute values are moderate (≲ a few × 10) for all $f_\Sigma$.

The extreme decrease of the inner solid density ρ_s in dust aggregates also affects the settling and/or growing times of dust particles. Since the maximum radius of dust particles a_s for ρ_s ≲ 10⁻⁴ g cm⁻³ is estimated to be much larger than l^gas, the drag force between gas and dust particles begins to follow Stokes' law, rather than Epstein's law, on the way to the midplane, or very probably to z ≈ Z₁. In the same way as the analytical investigation of the settling and size growth of a representative dust particle (Nakagawa et al. 1986), we can estimate how these processes and the relevant timescales change due to extremely low ρ_s of dust aggregates.

The evolution of dust particles in such situations is divided into the following three stages: (1) Dust particles with ρ_s ≈ 1 g cm⁻³ stick together to form the dust aggregates such as BCCA clusters, which have the extremely low internal density ρ_s ≲ 10⁻⁴ g cm⁻³ and a different fractal dimension from three. (2) Dust aggregates further settle and grow in mutual collisions, and eventually their radii exceed the mean free path of gas particles l^gas. (3) Dust aggregates continue to settle and grow under Stokes' law of drag force, and their radii gradually approach a_s with the extremely low ρ_s. In the case that the drag force is in accordance with Epstein's law (a ≲ l^gas), the settling time of a representative dust particle is given by

$$\begin{equation} t_\mathrm{E}(z) = \frac{8\Sigma ^\mathrm{gas}(R)}{\pi ^2\Sigma ^\mathrm{dust}(R)p_\mathrm{s}}\frac{t_\mathrm{K}}{1+(a_\mathrm{init}/a_\mathrm{s})} \log \left(\frac{a(z)H(R)}{a_\mathrm{init}z}\right), \end{equation} \tag{ 35 }$$

where a(z) = a_init + a_s(1 − z/H(R)) is the radius of the dust aggregate at the height z from the midplane and Σ^gas(R) is basically the same as Σ(R) given by Equation (1). Note that in stage (1) the radius of dust particles a(z) cannot directly relate with z due to the extensive decrease of ρ_s and the resultant increase of a_s during the transition from particles to aggregates. If the dust aggregate can settle to the narrow zone of the large ionization rate around the midplane during stages (1)–(2), we can easily estimate from Equation (35) that the settling (growth) time of the dust aggregate becomes somewhat longer than that obtained in Section 4.1.1. This is mainly because the dependence of the coefficient A (see Equation (12)) on ρ_s and a(z) canceled each other and the resultant factor for the settling time depends on a(z) logarithmically, as can be seen from Equation (35). Thus, the settling time of the aggregate toward a(z) = l^gas increases at most by several times.

If a dust aggregate grows further under the Stokes' law of drag force (a ≳ l^gas), the settling time in stage (3) is given by

$$\begin{eqnarray} t_\mathrm{S}(z) &=& \frac{12\Sigma ^\mathrm{gas}(R)}{\pi ^2\Sigma ^\mathrm{dust}(R)p_\mathrm{s}}\frac{l^\mathrm{gas}}{a_\mathrm{s}} \frac{t_\mathrm{K}}{1+(a_\mathrm{init}/a_\mathrm{s})}\nonumber\\ && \times\left[\frac{a_\mathrm{s}}{l^\mathrm{gas}}-\frac{a_\mathrm{s}}{a(z)}-\frac{1}{1+(a_\mathrm{init}/a_\mathrm{s})}\log \left(\frac{a(z)z_l}{l^\mathrm{gas}z}\right)\right],\nonumber\\ \end{eqnarray} \tag{ 36 }$$

where z_l is the height from the midplane for a(z) = l^gas. Because the timescale that a(z) becomes a_s under the Stokes' law is close to t_E(z_l), the total settling time for the dust aggregate with ρ_s ≲ 10⁻⁴ g cm⁻³ cannot be different from the value with ρ_s ≈ 1 g cm⁻³ obtained in Section 4.1.1 by an order of magnitude. In conclusion, the porous structure of dust aggregates affects moderately both the ratio of the ionization rates and the settling timescale. In order to clarify the precise values of these quantities in the thin layers of various disk models, we need further (numerical) investigation of the settling and size growth of dust aggregates.

Zsom et al. (2010) modeled the growth of dust aggregates at 1 AU in the midplane using the results of laboratory experiments and investigated the maximum aggregate size and mass that can be reached by coalescence. The collision model used was based on a continuous transition of aggregates between compact and porous spheres. They applied the collision model of dust particles to the Monte Carlo method and found that dust particles do not follow the growth–fragmentation cycles: the size growth is stopped by the bouncing barrier rather than fragmentation; the maximum mass of the particle is limited to ∼1 g (chondrule-sized particle). This mass for a compact sphere model is larger than the dust mass so that the energy loss of the emitted electron/positron becomes efficient (see Figures 1 and 2). Because the energy loss of emitted photons can be negligible for dust aggregates of this size, the ionization rates by ²⁶Al and ⁴⁰K around the midplane increase by factors of about 50 and 19 (see Table 3), respectively, as long as the settling of dust particles toward the narrow zone close to the midplane can be possible.

Dust particles in the protoplanetary disk generally absorb electrons or ions on their surface (Umebayashi 1983). Thus, the electric charge state of dust aggregates affects the fractal growth of aggregates, and the electric repulsion between colliding aggregates can impose a serious limitation on dust evolution (Okuzumi 2009). Okuzumi et al. (2011a, 2011b) generally investigated collisional growth of charged aggregates and found that the aggregates under a certain condition undergo bimodal growth where a certain number of aggregates stop growing due to the electrostatic repulsion and the rest continue growing by colliding with each other. They found that regions where the radionuclides become efficient sources for ionizing the gas can be classified as either bimodal growth, mentioned above, or total-freezeout growth, where all aggregates stop growing with a nearly monodisperse distribution. Furthermore, the charge state of aggregates can delay or prevent the settling of aggregates. Thus, we need a detailed study to clarify its effects on the ionization rates by radionuclides.

In order to understand the physical and chemical evolution of the protoplanetary disks, we must clarify the ionization state in all the evolutionary stages of various disks. Because hard X-ray emission from the central star can penetrate and ionize the disk, we need to calculate the hard X-ray ionization rate as a universal function of the disk structure and X-ray luminosity and temperature, which is similar to the "soft" X-ray ionization rate obtained by Igea & Glassgold (1999). Furthermore, X-rays contribute the photoelectric emission from the dust surface, which can change the charge state of dust particles considerably. We also need a detailed study of this process. When the settling and size growth of dust particles proceed extensively, the charge state of gas and dust particles in the inner region of the disk is affected by the ionization of radionuclides. In addition, both settling and size growth of dust particles vary the total surface area of dust particles, which provide important places for sticking and reaction (especially recombination) of various chemical species in gas phase. Thus, we must investigate the chemical network of gaseous and surface reactions simultaneously, which is beyond the scope of the present series of papers.

The results obtained in Section 4 depend slightly on models of the protoplanetary disk. As long as the state of a disk is quiescent enough to proceed the settling of dust particles, an increase in the ionization rates around the midplane occurs in the orbital region where the surface density is much larger than the range of energetic particles emitted by the decay of radionuclides, Σ(R) ≳ 100 g cm⁻². Essentials are the physical situations of a gaseous disk, i.e., whether or not the settling and the size growth of dust particles can proceed to a considerable extent in an actual disk. The turbulent motion of the gas affects both settling and size growth of dust particles, which are the competitive processes to the change of the ionization rates by radionuclides. Thus, we must investigate the settling and size growth of dust particles in disks of different turbulent strengths and clarify the physical conditions under which turbulent motion becomes important to the change of the ionization rates by radionuclides. We will deal with these issues in Paper II.

We thank the anonymous referee for careful reading of the manuscript and his comments, which were greatly helpful in improving our paper. We also thank Eric Feigelson for suggesting the importance of hard X-ray emission from YSOs to the ionization of disk. One of the authors (T.U.) is grateful to Y. Nakagawa for his hospitality in CPS.

Here we simply assume that molecules are smooth rigid elastic spheres not surrounded by force fields. Molecules affect each other's motion only at collisions in this case.

Consider collisions occurring between pairs of molecules of masses and diameters m₁, m₂ and d₁, d₂ in a gas mixture at rest in a uniform steady state. The total number of collisions between two molecules per unit volume and time is given by

$$\begin{equation} N_{12} = n_1 n_2 \sigma _{12}\left(\frac{8kT}{\pi \mu _{12}}\right)^{1/2}, \end{equation} \tag{ A1 }$$

where n₁ and n₂ are the number densities of molecules 1 and 2, and σ₁₂ = πd₁₂² = π(d₁ + d₂)²/4 and μ₁₂ = m₁m₂/(m₁ + m₂) are the geometric cross section for mutual collisions and reduced mass, respectively (for details see, e.g., Chapman & Cowling 1990). Changing the suffix 2 to 1 in Equation (A1), we get the number of collisions between pair of molecules 1, $({1}/{2})N_{11}=({1}/{2}){n_1}^2\sigma _{11}\sqrt{8kT/\pi \mu _{11}}=2{n_1}^2\pi {d_1}^2\sqrt{kT/\pi m_1}$ simply because N₁₁ counts each collision between a pair of molecules twice over. On the other hand, it should be noted that the average number of collisions of any one molecule with similar molecules per unit time is N₁₁/n₁ since each collision affects both molecules at once.

The collision interval for molecule 1, or mean time between successive collisions with all kinds of constituents, is

$$\begin{equation} \tau _1 = n_1/(N_{11}+N_{12}+\cdots), \end{equation} \tag{ A2 }$$

where the number of terms in the bracket is equal to the number of constituents in the mixture. The mean free path is the mean distance l₁ traveled by a molecule 1 between successive collisions in a given time t. Dividing the total distance $n_1\overline{c}_1t$ traveled by molecules 1 in this time by the total number of their collisions n₁t/τ₁, we obtain

$$\begin{equation} l_1 = \overline{c}_1\tau _1 = n_1\overline{c}_1/(N_{11}+N_{12}+\cdots) = 1/(\overline{\sigma }_1n_1), \end{equation} \tag{ A3 }$$

where $\overline{c}_1=\sqrt{8kT/\pi m_1}$ is the mean value of the thermal speed for molecules 1 and $\overline{\sigma }_1$ corresponds to the average cross section for mutual collision of molecules 1.

For the gas in the protoplanetary disk, we can consider only H₂ molecules and He atoms (hereafter denoted by suffixes 1 and 2) as neutral constituents in the mixture. The mean free path of neutrals l_n is found to be

$$\begin{eqnarray} l_\mathrm{n} &=\frac{n_1\overline{c}_1t + n_2\overline{c}_2t}{n_1t/\tau _1+n_2t/\tau _2} =\frac{n_1\overline{c}_1 + n_2\overline{c}_2}{N_{11}+2N_{12}+N_{22}} \nonumber \\ & =\frac{n_1/\sqrt{m_1}+n_2/\sqrt{m_2}}{{n_1}^2\sigma _{11}\sqrt{2/m_1}+2n_1n_2\sigma _{12}\sqrt{1/\mu _{12}}+{n_2}^2\sigma _{22}\sqrt{2/m_2}}.\nonumber\\ \end{eqnarray} \tag{ A4 }$$

Using the number density of neutrals n_n = n₁ + n₂, the average cross section for mutual collisions of neutrals is given by $\overline{\sigma }_\mathrm{n}=1/(l_\mathrm{n}n_\mathrm{n})$.

The data used for numerical calculation, the average cross sections for mutual collisions, and mean free paths at Earth's orbit (see Equation (13)) in the MMSN obtained for a gas mixture with solar abundance (X = 0.73, Y = 0.25) are summarized in Table 4. The average collisional cross section of He atoms $\overline{\sigma }_2$ is seemingly large, because this value, defined as $\overline{\sigma }_2=1/(l_\mathrm{2}n_2)$, corresponds closely to the geometric cross section for mutual collisions times the ratio of the number densities n_n/n₂ (≈7).

Table 4. Average Cross Sections for Mutual Collisions $\overline{\sigma }_i$ and Mean Free Paths l_i at the Orbit R = 1 AU in the MMSN for Neutral Gas Constituents i in a Gas Mixture with Solar Abundance

Constituent		Diameter	Cross Section	Mean Free Path
	(i)	di (cm)	$\overline{\sigma }_i$ (cm2)	li (cm)
H2	(1)	2.2 × 10−8a	2.54 × 10−15	1.32
He	(2)	2.4 × 10−8	1.95 × 10−14	1.01
Neutral	(n)	...	2.24 × 10−15	1.28

Note. ^aThe diameter is estimated to be d₁ = 3r₁, where r₁ = 7.41 × 10⁻⁹ cm is the internuclear distance of the H₂ molecule.

Download table as:  ASCIITypeset image

Based on the assumptions described in Section 3.2, we calculate analytically an average traveling path length 〈s〉 of the particles emitted by the decay of nuclide X in dust particles.

We consider a spherical dust particle of radius a and take a coordinate system in which its origin O is in the center of the dust particle and position A, where an energetic particle is emitted, is on the z-axis. Letting position B on the surface of the dust particle be the end point of the path of the emitted particle, we obtain the length s of line segment AB as

$$\begin{equation} s = -r\cos \theta +\sqrt{a^2-r^2\sin ^2\theta }, \end{equation} \tag{ B1 }$$

where r and θ are the distance from the origin O to A and the angle between the line AB and the z-axis, respectively. Using the assumption that the direction of the emitted particle is isotropic around A, the length $\overline{s}$ averaged over the emitted solid angle Ω is given by

$$\begin{eqnarray} \overline{s}&= \frac{1}{4\pi }\int _\Omega s\,d\Omega = \frac{1}{4\pi }\int _0^\pi \int _0^{2\pi }s\,d\varphi \,\sin \theta \,d\theta\nonumber\\ &= \frac{a}{2}\left[1+\frac{1-(r/a)^2}{2r/a}\log \frac{1+r/a}{1-r/a}\right], \end{eqnarray} \tag{ B2 }$$

where φ is the azimuthal angle of line segment AB in the above coordinate system.

To obtain the average traveling path length 〈s〉 of the emitted particle, we must further average $\overline{s}$ over position A, which is assumed to be distributed uniformly throughout the volume V of the dust particle. We also calculate it analytically and obtain

$$\begin{eqnarray} \langle s \rangle &=& \int _V \bar{s}\,dV / \int _V dV\nonumber\\ &=& \frac{3}{a^3}\int _0^a\frac{a}{2}\left[1+\frac{1-(r/a)^2}{2r/a}\log \frac{1+r/a}{1-r/a}\right]r^2\,dr = \frac{3}{4}a.\nonumber\\ \end{eqnarray} \tag{ B3 }$$

Thus, we can approximately regard that the average path length 〈s〉 of the emitted particle in the dust particle is somewhat smaller than its radius a.

In the case of positron annihilation two photons are emitted in the exactly opposite senses and the traveling path lengths of these photons in dust particles are correlated with each other. Hence, we must average the length s more carefully.

Let the positions B and C on the surface of the dust particle be the end points of the paths of two photons emitted in the opposite senses. Then the lengths s_AB and s_AC of line segments AB and AC are given by

$$\begin{eqnarray} s_\mathrm{AB} &= -r\cos \theta +\sqrt{a^2-r^2\sin ^2\theta },\nonumber\\ s_\mathrm{AC} &= r\cos \theta +\sqrt{a^2-r^2\sin ^2\theta }, \end{eqnarray} \tag{ B4 }$$

where θ is the angle between the line AB and the z-axis as above. If we restrict θ to a value less than a right angle (0 ⩽ θ ⩽ π/2), s_AB normally represents the shorter traveling path length s₋ of the emitted photons and s_AC represents the longer one s₊. Averaging s_± over the solid angle Ω_θ with this restricted extent of θ, we obtain the mean lengths of short and long traveling paths of photons emitted simultaneously. The results are as follows:

$$\begin{eqnarray} \overline{s}_\mathrm{short} &=& \frac{1}{2\pi }\int _{\Omega _\theta } s_-\,d\Omega = \frac{1}{2\pi }\int _0^{\pi /2} \int _0^{2\pi }s_-\,d\varphi \,\sin \theta \,d\theta\nonumber\\ &=& {-}\frac{r}{2}+\frac{a}{2}\left[1+\frac{1-(r/a)^2}{2r/a}\log \frac{1+r/a}{1-r/a}\right], \end{eqnarray} \tag{ B5 }$$

$$\begin{eqnarray} \overline{s}_\mathrm{long} &=& \frac{1}{2\pi }\int _{\Omega _\theta } s_+\,d\Omega = \frac{1}{2\pi }\int _0^{\pi /2} \int _0^{2\pi }s_+\,d\varphi \,\sin \theta \,d\theta\nonumber\\ &=& \frac{r}{2}+\frac{a}{2}\left[1+\frac{1-(r/a)^2}{2r/a}\log \frac{1+r/a}{1-r/a}\right]. \end{eqnarray} \tag{ B6 }$$

Note the position A corresponds to the spot where the positron annihilation occurs and that its distance from the origin r in Equations (B5) and (B6) is roughly equal to the range 〈s(e⁺)〉 of this positron.