RAPID COAGULATION OF POROUS DUST AGGREGATES OUTSIDE THE SNOW LINE: A PATHWAY TO SUCCESSFUL ICY PLANETESIMAL FORMATION

We adopt the minimum-mass solar nebula (MMSN) model of Hayashi (1981) with a solar-mass central star. The radial profiles of the gas surface density Σ_g and disk temperature T are given by Σ_g = 152(r/5 AU)^−3/2 g cm⁻² and T = 125(r/5 AU)^−1/2 K, respectively, where r is the distance from the central star. In this study, we focus on dust evolution outside the snow line, which is located at r ≈ 3 AU in the adopted disk model. The vertical structure is assumed to be in hydrostatic equilibrium, and hence the vertical structure of the gas density ρ_g is given by $\rho _g = (\Sigma _g/\sqrt{2\pi }h_g)\exp (-z^2/2h_g^2)$, where h_g = c_s/Ω is the gas scale height, c_s is the isothermal sound speed, and Ω is the Kepler frequency. The isothermal sound speed is given by $c_s = \sqrt{k_{\rm B}T/m_g}$, where k_B is the Boltzmann constant and m_g is the mean molecular mass. We assume the mean molecular weight of 2.34, which gives m_g = 3.9 × 10⁻²⁴ g and c_s = 6.7 × 10⁴(r/5 AU)^−1/4 cm s⁻¹. The assumed stellar mass (1 M_☉) leads to $\Omega = \sqrt{{\sc G}M_\odot /r^3} = 1.8\times 10^{-8}(r/5\, {\rm AU})^{-3/2}\, {\rm rad\, s^{-1}}$ and h_g/r = 0.05(r/5 AU)^1/4, where G is the gravitational constant.

In reality, protoplanetary disks can be heavier than the MMSN. The gravitational stability criterion (Toomre 1964) Σ_g < Ωc_s/πG ≈ 5.6 × 10³(r/5 AU)^−7/4 g cm⁻² allows the surface density to be up to about 10 times higher than the MMSN value. The dependence of our result on the disk mass will be analytically discussed in Section 4.

Initial dust particles are modeled as compact spheres of equal size a₀ = 0.1 μm and equal internal density ρ₀ = 1.4 g cm⁻³, distributed in the disk with a constant dust-to-gas surface density ratio Σ_d/Σ_g = 0.01. The mass of each initial particle is m₀ = (4π/3)ρ₀a³₀ = 5.9 × 10⁻¹⁵ g. In the following, we will refer to the initial dust particles as "monomers." We define the radius of a porous aggregate as $a = [(5/6N)\sum _{i=1}^N\sum _{j=1}^N({\bm x}_i-{\bm x}_j)^2]^{1/2}$, where N is the number of the constituent monomers and ${\bm x}_k (k=1,2,\dots,N)$ is the position of the monomers (Mukai et al. 1992). This definition is in accordance with previous N-body experiments (Wada et al. 2008; Suyama et al. 2008; Okuzumi et al. 2009) on which our porosity model is based (see Section 2.3.1).

Disk turbulence affects the collision and sedimentation of dust particles. To include these effects, we consider gas turbulence in which the turnover time and the mean-squared random velocity of the largest turbulent eddies are given by t_L = Ω⁻¹ and δv²_g = α_Dc²_s, respectively, where α_D is the dimensionless parameter characterizing the strength of the turbulence. The assumption for t_L is based on theoretical anticipation for turbulence in Keplerian disks (Dubrulle & Valdettaro 1992; Fromang & Papaloizou 2006; Johansen et al. 2006). The diffusion coefficient for the gas is given by D_g = δv²_gt_L = α_Dc²_s/Ω. If the gas diffusion coefficient is of the same order as the turbulent viscosity, α_D is equivalent to the so-called alpha parameter of Shakura & Sunyaev (1973). However, for simplicity, we do not consider the viscous evolution of the gas disk. We adopt α_D = 10⁻³ in our numerical simulations. A higher value of α_D would cause catastrophic collisional fragmentation of aggregates, which is not considered in this study (see Section 5.3).

We solve the evolution of the radial size distribution of dust aggregates using the method developed by Brauer et al. (2008a). This method assumes the balance between sedimentation and turbulent diffusion of aggregates in the vertical direction. Thus, the vertical number density distribution of aggregates is given by a Gaussian $({\cal N}/\sqrt{2\pi }h_d)\exp (-z^2/2h_d^2)$, where ${\cal N}(r,m)$ is the column number density of aggregates per unit mass and h_d(r, m) is the scale height of aggregates at orbital radius r and with mass m (Dubrulle et al. 1995). This approach is valid if the coagulation timescale is longer than the settling/diffusion timescale, which is true except for very tiny particles with short collision times (Zsom et al. 2011).

The evolution of the radial size distribution ${\cal N}(r,m)$ is given by the vertically integrated advection–coagulation equation, which reads (Brauer et al. 2008a)

$$\begin{eqnarray} \frac{\partial {\cal N}(r,m)}{\partial t} &=& \frac{1}{2}\int _0^m\!\! K(r,m^{\prime },m-m^{\prime }){\cal N}(r,m^{\prime }){\cal N}(r,m-m^{\prime })dm^{\prime } \nonumber \\ && - {\cal N}(r,m)\int _0^\infty\!\! K(r,m,m^{\prime }){\cal N}(r,m^{\prime })dm^{\prime } \nonumber \\ && - \frac{1}{r}\frac{\partial }{\partial r}\left[ rv_r(r,m){\cal N}(r,m)\right], \end{eqnarray} \tag{ 1 }$$

where v_r is the radial drift velocity and K is the vertically integrated collision rate coefficient given by

$$\begin{equation} K(r,m_1,m_2) = \frac{\sigma _{\rm coll}}{2\pi h_{d,1} h_{d,2}} \int _{-\infty }^\infty \Delta v \exp \left(-\frac{z^2}{2h_{d,12}^2}\right)dz. \end{equation} \tag{ 2 }$$

Here, σ_coll is the collisional cross section, h_{d, 1} and h_{d, 2} are the scale heights of the colliding aggregates 1 and 2, Δv is the collision velocity, and h_{d, 12} = (h⁻²_{d, 1} + h⁻²_{d, 2})^−1/2. As mentioned in Section 1, we neglect electrostatic and gravitational interactions between colliding aggregates and assume perfect sticking upon collision. Thus, the collisional cross section is simply given by σ_coll = π(a₁ + a₂)², where a₁ and a₂ are the radii of the colliding aggregates. The validity of neglecting fragmentation will be discussed in Section 5.3.

The dust scale height h_d in the sedimentation–diffusion equilibrium has been analytically obtained by Youdin & Lithwick (2007). For turbulence of t_L = Ω⁻¹ and D_g = α_Dc²_s/Ω, it is given by

$$\begin{equation} h_d = h_g\left(1+\frac{\Omega t_s}{\alpha _D}\frac{1+2\Omega t_s}{1+\Omega t_s} \right)^{-1/2}, \end{equation} \tag{ 3 }$$

where t_s is the stopping time of the aggregates. We use this expression in this study.

For the stopping time, we use

$$\begin{equation} t_s = \left\lbrace \begin{array}{@{}ll}t_s^{\rm (Ep)} \equiv {\displaystyle \frac{3m}{4\rho _g v_{\rm th} A}}, & a< {\displaystyle \frac{9}{4}} \lambda _{\rm mfp}, \\ [3mm] t_s^{\rm (St)} \equiv {\displaystyle \frac{4a}{9\lambda _{\rm mfp}}} t_s^{\rm (Ep)}, & a> {\displaystyle \frac{9}{4}} \lambda _{\rm mfp}, \end{array} \right. \end{equation} \tag{ 4 }$$

where $v_{\rm th} = \sqrt{8/\pi }c_s$ and λ_mfp are the thermal velocity and mean free path of gas particles, respectively, and A is the projected area of the aggregate. The mean free path is related to the gas density as

$$\begin{equation} \lambda _{\rm mfp} = \frac{m_g}{\sigma _{\rm mol}\rho _g}, \end{equation} \tag{ 5 }$$

where σ_mol = 2 × 10⁻¹⁵ cm² is the collisional cross section of gas molecules. Our gas disk model gives λ_mfp = 120(r/5 AU)^11/4 cm at the midplane. Equation (4) satisfies the requirement that the stopping time must obey Epstein's law t_s = t^(Ep)_s at a ≪ λ_mfp and Stokes' law t_s = t^(St)_s at a ≫ λ_mfp, respectively. Since t^(St)_s∝at_s^(Ep), an aggregate growing in the Stokes regime decouples from the gas motion more quickly than in the Epstein regime. In reality, Stokes' law breaks down when the particle Reynolds number (the Reynolds number of flow around the particle) is much greater than unity, but we neglect this in our simulations for simplicity. We will discuss this point further in Section 5.1.

The radial drift velocity is taken as

$$\begin{equation} v_r = - \frac{2\Omega t_s}{1+(\Omega t_s)^2}\eta v_K, \end{equation} \tag{ 6 }$$

where

$$\begin{equation} 2\eta \equiv - \biggl ({\displaystyle \frac{c_s}{v_K}} \biggr) ^2 \frac{\partial \ln \big(\rho _g c_s^2\big)}{\partial \ln r} \end{equation} \tag{ 7 }$$

is the ratio of the pressure gradient force to the stellar gravity in the radial direction and v_K = rΩ is the Kepler velocity (Adachi et al. 1976; Weidenschilling 1977; Nakagawa et al. 1986). The radial drift speed has a maximum ηv_K, which is realized when Ωt_s = 1. In our disk model, η scales with r as η = 4.0 × 10⁻³(r/5 AU)^1/2, and the maximum inward speed ηv_K = 54 m s⁻¹ is independent of r. Since η is proportional to the gas temperature, the maximum drift speed would be somewhat lower in colder disk models (Kusaka et al. 1970; Hirose & Turner 2011). Equation (6) neglects the frictional backreaction from dust to gas assuming that the local dust-to-gas mass ratio is much lower than unity or the stopping time of aggregates dominating the dust mass is much longer than Ω⁻¹. We examine the validity of this assumption in Section 5.2.1.

In this paper, we also consider the collisional evolution of aggregate porosities. We treat the mean volume V = (4π/3)a³ of aggregates with orbital radius r and mass m as a time-dependent quantity. The evolutionary equation for V(r,m) is given by

$$\begin{eqnarray} \frac{\partial \left(V{\cal N}\right)}{\partial t} &=& \frac{1}{2}\int _0^m [V_{1+2}K](r,m^{\prime },m-m^{\prime }) \nonumber \\ && \times {\cal N}(r,m^{\prime }){\cal N}(r,m-m^{\prime })dm^{\prime } \nonumber \\ && -V(r,m){\cal N}(r,m)\int _0^\infty K(r,m,m^{\prime }){\cal N}(r,m^{\prime })dm^{\prime } \nonumber \\ && - \frac{1}{r}\frac{\partial }{\partial r}[rv_r(r,m)V(r,m){\cal N}(r,m)], \end{eqnarray} \tag{ 8 }$$

where

$$\begin{eqnarray} [V_{1+2}K](r,m_1,m_2) &=& \frac{\sigma _{\rm coll}}{2\pi h_{d,1} h_{d,2}} \int _{-\infty }^\infty V_{1+2}\Delta v \nonumber\\ &&\times\exp\left(-\frac{z^2}{2h_{d,12}^2}\right)dz \end{eqnarray} \tag{ 9 }$$

with V_{1 + 2} being the volume of merged aggregates (described in Section 2.3.1). Equation (8) is identical to the original evolutionary equation for V derived by Okuzumi et al. (2009, their Equation (16)) except that here we take the vertical integration of the equation and take into account the radial advection of dust. In deriving Equation (8), we have assumed that the dispersion of the volume is sufficiently narrow at every r and m (see Okuzumi et al. 2009). This "volume-averaging" approximation allows us to follow the porosity evolution of aggregates without solving higher-order moment equations for the volume, and hence with small computational costs. This approximation is valid unless the porosity distribution at fixed r and m is significantly broadened by, e.g., collisional fragmentation cascades (Okuzumi et al. 2009).

2.3.1. Porosity Change Recipe

The functional form of V_{1 + 2} determines the evolution of aggregate porosities in our simulation. In this study, we give V_{1 + 2} as a function of the volumes of the colliding aggregates, V₁ = V(r, m₁) and V₂ = V(r, m₂), and the impact energy E_imp = m₁m₂Δv²/[2(m₁ + m₂)]. Before introducing the final form of our porosity change recipe (Equation (15)), we briefly review recent N-body collision experiments on which our recipe is based.

Collisional compression depends on the ratio between E_imp and the "rolling energy" E_roll (Dominik & Tielens 1997; Blum & Wurm 2000; Wada et al. 2007). The rolling energy is defined as the energy needed for one monomer to roll over 90° on the surface of another monomer in contact (Dominik & Tielens 1997). When E_imp ≪ E_roll, two aggregates stick without visible restructuring (the so-called hit-and-stick collision; see Figure 1(b)). In this case, the volume of the merged aggregate is determined in a geometric manner, i.e., independently of E_imp. When E_imp ≳ E_roll, internal restructuring occurs through inelastic rolling among constituent monomers (Dominik & Tielens 1997; see also Figure 1(c)). In this case, the final volume V_{1 + 2} depends on E_imp as well as on V₁ and V₂.

Zoom In Zoom Out Reset image size

For hit-and-stick collisions (E_imp/E_roll → 0), Okuzumi et al. (2009) obtained an empirical formula for V_{1 + 2},

$$\begin{equation} V_{1+2} = V_{1+2,{\rm HS}} \equiv V_1 + V_2 + V_{\rm void}, \end{equation} \tag{ 10 }$$

where V₁ and V₂(⩽ V₁) are the volumes of the two colliding aggregates, and

$$\begin{equation} V_{\rm void} = \min \left\lbrace 0.99 - 1.03 \ln \biggl ({\displaystyle \frac{2}{V_1/V_2+1}} \biggr), 6.94 \right\rbrace V_2 \end{equation} \tag{ 11 }$$

is the volume of the voids created in the collision (see Figure 1(b)). For V₁ ≈ V₂, the void volume is approximately equal to V₁, and hence the volume of the new aggregate is approximately given by V_{1 + 2} ≈ 3V₁. This is equivalent to a fractal relation $V \propto m^{3/d_f}$, where d_f ≈ 2 (see Section 4.2.1 of Okuzumi et al. 2009).

In the limit of E_imp ≫ E_roll and for a head-on collision of equal-sized aggregates (V₁ = V₂), Suyama et al. (2008) showed that V_{1 + 2} obeys the relation

$$\begin{equation} E_{\rm imp} = -\int _{2^{6/5}V_{1}}^{V_{1+2}} P(V) dV. \end{equation} \tag{ 12 }$$

Here, P ≡ −dE_imp/dV is the dynamic compression strength of the merged aggregate given by (Wada et al. 2008)

$$\begin{equation} P(V) = 2 \biggl ({\displaystyle \frac{5}{3}} \biggr) ^{6}\frac{b E_{\rm roll}}{V_0} \biggl ({\displaystyle \frac{\rho _{\rm int}(V)}{\rho _0}} \biggr) ^{13/3}N_{1+2}^{2/3}, \end{equation} \tag{ 13 }$$

where b = 0.15 is a dimensionless fitting parameter, V₀ = m₀/ρ₀ = (4π/3)a³₀ is the monomer volume, N_{1 + 2} = 2m₁/m₀ is the number of monomers contained in the merged aggregate, and ρ_int = 2m₁/V is the internal density of the merged aggregate. If we substitute Equation (13) into Equation (12), we obtain the equation that explicitly gives V_{1 + 2} as a function of E_roll and V₁,

$$\begin{equation} V_{1+2} =\left[ {\displaystyle \frac{(3/5)^5 E_{\rm imp}}{N_{1+2}^5bE_{\rm roll}V_0^{10/3}}} +\big(2V_1^{5/6}\big)^{-4} \right]^{-3/10}. \end{equation} \tag{ 14 }$$

This equation basically expresses the energy balance in collisional compression, but some caution is needed in interpreting it. Firstly, the initial state for the compression is chosen to be V = 2^6/5V₁, although the volume just after contact is V = 3V₁ (see above). This is based on the fact that compaction from V = 3V₁ to V = 2^6/5V₁ occurs through partial compression of the new voids, which requires little energy (Suyama et al. 2008). Secondly, the dynamic compression strength P depends on mass N_{1 + 2} as well as on internal density ρ_int, meaning that P is not an intensive variable. This is due to the fact that dynamically compressed parts in the merged aggregate have a fractal structure with a fractal dimension of 2.5 (Wada et al. 2008). In fact, Equations (12)–(14) are more naturally described in terms of variables in the 2.5-dimensional space, V_f∝a^5/2, ρ_f∝N_{1 + 2}/V_f, and P_f = −dE_imp/dV_f (see Wada et al. 2008; Suyama et al. 2008). An important point here is that aggregates become stronger and stronger against dynamic compression as they grow because of the N^2/3_{1 + 2} factor in P.

Equations (10) and (14) express how the volume of the merged aggregate is determined in the limits of E_imp ≪ E_roll and E_imp ≫ E_roll, respectively. To properly take into account the intermediate cases (E_imp ∼ E_roll), we adopt an updated analytic formula given recently by Suyama et al. (2012). This reads

$$\begin{equation} V_{1+2} =\! \left\lbrace \begin{array}{@{}l}\left[ \left(1- {\displaystyle \frac{E_{\rm imp}}{3bE_{\rm roll}}} \right)V_{1+2,{\rm HS}}^{5/6} + {\displaystyle\frac{E_{\rm imp}}{3bE_{\rm roll}}} \left(V_{1}^{5/6}+V_{2}^{5/6}\right) \right]^{6/5} \\ \ms \quad \big({\rm if}\ V_{1+2,{\rm HS}}^{5/6}>V_1^{5/6}+V_2^{5/6}\ {\rm and}\ E_{\rm imp} < 3bE_{\rm roll}\big), \\ \ms \ms \left[ {\displaystyle\frac{(3/5)^5(E_{\rm imp}-3bE_{\rm roll})}{N_{1+2}^5bE_{\rm roll}V_0^{10/3}}} +\left(V_1^{5/6}+V_2^{5/6}\right)^{-4} \right]^{-3/10} \\ \ms \quad \big({\rm if}\ V_{1+2,{\rm HS}}^{5/6}>V_1^{5/6}+V_2^{5/6}\ {\rm and}\ E_{\rm imp} > 3bE_{\rm roll}\big), \\ \ms \ms \left[ {\displaystyle\frac{(3/5)^5 E_{\rm imp}}{N_{1+2}^5bE_{\rm roll}V_0^{10/3}}} +V_{1+2,{\rm HS}}^{-10/3} \right]^{-3/10} \\ \ms \quad \big({\rm if}\ V_{1+2,{\rm HS}}^{5/6}<V_1^{5/6}+V_2^{5/6}\big), \\ \end{array} \right. \end{equation} \tag{ 15 }$$

where N_{1 + 2} is now defined as (m₁ + m₂)/m₀. Note that this equation reduces to Equation (10) when E_imp ≪ E_roll, and to Equation (14) when V₁ = V₂ and E_imp ≫ E_roll. Suyama et al. (2012) derived Equation (15) by taking into account small energy losses in the partial compression of the new voids. In addition, unlike Equation (14), Equation (15) takes into account the cases where colliding aggregates have different volumes and masses ($V_1 \not= V_2$, $m_1 \not= m_2$). Suyama et al. (2012) confirmed that Equation (15) reproduces the results of numerical collision experiments within an error of 20% as long as the mass ratio m₂/m₁(⩽ 1) between the colliding aggregates is larger than 1/16.

We comment on three important caveats regarding our porosity recipe. Firstly, Equation (15) is still untested for cases where colliding aggregates have very different sizes (m₂/m₁ < 1/16). Therefore, the validity of using Equation (15) is at present only guaranteed for the case where "similar-sized" (m₂/m₁ ≳ 1/16) collisions dominate the growth of aggregates. We will carefully check this validity in Section 3.2. Secondly, Equation (15) ignores offset collisions, in which a considerable fraction of the impact energy is spent for stretching rather than compaction (Wada et al. 2007; Paszun & Dominik 2009). For this reason, Equation (15) underestimates the porosity increase upon collision. Thirdly, we do not consider non-collisional compression (e.g., static compression due to gas drag forces), which could contribute to the compaction of very large aggregates. We will discuss the second and third points in more detail in Section 5.4.

In addition to V, we need to know the projected area A of aggregates to calculate the stopping time t_s. Unfortunately, a naive relation A = πa² breaks down when the fractal dimension of the aggregate is less than 2, since πa² increases faster than mass for this case while A does not. A projected area growing faster than mass means a coupling to the gas becoming stronger and stronger as the aggregate grows, which is clearly unrealistic. To avoid this, we use an empirical formula by Okuzumi et al. (2009) that well reproduces the mean projected area ${\overline{A}}$ of aggregates with monomer number N = m/m₀ and radius a for both fractal and compact aggregates. With this formula, all aggregates in our simulations are guaranteed to decouple from the gas as they grow. We remark that this treatment is only relevant for fractal aggregates with d_f ≲ 2; for more compact aggregates, the empirical formula reduces to the usual relation A ≈ πa².

The rolling energy E_roll has so far not been measured for submicron-sized icy particles, but can be estimated in the following way. It is anticipated by microscopic friction theory (Dominik & Tielens 1995) that the critical rolling force F_roll ≡ E_roll/(πa₀/2) is a material constant (i.e., E_roll is proportional to the monomer radius a₀). A rolling force of F_roll = (1.15 ± 0.24) × 10⁻³ dyn has recently been measured for micron-sized ice particles (Gundlach et al. 2011). Given that F_roll is independent of a₀, the measured force implies the rolling energy of E_roll = (πa₀/2)F_roll ≈ 1.8 × 10⁻⁸ erg for a₀ = 0.1 μm. We use this value in our simulations.

2.3.2. Collision Velocity

We consider Brownian motion, radial and azimuthal drift, vertical settling, and turbulence as sources of the collision velocity, and give the collision velocity Δv as the root sum square of these contributions,

$$\begin{equation} \Delta v = \sqrt{(\Delta v_B)^2+(\Delta v_r)^2+(\Delta v_\phi)^2+(\Delta v_z)^2+(\Delta v_t)^2}, \end{equation} \tag{ 16 }$$

where Δv_B, Δv_r, Δv_ϕ, Δv_z, and Δv_t are the relative velocities induced by Brownian motion, radial drift, azimuthal drift, vertical settling, and turbulence, respectively.

The Brownian-motion-induced velocity is given by

$$\begin{equation} \Delta v_B = \sqrt{\frac{\pi m_1m_2}{8(m_1+m_2)k_{\rm B}T}}, \end{equation} \tag{ 17 }$$

where m₁ and m₂ are the masses of the two colliding aggregates.

The relative velocity due to radial drift is given by Δv_r = |v_r(t_{s, 1}) − v_r(t_{s, 2})|, where t_{s, 1} and t_{s, 2} are the stopping times of the colliding aggregates, and v_r is the radial velocity given by Equation (6). Similarly, the relative velocity due to differential azimuthal motion is given by Δv_ϕ = |v'_ϕ(t_{s, 1}) − v'_ϕ(t_{s, 2})|, where

$$\begin{equation} v^{\prime }_\phi = -\frac{\eta v_K}{1+(\Omega t_s)^2} \end{equation} \tag{ 18 }$$

is the deviation of the azimuthal velocity from the local Kepler velocity (Adachi et al. 1976; Weidenschilling 1977; Nakagawa et al. 1986). Here, we have neglected the backreaction from dust to gas as we did for the radial velocity (see Sections 2.3.2 and 5.2.1).

For the differential settling velocity, we assume Δv_z = |v_z(t_{s, 1}) − v_z(t_{s, 2})|, where

$$\begin{equation} v_z = - \frac{\Omega ^2 t_sz}{1+\Omega t_s}. \end{equation} \tag{ 19 }$$

Equation (19) reduces to the terminal settling velocity v_z = −Ω²t_sz in the strong coupling limit Ωt_s ≪ 1, and to the amplitude of the vertical oscillation velocity at Ωt_s ≫ 1 (Brauer et al. 2008a).

For the turbulence-driven relative velocity, we use an analytic formula for Kolmogorov turbulence derived by Ormel & Cuzzi (2007, their Equation (16)). This analytic formula has three limiting expressions (Equations (27)–(29) of Ormel & Cuzzi 2007):

$$\begin{equation} \Delta v_t \approx \left\lbrace \begin{array}{@{}ll}\delta v_g {\rm Re}_t^{1/4}\Omega |t_{s,1} - t_{s,2}|, & t_{s,1} \ll t_\eta, \\ [3mm] (1.4\dots 1.7)\times \delta v_g\sqrt{\Omega t_{s,1}}, & t_\eta \ll t_{s,1} \ll \Omega ^{-1}, \\ \delta v_g \left({\displaystyle \frac{1}{1+\Omega t_{s,1}}} + {\displaystyle \frac{1}{1+\Omega t_{s,2}}} \right)^{1/2}, & \Omega t_{s,1} \gg 1, \end{array} \right. \end{equation} \tag{ 20 }$$

where Re_t is the turbulent Reynolds number, t_η = Re^−1/2_tt_L is the turnover time of the smallest eddies, and the numerical prefactor (1.4...1.7) in the second equality depends on the ratio between the stopping times, t_{s, 2}/t_{s, 1}. The turbulent Reynolds number is given by Re_t = D_g/ν_mol, where ν_mol = v_thλ_mfp/2 is the molecular viscosity. For t_{s, 1} ∼ t_{s, 2}, the maximum induced velocity is Δv_t ≈ δv_g, which is reached when Ωt_{s, 1} ≈ 1.

When two colliding aggregates belong to the Epstein regime and their stopping times are much shorter than t_η(≪ Ω⁻¹), the relative velocity driven by sedimentation and turbulence is approximately proportional to the difference between the mass-to-area ratios m/A of two colliding aggregates. In this case, as pointed out by Okuzumi et al. (2011a), the dispersion of the mass-to-area ratio becomes important for fractal aggregates of d_f ≲ 2, since the mean mass-to-area ratio of the aggregates approaches a constant, and hence the difference in $m/ {\overline{A}} (m)$ vanishes. To take into account the dispersion effect, we evaluate the differential mass-to-area ratio as $|\Delta (m/A)|^2 = |m_1/ {\overline{A}} _1-m_2/ {\overline{A}} _2|^2 + \epsilon ^2[(m_1/ {\overline{A}} _1)^2+(m_2/ {\overline{A}} _2)^2]$, where ${\overline{A}} _j = {\overline{A}} (m_j)$ (j = 1, 2) are the mean projected area of aggregates with mass m_j (see Section 2.3.1) and is the standard deviation of the mass-to-area ratio divided by the mean (for the derivation, see the Appendix of Okuzumi et al. 2011a). We assume = 0.1 in accordance with the numerical estimate by Okuzumi et al. (2011a).

We solve Equations (1) and (8) numerically with an explicit time-integration scheme and a fixed-bin method. The radial domain is taken to be outside the snow line, 3 AU ⩽ r ⩽ 150 AU, discretized into 100 rings with an equal logarithmic width Δln (r[AU]) = (ln 150 − ln 3)/100. The advection terms are calculated by the spatially first-order upwind scheme. We impose the outflow and zero-flux boundary conditions at the innermost and outermost radii (r = 3 AU and 150 AU), respectively; thus, the total dust mass inside the domain is a decreasing function of time. Our numerical results are unaffected by the choice of the boundary condition at the outermost radius, since dust growth at this location is too slow to cause appreciable radial drift within the calculated time. The coagulation terms are calculated by the method given by Okuzumi et al. (2009). Specifically, at the center of each radial ring we divide the mass coordinate into linearly spaced bins m_k = km₀(k = 1, 2, ..., N_bd) for m ⩽ N_bdm₀ and logarithmically spaced bins $m_k = m_{k-1}10^{1/N_{{\rm bd}}} (k=N_{{\rm bd}}\,{+}\,1,\dots)$ for m > N_bdm₀, where N_bd is an integer. We adopt N_bd = 40; as shown by Okuzumi et al. (2009), the calculation results converge well as long as N_bd ⩾ 40. The time increment Δt is adjusted at every time step so that the fractional decreases in ${\cal N}$ and $V{\cal N}$ fall below 0.5 (i.e., $\Delta t < - 0.5(\partial \ln {\cal N}/\partial t)^{-1}$ and $\Delta t < - 0.5(\partial \ln V{\cal N}/\partial t)^{-1}$) at all bins.

To begin with, we show the result of compact aggregation. In this simulation, we fixed the internal density ρ_int ≡ m/V of the aggregates to the material density ρ₀ = 1.4 g cm⁻³, and solved only the evolutionary equation for the radial size distribution ${\cal N}(r,m)$ (Equation (1)), as in previous studies (e.g., Brauer et al. 2008a). Figure 2 shows the snapshots of the radial size distribution at different times. Here, the distribution is represented by the dust surface density per decade of aggregate mass, $\Delta \Sigma _d/\Delta \log m \equiv \ln (10) m^2 {\cal N}(r,m)$. At each orbital radius, dust growth proceeds without significant radial drift until the stopping time of the aggregates reaches Ωt_s ∼ 0.1 (the dotted lines in Figure 2). However, as the aggregates grow, the radial drift becomes faster and faster, and further growth becomes limited only along the line Ωt_s ∼ 0.1 on the r–m plane. Figure 3 shows the evolution of the total dust surface density $\Sigma _d \equiv \int m{\cal N} dm = \int (\Delta \Sigma _d/\Delta \log m)d\log m$. We see that a significant amount of dust has been lost from the planet-forming region r ≲ 30 AU within 10⁵ yr. In this region, the dust surface density⁴ scales as r⁻¹, and hence the dust-to-gas surface density ratio ∝r⁻¹/Σ_g∝r^1/2 decreases toward the central star.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 4 shows the evolution of the dust size distribution observed at r = 5 AU. Here, in order to characterize the typical aggregate size at each evolutionary stage, we introduce the weighted average mass 〈m〉_m defined by

$$\begin{equation} \langle m \rangle _m \equiv \frac{\int m^2{\cal N}dm}{\int m{\cal N} dm} = \frac{1}{\Sigma _d}\int m \frac{\Delta \Sigma _d}{\Delta \log m}d\log m. \end{equation} \tag{ 21 }$$

The weighted average mass approximately corresponds to the aggregate mass at the peak of ΔΣ_d/Δlog m (see, e.g., Ormel et al. 2007; Okuzumi et al. 2011a). In Figure 4, the weighted average mass at each time is indicated by the short vertical line. At r = 5 AU, the growth–drift equilibrium is reached at t ≈ 4000 yr, and the typical size of the aggregates is 〈m〉_m ≈ 500 g in mass (≈4 cm in radius, ≈0.07Ω⁻¹ in stopping time). Note that the final aggregate radius is much smaller than the mean free path λ_mfp of gas molecules (the dashed arrow in Figure 4), which means that the gas drag onto the aggregates is determined by Epstein's law. As we will see in the following, porosity evolution allows aggregates to reach the Stokes drag regime at much smaller Ωt_s.

Zoom In Zoom Out Reset image size

Now we show how porosity evolution affects dust evolution. Here, we solve the evolutionary equation for V(r, m) (Equation (8)) simultaneously with that for ${\cal N}(r,m)$ (Equation (1)). The result is shown in Figure 5, which displays the snapshots of the aggregate size distribution ΔΣ_d/Δlog m and internal density ρ_int = m/V at different times t as a function of r and m. The evolution of the total dust surface density Σ_d is shown in Figure 6.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The left four panels of Figure 5 show how the radial size distribution evolves in the porous aggregation. At t < 10³ yr, the evolution is qualitatively similar to that in compact aggregation (Section 3.1). However, in later stages, the evolution is significantly different. We observe that aggregates in the inner region of the disk (r < 10 AU) undergo rapid growth and eventually overcome the radial drift barrier lying at Ωt_s ∼ 1 (dashed lines in Figure 5) within t ∼ 10⁴ yr. At this stage, the radial profile of the total dust surface Σ_d is hardly changed from the initial profile, as seen in Figure 6. In the outer region (r > 10 AU), aggregates do drift inward before they reach Ωt_s ∼ 1 as in the compact aggregation model. However, unlike in the compact aggregation, the inward drift results in the pileup of dust materials in the inner region (r ≈ 4–9 AU) rather than the loss of them from outside the snow line (see Figure 6). This occurs because most of the drifting aggregates are captured by aggregates that have already overcome the drift barrier. As a result, the dust-to-gas mass ratio in the inner regions is enhanced by a factor of several in 10⁵ yr.

The right four panels of Figure 5 show the evolution of the internal density ρ_int = m/V as a function of r and m. The first thing to note is that the dust particles grow into low-density objects at every location until their internal density reaches ρ_int ∼ 10⁻⁵–10⁻³ g cm⁻³. In this stage, the internal density decreases as ρ_int ≈ (m/m₀)^−1/2ρ₀, meaning that the dust particles grow into fractal aggregates with the fractal dimension d_f ≈ 2 (Okuzumi et al. 2009). The fractal growth generally occurs in early growth stages where the impact energy is too low to cause collisional compression, i.e., E_imp ≪ E_roll (e.g., Blum 2004; Ormel et al. 2007; Zsom et al. 2010). At m ∼ 10⁻⁴–10⁻⁶ g, the fractal growth stage terminates, followed by the stage where collisional compression becomes effective (E_imp ≫ E_roll). In this late stage, the internal density decreases more slowly or is kept at a constant value depending on the orbital radius. We will examine the density evolution in more detail in Section 3.2.2.

Figure 7 shows the evolution of the mass distribution function at r = 5 AU during t ≈ 1200–2500 yr. The evolution of the weighted average mass 〈m〉_m is shown in Figure 8. It is seen that the acceleration of the growth begins when the aggregate size a exceeds the mean free path of gas molecules, λ_mfp (the dashed arrow in Figure 7). This suggests that the acceleration is due to the change in the aerodynamical property of the aggregates. At a ≈ λ_mfp, the gas drag onto the aggregates begins to obey Stokes' law. In the Stokes regime, the stopping time t_s of aggregates increases rapidly with size (see Section 2.2). This causes the quick growth of the aggregates since the relative velocity between aggregates increases with t_s (as long as Ωt_s < 1). As a result of the growth acceleration, the aggregates grow from a ≈ λ_mfp to Ωt_s ≈ 1 within 300 yr, which is short enough for them to break through the radial drift barrier.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The decrease in the internal density plays an important role in the growth acceleration. More precisely, the low internal density allows the aggregates to reach a ≈ λ_mfp at early growth stages, i.e., at small Ωt_s. In fact, the growth acceleration was not observed in the compact aggregation, since the aggregate size is smaller than the mean free path at all Ωt_s < 1 (see Figure 4). A more rigorous explanation for this will be given in Section 4.

3.2.1. Projectile Mass Distribution

As noted in Section 2.3.1, our porosity change model has only been tested for collisions between similar-sized aggregates. To check the validity of using this model, we introduce the projectile mass distribution function (Okuzumi et al. 2009):

$$\begin{equation} C_{m_t}(m_p) \equiv \frac{m_p K(m_p,m_t){\cal N}(m_p)}{\int _0^{m_t} m_p^{\prime } K(m_p^{\prime },m_t){\cal N}(m_p^{\prime })dm_p^{\prime }}, \quad m_p \leqslant m_t, \end{equation} \tag{ 22 }$$

which is normalized so that $\int _0^{m_t} C_{m_t}(m_p) dm_p = 1$. The denominator of $C_{m_t}(m_p)$ is equal to the growth rate t⁻¹_grow ≡ dln m_t/dt of a target having mass m_t (see Okuzumi et al. 2009). Hence, the quantity $C_{m_t}(m_p) dm_p$ measures the contribution of projectiles within mass range [m_p, m_p + dm_p] to the growth of the target.

Figure 9 shows the projectile mass distribution per unit ln m_p, $m_pC_{m_t}(m_p)$, for targets with mass m_t = 〈m〉_m at r = 5 AU and at different t. We see that the growth of the m_t = 〈m〉_m target is dominated by projectiles within a mass range 0.1m_t ≲ m_p ⩽ m_t. In fact, the projectile mass distribution integrated over 0.1m_t ⩽ m_p ⩽ m_t exceeds 50% for all the cases presented in Figure 9. This means that the growth of aggregates is indeed dominated by collisions with similar-sized ones as required by the limitation of our porosity model. This is basically a consequence of the fact that the aggregate mass distribution ΔΣ_d/Δlog m is peaked around the target mass m ∼ 〈m〉_m (see Figure 7). The mass ratio m_p/m_t at the peak of $m_pC_{m_t}(m_p)$ reflects the size dependence of the turbulence-driven relative velocity Δv_t, which is the main source of the collision velocity in our simulation. At t ≲ 2000 yr (〈m〉_m ≲ 10³ g), the dominant projectile mass is lower than m_t(= 〈m〉_m), since both the target and projectiles are tightly coupled to turbulence (i.e., t_s(m_t), t_s(m_p) < t_η) and hence Δv_t vanishes at equal-sized collisions (see the first expression of Equation (20)). At t ≳ 2000 yr (〈m〉_m ≳ 10³ yr), the target decouples from small turbulent eddies (t_s(m_t) > t_η). This results in a shift of the dominant collision mode to m_p ≈ m_t because Δv_t no longer vanishes at equal-sized collisions (see the second expression of Equation (20)).

Zoom In Zoom Out Reset image size

3.2.2. Density Evolution History

To see the density evolution history in detail, we plot in Figure 10 the temporal evolution of the weighted average mass 〈m〉_m and the internal density of aggregates with mass m = 〈m〉_m at orbital radii r = 5 AU and 20 AU.

Zoom In Zoom Out Reset image size

As mentioned above, dust particles initially grow into fractal aggregates of d_f ≈ 2 until the impact energy E_imp becomes comparable to the rolling energy E_roll. With this information, one can analytically estimate the aggregate size at which the fractal growth terminates. Our simulation shows that the collision velocity between the fractal aggregates is approximately given by the turbulence-driven velocity in the strong-coupling limit (Equation (20) with t_s ≪ t_η). Assuming that the collisions mainly occur between aggregates of similar sizes (see Section 3.2.1), the reduced mass and the collision velocity are roughly given by m/2 and δv_gRe^1/4_tΩt_s, respectively. In addition, we use the fact that fractal aggregates with d_f ≈ 2 have a mass-to-area ratio which is comparable to their constituent monomers. This means that the stopping time of the aggregates is as short as the monomers and is hence given by Epstein's law. Thus, the impact energy is approximated as

$$\begin{equation} E_{\rm imp} \approx \frac{m}{4}\Delta v_t^2 \approx \frac{3}{8}m \biggl ({\displaystyle \frac{\delta v_g {\rm Re}_t^{1/4}\Omega }{\rho _g v_{\rm th}}} \biggr) ^2 \biggl ({\displaystyle \frac{3m}{4A}} \biggr) ^2. \end{equation} \tag{ 23 }$$

Furthermore, using the definitions for ρ_g, v_th, and Re_t, we have ρ_gv_th = (2/π)Σ_gΩ and Re_t = α_DΣ_gσ_mol/(2m_g) for the midplane. Substituting them into Equation (23) and using $\delta v_g =\sqrt{\alpha _D} c_s$ and m/A ≈ m₀/(πa²₀) = 4ρ₀a₀/3, we obtain

$$\begin{equation} E_{\rm imp} \approx \frac{3\pi ^2}{32\sqrt{2}}\alpha _D^{3/2} m c_s^2 \biggl ({\displaystyle \frac{\Sigma _g \sigma _{\rm mol}}{m_g}} \biggr) ^{1/2} \biggl ({\displaystyle \frac{\rho _0a_0}{\Sigma _g}} \biggr) ^2. \end{equation} \tag{ 24 }$$

Thus, the impact energy is proportional to the mass. We define the rolling mass m_roll by the condition E_imp = E_roll. Using Equation (24), the rolling mass is evaluated as

$$\begin{eqnarray} m_{\rm roll} &\approx & \frac{32\sqrt{2}}{3\pi ^2} \frac{E_{\rm roll}}{c_s^2\alpha _D^{3/2}} \biggl ({\displaystyle \frac{m_g}{\Sigma _g \sigma _{\rm mol}}} \biggr) ^{1/2} \biggl ({\displaystyle \frac{\Sigma _g}{\rho _0a_0}} \biggr) ^2 \nonumber \\ &\sim & 10^{-4} \,{\rm g} \biggl ({\displaystyle \frac{\alpha _D}{10^{-3}}} \biggr) ^{-3/2} \biggl ({\displaystyle \frac{T}{100 \,{\rm K}}} \biggr) ^{-1} \biggl ({\displaystyle \frac{\Sigma _g}{100 \,{\rm g\, cm^{-2}}}} \biggr) ^{3/2} \nonumber \\ &&\times \biggl ({\displaystyle \frac{F_{\rm roll}}{10^{-3} \,{\rm dyn}}} \biggr) \biggl ({\displaystyle \frac{\rho _0}{1\, {\rm g\, cm^{-3}}}} \biggr) ^{-2} \biggl ({\displaystyle \frac{a_0}{0.1 \,\mbox{$\mu $m}}} \biggr) ^{-1}, \end{eqnarray} \tag{ 25 }$$

where we have used that E_roll = (πa₀/2)F_roll (see Section 2.3.1). Using the relations a ≈ (m/m₀)^1/2a₀ and ρ_int ≈ (m/m₀)^−1/2ρ₀ for d_f ≈ 2 aggregates, the corresponding radius and internal density are found to be

$$\begin{equation} a_{\rm roll} \sim 1 \,{\rm cm} \biggl ({\displaystyle \frac{m_{\rm roll}}{10^{-4} \,{\rm g}}} \biggr) ^{1/2}, \end{equation} \tag{ 26 }$$

$$\begin{equation} \rho _{\rm int,\, roll} \sim 10^{-5}\, {\rm g\, cm^{-3}} \biggl ({\displaystyle \frac{m_{\rm roll}}{10^{-4} \,{\rm g}}} \biggr) ^{-1/2}. \end{equation} \tag{ 27 }$$

The triangles in Figure 10 mark the rolling mass at r = 5 AU and 20 AU predicted by Equation (25). The analytic prediction well explains when the decrease in ρ_int terminates.

The density evolution is more complicated at m > m_roll, where collisional compression is no longer negligible (i.e., E_imp > E_roll). At r = 5 AU, the internal density is approximately constant until the stopping time reaches Ωt_s = 1, and then decreases as ρ_int∝m^−1/5. At r = 20 AU, by contrast, the density is kept nearly constant until m ∼ 10² g (a ∼ 10² cm), and then decreases as ρ_int∝m^−1/8.

As shown below, the density histories mentioned above can be directly derived from the porosity change recipe we adopted. Let us assume again that aggregates grow mainly through collisions with similar-sized ones (m₁ ≈ m₂ and V₁ ≈ V₂). In this case, the evolution of ρ_int at E_imp ≫ E_roll is approximately given by Equation (14). Furthermore, we neglect the term (2V^5/6₁)⁻⁴ in Equation (14) assuming that the impact energy is sufficiently large (which is true as long as Ωt_s < 1; see below). Given these assumptions, the internal density of aggregates after collision, ρ_int = 2m₁/V_{1 + 2}, is approximately given by

$$\begin{equation} \rho _{\rm int} \approx \biggl ({\displaystyle \frac{3}{5}} \biggr) ^{3/2} \biggl ({\displaystyle \frac{E_{\rm imp}}{N_{1+2}bE_{\rm roll}}} \biggr) ^{3/10}N_{1+2}^{-1/5}\rho _0, \end{equation} \tag{ 28 }$$

where N_{1 + 2} = 2m₁/m₀. Since the impact energy E_imp ≈ m₁(Δv)²/4 is proportional to N_{1 + 2}(Δv)², Equation (28) implies that

$$\begin{equation} \rho _{\rm int} \propto (\Delta v)^{3/5} m^{-1/5}, \end{equation} \tag{ 29 }$$

where we have dropped the subscript for mass for clarity. Equation (29) gives the relation between ρ_int and m if we know how the impact velocity depends on them. Explicitly, if Δv∝m^βρ^γ_int, Equation (29) leads to

$$\begin{equation} \rho _{\rm int} \propto m^{(3\beta -1)/(5-3\gamma)}. \end{equation} \tag{ 30 }$$

In our simulation, the main source of the relative velocity is turbulence. The turbulence-driven velocity depends on t_s as Δv_t∝t_s at t_s ≪ t_η and $\Delta v_t \propto \sqrt{t_s}$ at t_η ≪ t_s ≪ t_L(= Ω⁻¹) (see Equation (20)). As found from Equation (4), the stopping time depends on ρ_int and m as t_s∝m/A∝m/a²∝m^1/3ρ^2/3_int in the Epstein regime (a ≪ λ_mfp) and as t_s∝ma/A∝m^2/3ρ^1/3_int in the Stokes regime (a ≫ λ_mfp). Using these relations with Equation (30), we find four regimes for density evolution,

$$\begin{equation} \rho _{\rm int} \propto \left\lbrace \begin{array}{@{}ll}m^{0}, & a\ll \lambda _{\rm mfp} \quad {\rm and}\ t_s \ll t_\eta, \\ m^{1/4}, & a\gg \lambda _{\rm mfp} \quad {\rm and}\ t_s \ll t_\eta, \\ m^{-1/8}, & a\ll \lambda _{\rm mfp} \quad {\rm and}\ t_\eta \ll t_s \ll t_L, \\ m^{0}, & a\gg \lambda _{\rm mfp} \quad {\rm and}\ t_\eta \ll t_s \ll t_L. \\ \end{array} \right. \end{equation} \tag{ 31 }$$

The circles, diamonds, and square in Figure 10 mark the size at which a = λ_mfp (i.e., t^(Ep)_s ∼ t_s^(St)), t_s = t_η, and Ωt_s = 1, respectively. At r = 5 AU, the sizes at which a = λ_mfp and t_s = t_η nearly overlap, and hence only two velocity regimes t_s = t^(Ep)_s ≪ t_η and t_η ≪ t_s = t^(St)_s ≪ t_L are effectively relevant. For both cases, Equation (31) predicts flat density evolution. At r = 20 AU, there is a stage in which t_s ≫ t_η and a ≪ λ_mfp, for which Equation (31) predicts ρ_int∝m^−1/8. These predictions are in agreement with what we see in Figure 10.

Equation (28) does not apply to the density evolution at Ωt_s > 1, where the collision velocity no longer increases and hence collisional compression becomes less and less efficient as the aggregates grow. However, if we go back to Equation (14) and assume that the impact energy E_imp is sufficiently small, we obtain V_{1 + 2} ≈ 2^6/5V₁, or equivalently V_{1 + 2}/m^6/5_{1 + 2} ≈ V₁/m^6/5₁, where m_{1 + 2} = 2m₁ is the aggregate mass after a collision. This implies that V/m^6/5 is kept constant during the growth, i.e., V∝m^6/5, and hence we have ρ_int = m/V∝m^−1/5. This is consistent with the density evolution at Ωt_s > 1 seen in the upper panel of Figure 10.

In this study, we have assumed that the stopping time t_s obeys Stokes' law whenever a ≳ λ_mfp. In reality, Stokes' law applies only when the particle Reynolds number (the Reynolds number of the gas flow around the particle) ${\rm Re}_p \equiv 2a|{\bm v}_d - {\bm v}_g|/\nu _{\rm mol}$ is less than unity, where $|{\bm v}_d - {\bm v}_g|$ is the gas–dust relative velocity. When Re_p ≳ 1, i.e., when the particle becomes so large and/or the gas–dust relative velocity becomes so high, the stopping time becomes dependent on the particle velocity (see, e.g., Weidenschilling 1977). In this subsection, we discuss how this effect affects our conclusion.

In general, the stopping time at a ≫ λ_mfp can be written as

$$\begin{equation} t_s = \frac{2m}{C_D\rho _g|{\bm v}_d - {\bm v}_g|A}, \end{equation} \tag{ 45 }$$

where C_D is a dimensionless coefficient that depends on Re_p. Stokes' law, which applies when Re_p ≪ 1, is given by C_D = 24/Re_p. In the opposite limit, Re_p ≫ 1, the drag coefficient C_D approaches a constant value (typically of order unity; e.g., C_D ≈ 0.5 for a sphere with 10³ ≲ Re_p ≲ 10⁵), which is known as Newton's friction law. Thus, in the Newton regime, the stopping time depends on the particle velocity, unlike in the Stokes regime. In this case, one has to calculate the stopping time and particle velocity simultaneously since the particle velocity in turn depends on the stopping time.

In the previous sections, we have ignored the Newton regime to avoid the above-mentioned complexity. However, it is easy to calculate the growth timescale in the Newton regime for given Ωt_s, for which the gas–dust relative velocity can be known in advance. Below, we show that the Newton drag sets the minimum value of $t_{\rm grow}|_{\Omega t_s = 1}$ (Equation (35)) for given orbital radius and internal density, which was not taken into account in Section 4. At the midplane, Equation (45) can be rewritten as $\Omega t_s = (2\sqrt{2\pi }/C_D)(c_s/|{\bm v}_d - {\bm v}_g|)m/(\Sigma _g A)$, where we have used that $\rho _g = \Sigma _g\Omega /(\sqrt{2\pi }c_s)$. When Ωt_s = 1, the gas–dust relative velocity is dominated by the dust radial velocity v_r = −ηv_K, so we can set $|{\bm v}_d - {\bm v}_g| \approx \eta v_K$. Thus, at the midplane, we obtain a relation

$$\begin{equation} \frac{(\rho _{\rm int}a)_{\Omega t_s=1}}{\Sigma _g} \approx \frac{3C_D}{8\sqrt{2\pi }}\frac{\eta v_K}{c_s} \approx 0.07 \biggl ({\displaystyle \frac{C_D}{0.5}} \biggr) \frac{\eta v_K}{c_s}, \end{equation} \tag{ 46 }$$

where we have used that m/A = 4ρ_inta/3. If C_D reaches a constant, ${(\rho _{\rm int}a)_{\Omega t_s=1}}/{\Sigma _g}$ no longer depends on aggregate properties. Putting this equation into Equation (35), we have

$$\begin{equation} t_{\rm grow}|_{\Omega t_s = 1} \approx 0.3 \biggl ({\displaystyle \frac{\Sigma _d/\Sigma _g}{0.01}} \biggr) ^{-1} \biggl ({\displaystyle \frac{C_D}{0.5}} \biggr) \biggl ({\displaystyle \frac{\eta v_K/c_s}{0.08}} \biggr) t_K. \end{equation} \tag{ 47 }$$

When C_D = 24/Re_p, Equation (47) reduces to the equation for the Stokes drag (Equation (42)), where $t_{\rm grow}|_{\Omega t_s=1}$ decreases with increasing aggregate size a. However, when Re_p becomes so large that C_D reaches a constant value, $t_{\rm grow}|_{\Omega t_s=1}$ no longer decreases with increasing a. Thus, we find that the Newton drag sets the minimum value of $t_{\rm grow}|_{\Omega t_s=1}$. For our disk model, in which Σ_d/Σ_g = 0.01 and ηv_K/c_s = 0.08(r/5 AU)^1/4, the minimum growth timescale is ≈0.2–0.3(C_D/0.5)t_K at r ≈ 3–10 AU.

Since the Newton drag regime was ignored in our model, the growth rate of aggregates was overestimated at high Re_p. As seen in the lower panel of Figure 11, the growth timescale $t_{\rm grow}|_{\Omega t_s=1}$ for the ρ_int = 10⁻⁵ g cm⁻³ aggregates falls below the minimum possible value given by Equation (47) at r ≲ 7 AU. This implies that dust growth is somewhat artificially accelerated in our simulation presented in Section 3.2. However, this artifact is not the reason why porous aggregates grow across the radial drift barrier in the simulation. Indeed, the drift timescale $t_{\rm grow}|_{\Omega t_s=1}$ is ≈40t_K at these orbital radii, and hence the minimum growth timescale still satisfies the condition for breaking through the drift barrier, Equation (38) (see Section 4). Thus, highly porous aggregates are still able to break through the radial drift barrier even if Newton's law at high particle Reynolds numbers is taken into account.

In summary, we have shown that Newton's friction law (C_D ≈ constant) at high particle Reynolds numbers sets a floor value for the growth timescale at Ωt_s = 1. In the numerical simulation presented in Section 3.2, the neglect of the Newton drag regime causes artificial acceleration of the growth of Ωt_s ≳ 1 aggregates. However, comparison with the drift timescale shows that the floor value of $t_{\rm grow}|_{\Omega t_s=1}$ is sufficiently small for dust to grow across Ωt_s = 1. Therefore, the deviation from Stokes' law at high particle Reynolds numbers has little effect on the successful breakthrough of the radial drift barrier observed in our simulation.

So far we have neglected the frictional backreaction from dust to gas when determining the velocities of dust aggregates (Equations (6) and (18)). Here, we discuss the validity of this assumption.

5.2.1. Effect on the Equilibrium Drift Velocity

Frictional backreaction generally modifies the equilibrium velocities of both gas and dust. The equilibrium velocities in the presence of the backreaction are derived by Tanaka et al. (2005) for arbitrary dust size distribution. The result shows that the radial and azimuthal velocities v_r and v_ϕ = v'_ϕ + v_K of dust particles with stopping time t_s are given by

$$\begin{equation} v_r = \frac{1}{1+(\Omega t_s)^2}v_{g,r} + \frac{2\Omega t_s}{1+(\Omega t_s)^2}v^{\prime }_{g,\phi }, \end{equation} \tag{ 48 }$$

$$\begin{equation} v^{\prime }_\phi = -\frac{\Omega t_s}{2[1+(\Omega t_s)^2]}v_{g,r} + \frac{1}{1+(\Omega t_s)^2}v^{\prime }_{g,\phi }, \end{equation} \tag{ 49 }$$

where

$$\begin{equation} v_{g,r} = \frac{2Y}{(1+X)^2+Y^2}\eta v_K, \end{equation} \tag{ 50 }$$

$$\begin{equation} v^{\prime }_{g,\phi } = -\frac{1+X}{(1+X)^2+Y^2}\eta v_K, \end{equation} \tag{ 51 }$$

are the radial and azimuthal components of the gas velocity relative to the local circular Keplerian motion, respectively, and

$$\begin{equation} X = \int \frac{\rho _d(m)}{\rho _g}\frac{1}{1+(\Omega t_s(m))^2} dm, \end{equation} \tag{ 52 }$$

$$\begin{equation} Y = \int \frac{\rho _d(m)}{\rho _g}\frac{\Omega t_s(m)}{1+(\Omega t_s(m))^2} dm, \end{equation} \tag{ 53 }$$

with ρ_d(m) being the spatial mass density of dust particles per unit aggregate mass.⁵ In the limit of X, Y → 0, the gas velocities approach v_{g, r} → 0 and v'_{g, ϕ} → −ηv_K, and hence Equations (48) and (49) reduce to Equations (6) and (18), respectively. Thus, the dimensionless quantities X and Y measure the significance of the frictional backreaction. As found from the integrands in Equations (6) and (18), the backreaction is non-negligible when the local dust-to-gas mass ratio exceeds unity and the aggregates dominating the dust mass couple tightly to the gas.

To test the effect of fractional backreaction, we have also simulated porous aggregation using Equations (48) and (49) instead of Equations (6) and (18) for the aggregate velocities. However, it is found that the effect of backreaction is so small that the resulting dust evolution is hardly distinguishable from that presented in Section 3. The upper panel of Figure 13 shows the temporal evolution of the gas velocities v_{g, r} and v'_{g, ϕ} observed in this simulation as a function of the weighted average mass 〈m〉_m. We see that the observed gas velocities deviate at most only by 9 m s⁻¹ ≈ 0.17ηv_K from the velocities when the backreaction is absent (dotted lines). As a result, the inward velocity −v_r of aggregates with m = 〈m〉_m is decreased only by 15% even when Ωt_s(〈m〉_m) ≈ 1 (see the black solid curve in the lower panel of Figure 13). The above result can be understood in the following way. As found from the definitions of X and Y (Equations (52) and (53)), the effect of the backreaction is significant only when the density of dust coupled to the gas (Ωt_s ≲ 1) is comparable to or higher than the gas density. When Ωt_s(〈m〉_m) ≲ 1, the density of the coupled dust at the midplane is ${\lesssim}\Sigma _d/h_d|_{\Omega t_s=1} \sim \Sigma _d/(h_g\sqrt{\alpha }_D) \sim (0.01/\sqrt{\alpha _D})\rho _{g,{\rm mid}} \sim 0.3\rho _{g,{\rm mid}} \lesssim \rho _{g,{\rm mid}}$, where ρ_{g, mid} is the midplane gas density and we have used that $h_d|_{\Omega t_s=1} \sim \sqrt{\alpha _D}h_g \sim 0.03h_g$ (Equation (3)) and Σ_d/Σ_g ≈ 0.01 (the latter is true as long as Ωt_s(〈m〉_m) ≲ 1). When Ωt_s(〈m〉_m) ≳ 1, the dust density does exceed the gas density at the midplane, but the most part of the dust mass is now carried by decoupled (Ωt_s > 1) aggregates, which do not affect the gas motion.⁶ Thus, the density of coupled dust is always lower than the gas density, and hence the backreaction effect is insignificant at all times.

Zoom In Zoom Out Reset image size

Furthermore, the effect on the differential drift velocity is even less significant because the decreases in the inward velocities nearly cancel out. As an example, the gray solid and dotted curves in the lower panel of Figure 13 show the differential radial velocity Δv_r between aggregates of stopping times t_s = t_s(〈m〉_m) and 0.3t_s(〈m〉_m) obtained from the simulations with and without the backreaction, respectively. We see that the maximum values of |Δv_r|, which are reached when Ωt_s(〈m〉_m) ≈ 0.7, differ only by 5%. Therefore, the frictional backreaction from dust to gas hardly affects the drift-induced collision velocity between dust aggregates.

5.2.2. Streaming Instability

In this study, we have assumed that all aggregate collisions lead to sticking. This assumption breaks down if the collisional velocity is so high that the collision involves fragmentation and erosion. If the mass loss due to fragmentation and erosion is significant, it acts as an obstacle to planetesimal formation (the so-called fragmentation barrier; e.g., Brauer et al. 2008a). Here, we discuss the validity and possible limitations of this assumption.

Recent N-body simulations predict that very fluffy aggregates made of 0.1 μm sized icy particles experience catastrophic disruption at collision velocities Δv ≳ 35 m s⁻¹ (Wada et al. 2009). If a large aggregate grows mainly through collisions with similar-sized ones (which is true in our simulations; see Figure 9), the collision velocity at Ωt_s ≈ 1 is dominated by the turbulence-driven velocity $\Delta v_t \approx \delta v_g \approx \sqrt{\alpha _D}c_s$ (Section 2.3.2). If the disk is optically thin and moderately turbulent (α_D = 10⁻³) as in our model, the collision velocity is ≈21 m s⁻¹ at r = 5 AU, so catastrophic disruption is likely insignificant for such collisions. However, if turbulence is as strong as α_D = 10⁻², the collision velocity at r = 5 AU and Ωt_s = 1 goes up to 67 m s⁻¹. In protoplanetary disks, strong turbulence with α_D ≳ 10⁻² can be driven by magnetorotational instability (MRI; e.g., Balbus & Hawley 1998). If such strong turbulence exists, fragmentation becomes no more negligible even for icy aggregates. Besides, the collision velocity can become higher than the above estimate when a large aggregate collides with much smaller ones, since the collision velocity is then dominated by the radial drift motion. For example, the differential radial drift velocity between an Ωt_s = 1 aggregate and a much smaller one is as high as ≈ηv_K ≈ 56 m s⁻¹ in optically thin disks. At such a high velocity, erosion by small aggregates can also slow down the growth of Ωt_s ≈ 1 aggregates, although net growth might be possible (see, e.g., Teiser & Wurm 2009; Teiser et al. 2011).

On the other hand, resupply of small dust particles by fragmentation/erosion has positive effects on dust growth. First, small dust particles stabilize MRI-driven turbulence because they efficiently capture ionized gas particles and thereby reduce the electric conductivity of the gas (e.g., Sano et al. 2000). This process generally leads to the reduction of the gas random velocity (and hence the reduction of turbulence-induced collision velocity), especially when the magnetic fields threading the disk are weak (Okuzumi & Hirose 2011). In addition, small fragments enhance the optical thickness of the disk, and thus reduce the temperature of the gas in the interior of the disk (given that turbulence is stabilized there). Since the radial drift velocity is proportional to the gas temperature, this leads to the reduction of the drift-induced collision velocity. In the limit of large optical depths, the gas temperature is reduced by a factor ≈(h/r)^1/4 ≈ 0.5 near the midplane (Kusaka et al. 1970), resulting in the reduction of the drift-induced collision velocity to 28 m s⁻¹. These effects may help the growth of large aggregates beyond the fragmentation barrier.

The size of monomers is another key factor. Although we have assumed monodisperse monomers of a₀ = 0.1 μm, the size of interstellar dust particles ranges from nanometers to microns. It is suggested both theoretically (Chokshi et al. 1993; Dominik & Tielens 1997) and experimentally (Blum & Wurm 2008) that the threshold velocity for sticking is roughly inversely proportional to a₀. Thus, inclusion of larger monomers generally leads to a decrease in the sticking efficiency. However, it is not obvious whether aggregates composed of multi-sized interstellar particles are mechanically weaker or stronger than aggregates considered in this study. For example, if the monomer size distribution dn₀/da₀ obeys that of interstellar dust particles, dn₀/dlog a₀∝a^−5/2₀ (Mathis et al. 1977), the total mass of the aggregates is dominated by the largest ones (m₀∝a³₀ and hence m₀dn₀/dlog a₀∝a^1/2₀). Nevertheless, the existence of smaller monomers can still be important, since the binding energy per contact E_break is proportional to a^4/3₀ (Chokshi et al. 1993; Dominik & Tielens 1997) and hence the total binding energy tends to be dominated by the smallest ones (E_breakdn₀/dlog a₀∝a^−7/6₀). The net effect of multi-sized monomers needs to be clarified by future numerical as well as laboratory experiments.

Another issue concerning the growth efficiency of icy aggregates arises from sintering. Sintering is the redistribution of ice molecules on solid surfaces due to vapor transport and other effects. In this process, ice molecules tend to fill dipped surfaces (i.e., surfaces with negative surface curvature) since the equilibrium vapor pressure decreases with decreasing surface curvature. In an aggregate composed of equal-sized icy monomers, this process leads to growth of the monomer contact areas (Sirono 2011b) and consequently to enhancement of the aggregate's mechanical strength such as F_roll. Significant growth of the contact areas could cause the reduction of the aggregate's sticking efficiency since the dissipation of the collision energy through internal rolling/sliding motion could then be suppressed (Sirono 1999). Furthermore, if the monomers have different sizes, sintering leads to the evaporation of smaller monomers (having higher positive curvature), which may result in the breakup of the aggregate (Sirono 2011a). Therefore, sintering can prevent the growth of icy aggregates near the snow line where sintering proceeds rapidly. Sirono (2011b) shows that the timescale of H₂O sintering falls below 10³ yr in the region between the snow line (3 AU) and 7 AU for the radial temperature adopted in our study. This is comparable to the timescale on which submicron-sized icy particles grow into macroscopic objects in this region (see Figure 7). However, if the disk is passive and optically thick (Kusaka et al. 1970), no icy materials (including H₂O and CO₂) undergo rapid sintering at r ≳ 4 AU (Sirono 2011b). Moreover, the required high optical depth can be provided by tiny fragments that would result from the sintering-induced fragmentation itself. Consistent treatment of the two competing effects is necessary to precisely know the location where sintering is really problematic.

To summarize, whether icy aggregates survive catastrophic fragmentation and erosion crucially depends on the environment of the protoplanetary disks as well as on the size distribution of the aggregates and constituent monomers. However, we emphasize that icy aggregates can survive within a realistic range of disk conditions as explained above. Indeed, the range is much wider than that for rocky aggregates for which catastrophic disruption occurs at collision velocities as low as a few m s⁻¹ (Blum & Wurm 2008; Wada et al. 2009; Güttler et al. 2010). In order to precisely predict in what conditions icy aggregates overcome the fragmentation barrier, we need to take into account the mass loss due to fragmentation/erosion and the reduction of collision velocities due to the resupply of small particles in a self-consistent way. This will be done in our future work.

Aggregates observed in our simulation have very low internal densities. This is a direct consequence of the porosity model we adopted (Equation (15)). Here, we discuss the validity and limitations of our porosity model.

As mentioned in Section 2.3.1, our porosity change recipe at E_imp ≳ E_roll is based on head-on collision experiments of similar-sized aggregates. In our simulation, dust growth is indeed dominated by collisions with similar-sized aggregates (see Section 3.2.1), so our result is unlikely affected by the limitation of the porosity model regarding the size ratio. By contrast, the neglect of offset collisions may cause underestimation of the porosity increase, since the impact energy is spent for stretching rather than compaction at offset collisions (Wada et al. 2007; Paszun & Dominik 2009). If this is the case, then the breakthrough of the radial drift barrier can occur even outside 10 AU.

On the other hand, the formation of very low density dust aggregates is apparently inconsistent with the existence of massive and much less porous aggregates in our solar system. For example, comets, presumably the most primitive dust "aggregates" in the solar system, are expected to have mean internal densities of ρ_int ∼ 0.1 g cm⁻³ (e.g., Greenberg & Hage 1990). Since our porosity model does not explain the formation of such large and less porous "aggregates," some compaction mechanisms have yet to be determined.

One possibility is static compression due to gas drag and self-gravity. Although static compression is ignored in our porosity model, it can contribute to the compaction of aggregates that are massive or decoupled from the gas motion. For relatively compact (ρ_int ∼ 0.1 g cm⁻³) dust cakes made of micron-sized SiO₂ particles, static compaction is observed to occur at static pressure >100 Pa (Blum & Schräpler 2004; Güttler et al. 2009). By contrast, the static compression strength has not yet been measured so far for icy aggregates with very low internal densities (ρ_int ≪ 0.1 g cm⁻³). However, for future reference, it will be useful to estimate here the static pressures due to gas drag and self-gravity.

The ram pressure, the gas drag force per unit area, is given by $P_{\rm ram} = C_D\rho _g|{\bm v}_d-{\bm v}_g|^2/2$, where C_D is the drag coefficient and $|{\bm v}_d-{\bm v}_g|$ is the gas–dust relative speed (see Section 5.1). At Ωt_s ≳ 1, the gas–dust relative speed is approximately equal to ηv_K. Thus, assuming Newton's drag law C_D ∼ 1 for Ωt_s ≳ 1 aggregates (Section 5.1), the ram pressure at Ωt_s ≳ 1 is estimated as

$$\begin{equation} P_{\rm ram} \sim \rho _g(\eta v_K)^2 \sim 10^{-5} \biggl ({\displaystyle \frac{\rho _g}{10^{-11}\, {\rm g\, cm^{-3}}}} \biggr) \biggl ({\displaystyle \frac{\eta v_K}{50 \,{\rm m\, s^{-1}}}} \biggr) ^2 {\rm Pa} \end{equation} \tag{ 54 }$$

independently of aggregate properties. Thus, if the static compression strength of our high porous aggregates is lower than 10⁻⁵ Pa, compression of the aggregates will occur at Ωt_s ≳ 1 due to ram pressure.

The static pressure due to self-gravity is estimated from dimensional analysis as

$$\begin{equation} P_{\rm grav} \sim \frac{Gm^2}{a^4} \sim 10^{-7} \biggl ({\displaystyle \frac{m}{10^{10} \,{\rm g}}} \biggr) ^{2/3} \biggl ({\displaystyle \frac{\rho _{\rm int}}{10^{-5}\, {\rm g\, cm^{-3}}}} \biggr) ^{4/3} {\rm Pa}. \end{equation} \tag{ 55 }$$

For m ∼ 10¹⁰ g and ρ_int ∼ 10⁻⁵ g cm⁻³, which correspond to the Ωt_s = 1 aggregates observed in our simulation (Figure 5), the gravitational pressure is much weaker than the ram pressure. However, since P_grav∝m^2/3, compression due to self-gravity becomes important for much heavier aggregates. For example, if ρ_int ∼ 10⁻⁵(m/10¹⁰ g)^−1/5 g cm⁻³ as is the case for the Ωt_s ≳ 1 aggregates observed in our simulation, P_grav exceeds P_ram at m ∼ 10¹⁷ g, which is comparable to the mass of comet Halley. Moreover, since P_grav∝ρ^4/3_int, gravitational compaction will proceed in a runaway manner unless the static compression strength increases more rapidly than P_grav. Thus, static compression due to self-gravity may be a key to fill the gap between our high porous aggregates and more compact planetesimal-mass bodies in the solar system.