PLANETESIMAL FORMATION AT THE BOUNDARY BETWEEN STEADY SUPER/SUB-KEPLERIAN FLOW CREATED BY INHOMOGENEOUS GROWTH OF MAGNETOROTATIONAL INSTABILITY

We consider a small region around the midplane which is rotating with the Keplerian frequency Ω at a distance r from a central star to study local dust motion and magnetohydrodynamics. The coordinates that we use are (x, y, z) where x is the radial distance from r, y is tangential distance, and z is vertical distance from the disk midplane.

We include centimeter to meter-size dust particles as super-particles. Total number of the super-particles is O(10⁷) and each super-particle represents O(10⁶)–O(10⁷) small dust particles. The equation of motion of the ith particle is given by

$$\begin{eqnarray} \frac{d \mathbf {v}_{i}}{d t}&=& {-}2\mathbf {\Omega }\times \mathbf {v}_{i}+3\Omega ^2x_{i}\hat{\mathbf {x}} -\frac{1}{\tau _{f}}\left(\mathbf {v}_{i}-\mathbf {u}\right) -\nabla \Phi, \end{eqnarray} \tag{ 1 }$$

where $\mathbf {u}$ is the gas velocity at the location of the ith particle, which is interpolated using gas velocities at the neighboring grid points. The third term on the right-hand side (r.h.s.) represents the specific gas drag force on the dust particle. We adopt the Epstein drag force,

$$\begin{eqnarray} &&\tau _f=\rho _{p} a/(\rho _g c_s), \end{eqnarray} \tag{ 2 }$$

where ρ_p and a are the internal density and radius of dust particles, ρ_g and c_s are the spatial density and sound velocity of surrounding gas. The Epstein law is valid for centimeter to meter-size dust particles at ∼5 AU if gas column density is less than twice that in the minimum solar nebula model (MMSN; Hayashi 1981). Even if we consider higher column density, the deviation in the drag force strength from the Epstein law would not be significant. The last term in Equation (1) is the self-gravitational force of the dust particles. We calculate the gravitational potential Φ from the interpolated dust spatial density ρ_d, solving the Poisson equation,

$$\begin{eqnarray} &&{\nabla }^2 \Phi =4\pi G{\rho }_d. \end{eqnarray} \tag{ 3 }$$

We calculate the self-gravity of the dust particles only in the situations where the GI is expected. We neglect the vertical gravity of the host star that causes settling of dust particles onto a thin layer (∼0.01H where H is the disk scale height; Schräpler & Henning 2004) near the midplane, because we are interested in the radial concentration of dust but not vertical settling. While this would not affect the growth of MRI around the disk midplane that we simulate (|z| < 0.25H), we do not have to resolve the thin dust layer, avoiding expensive computational cost.

For the disk gas, we use the isothermal resistive MHD equations,

$$\begin{eqnarray} \frac{\partial \mathbf {u}}{\partial t} +(\mathbf {u}\cdot \nabla )\mathbf {u}&=& {-}\frac{1}{\rho _g}\nabla \left(P+\frac{\mathbf {B}^2}{8\pi }\right)\nonumber\\ &&+\frac{1}{4\pi \rho _g}\left(\mathbf {B}\cdot \nabla \right)\mathbf {B} -2\mathbf {\Omega }\times \mathbf {u} +3\Omega ^2x\hat{\mathbf {x}} \nonumber \\ && -\beta c_{s}\Omega \hat{\mathbf {x}}-\frac{\epsilon }{\tau _f} \left(\mathbf {u}-\langle \mathbf {v} \rangle \right), \end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray} \frac{\partial \rho _g}{\partial t} +\nabla \cdot \left(\rho _g \mathbf {u}\right)&=&0, \end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray} \frac{\partial \mathbf {B}}{\partial t} &=& \nabla \times \left[\left(\mathbf {u}\times \mathbf {B}\right)- \eta \left(\nabla \times \mathbf {B}\right) \right], \end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray} P&=&c^2_{s}\rho _g, \end{eqnarray} \tag{ 7 }$$

where we assume constant c_s (an isothermal disk), $\langle \mathbf {v} \rangle$ is the mean velocity field of the dust particles in the grid cell, = ρ_d/ρ_g is the dust to gas ratio, and the term $-\beta c_{s}\Omega \hat{\mathbf {x}}$ expresses the global pressure gradient, which is separated from the local one. Let P₀ and δP be the global pressure and the local pressure variation due to the effect of MRI (P = P₀ + δP). Assuming that P₀∝r^q, the global pressure gradient is given by

$$\begin{eqnarray} &&{-}\frac{1}{\rho _g}\frac{\partial P_0}{\partial r} = -\frac{1}{\rho _g} \frac{P_0}{r} q = -\frac{H}{r} q c_{s}\Omega = -\beta c_{s} \Omega, \end{eqnarray} \tag{ 8 }$$

where H is the disk scale height defined by H = c_s/Ω. In our local model, β = qH/r is approximated to be constant. Note that both q and β are negative. Due to the radial pressure gradient, the gas rotation angular velocity deviates from Keplerian one as

$$\begin{eqnarray} \Omega _g &\simeq& \Omega \left(1 + \frac{1}{2}\left[\frac{H}{r}\right]^2 \frac{d \ln P}{d \ln r}\right)\nonumber\\ &=& \Omega \left(1 + \frac{1}{2}\frac{H}{r} \beta + \frac{1}{2}\left[\frac{H}{r}\right]^2 \frac{d \ln \delta P}{d \ln r}\right). \end{eqnarray} \tag{ 9 }$$

The last term on the r.h.s. of Equation (4) is the dust drag force (back-reaction of gas drag force on the dust particles). Except for this term, the equations of motion for disk gas and dust particles are the same as those in Paper II. The purpose of this paper is to study the effect of this term on dust concentration. When dust particles are accumulated and takes a large value of ≳ O(1), this term would influence the gas flow. The induction equation (Equation (6)) has the diffusion term by ohmic dissipation. Since we are interested in relatively inner regions (≲a few dozen AU), we neglect the ambipolar diffusion that could influence MRI growth in the outer regions of ≳ 100 AU (Chiang & Murray-Clay 2007). We treat the magnetic resistivity as a constant parameter for simplicity, though in reality it depends on the density of tiny grains (≲ μm).

We scale length, time, and velocity by H, 1/Ω, and c_s in simulations. We solve the MHD equations by combining the CIP scheme (Yabe & Aoki 1991) and MOC-CT method (Stone & Norman 1992a, 1992b). The dust density is allocated to the closest eight grids in the three-dimensional space using the cloud-in-cell model. This algorithm strictly conserves angular momentum transfer from dust to gas by using a similar method as Johansen et al. (2007). The boundary conditions in all directions are periodic. For the radial boundary, however, we take into account Keplerian differential rotation with the shearing-box model (Wisdom & Tremaine 1988; Hawley et al. 1995). While we are not interested in vertical sedimentation, we want to keep total dust mass in the whole simulation area, so we adopt the periodic boundary condition also for the vertical direction. The Poisson equation (3) is calculated by fast Fourier transform. This method requires the periodic boundary. For the radial direction, according to the sheared boundary, we shift the phase azimuthally in Fourier space after the Fourier transform in the periodic azimuthal direction, following Johansen et al. (2007). Our simulations are performed by a vector computer, NEC SX-9 at ISAS/JAXA.

The setup for the simulation in this paper is the same as CASE2 in Paper II. We assume non-uniform B_z to set a marginal MRI state where localized dead (stable) and active (unstable) regions coexist in the initial conditions. The initial magnetic field is $\mathbf {B}_0=(0, B_{0}\sin \theta, B_{0}\cos \theta)$, where θ = θ(x) is the angle between the magnetic field and the vertical axis. We assume a constant value of B₀ that is determined by plasma beta =400 to establish the initial equilibration. With the constant magnetic resistivity of η = 0.002H²Ω we adopt, a threshold vertical magnetic field for MRI is B_{z, crit} ∼ 0.2B₀ (Jin 1996). We set radially varying θ in which θ = 0° (B_z ≫ B_{z, crit}) in the central zone and θ = 85° (B_z < B_{z, crit}) in the side regions (see Figure 2(b) in Paper II), such that MRI grows only in the central zone. We set the radial width of the initially active region as L_u = 1.4H in all cases, and that of the dead regions as L_s = 4.0H in model-s40-* and L_s = 0.5H in model-s05-*, where the asterisk * represents other simulation parameters (see below).

We assume equal-size dust particles with τ_f = 1.0/Ω in model-s*-t10-* except for τ_f = 0.1/Ω in model-s40-t01-e010, neglecting their coalescence and fragmentation. At 5 AU in MMSN, τ_fΩ = 0.1 and 1.0 correspond to the dust sizes of approximately 3 and 30 cm, respectively. The super-particles are distributed uniformly such that the initial dust-to-gas ratio ₀ = 0.1 for model-s*-t*-e010 and ₀ = 0.01 for model-s40-t10-e001. Note that the dust-to-gas ratio in our simulation box corresponds to that near the midplane layer. If vertically global sedimentation of dust particles is taken into account, the dust-to-gas ratio in our simulation box should be larger than that averaged over the whole disk. For comparison, we also present the results of Paper II without the dust drag as model-*-*-test. The self-gravity of dust particles in Equation (1) is switched on only in the saturated state in which the GI is expected.

All of our simulations start with uniform gas density and pressure. The global pressure gradient is set to be β = −0.04. This assumed value of |β| is a few times smaller than that expected at ∼5 AU in MMSN. In order to compare the results with Johansen et al. (2007) and Paper II, we adopt the small value. As was argued in Paper II, the results would not be affected significantly by the value of |β|. The initial angular velocities of gas and particles are u_y/c_s = −(3/2)(x/H) + β/2 and v_y = −(3/2)(x/H), respectively. Because gas rotates slower, the particles migrate inward (negative direction of x) in the initial state. Initial disturbances are given randomly to the gas radial velocity with the amplitude of 0.001c_s. The size of our simulation box is (L_x, L_y, L_z) = ((2.5–9.5)H, 1.0H, 0.5H) and the resolution is dx = dy = dz = 0.01H. Eight particles are distributed in each grid initially and the total number is ∼O(10⁷). We have tested different numbers of distributed particles with different mass such that total mass is conserved and found that the results converge if the distributed number of particles for each grid is ≳ 8. The initial setup is summarized in Table 1.

Table 1. Setup of Individual Runs

Run	Ls	τfΩ	0	Self-gravity
model-s40-t10-e010	4.0H	1.0	0.10	Off and on
-t10-e001	4.0H	1.0	0.01	Off and on
-t10-test	4.0H	1.0	0.0	Off
-t01-e010	4.0H	0.1	0.10	Off
-t01-test	4.0H	0.1	0.0	Off
model-s05-t10-e010	0.55H	1.0	0.10	Off
-t10-test	0.55H	1.0	0.0	Off

Notes. L_s is the radial width of the initially dead region, τ_f is the friction time, and ₀ is the initial dust-to-gas density ratio; model-s40-t10-e010 and e001 are also restarted with the introduction of self-gravity of particles.

Download table as:  ASCIITypeset image

In this subsection, we discuss the results with τ_fΩ = 1.0, model-s40-t10-e010 with ₀ = 0.10, and model-s40-t10-e001 with ₀ = 0.01. Strong dust concentration was found in corresponding models without the dust drag in Paper II¡¡ (model-s40-t10-test), which is summarized in Figure 1. MRI is excited only in the initially active central region (−0.71 < x/H < 0.71). The MRI turbulence transfers mass and angular momentum of disk gas to establish a local rigid rotation in the initially active zone through the turbulent viscosity. Then, in the outer half of the active zone, gas flow is accelerated and migrates outward because of angular momentum gain. On the other hand, in the inner half, it is decelerated and migrates inward. As a result, gas mass is moved from the central zone to the side zones and gas pressure is lowered in the central zone. The effect of the pressure modulation extends by radial scale of ∼H. Equation (9) shows that

$$\begin{equation} \delta \tilde{u}_y = \frac{u_y - v_{\rm kep}}{c_s} \simeq \frac{\beta }{2}+ \frac{1}{2}\frac{H}{r} \frac{d \ln \delta P}{d \ln r} = - 0.02 + \frac{1}{2}\frac{d \ln \delta P}{d \ln (x/H)}. \end{equation} \tag{ 10 }$$

This equation shows that super-Keplerian regions are associated with locally positive pressure gradient with some offset due to global pressure gradient. We find that super-Keplerian regions are created in 0 ≲ x/H ≲ 2.0 and −2.8 ≲ x/H ≲ −2.5 (panels (a) and (b)). Although MRI turbulence has decayed in the result at tΩ = 70, the magnetic field has not been dissipated in the central zone (panel (c)). The MRI is suppressed by the disappearance of shear motion. If the rigid rotation is perturbed toward the original Keplerian motion, the retrieved shear motion causes MRI again to recover the rigid rotation. Thus, this profile is stable and strong.

Figure 2 presents the dust density in the saturated state in model-s40-t10-e010 with ₀ = 0.10 (panel (a)) and model-s40-t10-e001 with ₀ = 0.01 (panel (b)), in comparison with model-s40-t10-test without the dust drag (panel (c)). Because the dust drag depends on spatial density of the dust particles, the results depend on ₀. In all cases, after the particles are swept out of the active region by the MRI turbulence, they accumulate at the outer edge of the super-Keplerian zone at x/H ≃ 2.0 and −2.5. Particles leaving the simulation box from the small x (left hand) boundary reenter the simulation region from the large x (right hand) boundary after the shearing-box correction is taken into account. The dust that reentered from the right-hand boundary is halted at x/H ≃ 2.0, resulting in further increasing of the dust density. The number of locations of dust concentration is fewer in the case with the dust drag, because the drag smoothes out small amplitude fluctuations of gas pressure. In the case with the drag, the individual dust concentrated areas are broader in model-s40-t10-e010 than in model-s40-t10-e001, which is discussed below. We also found that velocity dispersion is lower for a denser clump, which was not observed in the case without the dust drag. The effect of the reduced velocity dispersion will be discussed in Section 4.

Zoom In Zoom Out Reset image size

Figure 3 shows the time evolution of the maximum density of dust particles (ρ_d) in the simulation cells. The dust density is scaled by the gas density averaged in the whole simulation region (〈ρ_g〉). The results are compared with those without the dust drag (dashed lines; model-Ls40-t10-test) for ₀ = 0.10 and ₀ = 0.01. In the case without the dust drag, only concentration relative to the initial state is measured, so these lines are drawn by the evolution of concentration assuming ₀ = 0.10 or ₀ = 0.01. The maximum density continues to increase monotonically in this case. However, the growth of the maximum value is saturated in both cases with the dust drag. The saturation is faster and the rate of increase is slightly smaller in model-s40-t10-e010 than in model-s40-t10-e001.

Zoom In Zoom Out Reset image size

In order to explain the broadening of the dust-accumulated region, in addition to the three-dimensional simulations, we performed the two-dimensional (x–z) version of the model-s40-t10-e010, in which the effect of the dust drag on the radial migration of the dust particles is more clearly shown. In Figure 4(a), we plot the difference of gas angular velocity from Kepler angular velocity, $\delta \tilde{u}_y=(u_y-v_{\rm kep})/c_s$, near the sub- and super-Keplerian boundary in the two-dimensional simulation. The velocities are averaged over the vertical direction. At tΩ = 55.0 (dashed line), the dust particles are expected to assemble at x/H ≃ 1.75, where $\delta \tilde{u}_y=0$. At tΩ = 87.4 (solid line), however, the region with $\delta \tilde{u}_y\sim 0$ becomes broader (1.7 ≲ x/H ≲ 1.9), because more dust particles have migrated to this region and their drag makes the gas flow close to Keplerian. Consequently, the dust particles are more broadly distributed there (Figure 4(b)). Figure 4(c) schematically illustrates the effect of the dust drag. The radial velocity (${\propto} \delta \tilde{u}_y$) of a migrating dust particle becomes slower as the dust particle approaches the dust concentrated region of $\delta \tilde{u}_y=0$, like a "traffic jam." Due to the modulation by the dust drag, the radial width of the Keplerian region is expanded, and the dust particles stop their inward migration before they reach the location at which $\delta \tilde{u}_y=0$ originally. Thus, the maximum dust density is self-regulated as shown in Figure 3.

Zoom In Zoom Out Reset image size

In the three-dimensional simulation, similar results are obtained, although small amplitude oscillations remain. The similarity implies that geometry is not the main cause for the broadening of the dust concentration region.

In Paper II, for smaller particles with τ_fΩ = 0.1 (model-s40-t01-test), we found that the dust concentration is not significant because the smaller dust particles are affected more by the turbulent diffusion. In the results in Section 3.1, we found that the dust drag suppresses the local turbulence and velocity dispersion of the dust particles. It could enhance the dust accumulation by decaying the turbulence around the dust-accumulated area, as found in Johansen et al. (2007).

However, the time evolution of the maximum scaled density of dust particles (ρ_d/〈ρ_g〉) in model-s40-t01-e010 is not significantly different from that in model-s40-t01-test (Figure 5). In both cases, after ρ_d rapidly increases by the diffusing out from the initially active region by the MRI turbulence tΩ ∼ 25, it gradually increases and is saturated at tΩ ≳ 50–70. The velocity dispersion of the dust particles is actually reduced from that in the case without the dust drag, but the effect is not strong enough to enhance the dust concentration.

Zoom In Zoom Out Reset image size

In Paper II, we found that with the smaller initial dead region L_s = 0.55H, the viscosity in the saturated state is α ∼ 10⁻², which is much larger than α ∼ 10⁻⁴ for the runs with L_s = 4.0H in model-Ls40-*. We found in Paper I that MRI turbulence does not decay if the magnetic Elsasser number radially averaged over the simulation region is smaller than unity in the initial state. The Elsasser number is defined by

$$\begin{eqnarray} &&\Lambda _{\rm m,ave}=v^2_{{\rm A}z}/\eta \Omega, \end{eqnarray} \tag{ 11 }$$

where $v_{{\rm A}z} = B_z /\sqrt{4 \pi \rho _g}$ is the z component of Alfvén velocity and ρ_g is the spatial density of the disk gas. The run with L_s = 0.55H corresponds to Λ_{m, ave} ∼ 0.5.

The stronger remnant turbulence limited ρ_d/〈ρ_g〉 to values less than 100 even for τ_fΩ = 1.0 in the case without the dust drag. We performed the run with L_s = 0.55H, adding the dust drag (model-s05-t10-e010). Although the drag force reduces the turbulent diffusion in the local concentrated region, the maximum dust density is not enhanced from that in the case without the dust drag (Figure 6).

Zoom In Zoom Out Reset image size

The GI is expected to arise when the self-gravity of dust particles is dominant over their thermal fluctuation (velocity dispersion), in other words, when the radius (size) R of a dust clump with mass M is smaller than the Jeans (Bondi) radius,

$$\begin{eqnarray} &&R_{\rm J}=\frac{2GM}{\sigma ^2} =2\left(\frac{\sigma }{c_s}\right)^{-2} \frac{M}{M_{\ast }} \left(\frac{H}{R}\right)^{-3} H, \end{eqnarray} \tag{ 12 }$$

where M = M(R) and σ = σ(R) are the total mass and mean velocity dispersion of particles in the region at distance R from the center of the clump, M_* is the host star's mass, c_s = ΩH, and $\Omega = \sqrt{GM_{\ast }/r^3}$. The background shear is included in σ. In order to numerically resolve GI, we set the grid size in the x-direction such that dx < 0.5R_J.

The condition of GI is often described by linear analysis of a uniform axisymmetric disk (e.g., Safronov 1969; Goldreich & Ward 1973; Sekiya 1983), which is essentially equivalent to Toomre's condition (Toomre 1964),

$$\begin{equation} 1 < Q = \frac{\Omega \sigma }{\pi G \Sigma _d}, \end{equation} \tag{ 13 }$$

where Σ_d is the unperturbed solid column density. If we use M ∼ πΣ_dR² and σ ∼ RΩ, Toomre's condition is identical to R < R_J except for a numerical factor of two. Because significant radial inhomogeneity develops before GI occurs in our case, we use the condition R < R_J that can be locally applied, rather than Q < 1.

The mass of each super-particle m is given by ρ_d0/n₀, where ρ_d0 and n₀ are the spatial density and number density of particles in the initial conditions. Then, the dust clump mass is given by

$$\begin{eqnarray} &&M\!=\!m N_R\!=\!\frac{N_R}{n_0}\rho _{d0}\!=\!\sqrt{\frac{8\pi ^3}{9}} \frac{N_R}{N_{R0}}\! \left(\!\frac{H}{r}\!\right)^{-2} \epsilon _0 \frac{\Sigma _{g0} r^2}{M_{\ast }} \left(\frac{R}{H}\right)^{3} M_{\ast }.\nonumber\\ \end{eqnarray} \tag{ 14 }$$

where N_R is the total particle number in R, N_R0 = (4π/3)n₀R³ is its initial value, ρ_d0 = ₀ρ_g0, and $\rho _{g0} = \Sigma _{g0}/\sqrt{2\pi } H$ by the assumption of a vertically isothermal disk. From Equtaions (12) and (14), the scaled Jeans radius is given by

$$\begin{eqnarray} \frac{R_{\rm J}}{H} &=& \sqrt{\frac{32\pi ^3}{9}} \left(\frac{\sigma }{c_s}\right)^{-2} \frac{N_R}{N_{R0}} \left(\frac{H}{r}\right)^{-5} \epsilon _0 \frac{\Sigma _{g0} r^2}{M_\ast } \left(\frac{R}{H} \right)^{3},\nonumber\\ \end{eqnarray} \tag{ 15 }$$

The initial dust-to-gas density ratio ₀ is given. The simulation results give dust enhancement in a clump N_R/N_R0, the clump size R/H, and velocity dispersion in the clump σ/c_s (the equations of motions we use are normalized by H and c_s). From the results and simulation parameters, we can evaluate R_J for any given values of Σ_g0r²/M_* and H/r that are specified by disk model through Equation (15).

The onset of GI depends on N_R and σ as shown by Equation (15). The dust drag suppresses N_R and reduces σ (Section 3). Here, using the arguments in Section 4.1, we examine the possibility of the GI in model-s40-t10-e010, which is one of the most promising runs for the GI.

Figure 7 shows N_R/N_R0 from the densest grid point (panel (a)) and the corresponding σ/c_s (panel (b)) in the results of model-s40-t10-e010 and model-s40-t10-test at tΩ = 58.0. Panel (c) shows Jeans radius R_J calculated for M_* = M_☉, r = 5 AU, H/r = 0.055, and Σ_g0 = 150 g cm⁻³(∼ Σ_MMSN at r = 5 AU, where Σ_MMSN is gas column density in MMSN). While the particle concentration is lowered by the dust drag only slightly (panel (a)), the velocity dispersion is significantly reduced (panel (b)). Since the positive effect for the GI (reduction in the velocity dispersion) dominates over the negative one (suppression in particle concentration), R_J calculated by Equation (15) is higher and the condition for the GI (R_J > R) is satisfied in the case with the dust drag (panel (c)).

Zoom In Zoom Out Reset image size

We carried out additional simulations including the self-gravity force of dust particles to demonstrate the formation of gravitationally bound clumps that may lead to planetesimals, in model-s40-t10-e010. We set M_* = M_☉, r = 5 AU, H/r = 0.055, and Σ_g0 = 280 g cm⁻³ ∼ 2Σ_MMSN(r = 5 AU). To reduce simulation cost, we introduced the self-gravity at tΩ = 96 when the dust concentration becomes saturated and R_J is much larger than the grid size dx.

Figure 8 shows the time evolution of the gravitational collapse. Shortly after introduction of the self-gravity, the elongated high density region is kinked (tΩ = 99) and it is separated into several clumps (tΩ = 100). The clumps grow by accreting surrounding dust particles and other clumps (tΩ = 120–140).

Zoom In Zoom Out Reset image size

To confirm that the clumps are gravitationally bound and estimate the mass of formed planetesimals, we define the range of a clump by its Hill radius. The time evolution of the clump is shown in Figure 9. Figure 9(a) shows the Hill radius, where the grid size is represented by a dashed line. Immediately after the introduction of the self-gravity at tΩ = 96, the Hill radius exceeds the grid size. After that, the clump is numerically resolved. Figure 9(b) shows that the velocity dispersion of the particles in the clump is always smaller than the surface escape velocity v_esc of the clump, which implies that the clump is gravitationally bound. If the gravitational collapse continues, it may form a planetesimal, although this simulation does not have resolution to follow the subsequent collapse.

Zoom In Zoom Out Reset image size

Figure 9(c) shows the temporal development of the clump mass M, which may correspond to the planetesimal mass. The abrupt jumps at tΩ ∼ 102 and ∼104 are caused by collisional merging with other clumps. Since the destruction process is not properly included in our simulation, such rapid growth may be unrealistic. A conservative estimate for the planetesimal mass may be the mass before the abrupt jumps, that is, ∼4 times Ceres' mass. However, note that this mass is close to the resolution of our simulation (Figure 9(a)) and the clump mass may be smaller in a higher-resolution simulation (Johansen et al. 2011).

We also performed the simulation with the self-gravity in model-s40-t10-e001, in which the Jeans radius is slightly larger than our grid size only for a short interval. The GI is not found in this case, but it might be seen in a high-resolution simulation.

In the simulation with addition of the self-gravity in Section 4.3, we assumed Σ_g0 that corresponds to ∼2Σ_MMSN at r = 5 AU. On the other hand, simulations before adding self-gravity are scaled by Σ_g0r²/M_* and H/r.

Here, fixing r = 5 AU and H/r = 0.055, we apply the results of individual runs for various Σ_g0 to derive a sufficient condition for gas column density to cause the GI. The condition for the GI is R < R_J(R). Since R_J∝Σ_g0 (Equation (15)), the condition is more easily satisfied for larger Σ_g0. In the saturated state in model-s40-t10-e010, the condition is satisfied even at Σ_g0/Σ_MMSN ∼ 1, while we showed the results with Σ_g0 = 2Σ_MMSN in Section 4.3. For smaller ₀ (model-s40-t10-e001), the critical column density is Σ_g0/Σ_MMSN ∼ 3. The smaller dust particles with τ_fΩ = 0.1 have no chance to excite the GI even in the weak residual turbulence (model-s40-t01-e010) as long as Σ_g0/Σ_MMSN < 20. In the stronger remnant turbulence (L_s = 0.5H), the GI is not expected unless Σ_g0/Σ_MMSN > 10 (model-s05-t10-e010).

Note that R_J is a function of ₀Σ_g0 (Equation (15)). That is, the possibility of the GI depends on the total column density of dust particles, but not on ₀. Thus, for example, R_J should be similar between the result with Σ_g0/Σ_MMSN = 1 in model-Ls40-t10-e010 (₀ = 0.10) and Σ_g0/Σ_MMSN = 10 in model-Ls40-t10-e001 (₀ = 0.01) at the same r. From the results of simulations in this paper, it is inferred that the GI may occur when ρ_{d0, crit} ≳ 0.03ρ_{g, MMSN}.