Dynamics of Porous Dust Aggregates and Gravitational Instability of Their Disk

We adopt the following surface densities of gas and dust (Weidenschilling 1977a; Hayashi 1981):

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{{\rm{g}}}=1700{f}_{{\rm{g}}}{\left(\displaystyle \frac{a}{1\mathrm{au}}\right)}^{-{\beta }_{{\rm{g}}}}\,{\rm{g}}\,{\mathrm{cm}}^{-2},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{{\rm{d}}}=\gamma {{\rm{\Sigma }}}_{{\rm{g}}},\end{eqnarray} \tag{ 2 }$$

where a is the distance from the central star, ${\beta }_{{\rm{g}}}$ is the power-law index, f_g is the ratio in the MMSN model, and γ is the dust-to-gas mass ratio. As the fiducial model, we consider the MMSN model (${\beta }_{{\rm{g}}}=3/2$ and $\gamma =0.018$) beyond the snowline (Hayashi 1981; Hayashi et al. 1985). We adopt the temperature profile

$$\begin{eqnarray}&&T={T}_{1}{\left(\displaystyle \frac{a}{1\mathrm{au}}\right)}^{-{\beta }_{{\rm{t}}}}\,{\rm{K}},\end{eqnarray} \tag{ 3 }$$

where T₁ is the temperature at $1\,\mathrm{au}$ and ${\beta }_{{\rm{t}}}$ is the power-law index. In the fiducial model, we adopt ${T}_{1}=120$ and ${\beta }_{{\rm{t}}}=3/7$ (Chiang & Youdin 2010). The isothermal sound velocity is calculated from the temperature as ${c}_{{\rm{s}}}=\sqrt{{k}_{{\rm{B}}}T/{m}_{{\rm{g}}}}$, where k_B is the Boltzmann constant and ${m}_{{\rm{g}}}=3.9\times {10}^{-24}\,{\rm{g}}$ is the mean molecular mass. The gas density at the disk midplane is ${\rho }_{{\rm{g}}}={{\rm{\Sigma }}}_{{\rm{g}}}/(\sqrt{2\pi }{c}_{{\rm{s}}}/{\rm{\Omega }})$, where ${\rm{\Omega }}=\sqrt{{{GM}}_{* }/{a}^{3}}$ is the Keplerian angular frequency and ${M}_{* }$ is the mass of the central star. Throughout this paper, we adopt ${M}_{* }={M}_{\odot }$, where ${M}_{\odot }$ is the solar mass. The mean free path of gas molecules is $l={m}_{{\rm{g}}}/{\sigma }_{{\rm{g}}}{\rho }_{{\rm{g}}}$, where ${\sigma }_{{\rm{g}}}=2\times {10}^{-15}\,{\mathrm{cm}}^{2}$ is the collisional cross section of gas molecules. The nondimensional radial pressure gradient is given as (Nakagawa et al. 1986)

$$\begin{eqnarray}\eta & = & -\,\displaystyle \frac{1}{2}{\left(\displaystyle \frac{{c}_{{\rm{s}}}}{a{\rm{\Omega }}}\right)}^{2}\displaystyle \frac{\partial \mathrm{log}({\rho }_{{\rm{g}}}{c}_{{\rm{s}}}^{2})}{\partial \mathrm{log}a}=(2.39{\beta }_{{\rm{g}}}+1.20{\beta }_{{\rm{t}}}+3.59)\\ & & \times {10}^{-4}\left(\displaystyle \frac{{T}_{1}}{120}\right){\left(\displaystyle \frac{a}{1\mathrm{au}}\right)}^{1-{\beta }_{{\rm{t}}}}.\end{eqnarray} \tag{ 4 }$$

For the fiducial model (${T}_{1}=120$, ${\beta }_{{\rm{g}}}=3/2$ and ${\beta }_{{\rm{t}}}=3/7$), η at 1 au is $0.77\times {10}^{-3}$.

2.2.1. Physical Properties

For simplicity, we will consider equal-mass dust aggregates with mass m_d. We assume that a dust aggregate consists of N monomers, where $N={m}_{{\rm{d}}}/{m}_{0}$. The monomer mass, density, and radius are m₀, ${\rho }_{0}$, and r₀, respectively, which satisfy ${m}_{0}=4\pi /3{\rho }_{0}{r}_{0}^{3}$. The gyration radius of a dust aggregate is given as (Mukai et al. 1992; Wada et al. 2008)

$$\begin{eqnarray}&&{r}_{{\rm{g}}}=\sqrt{\displaystyle \frac{1}{N}\displaystyle \sum _{k}{({{\boldsymbol{x}}}_{k}-{\boldsymbol{X}})}^{2}},\end{eqnarray} \tag{ 5 }$$

where ${{\boldsymbol{x}}}_{k}$ is the position of monomer k and ${\boldsymbol{X}}={\sum }_{k}{{\boldsymbol{x}}}_{k}/N$ is the position of the center of mass. The characteristic radius is defined by (Mukai et al. 1992)

$$\begin{eqnarray}&&{r}_{{\rm{c}}}=\sqrt{\displaystyle \frac{5}{3}}{r}_{{\rm{g}}}.\end{eqnarray} \tag{ 6 }$$

If the dust aggregates grow by ballistic cluster–particle aggregation (BCPA), each dust aggregate can be considered as a uniform sphere, where r_c corresponds to the physical radius. On the other hand, if the dust aggregates grow by ballistic cluster–cluster aggregation (BCCA), the dust aggregates are inhomogeneous with a fractal dimension of less than 3, and r_c corresponds to the maximum distance from the approximate center of mass (Okuzumi et al. 2009). We define the collisional cross section by using r_c. The characteristic volume and the mean internal density are calculated from r_c as follows:

$$\begin{eqnarray}&&{V}_{{\rm{c}}}=\displaystyle \frac{4\pi {r}_{{\rm{c}}}^{3}}{3},\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{\rho }_{\mathrm{int}}=\displaystyle \frac{{m}_{{\rm{d}}}}{{V}_{c}}.\end{eqnarray} \tag{ 8 }$$

The density of a dust aggregate with fractal dimension D is

$$\begin{eqnarray}&&{\rho }_{\mathrm{int}}={\left(\displaystyle \frac{{m}_{{\rm{d}}}}{{m}_{0}}\right)}^{1-3/D}{\rho }_{0}.\end{eqnarray} \tag{ 9 }$$

The fractal dimension of a BCCA cluster is $D\simeq 1.9$ (Mukai et al. 1992; Okuzumi et al. 2009), and thus the density of a BCCA cluster is ${\rho }_{\mathrm{BCCA}}={({m}_{{\rm{d}}}/{m}_{0})}^{-0.58}{\rho }_{0}$.

The interaction of dust aggregates with gas is characterized by the projected area A. Generally, the projected area A can be smaller than $\pi {r}_{{\rm{c}}}^{2}$. Okuzumi et al. (2009) obtained the projected area formula for both BCPA and BCCA clusters

$$\begin{eqnarray}&&A={\left(\displaystyle \frac{1}{{A}_{\mathrm{BCCA}}}+\displaystyle \frac{1}{\pi {r}_{{\rm{c}}}^{2}}-\displaystyle \frac{1}{\pi {r}_{{\rm{c}},\mathrm{BCCA}}^{2}}\right)}^{-1},\end{eqnarray} \tag{ 10 }$$

where A_BCCA is the projected area for the BCCA cluster (Minato et al. 2006)

$$\begin{eqnarray}\displaystyle \frac{{A}_{\mathrm{BCCA}}}{\pi {r}_{0}^{2}}\simeq \left\{\begin{array}{ll}12.5{N}^{0.685}\exp (-2.53/{N}^{0.0920}) & (N\lt 16)\\ (0.352N+0.566{N}^{0.862}), & (N\geqslant 16)\end{array}\right.,\end{eqnarray} \tag{ 11 }$$

and ${r}_{{\rm{c}},\mathrm{BCCA}}$ is the characteristic radius of the corresponding BCCA cluster,

$$\begin{eqnarray}&&\pi {r}_{{\rm{c}},\mathrm{BCCA}}^{2}\simeq {N}^{2/D}\pi {r}_{0}^{2}.\end{eqnarray} \tag{ 12 }$$

We define the area-equivalent radius from A as

$$\begin{eqnarray}&&{r}_{{\rm{A}}}=\sqrt{\displaystyle \frac{A}{\pi }}.\end{eqnarray} \tag{ 13 }$$

The area-equivalent radius is shown in Figure 1. If the density is sufficiently larger than the BCCA cluster density ${\rho }_{\mathrm{BCCA}}$, r_A is approximated by r_c. Since we investigate the evolution due to self-gravity compression, we consider the following range of parameters: ${m}_{{\rm{d}}}\gt {10}^{10}\,{\rm{g}}$ and ${\rho }_{\mathrm{int}}\gt {10}^{-5}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ (Kataoka et al. 2013b). In this parameter region, the density is much larger than ${\rho }_{\mathrm{BCCA}}$, and the corresponding fractal dimension is larger than 2. Thus, throughout this paper we assume ${r}_{{\rm{A}}}\simeq {r}_{{\rm{c}}}$. In short, we consider a dust aggregate to be a sphere with radius r_c, which has the usual mass–radius relation ${m}_{{\rm{d}}}=(4\pi /3){\rho }_{\mathrm{int}}{r}_{{\rm{c}}}^{3}$.

Zoom In Zoom Out Reset image size

2.2.2. Evolution

We consider the evolution of dust due to coagulation and self-gravity compression (Kataoka et al. 2013b). The compressive strength of a porous dust aggregate is (Kataoka et al. 2013a)

$$\begin{eqnarray}&&{P}_{\mathrm{comp}}\simeq \displaystyle \frac{{E}_{\mathrm{roll}}}{{r}_{0}^{3}}{\left(\displaystyle \frac{{\rho }_{\mathrm{int}}}{{\rho }_{0}}\right)}^{3},\end{eqnarray} \tag{ 14 }$$

where E_roll is the rolling energy (Dominik & Tielens 1997; Wada et al. 2007). Equilibrium is obtained when ${P}_{\mathrm{comp}}=P$ (Kataoka et al. 2013b), and from this we can calculate the equilibrium density. The self-gravity pressure is

$$\begin{eqnarray}&&{P}_{\mathrm{grav}}\simeq \displaystyle \frac{{{Gm}}_{{\rm{d}}}^{2}}{\pi {r}_{{\rm{c}}}^{4}},\end{eqnarray} \tag{ 15 }$$

and thus we obtain the equilibrium density (Kataoka et al. 2013b)

$$\begin{eqnarray}&&{\rho }_{\mathrm{eq}}\simeq 1.58{\left(\displaystyle \frac{{r}_{0}^{3}{\rho }_{0}^{3}G}{{E}_{\mathrm{roll}}}\right)}^{3/5}{m}_{{\rm{d}}}^{2/5}.\end{eqnarray} \tag{ 16 }$$

As the mass of the dust aggregate increases, its self-gravity pressure increases; this compresses it and increases its equilibrium density. In the fiducial model, we adopt ${E}_{\mathrm{roll}}=4.74\times {10}^{-9}\,\mathrm{erg}$, ${\rho }_{0}=1.0\,{\rm{g}}\,{\mathrm{cm}}^{-3}$, and ${r}_{0}=0.1\,\mu {\rm{m}}$.

Gravitational scattering among dust aggregates results in an increase in both ${\sigma }_{e}$ and ${\sigma }_{i}$ (Ida 1990; Stewart & Ida 2000). In Paper I, we adopted the simple stirring rate formula described by the Chandrasekhar relaxation time (Ida 1990), which is valid in the dispersion-dominated regime. Ohtsuki et al. (2002) examined the evolution of ${\sigma }_{e}$ and ${\sigma }_{i}$ using three-body integration, and they derived a semianalytic formula that is valid in both the dispersion-dominated and the shear-dominated regime. In this paper, we adopt their stirring rates:

$$\begin{eqnarray}&&{\left(\displaystyle \frac{d{\sigma }_{e}^{2}}{{dt}}\right)}_{\mathrm{grav}}=\displaystyle \frac{{a}^{2}{\rm{\Omega }}{h}^{4}{{\rm{\Sigma }}}_{{\rm{d}}}}{4{m}_{{\rm{d}}}}{P}_{\mathrm{VS}}\end{eqnarray} \tag{ 17 }$$

and

$$\begin{eqnarray}&&{\left(\displaystyle \frac{d{\sigma }_{i}^{2}}{{dt}}\right)}_{\mathrm{grav}}=\displaystyle \frac{{a}^{2}{\rm{\Omega }}{h}^{4}{{\rm{\Sigma }}}_{{\rm{d}}}}{4{m}_{{\rm{d}}}}{Q}_{\mathrm{VS}},\end{eqnarray} \tag{ 18 }$$

where h is the reduced Hill radius $h={(2{m}_{{\rm{d}}}/3{M}_{* })}^{1/3}$ and

$$\begin{eqnarray}{P}_{\mathrm{VS}} & = & 73\displaystyle \frac{\mathrm{log}(10{{\rm{\Lambda }}}^{2}/{\tilde{\sigma }}_{e,{\rm{r}}}^{2}+1)}{10{{\rm{\Lambda }}}^{2}/{\tilde{\sigma }}_{e,{\rm{r}}}^{2}}+\mathrm{log}({{\rm{\Lambda }}}^{2}+1)\\ & & \times \displaystyle \frac{72}{\pi {\tilde{\sigma }}_{e,{\rm{r}}}{\tilde{\sigma }}_{i,{\rm{r}}}}{\displaystyle \int }_{0}^{1}\displaystyle \frac{5{K}_{\lambda }-12(1-{\lambda }^{2}){E}_{\lambda }/(1+3{\lambda }^{2})}{\beta +(1/\beta -\beta ){\lambda }^{2}}d\lambda ,\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}{Q}_{\mathrm{VS}} & = & (4{\tilde{\sigma }}_{i,{\rm{r}}}^{2}+0.2{\tilde{\sigma }}_{e,{\rm{r}}}^{3}{\tilde{\sigma }}_{i,{\rm{r}}})\displaystyle \frac{\mathrm{log}(10{{\rm{\Lambda }}}^{2}{\tilde{\sigma }}_{e,{\rm{r}}}+1)}{10{{\rm{\Lambda }}}^{2}{\tilde{\sigma }}_{e,{\rm{r}}}}\\ & & +\,\mathrm{log}({{\rm{\Lambda }}}^{2}+1)\displaystyle \frac{72}{\pi {\tilde{\sigma }}_{e,{\rm{r}}}{\tilde{\sigma }}_{i,{\rm{r}}}}\\ & & \times {\displaystyle \int }_{0}^{1}\displaystyle \frac{{K}_{\lambda }-12{\lambda }^{2}{E}_{\lambda }/(1+3{\lambda }^{2})}{\beta +(1/\beta -\beta ){\lambda }^{2}}d\lambda ,\end{eqnarray} \tag{ 20 }$$

where $\beta ={\tilde{\sigma }}_{i,{\rm{r}}}/{\tilde{\sigma }}_{e,{\rm{r}}}$, ${K}_{\lambda }=K(\sqrt{3(1-{\lambda }^{2})}/2)$, ${E}_{\lambda }\,=E(\sqrt{3(1-{\lambda }^{2})}/2)$, and ${\rm{\Lambda }}={\tilde{\sigma }}_{i,{\rm{r}}}({\tilde{\sigma }}_{e,{\rm{r}}}^{2}+{\tilde{\sigma }}_{i,{\rm{r}}}^{2})/12$. The first and second terms correspond to the low-velocity and high-velocity limits, respectively. The reduced relative eccentricity ${\tilde{\sigma }}_{e,{\rm{r}}}$ and the inclination ${\tilde{\sigma }}_{i,{\rm{r}}}$ are ${\tilde{\sigma }}_{e,{\rm{r}}}=\sqrt{2}{\sigma }_{e}/h$ and ${\tilde{\sigma }}_{i,{\rm{r}}}=\sqrt{2}{\sigma }_{i}/h$, respectively (Nakazawa et al. 1989). The functions $K{(k)={\int }_{0}^{\pi /2}(1-{k}^{2}{\sin }^{2}\theta )}^{-1/2}d\theta$ and $E{(k)={\int }_{0}^{\pi /2}(1-{k}^{2}{\sin }^{2}\theta )}^{1/2}d\theta$ are the complete elliptic integrals of the first and second kind, respectively, where k is the modulus.

We assume that the average change in ${\sigma }_{e}$ for a dust collision is ${\sigma }_{e}^{2}\to {C}_{\mathrm{col}}{\sigma }_{e}^{2}$. We adopt ${C}_{\mathrm{col}}=1/2$, which corresponds to perfect accretion (Inaba et al. 2001). We discuss the effect of the imperfect accretion on the dust aggregate growth in Section 5.2. Using the nondimensional collision rate P_col (Nakazawa et al. 1989), the evolution equations for ${\sigma }_{e}$ and ${\sigma }_{i}$ are written as

$$\begin{eqnarray}&&{\left(\displaystyle \frac{d{\sigma }_{e}^{2}}{{dt}}\right)}_{\mathrm{col}}=-{C}_{\mathrm{col}}{P}_{\mathrm{col}}{h}^{2}{a}^{2}\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{d}}}}{{m}_{{\rm{d}}}}{\rm{\Omega }}{\sigma }_{e}^{2},\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&{\left(\displaystyle \frac{d{\sigma }_{i}^{2}}{{dt}}\right)}_{\mathrm{col}}=-{C}_{\mathrm{col}}{P}_{\mathrm{col}}{h}^{2}{a}^{2}\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{d}}}}{{m}_{{\rm{d}}}}{\rm{\Omega }}{\sigma }_{i}^{2}.\end{eqnarray} \tag{ 22 }$$

For the low-velocity regime, where ${\tilde{\sigma }}_{e},{\tilde{\sigma }}_{i}\lt 0.2$, P_col is independent of ${\tilde{\sigma }}_{e}$ and ${\tilde{\sigma }}_{i}$ (Ida & Nakazawa 1989; Inaba et al. 2001), and

$$\begin{eqnarray}&&{P}_{\mathrm{col}}\simeq {P}_{\mathrm{low}}=11.3\sqrt{\tilde{r}},\end{eqnarray} \tag{ 23 }$$

where ${\tilde{\sigma }}_{e}={\sigma }_{e}/h$ and ${\tilde{\sigma }}_{i}={\sigma }_{i}/h$ are the reduced eccentricity and inclination, respectively, and $\tilde{r}=2{r}_{{\rm{c}}}/{ha}$. For the medium-velocity regime, where $0.2\lt {\tilde{\sigma }}_{e},{\tilde{\sigma }}_{i}\lt 2$ (Ida & Nakazawa 1989; Inaba et al. 2001), P_col depends on ${\tilde{\sigma }}_{i,{\rm{r}}}$:

$$\begin{eqnarray}&&{P}_{\mathrm{col}}\simeq {P}_{\mathrm{med}}=\displaystyle \frac{{\tilde{r}}^{2}}{4\pi {\tilde{\sigma }}_{i,{\rm{r}}}}\left(17.3+\displaystyle \frac{232}{\tilde{r}}\right).\end{eqnarray} \tag{ 24 }$$

For the high-velocity regime, where $2\lt {\tilde{\sigma }}_{e},{\tilde{\sigma }}_{i}$, P_col is (Greenzweig & Lissauer 1992)

$$\begin{eqnarray}{P}_{\mathrm{col}}\simeq {P}_{\mathrm{high}} & = & \displaystyle \frac{{\tilde{r}}^{2}}{2\pi }\left(8{\displaystyle \int }_{0}^{1}d\lambda \displaystyle \frac{{\beta }^{2}{E}_{\lambda }}{{({\beta }^{2}+(1-{\beta }^{2}){\lambda }^{2})}^{2}}\right.\\ & & +\left.\displaystyle \frac{48}{\tilde{r}{\tilde{\sigma }}_{e,{\rm{r}}}^{2}}{\displaystyle \int }_{0}^{1}d\lambda \displaystyle \frac{{K}_{\lambda }}{{\beta }^{2}+(1-{\beta }^{2}){\lambda }^{2}}\right),\end{eqnarray} \tag{ 25 }$$

where the second term indicates the effect of gravitational focusing. Inaba et al. (2001) proposed the following formula for the general nondimensional collision rate:

$$\begin{eqnarray}&&{P}_{\mathrm{col}}=\min ({P}_{\mathrm{med}},{({P}_{\mathrm{high}}^{-2}+{P}_{\mathrm{low}}^{-2})}^{-1/2}).\end{eqnarray} \tag{ 26 }$$

We set ${P}_{\mathrm{col}}={P}_{\mathrm{high}}$ for a large random velocity, such as ${\tilde{\sigma }}_{i}\gt 10$. Note that this formula does not take into account the gas drag effect, which may alter the collision rate. In the present model, when the gas drag is strong, turbulent stirring causes the random velocity to be high. In this case, the collision rate formula is described by the geometrical cross section.

3.3.1. Gas Drag

Because of the hydrodynamic gas drag, ${\sigma }_{e}$ and ${\sigma }_{i}$ decrease with time as follows (Adachi et al. 1976; Inaba et al. 2001):

$$\begin{eqnarray}&&{\left(\displaystyle \frac{d{\sigma }_{e}^{2}}{{dt}}\right)}_{\mathrm{gas},\mathrm{drag}}=-\displaystyle \frac{2}{{t}_{{\rm{s}},e}}{\sigma }_{e}^{2},\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&{\left(\displaystyle \frac{d{\sigma }_{i}^{2}}{{dt}}\right)}_{\mathrm{gas},\mathrm{drag}}=-\displaystyle \frac{2}{{t}_{{\rm{s}},i}}{\sigma }_{i}^{2},\end{eqnarray} \tag{ 28 }$$

where ${t}_{{\rm{s}},e}$ and ${t}_{{\rm{s}},i}$ are the damping timescales of ${\sigma }_{e}$ and ${\sigma }_{i}$,

$$\begin{eqnarray}&&{t}_{{\rm{s}},e}=\displaystyle \frac{2{m}_{{\rm{d}}}}{\pi {C}_{{\rm{D}}}{r}_{{\rm{c}}}^{2}{\rho }_{{\rm{g}}}{v}_{{\rm{K}}}{\left(\tfrac{9{E}^{2}}{4\pi }{\sigma }_{e}^{2}+\tfrac{1}{\pi }{\sigma }_{i}^{2}+\tfrac{9}{4}{\eta }^{2}\right)}^{1/2}},\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&{t}_{{\rm{s}},i}=\displaystyle \frac{4{m}_{{\rm{d}}}}{\pi {C}_{{\rm{D}}}{r}_{{\rm{c}}}^{2}{\rho }_{{\rm{g}}}{v}_{{\rm{K}}}{\left(\tfrac{{E}^{2}}{\pi }{\sigma }_{e}^{2}+\tfrac{4}{\pi }{\sigma }_{i}^{2}+{\eta }^{2}\right)}^{1/2}},\end{eqnarray} \tag{ 30 }$$

where $E=E(\sqrt{3/4})$ is the elliptic integral of the second kind of argument $\sqrt{3/4}$, ${v}_{{\rm{K}}}=a{\rm{\Omega }}$ is the Keplerian velocity, and C_D is the nondimensional drag coefficient. We also adopted the corrections discussed by Kary et al. (1993) and Inaba et al. (2001).

The relative velocity between gas and dust is (Adachi et al. 1976)

$$\begin{eqnarray}&&u\simeq {({v}_{\mathrm{ran}}^{2}+{\eta }^{2}{v}_{{\rm{K}}}^{2})}^{1/2},\end{eqnarray} \tag{ 31 }$$

where ${v}_{\mathrm{ran}}={\left(\tfrac{5}{8}{\sigma }_{e}^{2}+\tfrac{1}{2}{\sigma }_{i}^{2}\right)}^{1/2}{v}_{{\rm{K}}}$ is the random velocity. Using the relative velocity given by Equation (31), we define the stopping time

$$\begin{eqnarray}&&{t}_{{\rm{s}}}=\displaystyle \frac{2{m}_{{\rm{d}}}}{\pi {C}_{{\rm{D}}}{r}_{{\rm{c}}}^{2}{\rho }_{{\rm{g}}}u},\end{eqnarray} \tag{ 32 }$$

which is nearly equal to ${t}_{{\rm{s}},e}$ and ${t}_{{\rm{s}},i}$.

The gas drag law changes with r_c (e.g., Adachi et al. 1976). If ${r}_{{\rm{c}}}\gtrsim l$, we use the Stokes drag or the Newton drag. For a low Reynolds number ($\mathrm{Re}\,\ll \,{10}^{3}$), the drag coefficient is approximated by ${C}_{{\rm{D}}}\simeq 24/\mathrm{Re}$ (Stokes drag), where $\mathrm{Re}\,=\,2{r}_{{\rm{c}}}u/\nu$. The viscosity ν is given by $\nu ={v}_{\mathrm{th}}l/2$, where ${v}_{\mathrm{th}}=\sqrt{8/\pi }{c}_{{\rm{s}}}$ is the thermal velocity. For a high Reynolds number (${10}^{3}\lt \mathrm{Re}\,\lt \,2\times {10}^{5}$), the drag coefficient is almost constant, ${C}_{{\rm{D}}}\simeq 0.4\mbox{--}0.5$ (Newton drag). If ${r}_{{\rm{c}}}\lesssim l$, we use the Epstein drag. Thus, we adopt the drag coefficient formula (Brown & Lawler 2003)

$$\begin{eqnarray}&&{C}_{{\rm{D}}}\\ =\,\left\{\begin{array}{ll}\displaystyle \frac{8{v}_{\mathrm{th}}}{3u} & ({r}_{{\rm{c}}}\lt 9l/4)\\ \displaystyle \frac{0.407}{1+8710/\mathrm{Re}}+\displaystyle \frac{24}{\mathrm{Re}}(1+0.150{\mathrm{Re}}^{0.681}) & ({r}_{{\rm{c}}}\gt 9l/4)\end{array}\right..\end{eqnarray} \tag{ 33 }$$

3.3.2. Turbulent Stirring

The random velocity of dust aggregates increases owing to the gas drag from the turbulent velocity field as (Paper I)

$$\begin{eqnarray}&&\displaystyle \frac{{{dv}}_{\mathrm{ran}}^{2}}{{dt}}=\displaystyle \frac{2{\tau }_{e}{v}_{{\rm{t}}}^{2}{\rm{\Omega }}}{S({\tau }_{e}+S)},\end{eqnarray} \tag{ 34 }$$

where $S={\rm{\Omega }}{t}_{{\rm{s}}}$ is the Stokes number, ${v}_{{\rm{t}}}=\sqrt{\alpha }{c}_{{\rm{s}}}$ is the magnitude of the turbulent velocity, α is the dimensionless turbulence strength (Cuzzi et al. 2001), and ${\tau }_{e}{{\rm{\Omega }}}^{-1}$ is the eddy turnover time. In the fiducial model, we adopt ${\tau }_{e}=1$ (Youdin 2011). For the isotropic turbulent velocity, the heating rate of ${\sigma }_{e}$ would be twice as large as that of ${\sigma }_{i}$, and thus we adopt the following formulae (Kobayashi et al. 2016):

$$\begin{eqnarray}&&{\left(\displaystyle \frac{d{\sigma }_{e}^{2}}{{dt}}\right)}_{\mathrm{turb},\mathrm{stir}}=\displaystyle \frac{4{\tau }_{e}{v}_{{\rm{t}}}^{2}{\rm{\Omega }}}{3{v}_{{\rm{K}}}^{2}S({\tau }_{e}+S)},\end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray}&&{\left(\displaystyle \frac{d{\sigma }_{i}^{2}}{{dt}}\right)}_{\mathrm{turb},\mathrm{stir}}=\displaystyle \frac{2{\tau }_{e}{v}_{{\rm{t}}}^{2}{\rm{\Omega }}}{3{v}_{{\rm{K}}}^{2}S({\tau }_{e}+S)}.\end{eqnarray} \tag{ 36 }$$

3.3.3. Turbulent Scattering

The fluctuations in gas density due to turbulence result in gravitational scattering of the dust aggregates. Okuzumi & Ormel (2013) derived the stirring rate of ${\sigma }_{e}$ due to this effect:

$$\begin{eqnarray}&&{\left(\displaystyle \frac{d{\sigma }_{e}^{2}}{{dt}}\right)}_{\mathrm{turb},\mathrm{grav}}={C}_{\mathrm{turb}}\alpha {\left(\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{g}}}{a}^{2}}{{M}_{* }}\right)}^{2}{\rm{\Omega }},\end{eqnarray} \tag{ 37 }$$

where C_turb is a nondimensional coefficient. From the semianalytical discussion, the nondimensional coefficient was obtained as (Okuzumi & Hirose 2011; Gressel et al. 2012; Okuzumi & Ormel 2013)

$$\begin{eqnarray}&&{C}_{\mathrm{turb}}=\displaystyle \frac{0.94{ \mathcal L }}{{(1+4.5{H}_{\mathrm{res},0}/H)}^{2}},\end{eqnarray} \tag{ 38 }$$

where H is the gas scale height, ${H}_{\mathrm{res},0}$ is the vertical half-width of the dead zone for the magnetorotational instability, and ${ \mathcal L }$ is a nondimensional saturation limiter that is less than or equal to unity. For simplicity, we set ${ \mathcal L }=1$. In the fiducial model, we adopt ${H}_{\mathrm{res},0}=H$, which leads to ${C}_{\mathrm{turb}}=3.1\times {10}^{-2}$. In Section 5, we discuss the effect of C_turb. The stirring rate of ${\sigma }_{i}$ is

$$\begin{eqnarray}&&{\left(\displaystyle \frac{d{\sigma }_{i}^{2}}{{dt}}\right)}_{\mathrm{turb},\mathrm{grav}}={\epsilon }_{i}^{2}\displaystyle \frac{d{\sigma }_{e}^{2}}{{dt}},\end{eqnarray} \tag{ 39 }$$

where ${\epsilon }_{i}$ is a nondimensional coefficient and is smaller than unity (Yang et al. 2012). We adopt ${\epsilon }_{i}=0.1$ (Kobayashi et al. 2016).

Considering the above processes, we obtain the evolution equations for ${\sigma }_{e}$ and ${\sigma }_{i}$:

$$\begin{eqnarray}\displaystyle \frac{d{\sigma }_{e}^{2}}{{dt}} & = & {\left(\displaystyle \frac{d{\sigma }_{e}^{2}}{{dt}}\right)}_{\mathrm{grav}}+{\left(\displaystyle \frac{d{\sigma }_{e}^{2}}{{dt}}\right)}_{\mathrm{col}}+{\left(\displaystyle \frac{d{\sigma }_{e}^{2}}{{dt}}\right)}_{\mathrm{gas},\mathrm{drag}}\\ & & +{\left(\displaystyle \frac{d{\sigma }_{e}^{2}}{{dt}}\right)}_{\mathrm{turb},\mathrm{stir}}+{\left(\displaystyle \frac{d{\sigma }_{e}^{2}}{{dt}}\right)}_{\mathrm{turb},\mathrm{grav}},\end{eqnarray} \tag{ 40 }$$

$$\begin{eqnarray}\displaystyle \frac{d{\sigma }_{i}^{2}}{{dt}} & = & {\left(\displaystyle \frac{d{\sigma }_{i}^{2}}{{dt}}\right)}_{\mathrm{grav}}+{\left(\displaystyle \frac{d{\sigma }_{i}^{2}}{{dt}}\right)}_{\mathrm{col}}+{\left(\displaystyle \frac{d{\sigma }_{i}^{2}}{{dt}}\right)}_{\mathrm{gas},\mathrm{drag}}\\ & & +{\left(\displaystyle \frac{d{\sigma }_{i}^{2}}{{dt}}\right)}_{\mathrm{turb},\mathrm{stir}}+{\left(\displaystyle \frac{d{\sigma }_{i}^{2}}{{dt}}\right)}_{\mathrm{turb},\mathrm{grav}}.\end{eqnarray} \tag{ 41 }$$

If the relaxation of ${\sigma }_{e}$ and ${\sigma }_{i}$ is sufficiently fast, we can obtain the equilibrium values of ${\sigma }_{e}$ and ${\sigma }_{i}$ from $d{\sigma }_{e}^{2}/{dt}=0$ and $d{\sigma }_{i}^{2}/{dt}=0$. In Section 4, we discuss the GI using the equilibrium values of ${\sigma }_{e}$ and ${\sigma }_{i};$ however, we note that the validity of this treatment is not trivial. The nonequilibrium effect is discussed in Section 5.

For a dust layer that is perfectly coupled with an incompressible fluid, the Roche density ${\rho }_{{\rm{R}}}$ is often used for the GI condition (Sekiya 1983; Yamoto & Sekiya 2004). In this case, a buckling mode develops, and dust-rich gas clumps form. Planetesimals may form in these dense clumps (Sekiya 1983; Cuzzi et al. 2008). However, in this paper, we focus on large dust aggregates that have a large Stokes number. In other words, dust aggregates can move relative to the gas, and dust can collapse through the gas, where the Roche criterion may not be applicable. In this case, Toomre's Q better describes the GI condition (Toomre 1964). We approximate the dust layer as a fluid with a finite velocity dispersion and calculate Q as

$$\begin{eqnarray}&&Q=\displaystyle \frac{{v}_{x}{\rm{\Omega }}}{{C}_{{\rm{T}}}G{{\rm{\Sigma }}}_{{\rm{d}}}},\end{eqnarray} \tag{ 42 }$$

where ${v}_{x}={({v}_{{\rm{K}}}^{2}{\sigma }_{e}^{2}/2)}^{1/2}$ is the radial component of the dust velocity dispersion and C_T is a nondimensional constant. The values of C_T are ${C}_{{\rm{T}}}=\pi$ for an ideal gas and ${C}_{{\rm{T}}}=3.36$ for collisionless particles (Toomre 1964); we adopted ${C}_{{\rm{T}}}=3.36$. In the majority of this paper, we adopt the condition for the GI as $Q\lt {Q}_{\mathrm{cr}}\simeq 2$ (Paper I).

For comparison, we also examine the Roche criterion. For the uniform dust layer perfectly coupled with the incompressible gas, the Roche density is (Sekiya 1983)

$$\begin{eqnarray}&&{\rho }_{{\rm{R}}}\simeq 0.6\displaystyle \frac{{M}_{* }}{{a}^{3}}.\end{eqnarray} \tag{ 43 }$$

We adopt this as the Roche density. Strictly speaking, the coefficient of the Roche density depends on the vertical structure and the equation of state of the dust layer. For instance, the coefficient for the Gaussian dust density distribution is 0.78 (Yamoto & Sekiya 2004). Furthermore, in this paper we do not consider the incompressible gas. Thus, Equation (43) is considered as an order-of-magnitude estimate of the Roche density for loosely coupling dust. From Equation (43) we define the nondimensional value Q_R (e.g., Youdin 2011),

$$\begin{eqnarray}&&{Q}_{{\rm{R}}}=\displaystyle \frac{{h}_{{\rm{d}}}}{{h}_{\mathrm{cr}}},\end{eqnarray} \tag{ 44 }$$

where ${h}_{{\rm{d}}}=a{\sigma }_{i}$ is the dust scale height and ${h}_{\mathrm{cr}}={{\rm{\Sigma }}}_{{\rm{d}}}/(\sqrt{\pi }{\rho }_{{\rm{R}}})$ is the critical scale height. If ${Q}_{{\rm{R}}}\lesssim 1$, the dust layer tends to be unstable in the sense of the Roche criterion.

We begin by considering the equilibrium state. We numerically calculate the equilibrium values of ${\sigma }_{e}$ and ${\sigma }_{i}$ by setting the right-hand sides of Equations (40) and (41) equal to zero. Then, from Equation (42), we calculate Q.

Figure 2 shows the equilibrium ${\sigma }_{e}$ and ${\sigma }_{i}$ for the fiducial model with $a=5\,\mathrm{au}$, $\alpha ={10}^{-3}$, and ${f}_{{\rm{g}}}=1$. In this model, ${\sigma }_{e}$ and ${\sigma }_{i}$ range from ${10}^{-5}$ to ${10}^{-3}$. Basically ${\sigma }_{e}$ and ${\sigma }_{i}$ have similar dependencies on m_d and ${\rho }_{\mathrm{int}}$. In most of the parameter regimes, ${\sigma }_{i}$ is smaller than ${\sigma }_{e}$. The ratio of ${\sigma }_{i}$ to ${\sigma }_{e}$ is discussed in detail below. Around the region where ${m}_{{\rm{d}}}\sim {10}^{15}\,{\rm{g}}$ and ${\rho }_{\mathrm{int}}\sim {10}^{-4}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$, ${\sigma }_{e}$ has the smallest value of about ${10}^{-5}$. The GI is likely to occur around this region in the fiducial model.

Zoom In Zoom Out Reset image size

We investigate Q on the m_d–${\rho }_{\mathrm{int}}$ plane. Figure 3 shows the result for the fiducial model. The minimum value of Q is at ${m}_{{\rm{d}}}=6.3\times {10}^{15}\,{\rm{g}}$ and ${\rho }_{\mathrm{int}}=1.4\times {10}^{-4}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$, where ${Q}_{\min }\simeq 0.6$, which is sufficiently small to allow the GI. Note that $Q\lt {Q}_{\mathrm{cr}}$ when ${\rho }_{\mathrm{int}}=1\times {10}^{-7}$ to $1\,{\rm{g}}\,{\mathrm{cm}}^{-3}$. Therefore, when the density is in the realistic range, the GI condition is inevitably satisfied during the evolution of the dust.

Zoom In Zoom Out Reset image size

The main heating and cooling mechanisms of ${\sigma }_{e}$ and ${\sigma }_{i}$ are shown on the m_d–${\rho }_{\mathrm{int}}$ plane for the fiducial model in Figure 4. In the region of low mass and low density, turbulent stirring is the dominant heating mechanism because the gas drag is strong. In the region of high mass and high density, turbulent scattering is the dominant heating mechanism, and in the region of high mass and low density, gravitational scattering is the dominant heating mechanism. The heating rate due to gravitational scattering is proportional to the scattering cross section, which is about ${({{Gm}}_{{\rm{d}}}/{v}_{\mathrm{ran}}^{2})}^{2}$ (e.g., Ida 1990). In the high-mass region, the random velocity is approximately given by the escape velocity because collisional damping and gravitational scattering are dominant. Thus, Figure 2 shows that in this region the random velocity decreases with decreasing ${\rho }_{\mathrm{int}}$. Therefore, gravitational scattering is stronger for the lower density if we fix the mass. The region of turbulent scattering for ${\sigma }_{i}$ is smaller than that for ${\sigma }_{e}$ because the scattering rate for ${\sigma }_{i}$ is smaller than it is for ${\sigma }_{e}$, as shown in Equation (39).

Zoom In Zoom Out Reset image size

We compare the Roche criterion with the Toomre criterion. In the fiducial model, there is no region for ${Q}_{{\rm{R}}}\lt 1$. As shown in Figure 3, there is the region for ${Q}_{{\rm{R}}}\lt 2$, but it is smaller than that for $Q\lt 2$. Thus, the Roche criterion is harder to satisfy than the Toomre criterion. The region for ${Q}_{{\rm{R}}}\lt 2$ relatively extends toward high ${\rho }_{\mathrm{int}}$. In this region, the main heating source is turbulent scattering. The heating rate of the inclination due to turbulent scattering is small, which leads to a thin dust layer. Thus, there the dust layer is more likely to be unstable in the sense of the Roche criterion. In the following we use the Toomre criterion, which is more optimistic.

We examine the dust aggregate dynamics in detail to clarify the physical background of the GI. Figure 5(a) shows the ratio of the random velocity to the surface escape velocity, ${v}_{\mathrm{esc}}=\sqrt{2{{Gm}}_{{\rm{d}}}/{r}_{{\rm{c}}}}$. If this ratio is less than unity, gravitational focusing is effective. In the low-mass region, the ratio is sufficiently larger than unity, which means that gravitational focusing is negligible. In this case, the collisional cross section is well approximated by the geometrical cross section, and the condition for runaway growth is not satisfied (e.g., Ohtsuki et al. 1993; Kokubo & Ida 1996). In the high-mass region, the escape velocity is comparable to the random velocity, and thus the collisional cross section is enhanced by a factor of about 2.

Zoom In Zoom Out Reset image size

When calculating u (Equation (31)), if ${v}_{\mathrm{ran}}/(\eta {v}_{{\rm{K}}})$ is small, then v_ran is negligible. Figure 5(b) shows ${v}_{\mathrm{ran}}/(\eta {v}_{{\rm{K}}})$. Other than when the mass and density are both high, this ratio is less than unity, and we can safely adopt the approximation $u\simeq \eta {v}_{{\rm{K}}}$.

The ratio of the random velocity to the Hill velocity, ${v}_{{\rm{H}}}={r}_{{\rm{H}}}{\rm{\Omega }}$, determines the regime of gravitational scattering, that is, the shear-dominated or dispersion-dominated regime (Weidenschilling 1989). Figure 5(c) shows ${v}_{\mathrm{ran}}/{v}_{{\rm{H}}}$. In the entire region, v_ran is larger than v_H, which indicates that the random velocity is dispersion dominated.

Figure 5(d) shows the ratio ${\sigma }_{i}/{\sigma }_{e}$, which is determined by the heating and cooling mechanisms. Figure 4 shows the main heating and cooling mechanisms on the m_d–${\rho }_{\mathrm{int}}$ plane. For the region where the main heating and cooling mechanism is gravitational scattering and collisions, respectively, ${\sigma }_{i}/{\sigma }_{e}$ is about 0.4–0.6, because gravitational scattering results in ${\sigma }_{i}/{\sigma }_{e}\simeq 0.5$ (Ida & Makino 1992). In the region with turbulent drag and gas drag, ${\sigma }_{i}/{\sigma }_{e}$ is larger than unity, and the damping rate of ${\sigma }_{e}$ due to gas drag is larger than that of ${\sigma }_{i}$. Thus, ${\sigma }_{i}$ is larger than ${\sigma }_{e}$. The ratio ${\sigma }_{i}/{\sigma }_{e}$ is about 0.8 where turbulent drag and collisions are dominant. In the region where turbulent scattering and collisions prevail, ${\sigma }_{i}/{\sigma }_{e}$ is less than 0.2, and the heating rate of ${\sigma }_{i}$ due to turbulent scattering is smaller than that of ${\sigma }_{e}$. This leads to a smaller value for ${\sigma }_{i}/{\sigma }_{e}$ (Yang et al. 2012; Kobayashi et al. 2016).

Figure 5(e) shows the drag coefficient C_D. There is no Epstein regime in Figure 5(e). The Epstein regime appears when we consider the small dust aggregates with ${m}_{{\rm{d}}}\lt {10}^{7}{\rm{g}}$. For $\mathrm{Re}\,\lt \,1$ (${C}_{{\rm{D}}}\gt 27.6$), the Stokes law is a good approximation. In the region where the mass is small and the density is large, the gas drag obeys the Stokes law. On the other hand, for massive dust aggregates, the gas drag obeys Newton's law. Thus, in this region, C_D is almost constant, at about 0.4–0.5.

The Stokes number is shown in Figure 5(f). For the low-mass, low-density region, S is less than unity, which means that the dust is well coupled with the gas. In the GI region, S is sufficiently larger than unity that the coupling does not occur.

We now consider the main heating and cooling processes shown in Figure 4. In the low-mass, high-density region, gas drag is the main cooling mechanism, while in the high-mass, low-density region, collisions dominate. Gas drag obeys the Stokes law in the low-mass, high-density region, as shown in Figure 5(e). In this case, the damping rate due to gas drag is $\propto {C}_{{\rm{D}}}{r}_{{\rm{c}}}^{2}{/m}_{{\rm{d}}}\propto {r}_{{\rm{c}}}{/m}_{{\rm{d}}}\propto {m}_{{\rm{d}}}^{-2/3}{\rho }_{\mathrm{int}}^{-1/3}$. Similarly, neglecting gravitational focusing, the damping rate due to collisions is $\propto {r}_{{\rm{c}}}^{2}{/m}_{{\rm{d}}}\propto {m}_{{\rm{d}}}^{-1/3}{\rho }_{\mathrm{int}}^{-2/3}$. Thus, for larger m_d and smaller ${\rho }_{\mathrm{int}}$, damping by collisions is more important than that by gas drag.

The region where gas drag is dominant for ${\sigma }_{e}$ is larger than that for ${\sigma }_{i}$. As shown in Figure 5(b), in most areas, ${\sigma }_{e},{\sigma }_{i}\lt \eta$. In this case, the damping rate due to gas drag can be approximated as ${(d{\sigma }_{e}^{2}/{dt})}_{\mathrm{drag}}\simeq -(3/{t}_{{\rm{s}}}){\sigma }_{e}^{2}$ and ${(d{\sigma }_{i}^{2}/{dt})}_{\mathrm{drag}}\,\simeq -(1/{t}_{{\rm{s}}}){\sigma }_{i}^{2}$. The damping time for ${\sigma }_{e}$ is ${t}_{{\rm{s}}}/3$, which is shorter than that for ${\sigma }_{i}$. On the other hand, the damping times due to collisions for ${\sigma }_{e}$ and ${\sigma }_{i}$ are the same. Thus, the gas drag region for ${\sigma }_{e}$ is larger than that for ${\sigma }_{i}$.

By a similar way to that used in Paper I, we derive the critical α for the existence of the GI region. First, we focus on the lower mass boundary of the GI region around ${\rho }_{\mathrm{int}}\simeq {10}^{-6}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$. In the low-density region, the main heating and cooling mechanisms are turbulent stirring and collisions, respectively. In this case, the equilibrium ${\sigma }_{e}$ is approximated by ${(d{\sigma }_{e}^{2}/{dt})}_{\mathrm{turb},\mathrm{stir}}+{(d{\sigma }_{e}^{2}/{dt})}_{\mathrm{col}}\simeq 0$. We neglect the second term in Equation (25) since ${v}_{\mathrm{ran}}\gt {v}_{\mathrm{esc}}$ and gravitational focusing is not effective, as shown in Figure 5(a). Furthermore, as shown in Figure 5(d), ${\sigma }_{i}/{\sigma }_{e}$ is almost constant in the low-density region. Thus, assuming ${\sigma }_{i}/{\sigma }_{e}=0.71$, we evaluate the first term of Equation (25) and obtain ${P}_{\mathrm{col}}\simeq 2.13{\tilde{r}}_{{\rm{c}}}^{2}$. Substituting this into Equation (21), we obtain the approximated ${(d{\sigma }_{e}^{2}/{dt})}_{\mathrm{col}}$. Figure 5(f) shows $S\gg 1$, where we obtain ${(d{\sigma }_{e}^{2}/{dt})}_{\mathrm{turb},\mathrm{stir}}\simeq 4{\tau }_{e}{v}_{{\rm{t}}}^{2}{\rm{\Omega }}/3{v}_{{\rm{K}}}^{2}{S}^{2}$ from Equation (35). We find that ${v}_{\mathrm{ran}}\lt \eta {v}_{{\rm{K}}}$ in Figure 5(b). Thus, from Equation (31), we approximate $u\simeq \eta {v}_{{\rm{K}}}$. Substituting this into Equation (32), we calculate S and finally obtain the approximated ${(d{\sigma }_{e}^{2}/{dt})}_{\mathrm{turb},\mathrm{stir}}.$ From Figure 5(e), Newton's drag is a good approximation; thus, we assume ${C}_{{\rm{D}}}\simeq 0.5$. Based on these approximations, we calculate ${\sigma }_{e}$ and obtain

$$\begin{eqnarray}&&Q\simeq 3.24\times {10}^{-2}\displaystyle \frac{{\alpha }^{1/2}\eta {\tau }_{e}^{1/2}{C}_{{\rm{D}}}{M}_{* }{{\rm{\Sigma }}}_{{\rm{g}}}}{{C}_{\mathrm{col}}^{1/2}{a}^{2}{m}_{{\rm{d}}}^{1/6}{\rho }_{\mathrm{int}}^{1/3}{{\rm{\Sigma }}}_{{\rm{d}}}^{3/2}}.\end{eqnarray} \tag{ 45 }$$

From $Q\lt {Q}_{\mathrm{cr}}$, we obtain the inequality

$$\begin{eqnarray}&&{m}_{{\rm{d}}}\gtrsim {m}_{\mathrm{low}}=1.16\times {10}^{-9}\displaystyle \frac{{\alpha }^{3}{\eta }^{6}{\tau }_{e}^{3}{C}_{{\rm{D}}}^{6}{M}_{* }^{6}{{\rm{\Sigma }}}_{{\rm{g}}}^{6}}{{C}_{\mathrm{col}}^{3}{Q}_{\mathrm{cr}}^{6}{a}^{12}{\rho }_{\mathrm{int}}^{2}{{\rm{\Sigma }}}_{{\rm{d}}}^{9}}.\end{eqnarray} \tag{ 46 }$$

Next, we focus on the upper mass boundary of the GI region around ${\rho }_{\mathrm{int}}\simeq {10}^{-1}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$. In the high-mass, high-density region, the main heating and cooling mechanisms are turbulent scattering and collisions, respectively. We obtain the equilibrium value for ${\sigma }_{e}$ from ${(d{\sigma }_{e}^{2}/{dt})}_{\mathrm{turb},\mathrm{grav}}+{(d{\sigma }_{e}^{2}/{dt})}_{\mathrm{col}}=0$. In evaluating ${(d{\sigma }_{e}^{2}/{dt})}_{\mathrm{col}}$, we adopt the same approximations described above. From $Q\lt {Q}_{\mathrm{cr}}$, we obtain the upper limit of m_d:

$$\begin{eqnarray}&&{m}_{{\rm{d}}}\lesssim {m}_{\mathrm{high}}=4.04\times {10}^{5}\displaystyle \frac{{C}_{\mathrm{col}}^{3}{Q}_{\mathrm{cr}}^{6}{{\rm{\Sigma }}}_{{\rm{d}}}^{9}}{{\alpha }^{3}{C}_{\mathrm{turb}}^{3}{\rho }_{\mathrm{int}}^{2}{{\rm{\Sigma }}}_{{\rm{g}}}^{6}}.\end{eqnarray} \tag{ 47 }$$

In Figure 4, m_low and m_high are plotted, and they are in rough agreement with the numerical results.

For the dust aggregates to trigger the GI, it is necessary that ${m}_{\mathrm{low}}\lt {m}_{\mathrm{high}}$. If this inequality is not satisfied, the GI region does not exist. Thus, we derive the following condition for the GI:

$$\begin{eqnarray}&&\alpha \lt {\alpha }_{\mathrm{cr},1}=2.65\times {10}^{2}\displaystyle \frac{{C}_{\mathrm{col}}{Q}_{\mathrm{cr}}^{2}{a}^{2}{{\rm{\Sigma }}}_{{\rm{d}}}^{3}}{\eta {\tau }_{e}^{1/2}{C}_{\mathrm{turb}}^{1/2}{C}_{{\rm{D}}}{M}_{* }{{\rm{\Sigma }}}_{{\rm{g}}}^{2}}.\end{eqnarray} \tag{ 48 }$$

Using the disk model, we rewrite ${\alpha }_{\mathrm{cr},1}$ as

$$\begin{eqnarray}{\alpha }_{\mathrm{cr},1} & = & 8.75\times {10}^{-3}{\tau }_{e}^{-1/2}{f}_{{\rm{g}}}{\left(\displaystyle \frac{\gamma }{0.018}\right)}^{3}{\left(\displaystyle \frac{a}{1\mathrm{au}}\right)}^{{\beta }_{{\rm{t}}}-{\beta }_{{\rm{g}}}+1}\\ & & \times {\left(\displaystyle \frac{{T}_{1}}{120}\right)}^{-1}{\left(\displaystyle \frac{{M}_{* }}{{M}_{\odot }}\right)}^{-1}{\left(\displaystyle \frac{2.39{\beta }_{{\rm{g}}}+1.20{\beta }_{{\rm{t}}}+3.59}{7.70}\right)}^{-1}\\ & & \times {\left(\displaystyle \frac{{C}_{\mathrm{turb}}}{3.1\times {10}^{-2}}\right)}^{-1/2}.\end{eqnarray} \tag{ 49 }$$

As shown in Section 5.1, this condition agrees well with the numerical results.

In the fiducial model (${\beta }_{{\rm{g}}}=3/2$ and ${\beta }_{{\rm{t}}}=3/7$), we find a weak dependence on a: $\propto {a}^{-1/14}$. Note that for the optically thin MMSN model (${\beta }_{{\rm{t}}}=1/2$ and ${\beta }_{{\rm{g}}}=3/2$) (Hayashi 1981; Hayashi et al. 1985), this dependence vanishes completely. Thus, for realistic disk models, the dependence on a is generally weak.

First, assuming the equilibrium random velocity, we examine the dependence of the critical α for the GI on the disk parameters. We adopt the dust evolution described in Section 2.2.2 and check whether its track crosses the GI region. Figure 6(a) shows the dependence on f_g. The GI is more likely with larger f_g and smaller α. The preference for small α occurs because turbulence is the main heating mechanism. From Equation (37), when the gas surface density (f_g) is large, the heating rate due to turbulent scattering is also large. In our model, the dust surface density also increases with f_g, since the dust-to-gas density ratio γ is fixed. When the dust surface density is large, this leads to strong self-gravity and a high collision frequency. Among these competing effects, the increases in self-gravity and collision frequency are dominant. Thus, the GI is more likely to occur when f_g is large.

Zoom In Zoom Out Reset image size

We found that the condition for the existence of the GI region (Equation (49)) agrees well with the numerical results for ${f}_{{\rm{g}}}\gt 0.3$, while it overestimates α for ${f}_{{\rm{g}}}\lt 0.3$. To determine the reason for this, in Figure 7 we plotted the main heating and cooling mechanisms for $\alpha ={10}^{-4}$ and ${f}_{{\rm{g}}}=0.1$. For low f_g, turbulent scattering is insignificant because its heating rate is proportional to ${{\rm{\Sigma }}}_{{\rm{g}}}^{2}$. In deriving Equation (49), we assumed that the upper boundary of the GI region is determined by a balance between turbulent scattering and collisions. However, this assumption breaks down, and thus, in this parameter region, Equation (49) differs from the numerical results.

Zoom In Zoom Out Reset image size

As expected from Equation (49), the dependence on a is weak. From the numerical results shown in Figure 6(b), the existence of the GI region is independent of a. The value of α required to cross the GI region slightly decreases with a, but its dependence is very weak.

As shown in Figure 6(c), the dependence on the dust-to-gas ratio γ is strong. The GI easily takes place for larger values of γ. The critical α for the existence of the GI region, ${\alpha }_{\mathrm{cr},1}$ (Equation (49)), is proportional to ${\gamma }^{3}$, which agrees well with the numerical results. However, the dependence of the critical α on γ for crossing the GI region is different; it is approximately proportional to ${\gamma }^{2}$. This dependence will be discussed below. The difference from ${\alpha }_{\mathrm{cr},1}$ increases with γ for $\gamma \gtrsim 0.05$.

The critical α decreases with increasing T₁, as shown in Figure 6(d). Increased temperatures lead to suppression of the GI. From Equation (4), η is proportional to T₁. The typical difference between the velocity of dust and that of gas is determined from $\eta {v}_{{\rm{K}}}$ because ${v}_{\mathrm{ran}}\lt \eta {v}_{{\rm{K}}}$. Thus, a higher T₁ means that the gas drag is strong, and this causes strong turbulent stirring and inhibits the GI. When ${T}_{1}=280$, the critical α that is often adopted is 0.4 times that when ${T}_{1}=120$.

Figure 6(e) shows that the critical α decreases with increasing ${\tau }_{e}$ as ${\tau }_{e}^{-1/2}$. A small ${\tau }_{e}$ means that the duration of the fluctuations in the turbulent velocity is short. In this case, from Equation (35), the heating rate due to turbulent stirring is small. Thus, the GI tends to occur for smaller ${\tau }_{e}$.

The critical α decreases with increasing C_turb as ${C}_{\mathrm{turb}}^{-1/2}$, as shown in Figure 6(f). From Equation (37), the heating rate due to turbulent scattering decreases with increasing C_turb. Thus, for smaller C_turb, the critical α is larger.

Figure 8 shows the influence of the power-law indices ${\beta }_{{\rm{g}}}$ and ${\beta }_{{\rm{t}}}$. In the fiducial model, the dependence on a is weak. However, if we adopt the general power-law indices, the dependence on a can be strong. As discussed in Section 4.3, the dependence on a vanishes only if ${\beta }_{{\rm{g}}}-{\beta }_{{\rm{t}}}=1$. As ${\beta }_{{\rm{g}}}$ increases, the dependence on a becomes weak, while as ${\beta }_{{\rm{t}}}$ increases, the dependence on a becomes strong. This behavior is consistent with Equation (49).

Zoom In Zoom Out Reset image size

In this section, we clarify the effect of the nonequilibrium random velocity on the onset of the GI. In the previous section, we assumed the equilibrium ${\sigma }_{e}$ and ${\sigma }_{i}$. However, under some conditions, the relaxation time of ${\sigma }_{e}$ and ${\sigma }_{i}$ may be comparable to the growth time of dust aggregates. In this case, the equilibrium values change before ${\sigma }_{e}$ and ${\sigma }_{i}$ reach equilibrium. Thus, the actual ${\sigma }_{e}$ and ${\sigma }_{i}$ may deviate from the equilibrium values.

We will simultaneously consider the evolution of the mass and random velocity. The mass evolution equation is

$$\begin{eqnarray}&&\displaystyle \frac{{{dm}}_{{\rm{d}}}}{{dt}}={C}_{\mathrm{grow}}{P}_{\mathrm{col}}{h}^{2}{a}^{2}{{\rm{\Sigma }}}_{{\rm{d}}}{\rm{\Omega }},\end{eqnarray} \tag{ 50 }$$

where C_grow is the sticking probability. The perfect accretion corresponds to ${C}_{\mathrm{grow}}=1$. For ${C}_{\mathrm{grow}}\ne 1$, we may need to change C_col. Strictly speaking, C_col depends on the energy dissipation rate on collisions, such as the restitution coefficient and the friction coefficient on the dust aggregate surface. For simplicity we assume ${C}_{\mathrm{col}}=1/2$ for any C_grow. We solve the differential equations given as Equations (40), (41), and (50) for ${\sigma }_{e}$, ${\sigma }_{i}$, and m_d, respectively. For simplicity, we assume that dust aggregates always have the equilibrium internal density, ${\rho }_{\mathrm{int}}={\rho }_{\mathrm{eq}}$.

Figure 9 shows the evolution of Q for the fiducial model with an initial mass of $m={10}^{10}\,{\rm{g}}$ and initial equilibrium values for ${\sigma }_{e}$ and ${\sigma }_{i}$. We assume the perfect accretion (${C}_{\mathrm{grow}}=1$). As shown below, in the case of the perfect accretion, the difference between the equilibrium and nonequilibrium models is the most significant. Toomre's Q deviates slightly from the equilibrium value since the dust aggregates grow during the relaxation of the random velocity. Among the models shown in Figure 6, the deviation of the minimum Q is 9.7% ± 13.2%. Thus, the nonequilibrium effect is not significant for the onset of the GI. We also evaluated the influence of the initial values by considering values for ${\sigma }_{e}$ and ${\sigma }_{i}$ that were 5 times and 1/5 times their equilibrium values. As shown in Figure 9, the difference in Q that stems from the initial values quickly vanishes as the system evolves.

Zoom In Zoom Out Reset image size

We compare the various relevant timescales. The growth timescale is defined as

$$\begin{eqnarray}&&{t}_{\mathrm{grow}}={m}_{{\rm{d}}}/\displaystyle \frac{{{dm}}_{{\rm{d}}}}{{dt}},\end{eqnarray} \tag{ 51 }$$

and the timescale of gravitational scattering is defined as

$$\begin{eqnarray}&&{t}_{\mathrm{grav}}={\sigma }_{e}^{2}/\left|{\left(\displaystyle \frac{d{\sigma }_{e}^{2}}{{dt}}\right)}_{\mathrm{grav}}\right|.\end{eqnarray} \tag{ 52 }$$

The other dynamical timescales t_col, ${t}_{\mathrm{gas},\mathrm{drag}}$, ${t}_{\mathrm{turb},\mathrm{stir}}$, and ${t}_{\mathrm{turb},\mathrm{grav}}$ are defined in the same way. The left panel of Figure 10 shows these timescales. If the growth is sufficiently slow compared with the dynamical timescales, the eccentricity reaches the equilibrium value. Figure 10 shows that the growth timescale is comparable to the dynamical timescales. Therefore the equilibrium eccentricity, which depends on the mass, changes before the eccentricity converges to the equilibrium value. However, since the growth timescale is not very different from the dynamical timescales, the gap from the equilibrium value is not so large.

Zoom In Zoom Out Reset image size

Figure 10 shows two sharp peaks of t_grav. This is because P_VS is negative if ${\sigma }_{i}/{\sigma }_{e}$ is small (Ohtsuki et al. 2002). In the relatively high velocity cases, gravitational scattering tends to realize ${\sigma }_{i}/{\sigma }_{e}\simeq 0.5$ (Ida 1990). Therefore, low ${\sigma }_{i}/{\sigma }_{e}$ leads to negative P_VS, while high ${\sigma }_{i}/{\sigma }_{e}$ leads to negative Q_VS. The right panel of Figure 10 shows that if ${\sigma }_{i}/{\sigma }_{e}$ is about 0.3 around ${m}_{{\rm{d}}}\sim {10}^{16}\,{\rm{g}}$, then P_VS is negative. Around the points where ${P}_{\mathrm{VS}}\simeq 0$, t_grav becomes very large.

The effect of the imperfect accretion is shown in Figure 11. The initial ${\sigma }_{e}$ and ${\sigma }_{i}$ are set as the equilibrium values. The evolution timescale is the fastest for the perfect accretion model. As C_grow decreases, the evolution timescale becomes longer, which is inversely proportional to C_grow. The difference between the equilibrium and nonequilibrium models becomes less with decreasing C_grow. This is because for smaller C_grow, the growth time is longer compared to the dynamical timescales, and thus the eccentricity can converge to the equilibrium value before the mass changes.

Zoom In Zoom Out Reset image size

In Figures 6 and 8, we examine the effect of the nonequilibrium velocity on the GI. If the nonequilibrium effect is considered, the value of α necessary for crossing the GI region decreases slightly; however, the difference is small. The nonequilibrium effect of the random velocity is thus insignificant. In summary, the above results justify the equilibrium random velocity model.

The dust evolution track is characterized by the monomer parameters E_roll, r₀, and ${\rho }_{0}$. We first examine the effect of E_roll. In Figure 12, it can be seen that E_roll has little effect on the critical α for crossing the GI region. For small E_roll, the critical α only slightly decreases. This is due to the structure of the GI region. As shown in Figure 13, the GI region stretches as α decreases. In particular, the upper bound on the GI region increases rapidly with α. For $\alpha =2\times {10}^{-3}$, we can see the elongated GI region for ${\rho }_{\mathrm{int}}\gt 0.1\,{\rm{g}}\,{\mathrm{cm}}^{-3}$. If this region appears, the dust evolution inevitably crosses the GI region.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To examine the condition for the existence of an elongated GI region, in Figure 14 we show the main heating and cooling mechanisms for $\alpha =2\times {10}^{-3}$. The lower boundary of the elongated region is determined by the balance between turbulent stirring and gas drag. We calculate ${(d{\sigma }_{e}^{2}/{dt})}_{\mathrm{turb},\mathrm{stir}}$ in the same way as in Section 4.3 except for C_D. Figure 5(e) shows that Stokes drag ${C}_{{\rm{D}}}\simeq 24/\mathrm{Re}$ is a good approximation. Thus, we adopt Stokes drag. Using these approximations, we calculate Q. From $Q\lt 2$, we obtain

$$\begin{eqnarray}&&{m}_{{\rm{d}}}\gt {m}_{\mathrm{low},2}=2.78\times {10}^{-2}\displaystyle \frac{{\alpha }^{3/2}{c}_{{\rm{s}}}^{3/2}{\nu }^{3/2}{\tau }_{e}^{3/2}{M}_{* }^{3/2}{{\rm{\Sigma }}}_{{\rm{g}}}^{3/2}}{{Q}_{\mathrm{cr}}^{3}{a}^{9/2}{\rho }_{\mathrm{int}}^{1/2}{G}^{3/2}{{\rm{\Sigma }}}_{{\rm{d}}}^{3}}.\end{eqnarray} \tag{ 53 }$$

Similarly, the upper boundary of this region is determined by the balance between turbulent scattering and gas drag. Thus, we obtain the upper boundary as

$$\begin{eqnarray}&&{m}_{{\rm{d}}}\lt {m}_{\mathrm{high},2}=5.62\times {10}^{3}\displaystyle \frac{{Q}_{\mathrm{cr}}^{3}{\nu }^{3/2}{{\rm{\Sigma }}}_{{\rm{d}}}^{3}}{{\alpha }^{3/2}{C}_{\mathrm{turb}}^{3/2}{c}_{{\rm{s}}}^{3/2}{\rho }_{\mathrm{int}}^{1/2}{{\rm{\Sigma }}}_{{\rm{g}}}^{3/2}}.\end{eqnarray} \tag{ 54 }$$

The condition for the existence of an elongated region is ${m}_{\mathrm{low},2}\lt {m}_{\mathrm{high},2}$, from which we obtain the critical α for crossing the GI region:

$$\begin{eqnarray}&&\alpha \lt {\alpha }_{\mathrm{cr},2}=58.7\displaystyle \frac{{Q}_{\mathrm{cr}}^{2}{a}^{2}{v}_{{\rm{K}}}{{\rm{\Sigma }}}_{{\rm{d}}}^{2}}{{C}_{\mathrm{turb}}^{1/2}{c}_{{\rm{s}}}{\tau }_{e}^{1/2}{M}_{* }{{\rm{\Sigma }}}_{{\rm{g}}}},\end{eqnarray} \tag{ 55 }$$

which can be rewritten as

$$\begin{eqnarray}{\alpha }_{\mathrm{cr},2} & = & 3.77\times {10}^{-3}{f}_{{\rm{g}}}{\tau }_{e}^{-1/2}{\left(\displaystyle \frac{\gamma }{0.018}\right)}^{2}{\left(\displaystyle \frac{a}{1\mathrm{au}}\right)}^{{\beta }_{{\rm{t}}}/2-{\beta }_{{\rm{g}}}+3/2}\\ & & \times {\left(\displaystyle \frac{{C}_{\mathrm{turb}}}{3.1\times {10}^{-2}}\right)}^{-1/2}{\left(\displaystyle \frac{{M}_{* }}{{M}_{\odot }}\right)}^{-1/2}{\left(\displaystyle \frac{{T}_{1}}{120}\right)}^{-1/2}.\end{eqnarray} \tag{ 56 }$$

We show this condition in Figure 6. If $\alpha \lt {\alpha }_{\mathrm{cr},2}$, the dust evolution crosses the GI region. Thus, we propose the following condition for the onset of the GI:

$$\begin{eqnarray}&&\alpha \lt \min ({\alpha }_{\mathrm{cr},1},{\alpha }_{\mathrm{cr},2}).\end{eqnarray} \tag{ 57 }$$

Zoom In Zoom Out Reset image size

In deriving this condition, we did not consider the nonequilibrium effect of the random velocity. However, as shown in Figure 6, that effect is insignificant. For the outer disk ($a=10\mbox{--}20\,\mathrm{au}$), we need a safety factor, such as $\alpha \lt 0.5\min ({\alpha }_{\mathrm{cr},1},{\alpha }_{\mathrm{cr},2})$.

In Paper I, we proposed a similar condition, and its parameter dependence is exactly the same as that of ${\alpha }_{\mathrm{cr},1}$. Other than when γ is particularly large or particularly small, the difference between ${\alpha }_{\mathrm{cr},1}$ and ${\alpha }_{\mathrm{cr},2}$ is very small. Thus, the simple condition proposed in Paper I ($0.5\times {\alpha }_{\mathrm{cr},1}$) is also a good estimate of the condition for crossing the GI region in most parameter regimes.