MODELING PLANETARY SYSTEM FORMATION WITH N-BODY SIMULATIONS: ROLE OF GAS DISK AND STATISTICS COMPARED TO OBSERVATIONS

Pre-main-sequence stars accrete gas through circumstellar disks. Although the origin of the disk viscosity is still controversial, the onset of MRI helps in the transportation of angular momentum through the disk (Balbus & Hawley 1991). In this paper, we adopt the ad hoc α-prescription (Shakura & Sunyaev 1973) to model the effect of mass transportation in the disk around a classical T-Tauri star. In such a disk, the kinematic viscosity is expressed as ν = αc_sh, where c_s and h are the sound speeds in the midplane disk and density scale height, respectively (Table 1). The surface density at stellar distance a is given as (e.g., Pringle 1981)

$$\begin{equation} \Sigma _g=\frac{\dot{M}}{3\pi \alpha c_s h }\left[1-\left(\frac{R_*}{a}\right)^{1/2}\right], \end{equation} \tag{ 1 }$$

where $\dot{M}$ is the gas accretion rate and R_* is the stellar radius. For a T-Tauri disk with stellar mass of 0.2 M_☉ ⩽ M_* < 2.0 M_☉, the observed accretion rate is (Natta et al. 2006; Vorobyov & Basu 2009)

$$\begin{equation} \dot{M}= 2.5 \times 10^{-8} \,M_\odot \,{\rm yr}^{-1} \left(\frac{M_*}{M_\odot }\right)^{1.3\pm 0.3}. \end{equation} \tag{ 2 }$$

In an optically thin disk, adopting the parameters in Table 1 and letting a ≫ R_*, we get Equation (3) from Equation (1),

$$\begin{equation} \Sigma _g=2400\, {\rm g \, cm^{-2}} \left(\frac{\alpha }{10^{-3}}\right)^{-1} \left(\frac{M_*}{M_\odot } \right)^{0.8} \left(\frac{a}{1\, \rm AU} \right)^{-1}. \end{equation} \tag{ 3 }$$

Table 1. Notations Used in This Paper, $\tilde{M}_*=M_*/M_\odot, \tilde{r}=r/1\,{\rm AU}$.

$T=280 \,\tilde{M}_* \tilde{r}^{-1/2}$	Disk temperature by stellar radiation
$c_s=1.2 \,{\rm km \, s}^{-1}\, \tilde{M}_*^{1/4} \tilde{r}^{-1/4}$	Sound speed in the ideal gas with adiabatic index γ = 1.4
$v_k=29.8 \,{\rm km \, s}^{-1}\, \tilde{M}_*^{1/2} \tilde{r}^{-1/2}$	Orbital speed of the Kepler motion
$h=r c_s/v_k =0.047\, \tilde{r}^{1/4} r$	Vertical scale height of the gas disk
ν = αcsh	Disk viscosity by α model
β = −ln Σg/ln a	Slope of the gas disk surface density

Download table as:  ASCIITypeset image

For an α-disk with constant α = 10⁻³, we derive the minimum mass solar nebular (MMSN), except with a different slope, β = −ln Σ_g/ln a (cf., β = 3/2 in MMSN; Hayashi 1981; Ida & Lin 2004a). We also assume that the protoplanetary disk has a sandwich-layered structure: inside a location a_crit, the protostellar disk is partly thermally ionized, while outside a_crit only the surface layer is ionized by stellar X-rays and diffuse cosmic rays, leaving in the central part of the disk a highly neutral and inactive "dead zone" (Gammie 1996). The viscosity in the MRI's active region could be one or two magnitudes larger than that of the dead zone (Sano et al. 2000). By assuming a constant accretion rate across the disk at a specific epoch in Equation (1), a positive density gradient is expected near a_crit, which helps to halt the embryos under type I migration (Kretke & Lin 2007; Kretke et al. 2009).

To model this effect, we let α_MRI and α_dead denote the α-values of the MRI's active and dead regions, respectively. The effective α for the disk is modeled as (Kretke & Lin 2007)

$$\begin{equation} \alpha _{\rm eff} (a)=\frac{\alpha _{\rm dead}-\alpha _{\rm MRI}}{2} \left[{\rm erf}\left(\frac{a-a_{\rm crit}}{0.1a_{\rm crit}}\right)+1\right]+\alpha _{\rm MRI}, \end{equation} \tag{ 4 }$$

where erf is the error function and 0.1a_crit is the width of the transition region. In this paper, we adopt α_MRI = 0.02, α_dead = 10⁻⁴ (Sano et al. 2000). The small value of α_dead is adopted to obtain a reasonable timescale of type II migration (see Equation (16)). Alternative choices of α_dead are also discussed in Section 4.1.

The location of the inner edge of the MRI's dead zone, a_crit, varies with disk temperature, kinematics, mass accretion rate, etc. Here, we adopt the expression from Kretke et al. (2009):

$$\begin{equation} a_{\rm crit}=0.16 \ {\rm AU}\left(\frac{\dot{M}}{10^{-8}M_{\odot }\,{\rm yr}^{-1}}\right)^{4/9} \left(\frac{M_{*}}{M_{\odot }}\right)^{1/3}\left(\frac{\alpha _{\rm MRI}}{0.02}\right)^{-1/5}. \end{equation} \tag{ 5 }$$

During the evolution of a T-Tauri star and its disk, $\dot{M}$ decreases from ⩽10⁻⁶ to ∼10⁻⁹ M_☉ yr⁻¹, according to the infrared excess observation (e.g., Gullbring et al. 1998). So, the location of a_crit migrates inward as time proceeds. To simplify the procedure, we let the gas disk deplete uniformly in a timescale τ_disk, so $\dot{M}=\dot{M}_0\exp (-t/\tau _{\rm disk})$. Assuming a disk mass of 0.02 M_☉ (Andrews & Williams 2005), we obtain $\dot{M}_0=0.02\,M_{\odot }/\tau _{\rm disk}$.

Substituting Equations (4) and (5) into Equation (3), we obtain the following surface density for the circumstellar disk:

$$\begin{eqnarray} &&\Sigma _{g}=\Sigma _{0} f _g \left(\frac{a}{\rm 1\, AU}\right)^{-1}\left(\frac{\alpha _{\rm eff}}{10^{-4}}\right)^{-1}\left(\frac{M_{*}}{M_{\odot }}\right)^{4/5} \exp \left({{-}\frac{t}{\tau _{\rm disk}}}\right),\nonumber\\ \end{eqnarray} \tag{ 6 }$$

where f_g is the gas enhancement factor. Σ₀ = 280 g cm⁻² is adopted in this work so that the total mass of the disk, up to 100 AU, is 0.02 M_☉ for f_g = 1, corresponding to the average disk mass in the Taurus–Auriga star formation region (Beckwith & Sargent 1996). For such a disk, Tommore's criterion for the onset of gravitational instability is expressed as $Q=\frac{c_s\Omega _K}{\pi G\Sigma _g} =340 (a/1\, \rm AU)^{-3/4}f_g^{-1}$. So Q > 1 holds, and gravitational instability will not occur up to 100 AU as long as f_g ⩽ 10.

We truncate the inner disk at the disk cavity due to the stellar magnetic field. Around the corotation radius, the stellar magnetic torque acts to extract angular momentum from the disk and spins down disk material. At the location where the stellar magnetic field completely dominates internal disk stresses, sub-Keplerian rotation leads to a free fall of the disk material onto the surface of the star in a funnel flow along magnetic field lines, which results in an inner-disk truncation (Königl 1991). The maximum distance of disk truncation is estimated at ∼9 stellar radii. Considering that the radius of the protostar is generally 2–3 times larger than its counterparts in the main-sequence stage, inner disk truncation would occur at <0.1 AU. Thus, we set the inner edge as 0.05 AU in this paper.

According to the core-accretion scenario of planet formation, planetary embryos form through cohesive collisions of heavy elements near the midplane of the circumstellar disk. Within the solid disk, an embryo will grow until it accretes all the dust material around its feeding zone, so that an isolation body is achieved (Kokubo & Ida 1998, 2002). During the later stage of planet formation, inward migration of planetesimals under gas drag and planetary embryos under type I migration may change the distribution of embryos, which results in an uncertainty of the initial embryos' masses. On the other hand, the MMSN model for the solid disk was derived according to the final solid masses of the present planets in the solar system (Hayashi 1981; Ida & Lin 2004a), which reflects the final distribution of the embryo, so we adopt the isolation masses with MMSN distributions, but with an enhancement factor f_d = f_g. In a gas disk, the instability of the isolation masses is suppressed until the gas disk is depleted. Assuming that the dynamical instability timescale of a group of isolated masses is the same as the disk dispersal timescale, the isolation masses are correlated with their mutual separations and are given by (Kokubo & Ida 2002; Zhou et al. 2007)

$$\begin{eqnarray} &&M_{{\rm iso}}=0.16 \,M_{\oplus } \left(f_d \eta _{\rm ice}\right)^{3/2}\left(\frac{a}{1\, \rm AU}\right)^{3/4}\left(\frac{M_*}{M_{\odot }}\right)^{-1/2} \left(\frac{k_{\rm iso}}{10} \right)^{3/2}\!,\nonumber\\ \end{eqnarray} \tag{ 7 }$$

where f_d is the enhancement factor of heavy elements over the MMSN model and η_ice = 1 inside a_ice and 4.2 outside a_ice, with a_ice being the location of the snow line, beyond which water is condensed as ice from disk gas (T ≃ 170 K). Here, k_iso is the separation of embryos scaled by R_H = (2 M_iso/3 M_*)^1/3a. For τ_disk = 1–10 Myr, k_iso = 8–10 at 1 AU and 7–9 at 10 AU for a disk with f_d = 1 (see Figure 11 of Zhou et al. 2007 and details therein).

Due to the variation of the disk accretion rate and stellar radiation at different epochs of the protostar, the location of the snow line varies (Garaud & Lin 2007). For a star with mass ⩽3 M_☉, we adopt the snow line location as (Kennedy & Kenyon 2008)

$$\begin{equation} a_{{\rm ice}}=2.7\, {\rm AU}\left(\frac{M_*}{M_{\odot }}\right)^{4/9}\left(\frac{\dot{M}}{10^{-8}M_{\odot } \,{\rm yr}^{-1}}\right)^{2/9}. \end{equation} \tag{ 8 }$$

To model the type I migration of the embryos, we adopt the expression of the migration timescale from Cresswell & Nelson (2006), which is also valid for eccentric orbits:

$$\begin{eqnarray} \tau _{\rm mig,I}&=&-\frac{1}{C_1}\frac{1}{2.7+1.1\beta } \left(\frac{M_*}{M_p}\right)\left(\frac{M_*}{\Sigma _g a^2}\right) \left(\frac{h}{r}\right)^2\nonumber \\ & &\times\left| \frac{1+\big(\frac{er}{1.3h}\big)^5}{1-\big(\frac{er}{1.1h}\big)^4} \right|\Omega ^{-1}, \end{eqnarray} \tag{ 9 }$$

where a negative (positive) value of τ_{mig, I} corresponds to inward (outward) migration, respectively; and M_p, r, e, and Ω are the mass, position, eccentricity, and angular velocity of the planet, respectively. Due to the MRI effect, β ≡ −∂lnΣ_g/∂lna can be negative near the location of maximum pressure, a_crit. Here, C₁ is a reduction factor. A large amount of literature has shown that to produce the observed planetary occurrence rate, C₁ ∈ [0.03, 0.3] is an appropriate range (e.g., Alibert et al. 2005; Ida & Lin 2008). To test its validity in our model, we set C₁ = 0.03, 0.1, 0.3 and execute some test runs. As shown in Figure 1(a), choosing both C₁ = 0.03 and 0.1 produces a very slow planet migration embedded in our disk model, which may render the formation of hot Jupiters very difficult. So, we set C₁ = 0.3 throughout this paper.

Zoom In Zoom Out Reset image size

The presence of the gas disk will damp the eccentricities of embedded embryos (Goldreich & Tremaine 1980). The e-damping timescale can be described as (Cresswell & Nelson 2006)

$$\begin{eqnarray} &&\tau _{\rm edap,I}=\frac{Q_e}{0.78}\left(\frac{M_*}{M_p}\right) \left(\frac{M_*}{a^2\Sigma _g}\right)\left(\frac{h}{r}\right)^4 \left[1+\frac{1}{4}\bigg(e\frac{r}{h}\bigg)^3\right]\Omega ^{-1},\nonumber\\ \end{eqnarray} \tag{ 10 }$$

where Q_e = 0.1 is a normalization factor to fit with the hydrodynamical simulations.

According to the conventional accretion scenario, giant planets form through three major stages (Perri & Cameron 1974; Mizuno 1980; Pollack et al. 1996): (1) Embryo growth stage: protoplanetary cores form and grow mainly via the bombardment of planetesimals before they attain isolation masses. (2) Quasi-hydrostatic sedimentation stage: the accretion of planetesimals tapers as their supply in the feeding zone is depleted. This induces quasi-hydrostatic sedimentation and the growth of the gaseous envelope due to the loss of entropy. (3) Runaway gas-accretion stage: when the mass of the gas envelope becomes comparable to that of the core, a runaway stage of gas accretion sets in continuously until the gas supply is exhausted by either the formation of a tidally induced gap near the protoplanet's orbit or the depletion of the entire disk.

Through quasi-static evolutionary simulations, Ikoma et al. (2000) derived the critical core mass where significant gas accretion occurs:

$$\begin{equation} M_{\rm crit} \sim 7\, M_\oplus \left(\frac{\dot{M}_{\rm core}}{10^{-7}\, M_\oplus \,{\rm yr}}\right)^{0.25} \left(\frac{\kappa }{1 \, {\rm cm^2\, g}^{-1}}\right)^{0.25}, \end{equation} \tag{ 11 }$$

where $\dot{M}_{\rm core}$ is the rate at which planetesimals are accreted onto the core and κ is the grain opacity (see also Rafikov 2006). Due to the uncertainty of $\dot{M}_{\rm core}$ and κ, we assume $\dot{M}_{\rm core} \propto f_g$ and adopt M_crit = 4 M_⊕f^0.25_g in this paper. A planet core beyond this critical mass will begin to accrete gas in a Kelvin–Helmholtz timescale, τ_KH ≃ 10⁹ yr(M_p/M_⊕)⁻³, so that dM_p/dt ≃ M_p/τ_KH (Ikoma et al. 2000). However, the accretion rate is also limited by the replenishing rate of the materials ($\dot{M}_{\rm disk}$), so the accretion rate of the planetary embryos is expressed as (Ida & Lin 2004)

$$\begin{equation} \frac{d\!M_p}{dt}=\min \left[ 10^{-9}\left(\frac{M_p}{M_{\oplus }}\right) ^4\,M_{\oplus }\,{\rm yr}^{-1},\dot{M}_{\rm disk} \right]. \end{equation} \tag{ 12 }$$

When the planet mass becomes sufficiently large, a tidally induced gap in the gas disk forms around its orbit (Lin & Papaloizou 1979). The critical mass for gap opening is determined by equating the timescale for type I torques to open a gap (in the absence of viscosity) with that for viscous diffusion to fill it. This gives (Armitage & Rice 2005)

$$\begin{equation} \frac{M_p}{M_*}\ge \alpha _{\rm eff}^{1/2}\left(\frac{h}{a}\right)^2. \end{equation} \tag{ 13 }$$

With α_eff ∼ 10⁻⁴ in the MRI's dead zone, the critical mass is given by

$$\begin{equation} M_{I,II}=7.5\left(\frac{a}{\rm 1\, AU}\right)^{1/2}\left(\frac{M_*}{M_{\bigodot }}\right)M_{\oplus }. \end{equation} \tag{ 14 }$$

After the giant planet opens a gap around the disk, it is embedded in the viscous disk and undergoes type II migration. As the mass of the planet grows and becomes comparable to the disk mass, migration slows and eventually stops. The migration speed with the slow-down effect can be expressed as (Alibert et al. 2005)

$$\begin{equation} \frac{da}{dt}=-\frac{3\nu }{2a}\times \min \left(1,\frac{2\Sigma _ga^2}{M_p}\right). \end{equation} \tag{ 15 }$$

So, the migration timescale is given as

$$\begin{eqnarray} \tau _{\rm mig,II}&=&0.6 \, {\rm Myr} \left(\frac{\alpha _{\rm eff}}{10^{-4}}\right)^{-1} \left(\frac{a}{\rm 1 \,AU}\right)\left(\frac{M_*}{M_{\odot }}\right)^{-1/2}\nonumber \\ & &\times \max \left(1,\frac{M_p}{2\Sigma _ga^2}\right). \end{eqnarray} \tag{ 16 }$$

During the migration of the giant planets, their eccentricities will be damped by the disk tide (Goldreich & Tremaine 1979, 1980; Ward 1988), unless the planet is very massive (∼20 M_J; Papaloizou et al. 2001). As the e-damping rate due to the gas disk is quite elusive for different mass regimes of the planets, we adopt an empirical formula for the e-damping timescale (Lee & Peale 2002):

$$\begin{equation} \tau _{\rm edap,II}= \tau _{\rm mig,II}/{K}, \end{equation} \tag{ 17 }$$

where K is a positive constant with values ranging from 10–100. To choose an appropriate value of K, we execute some test runs with K = 10, 30, and 100. As shown in Figure 1(b), the planet eccentricity is damped very quickly with K = 100. For K = 30, the eccentricity can be excited and effectively damped at the end of the simulation, when the gas disk has been nearly depleted. Therefore, we take K = 10 so that some planets in eccentric orbits can survive, in accordance with observations.

The gas accretion of a planet will stop when the gas material around the feeding zone (R_H = h) is exhausted; this gives the truncation mass of a gas giant:

$$\begin{equation} M_{\rm g,iso}(a) = 120\left(\frac{a}{1 \,\rm AU}\right)^{3/4}\left(\frac{M_*}{M_{\odot }}\right)M_{\oplus }. \end{equation} \tag{ 18 }$$

This limits the mass of a giant planet by pure gas accretion. However, collision and cohesive merging may increase the masses of giant planets up to several Jupiter masses.

Based on the model described above, we investigate the late-stage formation of the planetary system in this paper. The stellar mass is taken as 1 M_☉. The gas disk is derived from Equation (6). We further assume that embryos have obtained their isolated masses (see Equation (7)) with f_d = f_g. We ignore embryos with mass ⩽0.1 M_⊕ in the inner and outer orbits, so that there are in total 38–44 isolated embryos initially in circular and coplanar orbits in [0.5–13.5] AU for a system with a gas and a solid disk f_g = f_d = 1. The isolation masses and their radial extensions will be changed accordingly for different f_g = f_d.

The outer boundary of embryos is set according to the standard core growth, beyond which the core masses that can grow within 5 Myr are less than 0.1 M_⊕ (Kokubo & Ida 2002; Ida & Lin 2004a). Although this simplification may neglect their dynamical frictions with inner massive cores, this effect is similar to the damping effect by the gas-disk tide, which is more effective and already included in our simulations.

The angle elements (longitude of periastron and mean motion) of the embryos are randomly chosen. Figure 2 shows an example of the masses and initial locations of the embryos. According to Equation (7), a different τ_disk gives slightly different values of k_iso and M_iso (see Zhou et al. 2007).

Zoom In Zoom Out Reset image size

The acceleration of an embryo with mass M_i(i = 1,...N) is given as

$$\begin{eqnarray} \frac{d}{dt}\textbf {v}_i&=& -\frac{G(M_*+M_i)\textbf {r}_i}{r_i^3} + \sum _{j\ne i}^N G\!M_j \left[\frac{\textbf {r}_j-\textbf {r}_i}{|\textbf {r}_j-\textbf {r}_i|^3 } - \frac{\textbf {r}_j}{r_j^3}\right]\nonumber \\ & &\ -\frac{{\bf v}_i}{2 \tau _{\rm mig}}-\frac{(\textbf {v}_i \cdot \textbf {r}_i)\textbf {r}_i}{r_i^2\tau _{{\rm edap}}}, \end{eqnarray} \tag{ 19 }$$

where $\textbf {r}_i$ and $\textbf {v}_i$ are the position and velocity vectors of M_i relative to the star, respectively; and the third and fourth terms on the right-hand side are the accelerations that cause the eccentricity damping and migration, respectively. For embryos with M_i < M_{I, II} defined in Equation (14), τ_mig = τ_{mig, I} (Equation (9)) and τ_edap = τ_{edap, I} (Equation (10)) are used. For those with M_i > M_{I, II}, τ_mig = τ_{mig, II} (Equation (16)) and τ_edap = τ_{edap, II} (Equation (17)) are used instead. Note that, during the growth of an embryo, it may cross the critical mass M_{I, II}, thus passing from type I to type II migration, and the eccentricity-damping mode is switched.

In the absence of mutual planet perturbations (M_j = 0, j ≠ i), i.e., without the second term on the right-hand side of Equation (19), the secular variations of orbital elements of the embryo M_i under migration and eccentricity damping are derived from the classical perturbation theory:

$$\begin{eqnarray} \left\langle \frac{da}{dt}\right\rangle &=&{-}\frac{a}{\tau _{\rm mig}}-\frac{2e^2}{\tau _{\rm edap}(1+\sqrt{1-e^2})},\nonumber \\ \left\langle \frac{de}{dt}\right\rangle &=&{-}\frac{e(1-e^2)}{\tau _{\rm edap}(1+\sqrt{1-e^2})}, \\ \left\langle\frac{d\omega }{dt}\right\rangle&=&0,\nonumber \end{eqnarray} \tag{ 20 }$$

where ω is the argument of the periastron of the embryo orbit.

Assume that there are N planet cores with masses below M_{I, II} so that they undergo type I migration. Equation (9) gives $\tilde{a}/\dot{\tilde{a}}=\tau _{\rm mig,I}= k \tilde{M}_p^{-1} \tilde{a} \exp (t/\tau _{{\rm disk}})$ for a ≫ a_crit, where $\tilde{M}_p=M_p/M_\oplus$, $\tilde{a}=a/{\rm 1 \,AU,}$ and k ≈ 0.23 Myr if C₁ = 0.3 with the expression of Σ_g from Equation (6). An embryo with an initial semimajor axis $\tilde{a}_0$ will evolve with $\tilde{a}=\tilde{a}_0- \xi \tilde{M}_p$, where ξ = τ_disk[1 − exp (− t/τ_disk)]/k. So, the evolution of the relative separation between two neighboring embryos with the mass difference $\Delta \tilde{M}_p$ is given as

$$\begin{equation} \frac{\Delta \tilde{a}}{\tilde{a}} =\left(\frac{\Delta \tilde{a}_0}{\tilde{a_0}}\right) \frac{1-\xi \Delta \tilde{M_p}/\Delta \tilde{a}_0}{ 1-\xi \tilde{M}_p/\tilde{a}_0}, \end{equation} \tag{ 21 }$$

where $\Delta \tilde{a}_0$ is the initial separation. For a general case, $\Delta \tilde{M}/\Delta \tilde{a}_0 = \gamma \tilde{M}/ \tilde{a}_0$, and since $\xi \tilde{M}_p/\tilde{a}_0 < 1$ before a goes to a_crit, Equation (21) approximates to

$$\begin{equation} \frac{\Delta \tilde{a}}{\tilde{a}} \approx \left(\frac{\Delta \tilde{a}_0}{\tilde{a_0}}\right) [1+(1-\gamma) \xi \tilde{M}/\tilde{a}_0 ]. \end{equation} \tag{ 22 }$$

So, as long as γ < 1, the inward type I migration will increase the mutual distance scaled by its mutual Hill radii (R_H∝a). According to Zhou et al. (2007), such an increase will enhance the orbital crossing timescale of the system, for example, in the case of isolation masses in Equation (7) with γ = 3/4. However, due to the perturbation of giant planets and some other embryos, embryos in inner orbits will have their eccentricities excited, which may result in instability.

Similar to the previous analysis, giant planets undergoing type II migration with a timescale from Equation (16) may modify their mutual separations. As in the previous subsection, we assume that the planet mass is much smaller than the inner disk mass, and the disk has a constant α_eff. Thus, one has $\tilde{a}/\dot{\tilde{a}}=\tau _{\rm mig,II}= k^{\prime } \tilde{a}$ for a ≫ a_crit, where k' ≈ 0.6 Myr from Equation (16). So, $\tilde{a}=\tilde{a}_0-t/k^{\prime }$, and two neighboring planets will have constant separation $\Delta \tilde{a}_0$. When it is scaled by R_H∝a at time t,

$$\begin{equation} \frac{\Delta \tilde{a}}{\tilde{a}} =\left(\frac{\Delta \tilde{a}_0}{\tilde{a_0}}\right) \left(1-\frac{t}{k^{\prime }\tilde{a}_0 } \right)^{-1},\quad (t<k^{\prime }\tilde{a}_0). \end{equation} \tag{ 23 }$$

This means that the mutual separation of the embryos scaled by Hill radii will increase during inward type II migration. However, secular perturbations among them will excite their eccentricities during their convergent migration, thus destabilizing the system if their eccentricities are high enough.

The variations in disk mass and viscosity are modeled by two parameters, i.e., f_g in Equation (6) and α_dead in Equation (4). We set f_g = 0.3, 1, 3, and 10 and α_dead = 10⁻², 10⁻³, and 10⁻⁴ to study their individual contribution to the architecture of planetary systems. Five runs for each parameter are executed. We fix the surface density profile as well as the disk depletion timescale (τ_disk = 2 Myr) for these test runs.

In our model, α_dead mainly works on type I migration via Equation (9) and on type II migration via Equation (16). Since Σ_g∝(α_eff)⁻¹ (Equation (6)), τ_{mig, I}∝α_eff, i.e., planets in disks with smaller α_eff(∼ 10⁻⁴) migrate faster, which results in a smaller average semimajor axis (Figure 3(a)). For type II migration, if the planet mass is small (M_p ⩽ 2Σ_ga², which is about 90 M_⊕ at 5 AU and 170 M_⊕ at 10 AU), the migration speed decreases with α_eff, and thus larger α_eff (∼10⁻²) may also lead to a small average semimajor axis. Thus, disk viscosity may change the final location of planets. So, we take the most plausible value (α_dead = 10⁻⁴) in the following simulations (Sano et al. 2000). On the other hand, the average eccentricity is around 0.1, with large variations for all cases, which seems to be insensitive to α_dead. The reason is that small α_dead leads to a high density of the gas disk, resulting in two competing effects: (1) a fast type I migration so that a large frequency of embryo-scattering is expected, which enhances the eccentricities of embryos; and (2) a fast e-damping rate with a small e_mean. As the timescales of these two effects are comparable, e_mean is nearly independent of α_dead.

Zoom In Zoom Out Reset image size

There are three ways for f_g to influence orbital evolutions: the surface density of the gas disk (Equation (6)), initial isolation masses (Equation (7) with f_d = f_g), and the critical mass for the onset of runaway gas accretion (Equation (11)). Larger f_g may result in larger cores extended to an outer region. As we assume $\dot{M}_{\rm core} \propto f_d$, we have M_crit∝f^0.25_g. Compared to the masses of the initial embryos (∝f^3/2_g), it is much easier to form massive gas giants in a massive gas disk.

As shown in Figure 3(b), when f_g ⩾ 1, a_mean values are all around 3–4 AU. Although the migrations should be faster for planets in disks with larger f_g, the initial embryos extended to an outer region in these cases, which resulted in similar average locations. For a less massive disk (f_g = 0.3), these small embryos, initially in the inner region, never accreted gas, and thus they experienced type I migration throughout the evolution; therefore, they have smaller a_mean and e_mean. More effective eccentricity damping results in a smaller e_mean for f_g = 3 than for f_g = 1. When f_g = 10, the largest embryos can reach 80 M_⊕ via collisions in 1 Myr. Finally, planets left in these systems are mainly gas giants, and their eccentricities are damped much less effectively than those of small planets. This is why the e_mean of f_g = 10 is much larger than that of f_g = 3.

The role of the gas disk on planetary evolution is manifold: it causes the inward migration (type I or II) of protoplanets (Ward 1997; Lin & Papaloizou 1985) as well as the tidal damping of their orbital eccentricities. A long-survival gas disk may be helpful in the formation of giant planets, while planets in a short-lived disk may not have enough time to accrete gas, thus remaining either terrestrial or Neptunian planets. To investigate the effect of the disk depletion timescales (τ_disk), we set 11 values of τ_disk, evenly ranging from 0.5 to 5 Myr in a logarithmic scale (Haisch et al. 2001). For each τ_disk, we performed 20 runs of simulations by randomly choosing the orbital phase angles of the embryos. So, in total, we performed 220 simulations for f_g = 1. All simulations were stopped at t = 10 Myr.

An interesting problem of planet formation is the growth epoch of different types of planets. We define (somewhat arbitrarily) the following two types of planets.

Gas giant planets (GPs), including massive GPs (M_p ⩾ 30 M_⊕) and Neptune-sized GPs (10 M_⊕ ⩽ M_p < 30 M_⊕); and terrestrial planets (TPs), including super-Earth TPs (1 M_⊕ ⩽ M_p < 10 M_⊕) and sub-Earth TPs (M_p < 1 M_⊕). Figure 4 shows the distribution of planets' semimajor axes for the 220 runs of simulations at some epochs with f_g = 1. Only TPs are present at t = 10⁴ yr (Figure 4(a)). When t = 1 Myr, runaway gas accretion occurs; thus, a few GPs appear (Figure 4(b)). The most efficient growth of GPs occurs at 1–3 Myr; thus, the number of GPs increases very quickly (Figure 4(c)), until they become dominant members at the end of the simulations (Figure 4(d)).

Zoom In Zoom Out Reset image size

Figure 5 plots the correlations between the properties of the final systems and τ_disk. There are two regimes. For short-lived disks with τ_disk ⩽ 1 Myr, the forming planets have small average masses (70–200 M_⊕; Figure 5(c)), and the average semimajor axes of these systems are large (5–6 AU; Figure 5(d)), indicating that most of the surviving planets were formed in distant orbits, while migrations did not strongly affect their orbital architectures. The average eccentricities (∼0.15) are relatively large (Figure 5(a)), showing that the tidal damping effect is not very effective due to the short lifetime of the disk. However, planetary systems with a longer disk lifetime (τ_disk > 1 Myr) tend to have larger average planet masses (∼1 M_J), with lower average eccentricities (∼0.1) due to the long period of disk damping. The semimajor axes of these systems decrease from 4 to 2 AU, indicating that inward migration indeed plays an important role in sculpting their orbital configurations.

Zoom In Zoom Out Reset image size

Systems with τ_disk < 1.5 Myr do not have many differences from the surviving planetary systems, while N_mean decreases slightly as τ_disk surpasses 1.5 Myr (Figure 5(b)). The maximum of N_mean at ∼1.5 Myr is due to the mechanism for halting small planets near a_crit (see Section 5.1). In disks with smaller τ_disk, the effects due to gas disks are not so dominant, and the surviving number of planetary systems mainly depends on the interactions between embryos. In disks with a longer lifetime, giant planets experience sufficient type II migration and will be closer to the boundary of the MRI's region, where small planets were halted (see Section 5.1). Thus, small planets are more easily scattered out of the system or hit the host star, leaving fewer surviving planets.

One of the major effects that N-body simulations can describe is the eccentricity excitation due to mutual planetary perturbations. Such excitation may lead to some embryos being scattered out of the system, while others may be merged. Thus, the final number of planets should be greatly reduced from that of the initial embryos. Correlations between the number of surviving planets and the average mass, as well as the eccentricities of the planets in a planetary system, are shown in Figure 6. They are fitted by

$$\begin{eqnarray} e_{\rm mean}&=&0.65\times 0.67^{N_{\rm left}},\nonumber \\ M_{\rm mean}&=&1.28\, M_J \times 0.85^{N_{\rm left}}. \end{eqnarray} \tag{ 24 }$$

These correlations show that, in multi-planet systems like our solar system, planets have relatively lower average masses and eccentricities than those found in systems with one or only a few planets. Qualitatively, this law is easily understood. To achieve a longer orbital crossing time in a multi-planet system, either the eccentricities of the planets must be small or the planets must have large mutual separation scaled by their Hill radii; thus, smaller planetary masses are helpful in achieving a stable system (Zhou et al. 2007).

Zoom In Zoom Out Reset image size

Figures 7(a) and (b) show the semimajor axis distribution from both the observations and simulations. The observed distribution shows two peaks (Figure 7(a)). (1) At 1–3 AU: this is roughly the snow line (where water is frozen) of the system, where planet cores may be stalled under type I migration (Kretke & Lin 2007; Ida & Lin 2008), and subsequent accretion of gas makes them giant planets. (2) At 0.04–0.06 AU (or 3–5 days): this is roughly the inner edge of the gas disk, where planets under type II migration will be stopped (Lin et al. 1996). From our simulations, there are a large number of planets beyond 3 AU; thus, the lack of planets at >3 AU is due to observational bias.

Zoom In Zoom Out Reset image size

In addition, we show that there is an extra pileup of planets: (3) at around 0.2 AU (or around 30 days), which has not been revealed yet by observations. This location corresponds to the inner edge of the MRI's dead zone (a_crit), where small planets under type I migration may be halted. However, as the location of a_crit (see Equation (8)) moves inward, in the case of a short τ_disk, its migration speed is faster than that of the type I migration; therefore, it is difficult to halt these planets near a_crit. When τ_disk ∼ 1.5 Myr, $\frac{da}{dt}\mid _{\rm mig,I} \ \sim \ \frac{da_{\rm crit}}{dt}$, so if τ_disk > 1.5 Myr, the mechanism for halting small planets near a_crit works. After τ_disk, the migration rate is reduced due to gas depletion. In fact, a_crit ∼ 0.2 AU at t = 2τ_disk. This explains the accumulation of planets near 0.2 AU. Some small planets were scattered into 0.04 AU.

Our simulations also show that the majority (70%) of giant planets are located in 1–10 AU. Most of them are massive GPs and experienced type II migration. Due to the depletion of the gas disk, they could not migrate to the proximity of the star; hence, they halted outside 1 AU. Due to observational bias of radial velocity measurements, massive giant planets are revealed only at <4–5 AU.

Our simulations show that there are mainly two peaks for planet masses (Figure 7(c)).

• 1.  
  At 1–2 M_J. These are Jupiter-sized giant planets that have grown up with efficient gas accretion. Although our simulations reproduced this peak, there are not enough massive planets at M_p ∼ 5–8 M_J.
• 2.  
  At 0.1–3 M_⊕. There are a large number of terrestrial planets, which have grown up mainly by mutual collisions from the perturbations of outside giant planets. This peak has not been revealed by observations yet, and can be checked by future observations of terrestrial planets.One additional small peak is revealed by our simulations:
• 3.  
  Around 10 M_J. Massive embryos (M ∼ 80 M_⊕) formed after 1 Myr in disks with f_g = 10 and only experienced time-limited type II migration. More massive giant planets (>10 M_J) may be born through some other mechanisms, e.g., the gravitational instability scenario.

The desert at 10–20 M_⊕ consists of Neptune-sized planets. According to our model in Section 2.3, the Kelvin–Helmholtz timescales are quite short (∼10⁴ yr) for this mass regime, and runaway gas accretion quickly transforms them into giants. Only a few planets in the inner region with M_iso ∼ 10–20 M_⊕ survive in this desert. Also, merging from lower planet embryos may form Neptune-sized planets.

The distribution of eccentricities from simulations is similar to that from observations. In our simulations, the eccentricities of surviving planets vary from 0 to 0.84, while the observations show a maximum eccentricity of about 0.9 (Wright et al. 2009). We fit both distributions (simulations and observations) with an exponential decay function in the form of P(e)de = N(e)/N_tot ∼ exp (− Ae)de. Since ∫¹₀P(e)de = 1, we have

$$\begin{equation} P(e) =\frac{A \exp (-A e)}{1-\exp (-A)}, \end{equation} \tag{ 25 }$$

where A = 4.2 for observations and A = 7.8 for simulations of planets with V_r < 3 m s⁻¹, as shown in Figures 7(e) and (f). A larger value of A from simulations indicates a steeper slope, indicating that more planets with small eccentricities are still not detected, as shown in Figure 8(d).

Zoom In Zoom Out Reset image size

Figure 8 presents the correlation diagrams of a, e, and M_p from our simulations (right), with a comparison to the observational data (left). Our simulations reproduced all three correlation plots quite well. In particular, there is a gap in the M_p–e plot, indicating that a planet desert (0.005–0.08 M_J) depends on the eccentricity (Figures 8(c) and (d)). Figure 8(d) also shows that giant planets (M_p > 10 M_⊕) tend to have larger eccentricities on average. To show this more clearly, we reproduce the eccentricity–semimajor axis correlation plots and the eccentricity distribution plots for giant planets (M_p > 10 M_⊕) and terrestrial planets (M_p < 10 M_⊕) in Figure 9. Giant planets have average eccentricities ∼0.15 at all locations (except >10 AU), while terrestrial planets have average eccentricities ∼0.05. Although both of the e distributions can be fitted by exponential law in Equation (25), their coefficient A values are different: 7.78 for M_p > 10 M_⊕ and 21.0 for M_p < 10 M_⊕, indicating that small-mass planets tend to have smaller eccentricities.

Zoom In Zoom Out Reset image size

To understand the M_p–e desert, i.e., very few planets with masses 0.005–0.08 M_J have eccentricities larger than 0.3–0.4, we plot the e damping timescales τ_edap as a function of planet masses at 0.3, 1, 3, and 10 AU in Figure 10. For a planet, when M_p < M_{I, II}, τ_edap is calculated by Equation (10) using a mean eccentricity, e = 0.05 and Σ₀ = 280 g cm⁻². Otherwise, τ_edap is calculated by Equation (17), with viscosity α_dead = 10⁻⁴. According to Equation (10), small planets have longer τ_edap. The jumps of τ_edap occur at M_p = M_{I, II}. Note that τ_edap stays horizontal until the mass of a planet becomes comparable to the disk mass (M_p > 2Σ_pa²_p = 21 M_⊕(a/1 AU)). In the braking phase of type II migration, more massive planets have a longer e-damping timescale according to Equation (17). Assuming effective damping of eccentricity if τ_edap is less than ∼1 Myr, we obtain a mass regime 0.005–0.08 M_J, which is consistent with the M_p–e desert from both observations and our simulations in Figures 8(c) and (d). As shown in Figure 10, τ_edap values of massive GPs are in the horizontal region or beyond, so generally they have longer e-damping timescales than those of the TPs. Therefore, the mean eccentricity of TPs is smaller than that of GPs.

Zoom In Zoom Out Reset image size

During the entire evolution period we simulated (10 Myr), there were 606 collisions in the 20 runs of 2D simulations, while 471 collisions occurred in the corresponding runs of the 3D model. In the 2D runs, the number of embryos colliding with the host star or being ejected out of the system was 120, compared to 231 in 3D runs. This is consistent with the result of Chambers (2001). As seen in Figure 11(a), the cumulative distribution functions (CDFs) of collisions in two models show that most collisions in 3D runs occurred in relatively later stages. Thus, earlier collisions in the 2D model produced embryos with larger masses. To show the effects of earlier collisions on the final planetary architectures, we divided the collisions into three regions, as shown in Figure 11(b): regions I (M₀ < M_crit, M₁ < M_crit), II (M₀ < M_crit, M₁ > M_crit), and III (M₀ > M_crit, M₁ > M_crit), where M₀ represents the larger mass of the embryo before a two-body collision, M₁ represents the mass after the collision, and M_crit = 4 M_⊕ is the critical mass beyond which efficient gas accretion begins (see discussion below Equation (11)).

Zoom In Zoom Out Reset image size

Region I: Most collisions (2D: ∼75%; 3D: ∼85%) occurred in this region. These collisions influence only the Earth-like planets (TPs; see Section 4.2). Due to the shorter collision timescale in the 2D model, embryos have larger masses, on average, during earlier times and experience faster type I migrations and e damping, according to Equations (9) and (10). Therefore, TPs have a smaller mean eccentricity (cf., e_mean = 0.034 for 2D and e_mean = 0.047 for 3D) and semimajor axis (see Figure 11(c)). As a result of more collisions occurring in 2D simulations, the final systems are expected to have fewer and larger masses of TPs than those from 3D simulations. In fact, we find 11 TPs left among the 74 surviving planets in the 20 runs of 2D simulations, while 24 of the 98 surviving planets in the corresponding 3D simulations have a small average mass (Figure 11(c)).

Region II: About 8.5% (2D) and 8.4% (3D) of collisions occurred in this region. These collisions make the embryos with masses <M_crit grow beyond critical mass. However, this effect is limited by the small fraction of collisions and the slow gas accretion rate (Equation (12)) for embryos with masses slightly larger than M_crit = 4 M_⊕. Thus, the differences between 2D and 3D models due to collisions in this region can be ignored.

Region III: Roughly ∼16.5% (2D) and ∼6.6% (3D) of collisions occurred in this region. These embryos can accrete gas before collisions. As the speed of type I migrations is much faster than that of type II, the final locations of planets are determined mainly by their cores. As planet cores in the 2D model undergo more collisions, they have bigger masses and a fast migration speed, so their final average semimajor axis is smaller than that in the 3D model (Figure 11(d)). As the inner region has less gas to accrete, M_gas∝Σ_gaR_H, while R_H∝(M_p/3M_*)^1/3a and Σ∝a⁻¹, the final giant planets in 2D simulations have less mass (Figure 11(d)). The type I eccentricity damping rate is much smaller than that of type II (Figure 10), and from Equation (10), τ_{edap, I} ∼ a^1/2M⁻¹_p, so larger masses of planets in 3D simulations have shorter τ_{edap, I}; thus, the mean eccentricity (e_mean ≈ 0.086 in simulations) in the 3D model is smaller than that (e_mean ≈ 0.15) in the 2D model.

As most of our results in Section 5 are in statistics, we focus on the statistical differences between 2D and 3D models. The statistics of semimajor axes and masses of planets in the two models are presented in Figure 12. Some peaks and deserts in the 2D model are reproduced in the 3D model, such as the peaks at a_crit ∼ 2 AU, M_crit, and 1–2 M_J. However, some different characteristics still exist, especially the planet deserts at around 0.1 AU and 10 M_⊕ < M_p < 0.1 M_J in the 3D model. The lack of planets around 0.1 AU in the 3D model may be due to the prolonged collision timescale, and a much longer time of integration may give results qualitatively similar to the 2D model. The deficit of planets with masses 10 M_⊕ < M_p < 0.1 M_J in the 3D model makes the planet desert in Section 5.2 more obvious.

Zoom In Zoom Out Reset image size

Figure 13 shows the orbital inclination distribution for the surviving planets in the 20 runs of 3D simulations. As one can see, most of the planets remain at <2°; only a few planets with smaller masses have inclinations ∼10°.

Zoom In Zoom Out Reset image size

To summarize, we point out that the differences between 2D and 3D simulations are small, and most of the statistical results in the 2D model are qualitatively credible.