N-BODY SIMULATION OF PLANETESIMAL FORMATION THROUGH GRAVITATIONAL INSTABILITY AND COAGULATION. II. ACCRETION MODEL

According to the standard theory of planet formation, planets form from planetesimals, which are kilometer-sized solid bodies. However, the process of planetesimal formation is still controversial. The inward migration of meter-sized bodies due to gas drag is very rapid; its timescale is only 10² yr (Adachi et al. 1976; Weidenschilling 1977). The growth of dust to planetesimals due to particle adhering seems difficult in the standard disk model.

The gravitational instability scenario is an alternative scenario for this stage (Safronov 1969; Goldreich & Ward 1973). Dust particles settle into the midplane of a protoplanetary disk owing to the gravitational force from the central star if the gas flow is laminar. As the sedimentation of dust particles proceeds, the density at the midplane exceeds the critical density, and the dust layer finally becomes gravitationally unstable. The gravitationally unstable dust layer rapidly collapses by its self-gravity, forming km-sized planetesimals (Goldreich & Ward 1973; Coradini et al. 1981; Sekiya 1983). The rapid migration of meter-sized bodies can be avoided because the timescale of the gravitational instability is only on the order of a Kepler period. However, the gravitational instability scenario has a critical issue. As the sedimentation of dust grains toward the midplane proceeds, the vertical velocity shear increases and gives rise to the Kelvin–Helmholtz instability, which makes the dust layer turbulent. The turbulence prevents dust particles from settling. A lot of studies have been carried out on this issue (Weidenschilling 1980; Cuzzi et al. 1993; Champney et al. 1995; Weidenschilling 1995; Sekiya 1998; Dobrovolskis et al. 1999; Sekiya & Ishitsu 2000, 2001; Ishitsu & Sekiya 2002, 2003; Michikoshi & Inutsuka 2006). However, problems concerning the gravitational instability still remain unsettled.

Although many studies have been conducted on the condition of gravitational instability or its linear regime, there has been little study of the nonlinear regime of gravitational instability. Here we concentrate on the formation process of planetesimals assuming the gravitational instability. Tanga et al. (2004) performed N-body simulations of a gravitationally unstable particle disk. In their simulation, the optical depth is very small, and the particle disk is unstable, even initially. If we consider the formation of planetesimals through gravitational instability at a₀ ≃ 1 AU, the optical depth is high, and the initial velocity dispersion must be large owing to turbulent flow driven by the Kelvin–Helmholtz instability or magnetorotational instability (Balbus & Hawley 1991). In our analysis, we will take a close look at initially stable disks with large optical depth. Johansen et al. (2007) performed simulations using a code that solves the magnetohydrodynamic equations with a three-dimensional grid and includes particles with self-gravity. They mapped the particle density on the grid and solved the Poisson equation using a fast Fourier transform method. They included collisions among particles as damping the velocity dispersion of the particles within each grid. They found that particles collapse in an overdense region in the midplane and the gravitationally bound objects form with masses comparable to dwarf planets.

Michikoshi et al. (2007) performed a set of simulations in the self-gravitating collision-dominated particle disks without gas components (hereafter Paper I). The point-to-point Newtonian mutual interaction was calculated directly. They adopted the hard and soft sphere models as collision models, and examined parameter dependencies on the size of computational domain, restitution coefficient, optical depth, Hill radius of particles, and the duration of collision. In the hard sphere model, penetration between particles is not possible (e.g., Wisdom & Tremaine 1988; Salo 1991; Richardson 1994; Daisaka & Ida 1999; Daisaka et al. 2001). The velocity changes in an instant at collisions. In the soft sphere model, particles can overlap if they are in contact (e.g., Salo 1995). They found that the formation process of planetesimals is qualitatively independent of simulation parameters if initial Toomre's Q value, Q_init > 2. The formation process is divided into three stages: the formation of nonaxisymmetric wake-like structures (Salo 1995), the creation of aggregates, and the rapid collisional growth of the aggregates. The mass of the largest aggregates is larger than the mass predicted by the linear perturbation theory (hereafter linear mass M_linear; Goldreich & Ward 1973). However Michikoshi et al. (2007) could not find the saturation of the growth of the planetesimals in Paper I. The mass of the planetesimal formed in the simulation depends on the size of the computational domain. Almost all mass in the computational domain is absorbed by a single planetesimal, and thus the mass of the planetesimal is determined by the total mass in the computational domain.

In this paper, we use the alternative model of collisions, "accretion model." In physical terms, the accretion model corresponds to the more dissipative model than the hard and soft sphere models. In the accretion model, when two particles collide, if the binding condition is satisfied, the two particles are merged into one particle. Efficient energy dissipation occurs when two particles merge into one particle. In the soft and hard sphere models, since large aggregates consist of a lot of small particles and the number of particles does not decrease, the calculation is time-consuming. If we are not interested in internal states of planetesimals such as rotation or internal density profile, we can treat a large aggregate as one large particle. We expect that the final mass or spatial distribution of planetesimals can be estimated by the accretion model. The accretion model has an advantage from the viewpoint of calculation. The number of particles decreases as the calculation proceeds; thereby this model enables us to perform simulations of larger numbers of particles than that in the previous work. Taking this advantage, we can examine the large computational domains.

In this study, we neglect the effect of gas for the sake of simplicity. As many authors investigated, the effect of gas is very important. The solid particles drift radially due to gas drag (Adachi et al. 1976; Weidenschilling 1977), gas friction helps gravitational collapse (Ward 1976; Youdin 2005a, 2005b), gas turbulence prevents and helps concentration of solids (Weidenschilling & Cuzzi 1993; Barge & Sommeria 1995; Fromang & Nelson 2006; Johansen et al. 2006a), and clumps of solid particles form due to streaming instability (Youdin & Goodman 2005; Johansen et al. 2006b; Johansen et al. 2007). But once gravitational instability happens and large aggregates form, as a first step, we may be able to neglect gas drag because they are so large that they decouple from gas. Therefore our gas-free model may be applicable to the stages of gravitational instability and subsequent collisional growth.

In Section 2, we explain the simulation method. Our results are presented in Section 3. We compare the accretion model with the soft and hard sphere models and investigate the models of large computational domain. In Section 4, we summarize the results.

The method of calculation is the same as that used in Paper I except for the collision model, which is based on the method of simulation of planetary rings (e.g., Wisdom & Tremaine 1988; Salo 1991; Daisaka & Ida 1999; Ohtsuki & Emori 2000).

We introduce rotating Cartesian coordinates (x, y, z) called Hill coordinates, the x-axis pointing radially outward, the y-axis pointing in the direction of rotation, and the z-axis normal to the equatorial plane. The origin of the coordinates moves on a circular orbit with the semimajor axis a₀ and the Keplerian angular velocity Ω₀. Assuming |x_j|, |y_j|, |z_j| ≪ a₀, where (x_j, y_j, z_j) is the position of jth particle, we can write the equation of motion in nondimensional forms independent of the mass and the semimajor axis a₀, if we scale the time by Ω⁻¹₀, the length by Hill radius ha₀, and the mass by h³M_*, where h = (2m_p/3M_*)^1/3, m_p is the characteristic mass of particles, and M_* is the mass of the central star (e.g., Petit & Henon 1986; Nakazawa & Ida 1988):

$$\begin{equation} \begin{array}{c} \displaystyle{\frac{d^2 \tilde{x}_i}{d\tilde{t}^2}} = \displaystyle {2\frac{d \tilde{y}_i}{d\tilde{t}}} +3 \tilde{x}_i + \displaystyle {\sum _{j=1,j\ne i}^N \frac{\tilde{m}_j}{\tilde{r}_{ij}^3}(\tilde{x}_j- \tilde{x}_i)}, \\ \displaystyle {\frac{d^2 \tilde{y}_i}{d\tilde{t}^2}} = \displaystyle {-2\frac{d\tilde{x}_i}{d\tilde{t}}} + \displaystyle {\sum _{j=1,j\ne i}^N \frac{\tilde{m}_j}{\tilde{r}_{ij}^3}(\tilde{y}_j-\tilde{y}_i)}, \\ \displaystyle {\frac{d^2 \tilde{z}_i}{d\tilde{t}^2}} = -\tilde{z}_i + \displaystyle {\sum _{j=1,j\ne i}^N \frac{\tilde{m}_j}{\tilde{r}_{ij}^3}(\tilde{y}_j-\tilde{y}_i)}, \\ \end{array} \end{equation} \tag{ 1 }$$

where a tilde denotes corresponding nondimensional variables, $\tilde{r}_{ij}$ is the distance between particles i and j, and N is the number of particles.

We use the periodic boundary condition (Wisdom & Tremaine 1988). There are an active region and its copies. Inner and outer regions have different angular velocities. These regions slide upward and downward with shear velocity 3Ω₀L_x/2, where L_x is the length of the cell in the x-direction. We calculate gravitational forces from particles in these regions. We use the special-purpose computer GRAPE-6 for the calculation of self-gravity (Makino et al. 2003).

Particles used in the simulation are not realistic dust particles but superparticles that represent a group of small particles; thus, the parameter used in this paper is not exactly the physical restitution coefficient but corresponds to the coefficient of the dissipation due to collisions among particles. The collision models used in Paper I cannot treat more dissipative models where the tangential coefficients of restitution _t < 1 because we adopt _t = 1. The net restitution coefficient is = (²_n〈v²_n〉 + ²_t〈v²_t〉)/(〈v²_n〉 + 〈v²_t〉), where _n is the normal coefficient of restitution, v_n and v_t are normal and tangential velocities of a collision (Canup & Esposito 1995). The net restitution coefficient is always larger than 0.7 if we assume _t = 1 and 〈v²_n〉 = 〈v²_t〉. If we allow _t ≠ 1 in order to examine the model of < 0.7, we must consider rotations of individual particles. For the sake of simplicity, we ignore rotations of particles in this paper and assume _t = 1.

The way to treat a collision is to use the same method applied to the hard sphere model, but we consider the binding condition of colliding particles. If a pair overlaps, $\tilde{r}_{ij} < \tilde{r}_{\mathrm{p},i} + \tilde{r}_{\mathrm{p},j}$, where $\tilde{r}_{\mathrm{p},i}$ is a radius of particle i, and is approaching, $\tilde{\mathbf n}_{ij} \cdot \tilde{\mathbf v}_{ij}<0$, where $\tilde{\mathbf v}_{ij}$ is the relative velocity and $\tilde{\mathbf n}_{ij}$ is a unit vector along $\tilde{\mathbf r}_{ij}=(\tilde{x}_j, \tilde{y}_j,\tilde{z}_j) - (\tilde{x}_i,\tilde{y}_i,\tilde{z}_i)$, we assume that this pair collides, and changes its velocities by using the equations of collision:

$$\begin{equation} \tilde{\mathbf v}{_i^{\prime }}=\tilde{\mathbf v}{_i}-\frac{\tilde{m}_j}{\tilde{m}_i+\tilde{m}_j}(1+\epsilon)(\tilde{\mathbf n}_{ij}\cdot \tilde{\mathbf v}_{ij}) \tilde{\mathbf n}_{ij}, \end{equation} \tag{ 2 }$$

$$\begin{equation} \tilde{\mathbf v}_j^{\prime }=\tilde{\mathbf v}{_j}+\frac{\tilde{m}_i}{\tilde{m}_i+\tilde{m}_j}(1+\epsilon)(\tilde{\mathbf n}_{ij}\cdot \tilde{\mathbf v}_{ij}) \tilde{\mathbf n}_{ij}, \end{equation} \tag{ 3 }$$

where is a restitution coefficient. Using the updated velocities, we check whether the binding condition is satisfied by the Jacobi integral (e.g., Nakazawa & Ida 1988):

$$\begin{equation} \tilde{J}_{ij}= \frac{1}{2}\big(\dot{\tilde{x}_{ij}}^2+\dot{\tilde{y}_{ij}}^2+\dot{\tilde{z}_{ij}}^2\big)-\frac{3}{2}\tilde{x}_{ij}^2+\frac{1}{2}\tilde{z}_{ij}^2-\frac{3}{\tilde{r}_{ij}}+\frac{9}{2} < 0, \end{equation} \tag{ 4 }$$

where $(\dot{\tilde{x}}_{ij},\dot{\tilde{y}}_{ij},\dot{\tilde{z}}_{ij})$ is the relative velocity. If the pair is bound, two particles are merged into one particle conserving their momentum. The shape of the new particle is a sphere, the radius of which is determined by its mass.

We assume that all particles have the same mass m_p initially. The parameters of the present simulation are the distance from the central star a₀, the length of region L_x and L_y, the number of particles N, and the mass of particles m_p(r_p, ρ_p), where ρ_p is the density of a particle. The dynamical behavior is characterized by only two nondimensional parameters (Daisaka & Ida 1999); the initial optical depth τ and the ratio ζ = a₀h/2r_p. In the models of ζ>1, the centers of the particles are in their Hill sphere when two particles come into contact. The initial optical depth is given by Goldreich & Tremaine (e.g., 1982)

$$\begin{eqnarray} \tau &=&\frac{3\Sigma }{4 \rho _\mathrm{p} r_\mathrm{p}}=0.19 \left(\frac{\Sigma }{10 \;\mathrm{g\;cm}^{-2}} \right) \left(\frac{\rho _\mathrm{p}}{2 \;\mathrm{g\;cm}^{-3}} \right)^{-1}\nonumber\\ &&\times \left(\frac{r_\mathrm{p}}{20 \;\mathrm{cm}} \right)^{-1}, \end{eqnarray} \tag{ 5 }$$

where Σ is the surface density of particles. The ratio ζ is given by

$$\begin{equation} \zeta =105.528 \left(\frac{M_\odot }{M_* }\right)^{-1/3} \left(\frac{\rho _\mathrm{p}}{2 \;\mathrm{g}\;\mathrm{cm}^{-3}}\right)^{1/3} \left(\frac{a_0}{1 \;\mathrm{AU}}\right). \end{equation} \tag{ 6 }$$

Using the above two parameters, the other parameters can be estimated. The normalized most unstable wavelength of gravitational instability is (e.g., Sekiya 1983)

$$\begin{equation} \tilde{\lambda }_{\mathrm{most}}=12 \pi \tau \zeta ^2. \end{equation} \tag{ 7 }$$

The normalized size of computational domain $\tilde{L}$, the number of particles N, the normalized radius of particle $\tilde{r}_\mathrm{p}$, and Toomre's Q value are given by

$$\begin{equation} \tilde{L}_x = 12 \pi A_x \tau \zeta ^2, \end{equation} \tag{ 8 }$$

$$\begin{equation} \tilde{L}_y = 12 \pi A_y \tau \zeta ^2, \end{equation} \tag{ 9 }$$

$$\begin{equation} N=576 \pi A_x A_y \tau ^3 \zeta ^6, \end{equation} \tag{ 10 }$$

$$\begin{equation} \tilde{r}_\mathrm{p} = \frac{1}{2 \zeta }, \end{equation} \tag{ 11 }$$

$$\begin{equation} Q= \frac{\tilde{\sigma }_x}{6 \tau \zeta ^2}, \end{equation} \tag{ 12 }$$

where A_x is the ratio L_x/λ_most, A_y is the ratio L_y/λ_most, and $\tilde{\sigma }_x$ is the radial velocity dispersion scaled by a₀hΩ₀, respectively. In order to determine the size of the computational domain, we need to set A_x and A_y. The initial velocity dispersion is also a parameter. We determine it from the initial Toomre's Q_init value using Equation (12). As shown above, if we set τ, ζ, A_x, A_y, and Q_init, we can determine the initial state. The simulation parameters are summarized in Table 1. We call model 1 as the standard model. The number of particles is 1000–10,000.

Table 1. Initial Conditions and Results

Model		τ	ζ	Ax	Ay	Qinit	σx	〈Mpl/Mlinear〉	Npl
01	0.01	0.1	2.0	6.0	6.0	3.0	7.2	21.2	2
02	0.01	0.11	2.0	6.0	6.0	3.0	7.9	42.4	1
03	0.01	0.125	2.0	6.0	6.0	3.0	9.0	21.5	2
04	0.1	0.1	2.0	6.0	6.0	3.0	7.2	42.7	1
05	0.3	0.1	2.0	6.0	6.0	3.0	7.2	43.2	1
06	0.5	0.1	2.0	6.0	6.0	3.0	7.2	38.4	1
07	0.7	0.1	2.0	6.0	6.0	3.0	7.2	0.0	0
08	0.9	0.1	2.0	6.0	6.0	3.0	7.2	0.0	0
09	0.01	0.06	2.5	6.0	6.0	3.0	6.8	43.4	1
10	0.01	0.06	2.7	6.0	6.0	3.0	8.2	43.4	1
11	0.01	0.06	3.0	6.0	6.0	3.0	9.7	39.8	1
12	0.01	0.1	2.0	4.0	4.0	3.0	7.2	19.5	1
13	0.01	0.1	2.0	8.0	8.0	3.0	7.2	38.4	2
14	0.01	0.1	2.0	3.0	3.0	3.0	7.2	11.1	1
15	0.01	0.1	2.0	6.0	3.0	3.0	7.2	21.7	1
16	0.01	0.1	2.0	9.0	3.0	3.0	7.2	16.4	2
17	0.01	0.1	2.0	12.0	3.0	3.0	7.2	11.0	4
18	0.01	0.1	2.0	24.0	3.0	3.0	7.2	10.1	9
19	0.01	0.1	2.0	3.0	6.0	3.0	7.2	21.8	1
20	0.01	0.1	2.0	3.0	9.0	3.0	7.2	32.2	1
21	0.01	0.1	2.0	3.0	12.0	3.0	7.2	43.1	1
22	0.01	0.1	2.0	3.0	24.0	3.0	7.2	85.5	1
23	0.01	0.1	2.0	6.0	6.0	4.0	9.6	43.3	1
24	0.01	0.1	2.0	6.0	6.0	5.0	12.0	21.6	2

Notes. 〈M_pl〉, M_linear and N_pl are the mean planetesimal mass, the linear mass and the number of planetesimals at final state, respectively.

Download table as:  ASCIITypeset image

If ζ>1.5, two particles can be bound (Ohtsuki 1993; Salo 1995), so we need to adopt ζ>1.5 in order to investigate the formation process of planetesimals. In this paper, we use ζ = 2, which is much smaller than the realistic value ζ ≃ 100. Thus the physical size of a particle in our simulation is very large. It indicates that the effect of the self-gravity is relatively weaker than that of inelastic collisions. To investigate gravitational instability, the particle density should be higher than the Roche density. The particle density is ρ_p = 4ζ³ρ_R, where ρ_R = (9/4π)M_*/a³₀ is the Roche density. Thus when ζ = 2, ρ_p ≫ ρ_R. So we can study gravitational instability with this parameter. But this approximation may change the result when comparing between ζ = 2 and ζ = 100, although gravitational instability occurs in both models. In the future work, we will perform the simulation with a larger ζ value by using next-generation supercomputers.

In most of the models, we set Q_init = 3, which corresponds to $\tilde{\sigma }_x \simeq 7.2$ for the standard model. The characteristic velocity of the actual turbulence v_turb is about ηv_K, where v_K is the Kepler velocity, η = −12(c²_s/v²_K)(∂log P/∂log r), c_s is the sound velocity of gas, P is the pressure of gas, and r is the distance from the central star (Adachi et al. 1976; Weidenschilling 1977). At r = 1 AU, in the Hill unit, the characteristic velocity of the turbulence is $\tilde{v}_\mathrm{turb} \simeq 1.3 \times 10^7$ that is much larger than $\tilde{\sigma }_x$. If the turbulence weakens, Q decreases gradually from Q ≫ 1, and gravitational instability occurs finally. Therefore, we choose the initial condition that the disk is gravitationally stable (Q_init = 3).

We set the initial velocities and positions of particles using the above parameters. A Keplerian orbit has six parameters: the position of the guiding center (x_g, y_g), eccentricity e, inclination i, and the two phases of epicycle and vertical motions (e.g., Nakazawa & Ida 1988). The eccentricity e and inclination i of particles are assumed to obey the Rayleigh distribution with the ratio $\sqrt{\langle e^2 \rangle / \langle i^2 \rangle }=2$ (Ida & Makino 1992). We set the position of the guiding center and the two phases uniformly, avoiding overlapping.

The equation of motion is integrated with a second-order leapfrog scheme. We adopt the variable time step used by Daisaka & Ida (1999). The time step formula is given by $\Delta t = \eta \min _{i} |\mathbf {a}_i|/|\dot{\mathbf a}_i|$, where a_i and $\dot{\mathbf a}_i$ are the acceleration of particle i and its time derivative, and η is a nondimensional time step coefficient, respectively.

3.1. Time Evolution

Figure 1 shows snapshots in model 1 at t/T_K = 0.0, 0.4, 0.8, 1.2, 1.6, and 2.0 where T_K is the orbital period 2π/Ω₀. At t/T_K = 0.4, we cannot see any spatial structures and large bodies. The formation process of planetesimals is the same as those of the hard and soft sphere models (Paper I). The nonaxisymmetric wake-like structure forms at t/T_K = 0.8, 1.2, which appears if we consider self-gravity. The nonaxisymmetric wake-like structure is also seen in planetary rings (Salo 1995). The density in the wake-like structure is higher than the mean surface density Σ. For example, the surface density in the dense region is 1.52Σ at t/T_K = 1.6. Planetesimal seeds form in the dense parts of these wakes by fragmentation (t/T_K = 1.6, 2.0). Here we consider the formation of the wake-like structure and planetesimal seeds as the gravitational instability stage.

Zoom In Zoom Out Reset image size

Figure 2 shows snapshots in model 1 at t/T_K = 2.0, 3.0, 4.0, 5.0, 6.0, and 7.0. Once planetesimal seeds form, the nonaxisymmetric wake-like structures disappears (t/T_K = 3.0). Planetesimal seeds grow rapidly by mutual coalescence (t/T_K = 3.0, 4.0, 5.0). In this stage, the disk is gravitationally stable. The gravitation seems to merely enhance the coalescent growth. When almost all planetesimal seeds are absorbed by a few planetesimals, the growth slows down (t/T_K = 6.0, 7.0). This stage is the collisional growth stage.

Zoom In Zoom Out Reset image size

Figure 3 shows the time evolution of the number of large particles, the mass of the largest and second largest particles, the ratio of the mass of large particles to the total mass, and the velocity dispersion of field particles for model 1. The time evolution of the number of large particles is similar to those of the hard and soft sphere models (Paper I). The number of large particles has a maximum at t/T_K ≃ 1.7, and its value is about 1.2, which is slightly larger than that in the hard and soft sphere models. However, the number of large particles decreases rapidly after t/T_K ≃ 1.7. This decay is caused by the fast coalescence of large particles.

Zoom In Zoom Out Reset image size

The mass of the largest particle monotonically increases up to about 44 M_linear, which is twice larger than that of the hard and soft sphere models, where M_linear = πΣ(λ_most/2)² is the planetesimal mass estimated by the linear theory. In the accretion model, the growth of planetesimals is enhanced because sticking of dust grains is efficient. The mass of the second largest particle changes discontinuously, when the second largest particle is absorbed by the largest one.

The ratio of the mass of the large particles to the total mass is almost the same as that of the hard and soft sphere models. The ratio monotonically increases up to about 1, and most of the mass in the computational domain is finally absorbed by the large planetesimals.

The velocity dispersion decreases until t/T_K ≃ 1.3 because of the dissipation due to inelastic collisions. When the velocity dispersion becomes sufficiently small (Q ≃ 2), the gravitational instability occurs and a lot of planetesimal seeds form. As the field particles are scattered by the planetesimal seeds, the velocity dispersion increases.

The strength of self-gravity is measured by two dimensionless quantities: (e.g., Youdin 2005a),

$$\begin{equation} Q = \frac{\sigma _x \Omega _0}{\pi G \Sigma }, \end{equation} \tag{ 13 }$$

$$\begin{equation} Q_R = \frac{ \rho _\mathrm{R}}{ \rho } \simeq \frac{9}{4\pi } \frac{h_d \Omega _0^2}{G \Sigma }, \end{equation} \tag{ 14 }$$

where ρ is the mass density of the dust layer, and h_d is the thickness of the dust layer. The value Q_R is the ratio of the Roche density to the dust density. If these values are sufficiently small, the disk is gravitationally unstable. Figure 4 shows time evolution of Q and Q_R. The value Q_R is always larger than Toomre's Q value. The ratio Q/Q_R is approximately equal to 0.5. If the hydrostatic balance (σ_x ≃ Ω₀h_d) can be assumed, Q/Q_R = 4/9 ≃ 0.44.

Zoom In Zoom Out Reset image size

If the gravitational instability occurs, there are no remarkable qualitative differences in the planetesimal formation process among different collision models. The growth of particles in the accretion model is more efficient than in the hard and soft sphere models; thus, the maximum mass of planetesimals is slightly larger than that in the hard and soft sphere models.

3.2. Parameter Dependence

3.2.1. Physical Parameters

We summarize the parameter dependences of the simulation results on the restitution coefficient , the optical depth τ, the ratio of Hill radius to the particle diameter ζ, the initial Toomre's Q value. We focus on the mass of the largest particle, and Toomre's Q value.

Figure 5 shows the dependence on the restitution coefficient . We vary the restitution coefficient from 0.01 to 0.9 (models 1, 4, 5, 6, 7, 8). No gravitational instability occurs; in other words, Q never becomes less than about 2, for = 0.9, and no large body forms. There are two reasons why planetesimals do not form. The velocity dispersion increases monotonically if is larger than the critical value _c ≃ 0.7 (Goldreich & Tremaine 1982; Salo 1995; Daisaka & Ida 1999; Ohtsuki 1999). The velocity dispersion decreases due to inelastic collisions and increases due to gravitational scattering. If the restitution coefficient is large, the dissipation is relatively weak; thus, the velocity dispersion does not decrease. Because Q value does not reach the critical value, gravitational instability does not occur. Therefore, planetesimal seeds do not form, and the subsequent collisional growth does not happen. We cannot apply the formation process of planetesimal stated in our paper to the model with = 0.9. With a larger restitution coefficient, the reduction of the relative velocity on collisions is smaller. In addition, the collision velocity increases as the velocity dispersion increases. The coagulation is not efficient with the large restitution coefficient. Thus planetesimals do not form by coagulation within the simulation time.

Zoom In Zoom Out Reset image size

We study the effect of the optical depth τ by comparing results with τ = 0.1, 0.11, and 0.125 (models 1, 2, 3). The dependence on the optical depth τ is shown in Figure 6. The results are similar to those of the hard and soft sphere models (Paper I). The largest masses and Toomre's Q value are on the same order.

Zoom In Zoom Out Reset image size

We set ζ = 2.5, 2.75, and 3.0, and fix the other parameters in the standard model except for the optical depth τ = 0.06 (models 9, 10, 11). The dependence on ζ is shown in Figure 7. The result is also similar to that of the hard and soft sphere models (Paper I). The largest masses are M_max/M_linear = 43.3, 43.2, and 39.8, respectively. The largest masses are almost the same, which is about 40M_linear. Figure 7 shows the similar time evolution of Toomre's Q value because the masses of planetesimals are almost the same.

Zoom In Zoom Out Reset image size

Figure 8 shows the time evolution of Toomre's Q values for Q_init = 3.0, 4.0, and 5.0 (models 1, 23, 24). If the growth is due to gravitational instability, the results do not depend on the initial Toomre's Q value Q_init. Time evolutions are similar for these models. Toomre's Q value decreases through inelastic collisions between particles. As Q decreases, self-gravity becomes stronger. When Q becomes less than about 2, the disk becomes gravitationally unstable and the wake-like structure develops. Then Q increases due to scattering of field particles by planetesimal seeds.

Zoom In Zoom Out Reset image size

In Paper I, we showed that there is no remarkable difference between the hard and soft sphere models on the formation process of planetesimals. We confirm that the results of the accretion model are also essentially the same as these of the hard and soft sphere models. The dissipation rate due to collision is a crucial parameter, but the planetesimal formation process is independent of the collision model. Before the gravitational instability, collisions merely decrease the velocity dispersion. In the late stage, coalescence among planetesimal seeds occurs, but this process is independent of collision models.

3.2.2. Simulation Domain Size

In the hard and soft sphere models, the mass of the largest planetesimal depends on the size of the computational domain, and we could not find the saturation of growth (Paper I). The growth of planetesimals continues until most of the mass in the computational domain is absorbed by the largest planetesimal. We expect to find the saturation of growth because the accretion model enables us to perform a wider simulation with a large number of particles.

First we assume that the ratio of the length of the computational domain in the x-direction to the y-direction is unity; in other words, the shape of the computational domain is a square. We vary the size of the computational domain A = 4 to 8. Results are shown in Figure 9 (models 12, 1, 13). The final masses of the largest planetesimals and Toomre's Q values are M_max/M_linear = 19.5, 44.0, 66.5, Q = 13.6, 18.2, 21.9, respectively. As the size of the simulation domain increases, the final planetesimal mass and the final Toomre's Q value increase.

Zoom In Zoom Out Reset image size

Next, we vary the ratio of the length of the computational domain in the x- and y-directions. Now the shape of the computational domain is rectangular. We vary the length in the y-direction, A_y, from 3–24 fixing the length in the x-direction, A_x, (the left panel in Figure 10; models 14, 15, 16, 17, 18), and we vary A_x from 3–24 fixing A_y (the right panel in Figure 10; models 14, 19, 20, 21, 22). The masses of the largest planetesimal and Toomre's Q values are M_max/M_linear = 11.1, 21.7, 18.8, 15.9, 16.8, Q = 10.4, 13.6, 10.2, 12.2, 11.3 (models 14, 15, 16, 17, 18) and M_max/M_linear = 11.1, 21.8, 32.2, 43.1, 85.5, Q = 10.4, 17.1, 15.9, 21.3, 24.1 (models 14, 19, 20, 21, 22), respectively. In the model with fixed A_x, the mass of the largest planetesimal increases with A_y, while in the model with fixed A_y, the largest mass is roughly constant. These results indicate that the mass of the largest planetesimal is determined by A_y. If A_x is sufficiently large, multiple planetesimals form. Toomre's Q value increases with A_y because of the scattering of field particles by the large bodies.

Zoom In Zoom Out Reset image size

Here we estimate the final mass of a planetesimal, M_pl, in the computational domain in a similar way to estimate the isolation mass of protoplanets in the oligarchic growth picture (Kokubo & Ida 1998). We assume that particles whose orbital radii are within γr_H of the planetesimal are finally absorbed by the planetesimal, though they are not in the Hill sphere of the planetesimal initially, where γ is a nondimensional factor and r_H is the Hill radius of the planetesimal. Therefore, the planetesimal absorbs particles in the rectangle of 2γr_H in width and L_y in length. The Hill radius of the planetesimal is r_H = a₀(2M_pl/3M_*)^1/3, and we define

$$\begin{equation} M_\mathrm{H}=2 \gamma r_{\mathrm{H}} L_y \Sigma. \end{equation} \tag{ 15 }$$

The planetesimal with the mass M_pl can absorb the mass M_H, and M_H increases as the planetesimal grows. When M_pl ≃ M_H, the growth becomes slow. We estimate the mass of the planetesimal by this condition. In this case the number of planetesimals is about L_x/(2γr_H). We obtained the mass of the planetesimal

$$\begin{equation} \tilde{M}_\mathrm{pl}=\min (576 \sqrt{6} \zeta ^6 \tau ^3 (\gamma A_y)^{3/2}, \tilde{M}_\mathrm{total}), \end{equation} \tag{ 16 }$$

where M_total are the total mass in the computational domain. The number of the planetesimals is N_pl is

$$\begin{equation} N_\mathrm{pl}= \max \left(\pi \sqrt{\frac{3}{8}} \frac{A_x }{\gamma ^{3/2}A_y^{1/2}},1 \right). \end{equation} \tag{ 17 }$$

The mass of the planetesimals is equal to the total mass in the computational domain and the number of the planetesimal is unity if the following condition is satisfied:

$$\begin{equation} A_x < \frac{2 \sqrt{6}}{3 \pi } A_y^{1/2} \gamma ^{3/2}. \end{equation} \tag{ 18 }$$

The comparisons between the analytical estimation and the simulations are shown in Figures 11 and 12. If we assume γ = 2.5, the analytical estimation agrees with the simulation results. The parameter γ = 2.5 indicates that the mean orbital separation of planetesimals is 5r_H. This is narrower than 10r_H, which is the typical orbital separation of protoplanets caused by oligarchic growth (Kokubo & Ida 1998, 2000). This narrow orbital separation is possible for a dissipative system. We will discuss this issue later. For instance, Figure 13 shows the final state in model 18. In this Figure, there are 9 planetesimals in the computational domain whose length in the x-direction is 400. Thus the mean orbital separation is about 44. The Hill radius of the mean planetesimal (〈M_pl/M_linear〉 ≃ 10.1) is about 9.5. Therefore the separation is estimated to about 9.5 × 5 = 48, which is consistent with Figure 13.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Our analytical estimation gives good agreement with the numerical results in each model. Figure 11 shows the mean mass and the number of planetesimals as a function of the size of the computational domain in the model where the domain is square. The mean mass of planetesimals increases with the size of the computational domain. From the Equation (18), when A_x = A_y < 4.2, the number of planetesimals is equal to unity and the mass of the planetesimal is approximately equal to the total mass in the computational domain. The slope of the line of the mean mass of planetesimals changes at A_x = A_y = 4.2 and the number of the planetesimals is larger than unity if A_x = A_y > 4.2.

Figure 12 shows the result for the model of the rectangular computational domain. In the models of A_y = 3, the mean mass of planetesimals is constant and the number of planetesimals is larger than unity if A_x > 3.6. In the models of A_x = 3, the slopes of the lines change at A_y = 2.1. The growth of planetesimals cannot be saturated although A_y is large. In these models, the Hill radius of the planetesimal is longer than the size of the computational domain in the x-direction. The only one planetesimal forms, and it absorbs most mass in the computational domain. It should be remembered that the local calculation of planetesimal formation has the domain-size dependence.

If we assume that the surface density of dust is

$$\begin{equation} \Sigma = 10 \times f_\mathrm{d} \left(\frac{a_0}{1 \;\mathrm{AU}} \right)^{-3/2} \;\mathrm{g}\;\mathrm{cm}^{-2}, \end{equation} \tag{ 19 }$$

where f_d is a scaling factor, the linear mass M_linear is given by

$$\begin{equation} M_\mathrm{linear}=3.5\times 10^{18} f_\mathrm{d}^3 \left(\frac{a_0}{1 \;\mathrm{AU}}\right)^{3/2} \left(\frac{M_*}{M_\odot }\right)^{-2} \;\mathrm{g}, \end{equation} \tag{ 20 }$$

and from Equation (16) the estimated planetesimal mass M_pl is given by

$$\begin{eqnarray} M_\mathrm{pl}&=&9.0\times 10^{18} f_\mathrm{d}^3 \left(\frac{a_0}{1 \;\mathrm{AU}}\right)^{3/2} \left(\frac{M_*}{M_\odot }\right)^{-2} \left(\frac{\gamma }{2.5}\right)^{3/2} A_y^{3/2} \;\mathrm{g} \nonumber\\ &=& 2.6 \left(\frac{\gamma }{2.5}\right)^{3/2} A_y^{3/2} M_\mathrm{linear}, \end{eqnarray} \tag{ 21 }$$

where we assume A_x > 2.1A^1/2_y(γ/2.5)^3/2. The factor f_d = 1 roughly corresponds to the minimum mass solar nebula model inside the snow line (Hayashi 1981). From Equation (21), the estimated planetesimal mass becomes larger than the linear mass if the feeding zone is large, (A_yγ)^3/2 > 1.6.

The estimated and linear masses exhibit the same dependencies on f_d, a₀, and M_*, and the estimated planetesimal mass is on the same order of the linear mass when A_yγ ≃ 1. These results are explained as follows. Here we rewrite the Hill radius by the most unstable wavelength

$$\begin{equation} r_\mathrm{H} = \frac{\lambda _\mathrm{most}}{(3 \pi)^{1/3}} \left(\frac{M_\mathrm{pl}}{M_\mathrm{linear}} \right)^{1/3}, \end{equation} \tag{ 22 }$$

where the most unstable wavelength for Q = 1 is

$$\begin{equation} \lambda _\mathrm{most} = \frac{2\pi ^2 G \Sigma }{\Omega _0^2}. \end{equation} \tag{ 23 }$$

The most unstable wavelength of the gravitational instability and the Hill radius of the linear mass are on the same order. The condition Q = 1 corresponds to the balance between the self-gravity and the effect of the rotation. On the other hand, the Hill radius is determined by the balance between the self-gravity and the tidal force. Thus, r_H and λ_most are on the same order, and the planetesimal mass is on the order of the linear mass only if A_yγ ≃ 1. In general, because A_y ≫ 1, the planetesimal mass is larger than the linear mass.

Collisions among planetesimals are more frequent than the realistic case in the present simulations since the bulk density of super-particle is much smaller than the realistic value. Therefore the orbital separation of large planetesimals is different from that of the standard oligarchic growth derived by N-body simulations of planetesimals (Kokubo & Ida 1998, 2000). If the number density is low, the orbital separation is determined by the balance between the orbital repulsion and their growth (Kokubo & Ida 1998). But if the number density is high, and the collisional dissipation is effective, the orbital repulsion is relatively weak. Thus, in the collisional oligarchic growth, the orbital separation can be narrower than 10R_H (Goldreich et al. 2004).

We performed local N-body simulations of planetesimal formation through gravitational instability and collisional growth without a gas component using a shearing box method. The accretion model was adopted as the collision model. In the accretion model, when particles collide, if the binding condition is satisfied, two particles merge into one. The number of particles decreases as the calculation proceeds; thus, this model enables us to perform the long-term and large-scale calculations. We compared the results of the accretion model with those of the hard and soft sphere models used in Paper I.

The formation process of planetesimals is the same as that of the hard and soft sphere models. It is divided into three stages: the formation of nonaxisymmetric structures; the creation of planetesimal seeds; and the collisional growth of the planetesimals. The velocity dispersion of particles decreases gradually due to the inelastic collisions. The gravitational instability occurs when the velocity dispersion of particles becomes a critical value determined by Q ≃ 2. The nonaxisymmetric structures form and fragment into planetesimal seeds, the masses of which are on the order of the linear mass. In this paper, we consider the formation of nonaxisymmetric structures and planetesimal seeds as gravitational instability. Planetesimal seeds grow rapidly by coagulation and the large planetesimals form. Most of the mass in the computational domain is absorbed by a few large planetesimals.

We studied the effect of physical parameters on planetesimal formation: the restitution coefficient , the initial optical depth τ, the dependence of planetesimal formation on the ratio ζ = r_H/2r_p, and the initial Toomre's Q value. We found no qualitative differences among the collision models. In the accretion model, the merged particles cannot split; thus, the growth of particles is more efficient than in the hard and soft sphere models. If the restitution coefficient is smaller than the critical value, the time evolution is similar to the standard model. If τ is large, the collision frequency is high, and the dissipation of the kinetic energy is efficient. Thus, the gravitational instability occurs more quickly for larger τ. In this narrow parameter range ζ = 2.5, 2.75, 3.0, there is no remarkable dependence on ζ.

By the long-term and large-scale calculations using the accretion model, we found that the mean mass of the planetesimals depends on the size of the computational domain. We found that the mean mass of planetesimals is proportional to L^3/2_y. However, this mass is independent of L_x if L_x is sufficiently large. The mean mass of planetesimals is estimated by Equation (16). Large planetesimals sweep small planetesimals in the rotational direction. If the orbital separation of planetesimals is sufficiently large, the planetesimals cannot collide. The typical orbital separation is about 5r_H in the present simulation The dependence of the planetesimal mass on the size of the computational domain indicates that we cannot simply discuss the realistic mass of the planetesimal using the local calculation.

The effect of gas was neglected in our simulation in order to study the physical process of gravitational instability of dust particles and subsequent collisional evolution as a first step. Once gravitational instability happens and large aggregates form, we may be able to neglect gas drag because they are so large that they decouple from gas. So our gas-free model may be applicable to the stages of gravitational instability and subsequent collisional growth. The simple gas-free model provides thorough understanding of the gravitational instability and collisional growth of planetesimals. However, it is necessary to include gas drag so as to investigate the realistic formation process of planetesimals, since the stopping time of small dust particles is shorter than the Kepler time. We will investigate the gravitational instability of a particle disk embedded in a laminar gaseous disk in the next paper.

The authors thank the anonymous referee for helpful comments and suggestions. This research was partially supported by MEXT (Ministry of Education, Culture, Sports, Science and Technology), Japan, the grant-in-aid for Scientific Research on Priority Areas, "Development of Extra-Solar Planetary Science," and the Special Coordination Fund for Promoting Science and Technology, "GRAPE-DR Project." S.I. is supported by grants-in-aid (15740118, 16077202, and 18540238) from MEXT of Japan.