UNSTABLE PLANETARY SYSTEMS EMERGING OUT OF GAS DISKS

Out of over 360 planetary systems discovered so far, about 12% are known to be multiplanet systems (http://exoplanet.eu/). Also, recent observations have started revealing that many of the previously detected single planets are accompanied by an additional planet on a further orbit (e.g., Wittenmyer et al. 2007; Wright et al. 2007). It will become increasingly more important to understand the formation and evolution of multiplanet systems, which can explain the observed properties of extrasolar planetary systems.

Recent numerical N-body simulations of planetary systems without a gas disk demonstrated that dynamical instabilities occurring in the multiplanet systems, which are characterized by orbital crossings, collisions, and ejections of planets, could increase planetary eccentricities (e) efficiently (Rasio & Ford 1996; Weidenschilling & Marzari 1996). These studies successfully reproduced the observed eccentricity distribution of extrasolar planets (Ford & Rasio 2008; Chatterjee et al. 2008, hereafter C08; Jurić & Tremaine 2008, hereafter JT08).

Such N-body simulations also suggest that the planet–planet interactions alone cannot explain small semimajor axes (a) of the observed planets, if giant planets are formed beyond the ice line as expected from the standard planet formation theory. More specifically, starting with giant planet systems beyond 3 AU, C08 found that it is difficult to scatter planets within ∼1 AU. This is because planet–planet scattering conserves energy and therefore is not particularly efficient in shrinking the planetary orbits.

The disk–planet interactions, on the other hand, are known to decrease semimajor axes of planets efficiently (Ward 1997). The overall effect of such interactions on the orbital eccentricity is highly uncertain, and depends on the detailed disk structure, as well as planetary masses. Generally, disk–planet interactions lead to eccentricity damping, but for planets massive enough to open a clean gap in the disk, eccentricity can increase rapidly depending on the level of saturation of corotation resonances (Goldreich & Sari 2003; Moorhead & Adams 2008). However, hydrodynamic simulations show that the disk–planet interactions typically lead to e < 0.2 (e.g., D'Angelo et al. 2006).

The time to dynamical instability scales with the distance between planets (e.g., C08). Thus, all the N-body studies of planet–planet scattering assume initially dynamically active planetary systems (i.e., distance between planets being less than a few Hill radii). However, with the aid of a gas disk, planets which would not easily reach dynamical instability could experience strong interactions. For example, when planets are embedded in the inner cavity of a disk, the surrounding disk could drive the outer planet closer to the inner ones, triggering dynamical instability (e.g., Adams & Laughlin 2003; Moorhead & Adams 2005). The orbital eccentricity can also be increased if planets are trapped in mean-motion resonances (MMRs) during such a convergent migration (e.g., Snellgrove et al. 2001; Lee & Peale 2002). Alternatively, when a gas disk annulus is left between planets, the eccentricities of planets could still increase by repeated resonance crossings due to divergent migration (Chiang et al. 2002). Thus, whether the combined effects of disk–planet and planet–planet interactions lead to eccentricity excitation, or damping, should be studied carefully. One of the goals of our study is to verify the initial assumptions of N-body studies (i.e., planets are on nearly circular and coplanar orbits while the gas disk is still present), and figure out whether the eccentricity distribution is largely determined before, or after, the disk dissipation.

In this paper, we numerically study the evolution of three-planet systems in a dissipating gas disk, and constrain the "initial" conditions of planetary systems which can reproduce the a and e distributions simultaneously. We calculate the disk–planet interactions directly so that the disk and planetary orbits evolve self-consistently. Also, we take account of the effect of saturation of corotation resonances on eccentricity damping. We introduce the numerical methods in Section 2, and the initial conditions in Section 3. In Section 4, we show that the observed a–e scatter plot can be reproduced well for a reasonable range of disk masses. We also discuss the mass distribution and MMRs for representative cases. Finally, in Section 5, we compare this work with some recent observations, and summarize our results.

To simulate multiplanet systems in gas disks, we use a hybrid code which combines the symplectic N-body integrator SyMBA (Duncan et al. 1998) with a one-dimensional gas disk evolution code (Thommes 2005). SyMBA utilizes a variant of the so-called mixed-variable symplectic (MVS) method (Wisdom & Holman 1991), which treats the interaction between planets as a perturbation to the Keplerian motion around the central star, and handles close encounters between bodies with a force-splitting scheme which has the properties of an adaptive time step (Duncan et al. 1998). When the bodies are well separated, SyMBA has the speed of the MVS method, while during the close encounters, the time step for the relevant bodies is recursively subdivided.

On the other hand, the gaseous disk evolves both viscously and via gravitational interaction with planets, according to a general Navier–Stokes equation. Following the standard prescription by Lin & Papaloizou (1986), the gas disk is divided into radial bins, which represent disk annuli with azimuthally and vertically averaged disk properties like surface mass density, temperature, and viscosity. Viscous evolution of the disk is calculated by using the standard alpha viscosity prescription (Shakura & Sunyaev 1973), while the disk–planet interactions modify the disk evolution via the torque density formulated in Ward (1997; see also Menou & Goodman 2004). The calculated torque density is used in turn to determine the migration rates of planets.

In our simulations, a disk extends from 0.02 to 100 AU. The time step of simulations is typically 0.05 yr. We remove a planet from the simulation when the heliocentric distance of a planet becomes larger than 1000 AU, or smaller than 0.02 AU. These cases are referred to as "ejection", or "collision", respectively, throughout the paper. We assume that planets collide when the separation between two planets become less than the sum of the planetary radii. Planets are treated as sticky spheres, so that the collision conserves total mass and momentum. We refer to a collision between planets as "merger" in this paper.

2.1. Gas Accretion onto a Planet

The above code follows the evolution of planetary systems as planets gravitationally interact with each other, migrate, and open gaps in the disk. However, planets are also expected to clear the gas annuli between them as they grow by accreting gas from the surrounding disk.

Once a planet becomes massive enough to open a gap in the disk, the gas accretion rate is controlled by how quickly the disk can supply gas to the planet, rather than how quickly the planet can accrete gas (Bryden et al. 2000; Tanigawa & Watanabe 2002). Since all the planets in our simulations are more massive than Neptune, they are expected to have circumplanetary disks, from which they accrete. Although our code does not resolve such disks, we can mimic the accretion effect by adopting the results of hydro simulations. Tanigawa & Watanabe (2002) showed that, without a gap-opening effect, a planet accretes gas within a few Hill radii on the following timescale:

$$\begin{equation} \tau _{\rm subdisk}=30.76\; {\rm yr} \left(\frac{M_p}{M_J}\right)^{-1/3} \left(\frac{\Sigma (r)}{1.7\times 10^3\,{\rm g\,cm^{-2}}}\right)^{-1}. \end{equation} \tag{ 1 }$$

Here, M is the mass and Σ is the surface mass density. The subscripts p and J represent the planet and Jupiter, respectively. Due to the gap-opening effect, the true accretion timescale is likely longer than the above estimate (e.g., D'Angelo et al. 2003). However, since our code handles gap opening due to disk–planet interactions separately, we adopt the above accretion timescale for all of our planets.

2.2. Eccentricity Evolution

For simplicity, we focus on eccentricity evolution due to the first-order resonances. Also, to evaluate the effects of eccentricity excitation due to planet–planet interactions, we neglect the e excitation due to disk–planet interactions.

When a planet is too small to open clean gaps, the eccentricity evolution is dominated by co-orbital Lindblad resonances which damp eccentricity (e.g., Artymowicz 1993). In such a case, the eccentricity damping timescale can be expressed as follows (e.g., Kominami & Ida 2004):

$$\begin{equation} \tau _{\rm{edamp}}=-\frac{e_{p}}{\dot{e_{p}}}=\left(\frac{h_p}{r_p}\right)^4 \left(\frac{M_{*}^{2}}{M_{p}\bar{\Sigma _p} r_p^{2}}\right)\Omega (r_p)^{-1}. \end{equation} \tag{ 2 }$$

Here, h is the pressure scale height, r is the orbital radius, and Ω is the orbital frequency. The subscript * denotes the star. Note that the equation is evaluated at the location of a planet. Also, for the disk's surface mass density, we take the average of the pericenter, semimajor axis, and apocenter of a planetary orbit.

On the other hand, when a planet becomes large enough to open a clean gap, the effect of co-orbital Lindblad resonances diminishes, and the competing effects of e excitation due to external Lindblad resonances and e damping due to corotation resonances become significant (Goldreich & Tremaine 1980). In such a case, the eccentricity damping can be written as follows (Goldreich & Sari 2003):

$$\begin{equation} \tau _{\rm edamp}=-\frac{e_{p}}{\dot{e_{p}}}= \frac{1}{K_{e}}\left(\frac{w}{r}\right)^4 \left(\frac{M_{*}^{2}}{M_{p}\Sigma r^{2}}\right)\Omega ^{-1}, \end{equation} \tag{ 3 }$$

where w = r(3πα)^−1/3(r/h)^2/3(M_p/M_*)^2/3 is the gap width, which is determined by balancing the gap-opening tidal torque from the principal Lindblad resonances with the gap-closing viscous torque, and α is the standard viscosity parameter by Shakura & Sunyaev (1973). The above equation is evaluated at the gap edges, where the contribution from the resonances is the largest. The damping efficiency is governed by K_e, which is defined as below by following the approach of Goldreich & Sari (2003):

$$\begin{equation} K_{e}=\left[1.046\, F(p) - 1 \right]. \end{equation} \tag{ 4 }$$

In this equation, F(p) is the saturation function of corotation torques that is numerically evaluated by Ogilvie & Lubow (2002) and interpolated by Goldreich & Sari (2003) as

$$\begin{eqnarray} F(p) & \simeq & \frac{(1+0.65\, p^{3})^{5/6}}{(1+1.022\, p^{2})^{2}}, \end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray} p & \sim & \left(\frac{r}{h}\right)^{2/9}\left(\frac{M_{*}}{M_{p}}\right)^{1/9}\frac{e_{p}}{\alpha ^{1/9}}. \end{eqnarray} \tag{ 6 }$$

For K_e = 0.046 (and hence F(p) = 1), the corotation torques are unsaturated, and fully contribute to the eccentricity damping, while for the other extreme K_e = 0, the effects of corotation torques are negligible, and there is no damping. Negative values of K_e corresponds to e excitation, but we do not take the effect into account. Note, however, that hydrodynamic simulations show that the disk–planet interactions typically lead to e < 0.2 (e.g., D'Angelo et al. 2006), and therefore may not be able to explain planets with high eccentricities.

2.3. Disk Dissipation Timescale

We study evolution of three-planet systems in a gas disk by changing various parameters. Specifically, we study five different disk masses with 200 planetary systems, by changing the gas dissipation time τ_GD between 2 and 4 Myr. We run each set of runs with and without the effect of saturation of corotation resonances. Our initial conditions are summarized in Table 1.

We focus on three-planet systems, and define their initial properties following C08. We determine the planetary masses by adopting a simplified core accretion scenario. Specifically, we sample planetary core masses randomly from a distribution over 1–100 M_E (where M_E is the Earth mass) uniform in M^1/5_core, and assume that each core accretes all gas within the "feeding zone" that is extending over Δ = 8R_Hill, core, and centered on the core:

$$M_{p}=2\pi a\Delta \Sigma +M_{\rm core}.$$

Here, the core's Hill radius is defined as R_Hill, core = (M_core/(3M_*))^1/3a. The size of the feeding zones is a typical distance between planetary embryos (Kokubo & Ida 2002).

As in C08, the semimajor axes are chosen so that the distance between planets is scaled with K = 4.4 times the Hill radius of the ith planet:

$$a_{i+1}-a_{i}=K\, R_{{\rm Hill},i},$$

with a₁ = 3 AU. We fix the initial semimajor axis of the innermost planet following the common assumption that giant planets form beyond the "ice line," where the solid density is higher due to the condensation of icy and/or carbonaceous material (Lewis 1974; Lodders 2004). From this prescription, we obtain planets with mass ranging over 0.4–4 M_J, between 3 and 6.5 AU.

We make two independent sets of 100 three-planet systems. Since typical planetary eccentricities and inclinations just after planet formation are not known, we use two different ranges of initial e and i. For one of the sets (marked with d1 in Table 1), we choose initially moderate e and i, which are randomly drawn from a uniform distribution in the range of e = 0–0.1, and from a distribution uniform in cos i over the range of i = 0–10 deg, respectively. For the second set (marked with d2 in Table 1), on the other hand, we choose initially small e and i over the range of e = 0–0.05, and i = 0–0.03 deg, respectively. The former set is identical to Mass distribution 3 in C08. In both sets, phase angles are randomly chosen from a uniform distribution over 0–360 deg. Also, the gas disk dissipation time τ_GD is chosen randomly between 2 and 4 Myr for each three-planet system.

Since our planets are nearly fully grown, we use evolved gas disks as the initial gas disks, instead of using the primordial ones. Each initial disk is generated by evolving the minimum mass solar nebula (MMSN)-type disk with Σ = 10³(a/AU)^−3/2 g cm⁻² for 5, 6, 7, 8, and 9 Myr without planets, under the disk's viscosity α = 5 × 10⁻³. This corresponds to five different initial disk masses in the range of ∼0.3–1.6 M_J. All disks are extending from 0.02 to 100 AU.

These initial conditions are summarized in Table 1.

In Sections 4.1 and 4.2, we focus on the results of sets with initially moderate e and i (i.e., sets denoted by d1, see Table 1). Without a gas disk, C08 showed that N-body simulations of this type of three-planet systems reproduce the observed e distribution of exoplanets very well. On the other hand, by assuming that planets are formed beyond the ice-line (∼3 AU), their simulations showed that it is difficult to scatter planets within ∼1 AU. In Section 4.1, we discuss the effects of the initial disk mass on the a–e distribution, while in Section 4.2, we study the effects of the saturation of corotation resonances. In Section 4.3, we investigate the corresponding results for initially small e and i. We discuss the evolution for typical cases in Section 4.4, and the mass distribution in Section 4.5. Finally, in Section 4.6, we study the MMRs seen in our simulations.

4.1. Effect of the Initial Disk Mass

First, we investigate the effects of initial disk masses on the final distribution of a and e by neglecting the effects of saturation of corotation resonances (i.e., t5d1–t9d1). The a–e scatter plots of planets after 100 Myr for t6d1–t9d1 are plotted in Figure 1. Also plotted is the observed a–e distribution. As expected, more massive disks have more planets with smaller semimajor axes, due to efficient planet migration.

Zoom In Zoom Out Reset image size

The overall trend of the a–e scatter plot is reproduced fairly well, especially for t6d1 and t7d1. We can quantify this by using the Kolmogorov–Smirnov (K–S) test. We perform the K–S tests against the null hypothesis that two arbitrary distributions are drawn from the same underlying distribution, and quote the significance level probabilities in Table 3. We choose to reject the null hypothesis for P < 0.1. To compare the observed and simulated distributions, we select planets between 0.2 and 6 AU with mass ranging over 0.3–4 M_J. The lower limit of the semimajor axis corresponds to the orbital radius below which the tidal interaction with the central star becomes significant, while the upper limit is motivated by the current maximum orbital radius of a planet detectable by radial velocity observations. Our two dimensional K–S tests for the observed and simulated distributions indeed show that we cannot reject the null hypothesis for t6d1 and t7d1. The agreement with the observed distribution becomes worse for more, or less massive disk cases (t5d1, t8d1, and t9d1). Thus, our results indicate that there may be an optimum disk mass to reproduce the observed a–e distribution.

Now we look into the a and e distributions separately. Figure 2 compares the observed and simulated a and e distributions at 100 Myr for t6d1, t7d1, and t8d1. For the a distribution, the K–S test shows that we cannot reject the null hypothesis for t6d1 and t7d1 at 100 Myr, while for t8d1, the observed and simulated semimajor axis distributions are significantly different from each other. Thus, our results indicate that a disk mass less than ∼0.5 M_J does not lead to significant planet migration. Efficient planet migration requires that the outer disk mass be comparable to, or larger than, the planetary mass (Ivanov et al. 1999). Since the mass range of our planets is 0.4–4 M_J, it is expected that lower disk masses will lead to little planetary migration.

Zoom In Zoom Out Reset image size

For the e distribution, on the other hand, the K–S test shows that the null hypothesis cannot be rejected for t7d1 at 100 Myr. The distributions of t8d1 and t9d1 are dominated by planets with relatively small eccentricities (e ≲ 0.3), indicating that these cases are more dynamically stable compared to t7d1. On the other hand, t6d1 has a deficit of planets with e ≲ 0.1 and an overabundance of planets with e ∼ 0.2. A similar deficit of planets with low eccentricities is reported by JT08 for disk-less cases.

Table 2 summarizes the dynamical outcomes of our simulations, and lists the numbers of planets which are ejected out of the system, collide with the central star, and merge with another planet, before and after τ_GD. Rates of each of these events with respect to the number of multiple-planet systems are shown in parentheses. All of these rates increase for more massive disks, which implies that stronger convergent/divergent migration in the disk leads to more frequent dynamical instabilities. The most common outcome of such an instability is ejection, and 10%–20% of the planetary systems eject a planet before, and/or after τ_GD. Since the occurrence rates of dynamical instabilities are comparable before and after τ_GD, we check whether the observed a–e distribution is well-reproduced at τ_GD.

Table 2. Ejected, Collided, and Merged Planets

Set No.	Ejections			Collisions			Mergers
	Before τGD	After τGD	Total	Before τGD	After τGD	Total	Before τGD	After τGD	Total
t5d1	110 (36.7)	11 (23.4)	121	21 (7)	9 (19.1)	30	58 (19.3)	2 (4.26)	60
t6d1	62 (20.7)	35 (23.6)	97	6 (2)	11 (7.43)	17	42 (14)	7 (4.73)	49
t7d1	45 (15.3)	46 (21.6)	91	1 (0.333)	9 (4.23)	10	25 (8.33)	8 (3.76)	33
t8d1	35 (11.7)	29 (13.6)	64	0 (0)	7 (3.27)	7	15 (5)	2 (0.935)	17
t9d1	25 (8.33)	26 (10.2)	51	0 (0)	9 (3.53)	9	12 (4)	3 (1.18)	15
t5d1cr	118 (39.3)	16 (29.6)	134	18 (6)	3 (5.56)	21	54 (18)	4 (7.41)	58
t6d1cr	83 (27.7)	31 (27.9)	114	11 (3.67)	14 (12.6)	25	41 (13.7)	8 (7.21)	49
t7d1cr	39 (13)	46 (20.8)	85	1 (0.333)	10 (4.52)	11	19 (6.33)	9 (4.07)	28
t8d1cr	47 (15.7)	27 (12.1)	74	0 (0)	2 (0.893)	2	14 (4.67)	3 (1.34)	17
t9d1cr	31 (10.3)	36 (14.5)	67	1 (0.333)	2 (0.806)	3	11 (3.67)	3 (1.21)	14
t5d2	99 (33)	14 (30.4)	113	19 (6.33)	3 (6.52)	22	87 (29)	3 (6.52)	90
t6d2	61 (20.3)	31 (23.5)	92	8 (2.67)	5 (3.79)	13	61 (20.3)	3 (2.27)	64
t7d2	45 (15.3)	29 (16.2)	74	4 (1.33)	5 (2.79)	9	44 (14.7)	7 (3.91)	51
t8d2	33 (11)	19 (8.88)	52	2 (0.667)	5 (2.34)	7	37 (12.3)	4 (1.87)	41
t9d2	25 (8.33)	36 (16)	61	0 (0)	2 (0.889)	2	41 (13.7)	6 (2.67)	47
t5d2cr	95 (31.7)	12 (33.3)	107	18 (6)	1 (2.78)	19	89 (29.7)	1 (2.78)	90
t6d2cr	53 (17.7)	33 (23.6)	86	11 (3.67)	8 (5.71)	21	62 (20.7)	4 (2.86)	66
t7d2cr	41 (13.7)	30 (17.4)	71	5 (1.67)	2 (1.16)	7	52 (17.3)	6 (3.49)	58
t8d2cr	28 (9.33)	20 (10.1)	48	3 (1)	3 (1.51)	6	51 (17)	4 (2.01)	55
t9d2cr	27 (9)	17 (7.83)	44	0 (0)	4 (1.84)	4	43 (14.3)	2 (0.922)	45

Notes. Numbers of planets which are ejected from the system, collided with the central star, and merged with another planet, before and after the gas dissipation time τ_GD, as well as throughout the simulations. Larger number of collisions before τ_GD indicates more efficient planet migration, while the larger number of ejections indicates higher occurrences of dynamical instability. The percentages of planets in more than two-planet systems which are ejected/collided/merged are shown inside the brackets.

Download table as:  ASCIITypeset image

We perform the two-dimensional K–S tests for the observed a–e distribution and the corresponding distributions of our simulations (t5d1–t9d1) at τ_GD. The null hypothesis is rejected for all of these sets. These results indicate that planet–planet interactions after the gas disk dissipation must have played a significant role in determining the final a–e distribution. The one-dimensional K–S tests support this statement. For the a distribution, the null hypothesis is rejected for all the sets but t7d1 at τ_GD, while for the e distribution, the null hypothesis is rejected for all five sets at τ_GD.

Overall, our results imply that the a distribution is largely determined while a gas disk is still present, since the semimajor axis distribution does not change dramatically between τ_GD and 100 Myr. The e distribution is largely determined by planet–planet interactions after the disk dissipates, because the orbital eccentricities generally stay low while a gas disk is still present. Thus, we find that the initial assumption of nearly circular orbits in the previous N-body studies is reasonable.

4.2. Effect of the Saturation of Corotation Resonances

4.3. Effects of Different Initial e and i

We also study the evolution of three-planet systems with initially small e and i (sets denoted with d2, see Table 1). The overall results are similar to the cases with moderate e and i, for both with and without the saturation of corotation resonances.

For cases without the saturation of corotation resonances (t5d2–t9d2), the two-dimensional K–S tests against the observed and final distributions show that the null hypothesis cannot be rejected for t6d2 and t7d2. On the other hand, for cases with the saturation of corotation resonances (t5d2cr–t9d2cr), the corresponding tests show that the null hypothesis cannot be rejected for t7d2cr and t8d2cr. Thus, again, there seems to be an optimum disk mass to reproduce the a–e distribution.

One significant difference between the results for initially small or moderate e and i values is the rate of mergers. When initial systems are closer to coplanar, the merger rates tend to be higher. By comparing t5d1–t9d1 with t5d2–t9d2, we find that the merger rates before τ_GD are ∼1.5–3 times higher in t5d2–t9d2, while the merger rates after τ_GD are comparable. Similarly, the comparison of t5d1cr–t9d1cr with t5d2cr–t9d2cr shows that the merger rates before τ_GD are ∼1.5–5 times higher in the latter sets. We find that most of these mergers occur immediately after the start of the simulations, within 10⁴ yr. Ford et al. (2001) studied the dynamical evolution of equal-mass two-planet systems, and found that mergers tend to produce low eccentricity (e < 0.1) planets. However, in our simulations, we do not find any statistically significant change in the fraction of low e planets between d1 and d2 cases.

Figure 3 summarizes the a–e scatter plots for t7d1, t7d1cr, t7d2, and t7d2cr at 100 Myr. The two-dimensional K–S tests for the observed distribution covering the range 0.2 AU ⩽a ⩽ 6 AU, and a simulated distribution reject the null hypothesis for all the cases.

Zoom In Zoom Out Reset image size

4.4. Various Evolution Cases

Figures 4–8 show some examples of evolution of three-planet systems. Figure 4 is a dynamically "stable" system, in which there are no ejections, mergers, nor collisions for 100 Myr. Among the survived systems of our "successful" sets, which have the K–S probability of the a–e scatter plots of P > 0.1 (t6d1, t7d1, t6d1cr, t7d1cr, t6d2, t7d2, t7d2cr, and t8d2cr), the fraction of these "stable" cases increases as the initial disk mass decreases, and spans over ∼3%–19%.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figures 5 and 6 show the cases where dynamical instability occurs while the gas disk is still present. In the former figure, gravitational interactions between planets lead to the ejection of a planet, while in the latter, they lead to the collision with the central star. In these cases, the ejection/collision does not result in a high eccentricity of the remaining planets, probably due to the eccentricity damping in the disk. We find that about 40% of the planetary systems in the successful sets experience either collisions or ejections before the gas disk dissipates.

Figures 7 and 8 show the cases where dynamical instability occurs long after τ_GD. In the former figure, the strong interactions between planets lead to the ejection of one of the planets, and the collision of another planet with the central star. The planet left behind gains a large eccentricity. In the latter figure, such interactions lead to the collision of a planet with the central star, again, leaving the other planets on eccentric, and inclined orbits. About 40% of systems in the successful sets become dynamically unstable after the disk dissipation.

The a–e distribution for all the successful sets show the same trend as the representative case in Section 4.1—the observed a–e distribution is poorly reproduced at τ_GD, while at the end of the simulations (100 Myr), the observed distribution is well reproduced. Thus, although the rates of systems that become dynamically unstable are similar before and after τ_GD, we find that the eccentricity distribution is largely determined after the disk dissipation, as the previous N-body simulations implicitly assumed.

4.5. Mass Distribution

We also studied the a − M and e − M distributions for these sets. We find that the sets with the K–S probability of the a–e scatter plots being P > 0.1 also have a high K–S probability for a–M and e–M distributions (see Table 3).

In Figure 9, we plot a and e against M_psin i for a representative case of t7d1. The other cases (t6d1, t6d1cr, t7d1cr, t6d2, t7d2, t7d2cr, and t8d2cr) look similar to this. Here, we choose i randomly from 0–180 deg for the simulated planets. These plots mimic the expected a–M and e–M plots, assuming the simulated systems were observed, and the planetary orbits were inclined randomly with respect to the plane of the sky. From these plots, we expect that the future observations will find lower mass planets (M_psin i ≲ 0.3 M_J) on large semimajor axis (≳1 AU).

Zoom In Zoom Out Reset image size

4.6. Mean-motion Resonances

Although it is still too early to derive any statistical trend, some of the extrasolar planetary systems are observed to be in MMRs. Thus, it is interesting to investigate whether any of the simulated systems are in such configurations.

At the end of the simulations, our "successful" cases, t6d1, t7d1, t6d1cr, t7d1cr, t6d2, t7d2, t7d2cr, and t8d2cr have ∼20–70 multiplanet systems. For each of these systems, we estimate whether they are in a particular resonance by using the following resonance variable (Murray & Dermott 1999):

$$\begin{equation} \varphi = j_1\lambda _{o}+j_2\lambda _{i}+j_3\varpi _{o}+j_4\varpi _{i}, \, \end{equation} \tag{ 7 }$$

where λ and ϖ are the mean longitude and longitude of pericenter, respectively, and the subscripts i and o indicate inner and outer planets. Here, we focus on near-coplanar cases, and thus neglect the terms regarding the longitude of ascending node. When planets are in p + q:p MMR, we can define (j₁,  j₂,  j₃,  j₄) = (p + q, − p, − q,  0) or (p + q, − p,  0, − q).

We follow first- to third-order resonances (2:1, 3:2, 3:1, 5:3, 4:1, 5:2), as well as some higher order resonances (5:1, 7:3, 6:1, 7:2, 7:1, 8:1, 9:2, 9:1, 10:1). In Table 4, we summarize the numbers and kinds of MMRs seen in our systems at the end of the simulations. We find that most of these systems get trapped in MMRs while a gas disk is still present, either via migration or planet–planet scattering. A recent N-body study of three-planet systems (with no gas disk) by Raymond et al. (2008) showed that planet–planet scatterings can populate both low- and high-order MMRs. Our simulations confirm their result, and suggest that the combined effect of disk–planet and planet–planet interactions can generate planets in a variety of MMRs.

Table 4. Numbers of Planetary Systems in MMRs

Set No.	2:1	3:2	3:1	4:1	5:2	5:1	6:1	7:2	7:1	8:1	9:2	9:1	10:1
t6d1	6	0	16	5	4	1	2	2	0	0	0	0	0
t7d1	23	1	22	7	5	1	1	1	1	1	0	0	0
t6d1cr	5	0	6	5	2	0	0	1	0	1	0	0	0
t7d1cr	23	0	17	7	7	3	5	2	2	0	0	0	1
t6d2	10	0	17	7	4	0	1	0	0	0	0	0	0
t7d2	14	2	19	11	6	0	2	1	1	1	1	0	1
t7d2cr	14	2	19	11	6	0	2	1	1	1	1	0	1
t8d2cr	21	0	28	7	13	6	4	2	0	0	0	1	0

Notes. The numbers of systems which have each kind of MMRs. The shown sets are successful cases which have P > 0.1 for the K–S test of the observed and simulated a–e scatter plot. Most of these planets are trapped in MMRs while a gas disk is still present.

Download table as:  ASCIITypeset image

In our simulations, most multiplanet systems turn out to be in MMRs (∼70%–95%), while in Raymond et al. (2008), roughly 5%–10% of dynamically unstable systems ended up being in MMRs. This high rate of planets in MMRs seen in our simulations is clearly inconsistent with the observed systems. Adams et al. (2008) showed that the stochastic forcing effects of turbulence tend to pull planets out of MMRs. If this is the case, the number of planets in MMRs may be much lower.

We have studied the evolution of multiple-planet systems in a dissipating gas disk by means of a hybrid code that combines the N-body symplectic integrator SyMBA, and a one-dimensional gas disk evolution code. We simulate disk accretion, gap opening by planets, and planet migration in a self-consistent manner, and also take account of the effects of eccentricity damping by disk–planet interactions as well as gas accretion by planets. The main goal of this study is to investigate different plausible scenarios and understand how various initial conditions affect the final distributions of observable orbital properties. Specifically, we have investigated the evolution of three-planet systems in a gas disk by utilizing a range of planetary and orbital properties (e.g., M_p, a, e, and i), as well as disk masses, and taking account of the absence/presence of saturation of corotation resonances.

The initial conditions of our planetary systems are motivated by the standard planet formation theory. We assume that giant planets are formed via core accretion, beyond the ice line, on initially nearly circular, and coplanar orbits (see Section 3). Starting with three-planet systems that are initially fully embedded in gas disks, we have shown that the observed a–e distribution can be well reproduced by such models (see Section 4).

It is interesting to further compare the orbital properties predicted by our simulations with current observations. Recently, Wright et al. (2009) pointed out possible differences in distributions of orbital parameters between single- and multiple-planet systems. They found that the highly eccentric orbits (e > 0.6) are predominantly associated with apparently single-planet systems. The results of our simulations agree with this, and show that the high eccentricities (e ≳ 0.5–0.6) predominantly occur among single-planet systems at the end of the simulations. This implies that currently observed apparently single-planet systems may have experienced strong dynamical instabilities in the past.

They also pointed out that planets with masses >1 M_J have a broader eccentricity distribution, compared to less massive ones. Although planetary masses in our simulations span only the relatively narrow range 0.4–4 M_J, we do find a broader eccentricity distribution for planets with masses above 1 M_J. Such a trend may be naturally explained as a result of strong gravitational interactions between planets, which tend to remove the less massive planet out of the system via the ejection, or the collision, leaving the more massive one on an eccentric orbit.

Another striking trend that Wright et al. (2009) pointed out is the difference in semimajor axis distributions for single- and multiple-planet systems. The a distribution of single-planet systems is characterized by the so-called three-day peak and the jump in planetary abundance beyond 1 AU, while the corresponding distribution for multiple-planet systems is rather uniform. Such a difference may be well explained if the evolution of single-planet systems has been dominated by planet–planet scatterings as opposed to planet migration, while that of multiple-planet systems has been dominated by planet migration. Although we do not see such a difference in orbital distributions between single- and multiple-planet systems in our simulations, this may be partly because we do not have good numerical resolutions at small orbital radii, where the difference between these two distributions becomes prominent. Interestingly, we find a group of planets with a range of e that are located around 0.1 AU just before they are ejected out of the systems, or collide with the central star. Although our code does not include the tidal evolution directly, separate calculations of the tidal evolution of these systems indicate that they could be candidates for the observed close-in, tidally affected planets. Nagasawa et al. (2008) pointed out that such close-in planets can be formed even without migration in gas disks. They studied the orbital evolution of systems with three Jupiter-mass planets by including the tidal circularization, and showed that the Kozai mechanism combined with tidal orbital circularization (i.e., Kozai migration) can form these close-in planets for ∼30% of their simulations.

Recent observations of the Rossiter–McLaughlin effect started revealing an interesting subset of close-in exoplanets which have clearly misaligned orbits, or even retrograde ones (e.g., Winn et al. 2009; Narita et al. 2009). Out of at least 17 transiting systems with the observed projected stellar obliquity, eight have a clearly misaligned orbit (>10 deg), among which three are on retrograde orbits. We find that, from both simulations of three-planet systems with and without gas disks, the fractions of clearly inclined planets are ∼15% and ∼65%, respectively, while the rate of planets on retrograde orbits is less than 1% for both types of simulations. Thus, our results indicate that mechanisms other than planet–planet scatterings may be responsible for the observed high rate of close-in, retrograde planets. Note, however, Nagasawa et al. (2008) showed that inclinations could be more broadly distributed if close-in planets are formed via Kozai migration. If this is the case, the inclination distribution of close-in planets may be significantly different from that of planets on further orbits.

Our investigations are affected by several significant uncertainties. First of all, the initial conditions for these kinds of simulations are highly uncertain. For example, we assume that nearly fully grown giant planets are initially embedded in a gas disk. In reality, however, planets would start opening gaps as they grow, and therefore they should not be fully embedded in a gas disk initially. However, the planets in our simulations open gaps in a time on the order of the orbital period (i.e., less than several tens of years), which is shorter than, or at most comparable to, both the dynamical and migration timescales. Therefore, we do not expect a huge difference in the outcome due to this approximation.

Also, our choice of the initial planetary properties like mass and semimajor axis, although motivated by the core accretion scenario, is rather arbitrary. To better approximate the initial conditions, we could have performed planet formation simulations as in Thommes et al. (2008). However, such simulations are computationally very expensive for executing statistical studies.

Second, planet migration in our models does not take into account the effects of corotation resonances. Since corotation torques are sensitive to sharp gradients in the surface mass density, they may have a significant effect on some of the planets simulated here, which are massive enough to open a gap in the disk, but not massive enough to open a clean gap and saturate corotation resonances. The so-called Type III migration due to the corotation torques tends to accelerate the inward planet migration (Masset & Papaloizou 2003), and thus we are likely to underestimate planet migration rates for intermediate-mass (Saturn-like) planets, which do not open a clean gap.

Regarding eccentricity evolution, we focus on the effects of first-order resonances. Although this is a fair approximation for planets with small eccentricities (Goldreich & Tremaine 1980), higher-order resonances may become important if e is excited to a large value in the disk. Taking account of higher-order resonances, Moorhead & Adams (2008) semianalytically showed that the eccentricity evolution timescales are decreased by a factor of a few, but the overall trend of eccentricity evolution (e.g., damping or driving of e) does not change. Furthermore, they suggested that eccentricity decreases for 0.1 ≲ e ≲ 0.5 as long as corotation resonances are unsaturated, while eccentricity can rapidly increase when corotation resonances are fully saturated. On the other hand, current hydrodynamic simulations show that disk–planet interactions excite eccentricity up to only ∼0.2 (e.g., D'Angelo et al. 2006). As already mentioned, our simulations do not take into account the possibility of eccentricity excitation due to disk–planet interactions. Thus, our simulations may underestimate the eccentricity excitation rate for massive planets, which open a clean gap, and have a sufficient eccentricity to saturate corotation resonances. Clearly, this is an issue which requires a further investigation.

Finally, the gas disk in our simulations is removed exponentially once the randomly selected disk dissipation time is reached. Although such a disk removal is included to mimic the effect of photoevaporation, we did not model its physics directly. This is because the photoevaporation rate is difficult to estimate accurately, due to its sensitivity to the stellar flux, which in turn depends on the disk accretion rate, as well as the stellar environment (Matsuyama et al. 2003).

Bearing these in mind, our simulations successfully reproduced the general trends of the observed properties, starting with the initial conditions expected from the core accretion scenario. We summarize our findings below.

• 1.  
  Although the frequency of occurrence of dynamical instabilities before and after the disk dissipation is comparable, the e distribution is largely determined by planet–planet interactions after τ_GD.
• 2.  
  The a distribution is largely determined by disk–planet interactions. To explain the current a distribution, the disk mass has to be comparable to, or larger than, the planetary mass.
• 3.  
  There may be an optimum disk mass to reproduce the observed a–e distribution.
• 4.  
  Dynamical instabilities in a gas disk which involve ejections/collisions/mergers tend not to lead to large eccentricities of remaining planets.
• 5.  
  Initially nearly coplanar systems tend to have a higher merger rate between planets.
• 6.  
  For the range of disk masses we use, we do not see a significant difference in outcomes when taking account of the saturation of corotation resonances, and hence the reduction in eccentricity damping.
• 7.  
  We find that the combined effects of disk–planet and planet–planet interactions lead to both low- and high-order MMRs. Our results also indicate that, without perturbing effects (e.g., turbulence), too many planets may be trapped in MMRs.
• 8.  
  Starting with relatively well aligned, prograde orbits, we find that planet–planet interactions are not efficient in producing retrograde orbits.

This work was supported by NSF Grant AST-0507727 (to F.A.R.) and by a Spitzer Theory grant, as well as a grant from the Natural Sciences and Engineering Research Council of Canada (to E.W.T.). We thank the anonymous referee for careful reading, and helpful suggestions, which have improved this paper significantly. S.M. thanks Ralph Pudritz for making SHARCNET available for many of the simulations shown in this work.