PLANET–PLANET SCATTERING IN PLANETESIMAL DISKS. II. PREDICTIONS FOR OUTER EXTRASOLAR PLANETARY SYSTEMS

We adopt a simple and well-defined version of the conceptual model shown in Figure 1 as the initial conditions for our scattering simulations. We assume that three planets form randomly separated by 4–5 mutual Hill radii R_H,m,

$$\begin{equation} R_{H,m} = \frac{1}{2} \, (a_1 + a_2) \, \left(\frac{M_1 +M_2}{3\, M_\star }\right)^{1/3}. \end{equation} \tag{ 1 }$$

Here, a is the orbital semimajor axis, M is the planetary mass, M_⋆ is the stellar mass (fixed at 1 M_☉), and subscripts 1 and 2 refer to the inner and outer planet, respectively. The spacing of 4–5 R_H,m, which is common across all but one of our sets of runs, was chosen to yield systems that are unstable on a 10⁵–10⁶ year timescale (Chambers et al. 1996; Marzari & Weidenschilling 2002; Zhou et al. 2007; Chatterjee et al. 2008). The exception is our ensemble of runs with three 3 M_J planets, for which the number of unstable cases was so small that we adjusted the spacing to be 3.5–4 R_H,m and re-ran the simulations (we still use the set of more widely spaced, mostly stable simulations in our analysis of MMRs in Section 4).

The outermost planet of the three planets was placed two (linear) Hill radii,

$$\begin{equation} R_H = a \left(\frac{M}{3M_*} \right)^{1/3}, \end{equation} \tag{ 2 }$$

interior to 10 AU. We then chose a single random variable, uniformly distributed in the range between 4 and 5, to define the interplanetary spacing. The two additional planets were placed on the appropriate orbits, moving toward the (solar-mass) star. Planets were given zero eccentricity and randomly chosen mutual inclinations of less than 1°. We performed 10 sets of simulations, varying the planetary mass distribution in each set. For our two largest sets (1000 simulations each), we randomly selected planet masses according to the observed distribution of exoplanet masses (Butler et al. 2006),

$$\begin{equation} \frac{{d}N}{{d}M} \propto M^{-1.1}. \end{equation} \tag{ 3 }$$

In the Mixed1 set, we restricted the planet mass M_p to be between a Saturn mass M_Sat and three Jupiter masses M_Jup. For our Mixed2 set, the minimum planet mass was decreased to 10 M_⊕.¹⁰

Although, strictly, we have little knowledge of the actual mass function of extrasolar planets at orbital radii beyond 5 AU, we regard the runs set up with the mass function observed at smaller radii as the most realistic. To gain additional insight into the behavior of more idealized systems, we performed four additional sets of simulations with equal-mass planets (500 simulations each), and four sets that included radial mass gradients (250 simulations each). These sets are named with the appropriate planet masses, starting with the closest to the star outward; for example, in the 3J-J-S set the three Jupiter-mass (3 M_J) planet is closest to the star, followed by a Jupiter-mass planet and a Saturn-mass planet. The eight sets are: 3J-3J-3J, 3J-3J-3J, S-S-S, N-N-N, 3J-J-S, S-J-3J, J-S-N, and N-S-J, where 3J refers to a planet mass of 3 M_J, J is M_J, S is a Saturn mass M_S, and N is 30 M_⊕ (the "N" is intended to refer to Neptune, although the mass is augmented by just under a factor of 2 from Neptune's true mass of 17 M_⊕ to maintain roughly a factor of 3 between the different planet masses).

Each of our ∼5000 simulations was run twice: once with just the three planets and once also including an external planetesimal disk. The simulations without disks were presented in Raymond et al. (2008a, 2009b) and are used in this paper only as a comparison sample. Simulations with disks included 1000 planetesimal particles distributed between 10 and 20 AU following a radial surface density profile Σ ∝ r⁻¹, roughly consistent with submillimeter observations of outer disks around young stars (Andrews & Williams 2007). In all cases, the total mass of the disk was 50 M_⊕. We note that the inner edge of the disk lies interior to the radius where a test particle in the restricted three-body problem would be stable, so the disk is in immediate dynamical contact with the outer planet.

Each simulation was integrated for 100 Myr using the hybrid integrator in the Mercury simulation package (Chambers 1999) with a 20 day time step. There is the potential for significant numerical error in the integration for objects with small perihelion distances, which is usually manifested in terms of a secular increase in energy for close-in particles (Rauch & Holman 1999; Levison & Duncan 2000). With our time step of 20 days, the conventional wisdom that at least 10 steps are required per orbit implies that we would expect to be able to resolve the orbit of a body at 0.67 AU. In fact, a more detailed test using a test particle being forced into the star by a giant planet via the Kozai mechanism shows that the integration error remains less than 10⁻⁴ down to a perihelion distance of about 0.4 AU.

We gauged the fidelity of the outcome of each of our simulations based on the simple energy criterion dE/E < E_X, where E_X is an energy threshold. For the simulations without disks, we used the value of Barnes & Quinn (2004), who showed that E_X = 10⁻⁴ is adequate to test for dynamical stability of multi-planet systems. For the simulations with disks, the typical simulation error was actually dE/E ∼ 10⁻⁴, even for cases in which the planets were stable and experienced very little orbital evolution. By sifting through a large number of examples and looking at the energy conservation of individual bodies, it appeared that an error threshold of E_X ≈ 10⁻³ was adequate to yield reliable results. We therefore set a threshold E_X = 5 × 10⁻⁴. Runs that failed to meet our energy criterion were either rerun or rejected. For the Mixed1 and Mixed 2 sets, all runs with dE/E > E_X were re-run with smaller time steps: 5 days for simulations with no planetesimal disks and 10 days for simulations with disks. Tests showed that limits of 5 (10) days sufficed to accurately resolve orbits down to perihelion distances of 0.15 (0.25) AU. In each set, some of the re-run simulations (typically 15–35) still did not meet our energy criterion and were discarded from the sample. For the remaining sets of simulations, our finite computational resources did not allow us to re-run those simulations with disks that failed the energy test. In these cases, we discarded the examples with poor energy conservation based on the initial 20 day time step runs. We do not believe that this introduced any substantial bias in our sample because certain expected outcomes occurred reliably; for example, the eccentricity distributions for the J-J-J and 3J-3J-3J simulations are virtually identical with and without disks.

The details of our planetary initial conditions are motivated primarily by considerations of simplicity, and by the desire to pick planetary separations that are not grossly unstable. Our spacing of 4–5 R_H,m would, however, be broadly consistent with models in which planets became temporarily trapped in MMRs during the gaseous disk phase (Snellgrove et al. 2001; Lee & Peale 2002; Thommes et al. 2008a). Figure 2 shows the location of the 2:1, 3:2, and 5:3 MMRs as a function of the summed mass of the planets involved. These are resonances that might be commonly populated during the gas disk phase, with the 5:3 MMR being generally unstable once the gas dissipates (Kley et al. 2004; Pierens & Nelson 2008). For a range of planet masses between 0.5 M_J and 5 M_J, at least one of these resonances lies in the range probed by our simulations. We observe, of course, that a model based specifically on the idea of resonance trapping followed by release would differ in detail from ours, since in this case the planetary spacing in units of mutual Hill radii would correlate with the planetary masses.

Zoom In Zoom Out Reset image size

Our choices of disk mass and radial extent are also somewhat arbitrary. A disk mass of 50 M_⊕ is similar to that used in successful outer solar system models, such as the Nice model (Gomes et al. 2005), and it is also comparable to the summed core masses that one would expect for three giant planets. It is thus roughly consistent with the idea that planetesimals in a disk with a continuous surface density distribution accreted to form cores inside some characteristic radius, while failing to do so further out. Our inner disk radius of 10 AU, on the other hand, is no more than a guess. It could typically be larger, in which case all of the disk-driven effects discussed in this paper would only be manifest for planets further out.

Given the spacing of 4–5 mutual Hill radii, we expect that the typical instability timescale in our simulations should be 10⁵–10⁶ years (Marzari & Weidenschilling 2002; Chatterjee et al. 2008). Indeed, the median time for the first close encounter between planets in the cases without disks ranged from 58,000 years (S-S-S) to 3.1 Myr (3J-3J-3J). Instabilities typically set in earlier for systems containing less massive planets (see Chambers et al. 1996), and when the planet mass increases outward (i.e., N-S-J simulations went unstable more quickly than J-S-N simulations).

Dynamically unstable systems tend to stabilize by destroying one or more planets, usually via collision with another planet or hyperbolic ejection. Systems with equal-mass planets tend to undergo far more close encounters over a much longer timespan than systems with significant mass differences between planets. The planets that survive in systems with equal-mass planets have much larger eccentricities than those in unequal-mass systems (Ford et al. 2003; Raymond et al. 2008a). The duration of instability is also mass dependent: systems with low but equal-mass planets remain unstable for longer periods than systems in which the planet masses are scaled up.

Figure 3 shows the cumulative distribution of (1) the timescale from the start of each simulation to the first close encounter between planets (shown in black) and (2) the duration of the system instability, i.e., the time from the first to the last close encounter. We plot all 10 sets of runs, both with and without planetesimal disks. For this unstable subset of runs, the 100 Myr duration of the simulations is sufficient to observe the onset and full duration of the instability in almost all cases.

Zoom In Zoom Out Reset image size

With the sole exception of the N-S-J runs, the onset of instabilities in systems with disks occurs sooner than in those without disks (Figure 3). One important effect of the presence of the planetesimal disks is clearly to stabilize systems that might otherwise be unstable on long timescales. This effect is most pronounced in systems with low-mass outer planets (e.g., the J-S-N cases). The reason for this is likely related to the mass dependence of the migration rate of a planet through a planetesimal disk, which is a strongly decreasing function of planet mass once the mass exceeds the mass of planetesimals within a few Hill radii (Ida et al. 2000; Kirsh et al. 2009). This migration can cause two planets (usually the outer two) to cross a low-order mutual MMR (see Figure 2). MMR crossings lead to an increase in planetary eccentricities (e.g., Chiang et al. 2002; Tsiganis et al. 2005), and trigger an instability that leads to close encounters.

Disks also reduce the duration of the unstable phase, but only in systems containing lower-mass outer planets. The most dramatic example are the N-N-N simulations, for which the typical period of close encounters was reduced by a factor of more than 100 (Figure 3). This occurs because the orbits of high-eccentricity planets are circularized via dynamical friction with the planetesimal disk (e.g., Thommes et al. 1999), and this process is faster and more efficient for lower-mass planets. As scattered planets are re-circularized by the disk, their periastron distances are increased, they are removed from the region of scattering and future encounters between planets do not occur. In contrast, the planetary dynamics in systems with a massive planet adjacent to the planetesimal disk are effectively shielded from the effects of the planetesimal disk such that the duration of instabilities is the same as for simulations with no disks.

In planetary systems with no dissipation, once an instability occurs the endpoint is almost inevitably a violent event: a planet–planet collision, a planet–star collision, or a hyperbolic ejection from the system.¹¹ However, in systems with dissipation from gaseous or planetesimal disks, the situation changes and the number of potential outcomes increases. The key new processes that arise from the presence of a planetesimal disk occur primarily for planets whose masses are similar to the planetesimal disk mass. Given our disk mass of 50 M_⊕, it is the simulations that include at least one roughly Saturn-mass or smaller planet that are affected. Interactions between the planetesimal disk and lower-mass planets can lead to angular momentum exchange and radial migration (Fernandez & Ip 1984). In addition, dynamical friction from the planetesimal disk can re-circularize the orbit of a scattered planet. Both of these processes feature prominently in recent models of the dynamical evolution of the solar system's giant planets (Malhotra 1995; Thommes et al. 1999; Tsiganis et al. 2005; Ford & Chiang 2007).

Figure 4 shows the range of outcomes of the unstable simulations for each giant planet configuration. The number of planets surviving per system increases for simulations with disks, but only for cases that contain a low-mass planet. The most dramatic difference is seen in the N-N-N simulations, for which the number of planets per unstable system increased from 1.98 to 2.87, meaning that only one out of every 7.5 unstable simulations with disks destroyed a planet. The pie charts in Figure 4 show what happened to the destroyed planets in each case. For almost all cases, the dominant destruction mechanism was ejection from the system, which accounted for at least 3/4 of the destroyed planets in all configurations except S-S-S and N-N-N. In the N-N-N simulations without disks, ejection and collisions played a comparable role, which is not surprising given that the ratio of the escape speed from the planet to the escape speed from the planetary system (sometimes called the "Safronov number," Θ), is the smallest for any system, only about 2 (Ford et al. 2001; Goldreich et al. 2004). In the N-N-N simulations with disks, gravitational kicks during close encounters were not strong enough to eject a single planet in any simulation. The reason for this is twofold. First, dynamical friction efficiently circularizes the orbits of low-mass planets such that the typical encounter speed is low. Second, dynamical friction is also effective at effectively capturing scattered planets in the disk since the total disk mass actually exceeds the planet mass (30 M_⊕) in this case.

Zoom In Zoom Out Reset image size

We have seen that the planetesimal disk can have a strong influence on the planets' dynamical evolution. The planets have a back-reaction on the disk itself, pumping up planetesimal eccentricities via dynamical friction and removing planetesimals from the system via dynamical ejection or collision. The number and strength of planet–planetesimal interactions depend most strongly on the planetary orbits. If the planets' orbits remain circular then the planets can interact with only a fraction of the disk and encounter velocities depend mainly on the planetesimal eccentricities. If, however, planetary orbits are eccentric then the planets can interact with a larger fraction of the disk and encounter speeds will be larger, resulting in more frequent ejection of planetesimals. Thus, the perturbations felt by the planetesimal disk will correlate with the planetary orbits and therefore the dynamical stability of the planets.

Figure 5 shows the timescale for the destruction of a given fraction of the planetesimal disk for all of the Mixed1 simulations, where each curve corresponds to a single simulation. Here, we have subdivided the "stable" simulations, which never experience close encounters, into moderately stable and stable subsets. We define a moderately stable system as one which never experienced a close encounter between planets but for which the final mass-weighted planetary eccentricity is larger than 0.025. These systems have thus experienced significant planet–planet perturbations despite the absence of close planetary encounters.

Zoom In Zoom Out Reset image size

The unstable Mixed1 systems destroy the majority of their planetesimal disks on a 10⁵–10⁶ year timescale (Figure 5), although late-onset instabilities can be seen as nearly vertical lines at later times. There is often a delay for the destruction of the last 10% of planetesimals because in most cases these have large inclinations such that the close encounters with the planets needed to reach zero energy are less frequent. In fact, many unstable cases do not destroy the entire planetesimal disk. Moderately stable systems can have a range in the timescale and amount of planetesimal destruction depending on the planet masses and evolution. For relatively massive planets, the sudden eccentricity increase that accompanies a resonance crossing leads to a corresponding jump in the eccentricities of many or most disk particles, leading to orbit crossings, close encounters and ejection. For less massive planets, secular perturbations are weaker so less of the disk is destabilized during such events. One configuration that is quite efficient at clearing out planetesimals is a lower-mass outer planet and a high-mass middle planet. The outer planet's low mass allows it to migrate outward somewhat due to planetesimal scattering. The scattered planetesimals, in turn, are quickly ejected by encounters with the high-mass middle planet.

The stable Mixed1 systems, which are relatively high mass, destroy only a relatively small portion of the planetesimal disk. The outer planet interacts directly with planetesimals within the stability boundary and also may excite planetesimals that happen to be in resonances onto planet-crossing orbits. These directly interacting planetesimals can be ejected by the outer planet or can be "passed" inward to interact with the inner planets, which typically eject the planetesimals. Thus, the inner planetesimal disk is quickly cleared out and interactions between the planets and the planetesimal disk become infrequent, with 50%–75% of the disk remaining.

The evolution of the planetesimal disks in our simulations offers a rich data set to understand the connection between planetary dynamics and disk structure. In a future paper, we will correlate the planetary and planetesimal components of each system and calculate the infrared detectability of these disks for each case. For the remainder of this paper, we restrict ourselves to the planetary dynamics.

Figure 6 shows the cumulative eccentricity distributions for the innermost and outermost planets in all of our unstable simulations, with and without disks. We also plot the observed exoplanet eccentricity distribution, excluding planets inside 0.1 AU which are likely to have had their orbits altered by tides (Jackson et al. 2008). We note that, since our simulations start with planets at fairly large radii, we do not populate the full radial range across which the observed distribution is determined. Currently, however, the evidence for any radial dependence to the eccentricity distribution is marginal (Ford & Rasio 2008), so the comparison between the observed and simulated distributions is justifiable.

Zoom In Zoom Out Reset image size

Although the result is now well established (Chatterjee et al. 2008; Jurić & Tremaine 2008), it is still startling to see from Figure 6 how easily planet–planet scattering can reproduce the observed eccentricity distribution. With just the simple assumption that planets have the observed mass distribution and form on marginally unstable orbits (i.e., the Mixed1 simulations), we can produce a statistical match to the observed distribution. The planet–planet scattering model appears to be a very robust and simple way to explain the observed systems.

As expected, the difference between planetary eccentricities for simulations with and without disks is strongly mass dependent and clearly seen in a comparison between the Mixed1 and Mixed2 simulations. The effect of the planetesimal disk is clear for M_p ≲ M_S, and for the lowest-mass case, N-N-N, orbital re-circularization from planet–planetesimal disk interactions is so strong that there are no high-eccentricity planets at all—the most eccentric planet in all the N-N-N simulations with disks is just 0.06. As seen in previous work, scattering among equal-mass planets produces larger eccentricities than scattering among planets with different masses (Ford et al. 2003; Raymond et al. 2008a), and this is clearly seen in Figure 6.

Variations between the eccentricities of inner and outer planets are caused by two different effects: scattering from other planets and eccentricity damping from the planetesimal disk. For the Mixed1 simulations, there is very little difference between the inner and outer eccentricity distributions because there is no planetary mass gradient and the planet masses are too large to be significantly affected by the planetesimal disk. For the Mixed2 simulations with disks, outer planets have lower eccentricities. This offset is not seen for the simulations without disks. Thus, the lower-eccentricity outer planets are due to damping from the planetesimal disk. For simulations with equal-mass planets, there is a modest inner versus outer eccentricity difference because of disk damping but only for the S-S-S cases, as the Jupiter-mass planets (J-J-J) barely feel the disk and 30 M_⊕ planets (N-N-N) are all quickly circularized. For systems with radial gradients in planet mass, there is an inner versus outer eccentricity difference caused by the simple mass dependence of the recoil velocity during a planetary encounter. During instabilities involving different-mass planets, this leads to more massive planets having smaller eccentricities than less massive planets. For systems with negative mass gradients (3J-J-S and J-S-N), this causes the outer planets to have higher eccentricities than the inner planets, and the effect is reversed for systems with positive mass gradients (S-J-3J and N-S-J). For the lower-mass cases with mass gradients (J-S-N and N-S-J), the effects of preferential eccentricity damping of outer planets by the planetesimal disk are also clearly seen.

In Figure 7, we show the eccentricity of the innermost (and thus most easily detectable) planet as a function of the final total mass of the planetary system, using just the Mixed1 and Mixed2 simulations. In the presence of disks, there is a strong correlation between these quantities (Raymond et al. 2009a). For systems with total masses of less than ∼1 M_J, circular or near-circular orbits dominate. For higher system masses, on the other hand, planetary eccentricities display a large range, from near zero to >0.8. This correlation is much weaker if we plot, instead, planetary eccentricity against individual planet mass for all surviving planets (Figure 8). As is obvious from the cumulative distributions, disks do act to circularize low-mass planets more than high-mass planets, but there is no sharp transition between the two regimes. We interpret this difference as being due to the fact that it is the total mass of the planetary system that determines the disk lifetime and the ability of the disk to damp eccentricities, rather than any individual planet mass (see discussion in Section 3.6).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The semimajor axis, beyond which our model predicts a transition to low-eccentricity orbits for M < M_J planets, depends upon the characteristic inner radius of the primordial planetesimal disk (note that the mass at which the transition occurs will also vary with this quantity). This is a free parameter, which we have fixed at 10 AU. If, in reality, planetesimal disks do typically extend in to about 10 AU, our results suggest that the transition to the low-eccentricity regime may be observable with feasible extensions of existing radial velocity surveys (if the inner extent of planetesimal disks is at substantially larger radii astrometry or direct imaging may offer better prospects). To quantify this, we have calculated the radial velocity amplitude K,

$$\begin{equation} K = \frac{1}{\sqrt{1-e^2}} \left(\frac{M}{M_*} \right) \sqrt{\frac{GM_*}{a}} \sin i, \end{equation} \tag{ 4 }$$

for the innermost surviving planet for all of the Mixed1 and Mixed2 runs. The distribution of e as a function of K is shown in Figure 9, assuming M_* = M_☉ and (for simplicity) that sin i = 1. Since the detectability of a planet via radial velocity measurements depends upon the period as well as the amplitude, the points are further coded to indicate relatively short-period (P < 2000 days), intermediate-period (2000 days < P < 4000 days), and long-period (P > 4000 days) planets.

Zoom In Zoom Out Reset image size

A measurement of e(K) for the innermost planet provides no information as to the total mass of the planetary system, and as a consequence the transition to lower eccentricities at lower K resembles that shown in Figure 8 (e as a function of individual planet mass) rather than the sharper transition seen in Figure 7 (e as a function of the final total system mass). Nonetheless, we see a clear trend to lower eccentricities that sets in for radial velocity amplitudes K ≲ 10 m s⁻¹. Although there are high e outliers at all values of K, for K = 1–5 m s⁻¹ the typical eccentricity lies in the 10⁻²–0.1 range, whereas for K > 10 m s⁻¹ values e > 0.1 are obtained. Comparing the results with and without disks makes clear that this difference is caused by the damping effects of planetesimals.

Detection of planets with radial velocity amplitudes K < 10 m s⁻¹, although not easy, is already possible (almost 40 such systems are currently known). The primary obstacle to observational study of the outer planet region, where the dynamical effects of planetesimal disks should become evident, is rather the long periods of planets. As is clear from Figure 9, essentially all of the low-mass planets, whose orbits have been partially circularized by interaction with planetesimals, have periods in excess of 4000 days. High precision radial velocity measurements, capable of finding planets with K ≈ 5 m s⁻¹ and periods in excess of 10 years, are therefore required in order to start probing the dynamics discussed in this paper.

Inspection of Figure 9 suggests a second observational signature of planetesimal disk dynamics: the existence of extremely low eccentricity planets. In the presence of disks, a small but measurable fraction of even massive planets (with K = 10–30 m s⁻¹) end up with e < 10⁻², whereas such systems are rare in the absence of disks. Extremely precise measurements of eccentricity, capable of detecting an excess (above the predictions of pure scattering models) of truly circular orbits, would therefore be valuable even if such observations were only available for massive planets. The interpretation of an excess of circular orbits, however, would likely be more ambiguous than a measurement of the full e(K) distribution, since it could reflect small amounts of gas damping or an admixture of systems that only formed one planet, as well as the dynamical influence of planetesimals.

Planet scattering can lead to large mutual inclinations between the surviving planetary orbits, as well as misalignment with respect to the initial plane of the planets (Chatterjee et al. 2008; Jurić & Tremaine 2008). Mutual inclinations between the orbital planes of planets in multi-planet systems can be detected via astrometry (e.g., Bean & Seifahrt 2009), while any misalignment between the orbital plane and the plane perpendicular to the stellar spin axis can be detected directly for transiting planets via the Rossiter–McLaughlin effect (Gaudi & Winn 2007; Winn et al. 2005). Another interesting consideration is that the initial planetary plane coincides with the plane of the planetesimal disk and presumably also with the plane of the dust disk that should be produced by collisional grinding of planetesimals (e.g., Wyatt 2008).

Figure 10 (dashed lines) shows the distribution of inclinations with respect to the initial orbital plane for the unstable systems in all 10 sets of simulations, both with and without planetesimal disks. As was the case for eccentricities, larger inclinations are generated in systems with equal-mass planets than in systems with large mass ratios. Inclinations in the few to ∼15° range are common in most sets of simulations, but only the equal-mass systems are able to generate larger inclinations. Unlike the eccentricity distributions for which higher-mass planets have higher eccentricities, all four equal-mass systems (without disks) have virtually identical inclination distributions.

Zoom In Zoom Out Reset image size

As expected, the simulations with disks yielded smaller inclinations than the simulations without disks, and the effect was stronger for lower-mass planets. Figure 11 shows the inclination of the innermost planet, which as in the case of eccentricity is controlled by the total planetary mass rather than the individual planet masses. The same sharp break at M_tot ≈ 1 M_J is seen to divide planets with small inclinations versus those with a large range in inclination.

Zoom In Zoom Out Reset image size

The mutual inclination Δi between two planetary orbits with inclinations i₁ and i₂ is given by

$$\begin{equation} \cos {\Delta i} = \cos {i_1} \cos {i_2} + \sin {i_1} \sin {i_2} \cos {\left(\Omega _1-\Omega _2\right)}, \end{equation} \tag{ 5 }$$

where Ω₁ and Ω₂ refer to the longitudes of ascending node. Figure 10 shows the distributions of the mutual inclination Δi between the innermost and the adjacent planet in unstable systems for which two or more planets survived (solid lines). The Δi values are consistently larger than the inclinations with respect to the initial orbital plane. During a close planetary encounter, if one planet is scattered in the $\hat{z}$-direction, then the other planet must receive a kick in the $-\hat{z}$-direction. For a single encounter that dominates the planets' motion in the vertical direction, such a kick would also cause the two planets' longitudes of ascending node to be offset by roughly 180°, depending on their eccentricities. Thus, the instantaneous mutual inclination will be close to the sum of the two planets' inclinations with respect to their initial plane (Equation (5)). Although the alignment of the two planets' nodes will change in time due to secular perturbations, angular momentum conservation requires that the mutual inclination remains relatively large.

As for the case of inclination with respect to the initial orbital plane, equal-mass planets yield much larger Δi values than planets with large mass ratios. However, in this case the lower-mass equal-mass systems have higher mutual inclinations than the higher-mass cases. Ten percent of the unstable planets in the equal-mass systems without disks had inclinations larger than 35° (3J-3J-3J), 39° (J-J-J), 42° (S-S-S), and 49° (N-N-N), compared with 26° for the Mixed1 simulations. The larger inclinations in lower-mass equal-mass planetary systems are caused by the ease with which the high-mass planets can eject each other. The higher-mass systems for which two planets survive have undergone far fewer close encounters than the lower-mass systems with two surviving planets: the median number of encounters in the unstable 3J-3J-3J, J-J-J, S-S-S, and N-N-N systems for which two plants survived was 81, 175, 542, and 1871, respectively. However, the mass-weighted eccentricities of the two-planet equal-mass systems were slightly higher for the higher-mass planets. Thus, large planetary eccentricities are linked with the strength of encounters between planets while large inclinations are linked with a large number of scattering events. Note that the number and strength of close encounters in systems with significant mass ratios are far less, thus explaining their lower eccentricities and inclinations (see Raymond et al. 2009b).

As seen in Figure 3, the duration of the close encounter phase (as well as the total number of encounters) is reduced for the simulations with planetesimal disks. Therefore, we expect lower inclinations for simulations with planetesimals—this is indeed clearly seen in Figure 10. The extreme damping for the N-N-N simulations is again seen in terms of their very low inclinations and mutual inclinations in systems with planetesimal disks.

What sets the boundary between typically eccentric (and inclined) and typically circular (and coplanar) planets, which we have determined empirically lies at a total system mass of about M_tot ≈ 0.7 M_J? Plausibly, this transition mass may be set simply by the reservoir of disk angular momentum that is able to interact with the planet and circularize its orbit.

To illustrate the point, we construct an (over)simple model of the circularization process. Let us assume that scattering (with or without disks) typically yields a single planet with some characteristic semimajor axis a, eccentricity e, and inclination i. The orbital angular momentum deficit (AMD) quantifies the difference between the orbital angular momentum for a circular orbit as compared to an eccentric and inclined one at the same semimajor axis. The AMD is given by

$$\begin{equation} {\rm AMD} = M \sqrt{GM_*a} \left(1 - \cos i \sqrt{1-e^2} \right). \end{equation} \tag{ 6 }$$

To fully circularize the planet, we need to add this much angular momentum via planetesimal interactions. If these interactions occur primarily at or near apocenter in a disk with surface density profile,

$$\begin{equation} \Sigma = C r^{-1}, \end{equation} \tag{ 7 }$$

with C being a constant, the mass of the disk within a radial zone of half-width n Hill radii (at the apocenter distance) is

$$\begin{equation} \Delta M_{\rm disk} = 4 \pi n C \left(\frac{M}{3 M_*} \right)^{1/3} a (1+e). \end{equation} \tag{ 8 }$$

Assuming the planetesimals to have circular orbits, the total angular momentum of the disk that interacts with the planet is, approximately,

$$\begin{equation} \Delta L_{\rm disk} \simeq \Delta M_{\rm disk} \sqrt{GM_* a(1+e)}. \end{equation} \tag{ 9 }$$

By equating the available angular momentum to the AMD, we find that the disk can circularize planets for masses,

$$\begin{equation} M \lesssim \frac{ (4 \pi n C)^{3/2} }{ (3 M_*)^{1/2} } \frac{ (1+e)^{9/4} }{ (1 - \cos i \sqrt{1-e^2})^{3/2} } a^{3/2}. \end{equation} \tag{ 10 }$$

This is an analog of the usual isolation mass, except here in the case of circularization rather than accretion. One should note that the factor involving e is not of order unity—typically it is quite large (∼10²).

This analysis is quite crude, but it does predict transition masses that are of the same order of magnitude as those observed. For example, substituting our disk parameters and adopting planetary orbital elements of a = 8 AU, e = 0.4, and i = 10°, we obtain (for n = 2) a transition mass M ≈ 0.6 M_J. This suggests that we should interpret the transition mass as being a consequence of the finite angular momentum, within the disk, available to circularize planets. So, why is it that the eccentricity is controlled by the total planet mass rather than the individual planet mass (e.g., Figures 7 and 8)? The answer appears to be that it is the total mass which regulates the eccentricity excitation, which in turn determines the amount of damping by scattering disk particles.

Planet–planet scattering tends to segregate systems by mass, with more massive planets closer-in and less massive planets farther out (Chatterjee et al. 2008). We expect this effect to be magnified in systems with planetesimal disks because the disk can "trap" low-mass planets that would otherwise be ejected from the system. To look at the radial mass evolution of our simulations, we restrict ourselves to the Mixed1 and Mixed2 cases which did not start with pre-defined radial mass gradients.

Figure 12 shows the median planet mass as a function of orbital distance for the unstable Mixed1 and Mixed2 systems, with and without disks (also shown is the 10%–90% range of planet masses in a given radial bin). For the Mixed1 systems, there is a clear segregation of high-mass planets to the inner system and low-mass planets in the outer system, with very little difference for the simulations with and without planetesimal disks. The most massive planets tend to reside at 4–9 AU, close to the starting initial orbits.

Zoom In Zoom Out Reset image size

Radial mass segregation is also clear in the Mixed2 simulations but there are additional subtleties. As for the Mixed1 cases, high-mass planets are confined within 10 AU and low-mass planets outside 10 AU. There is a large peak in the median Mixed2 planet mass between 8 and 9 AU, both with and without planetesimal disks. There are actually 2–3 times more planets in this 1 AU-wide bin than in adjacent bins for the Mixed2 simulations with and without disks. Perhaps surprisingly, any differences in mass segregation caused by disks appear to be modest even for the Mixed2 simulations.

Scattering tends to spread planetary systems out. We can quantify this by looking at the planetary separations at the end of the simulations in units of mutual Hill radii R_H,m (recall that all simulations started separated by 4–5 R_H,m). Figure 13 shows the separation of the two inner planets as a function of the total mass in surviving planets for the Mixed1 and Mixed2 simulations, with and without disks. There are several interesting pieces of information in Figure 13. First, the typical interplanetary spacing is ∼4–30 R_H,m, showing that most systems have indeed spread out. Second, the interplanetary spacing decreases significantly for larger system masses, although there remains a large spread for any given mass. This general trend can be explained by the much larger number of scattering events undergone by lower-mass systems before the destruction (usually by ejection) of a planet.¹² The spread in separation is due in part to stochastic variations in scattering events from simulation to simulation and in part to variations in planetary mass ratios—the larger number of scattering events for equal-mass systems causes them to be more widely spaced than systems with larger mass ratios. Third, there is an abrupt break between the separations of lower-mass systems with and without planetesimal disks. For M_tot ≲ 0.7 M_J, systems with disks are much more compact than systems without disks. This is a result of the damping of the planetesimal disk: low-mass planets on wide orbits have their eccentricities efficiently damped and their perihleion distances increased such that they avoid additional close encounters with the inner planet. This is evidenced by the fact that the break between the simulations with and without disks occurs at the same total planet mass as for the eccentricity and inclination distributions (see Figures 7 and 11).

Zoom In Zoom Out Reset image size

For systems in which three planets survive, the outer two planets tend to be more widely spaced than the inner two. Figure 14 shows the separation of each pair of planets for unstable systems which preserved all three planets. In the bulk of cases, the outer pair of planets are indeed more widely spaced than the inner pair. However, there do exist many cases for which the outer planets are more compact, as well as a few cases in the Hill unstable region which are likely to be unstable on slightly longer timescales. Systems with a more compact outer pair of planets tend to be those with lower-mass middle planets. As the outer planet scatters planetesimals inward, a massive middle planet can eject them from the system, causing the planet to move inward and leading to a compact inner system but a spread out outer system. In contrast, a low-mass middle planet cannot generally eject the scattered planetesimals and so simply scatters them toward the inner planet, causing the middle planet to move outward, yielding a compact outer system and a more spread out inner system.

Zoom In Zoom Out Reset image size

Planet–planet scattering creates a population of high-eccentricity planets with large apocenter distances (Veras et al. 2009; Scharf & Menou 2009). For the vast majority of cases, these planets, which can have semimajor axes as large as ≳10, 000 AU, represent a transient phase on the path to dynamical ejection.¹³  These planets are of interest because they represent the source of the so-called free-floating planets. In some cases, their large separations mean that they could be good targets for direct detection, although they probably only exist for the first 10–100 Myr of a star or star cluster's lifetime.

In our simulations, we imposed an ejection radius of 100 AU, so we were unable to probe very distant planets. Nonetheless, we can study planets that ended up on transient, or in some cases stable, orbits beyond the initial radial extent of our initial conditions. To address the issue of planets on widely separated orbits, we kept track of the time spent by each planet in each simulation with an aphelion distance in radial bins of 25–50 AU, 50–75 AU, and 75–100 AU. We also characterized planets in each of those bins as either being in a transitional state (usually on the path to ejection) or on stable orbits. Inspection of the final eccentricity as a function of semimajor axis distributions, plotted in Figure 15, shows immediately that a large number of scattered planets survived on stable orbits with large aphelion distances.

Zoom In Zoom Out Reset image size

Table 2 shows the statistics of planets at large orbital separations in our simulations. The table includes information about both stable and transitional planets with aphelion distances Q > 25 AU and Q > 50 AU. First of all, we see that scattering among equal-mass planets is far more efficient at producing stable planets with large orbital radii than scattering among planets with mass gradients. This is simply because of the much larger number of scattering events that occur in equal-mass systems. Next, we see that this trend does not hold for the fraction of systems which experience a transitional, high-Q phase. Lower-mass planets have a higher probability of experiencing this transitional phase. We interpret this as an ejection timescale issue—the number of encounters and duration of the instability phase increases dramatically for lower-mass systems (see Section 3.1) such that the typical ejection is very drawn out. In contrast, for higher-mass planets a smaller number of scattering events is needed to eject a planet. Our analysis only includes outputs every 10⁵ years; thus, ejections which occur in less than that time are not registered or counted in Table 2. For equal-mass and mixed systems, the mean time during our 100 Myr simulations for which a planet could be observed with Q > 25 [50] AU correlates with the fraction of transitional systems, with typical observable probabilities of a couple of percent (the probability is here calculated as 〈t〉 divided by the 100 Myr simulation length). For the mass gradient simulations, 〈t〉 is longer for negative mass gradients (for 3J-J-S and J-S-N). This is because the most massive planet is closer-in in those cases, deeper in the star's potential well, such that a larger number of encounters is needed to reach the system's escape velocity.

As expected, the higher-mass systems are barely affected by the presence of planetesimals disks in terms of the fraction of orbits which are at large distances and the fraction that appear in the transitional phase. However, the duration of the transitional ejection phase is lengthened by a factor of a few for the equal-mass cases, meaning that the time required to eject those planets is increased due to the damping from the disk. For the lower-mass planets, there are fewer planets at large orbital distances in terms of both stable and transitional configurations. This is simply due to the disk's ability to "trap" planets by orbital circularization. This circularization occurs within the disk, so trapped planets are generally confined to orbits within the outer edge of the disk, at 20 AU.

Can observations of long-period planets tell us anything about the inner planets that spawned them? Assuming that outer planetesimal disks are ubiquitous, equal-mass and higher-mass planetary systems are more likely to populate distant stable orbits than lower-mass systems or those with large mass ratios between planets. However, a wide range of configurations produces long-lived transient planets on large orbits that are generally on their way toward dynamical ejection. We suspect that combining observations of long-period planets with infrared observations of the outer planetesimal disk may constrain the problem, but that is beyond the scope of the current paper.

A more dynamically sophisticated way to interpret the radial distribution of planets in multiple planet systems is to measure their separation in terms of the minimum required for Hill stability. This approach incorporates eccentricity and inclination information consistently and, if there are only two planets in the system, connects to the formal definition of Hill stability that has been proven for that case (Marchal & Bozis 1982; Gladman 1993). Observationally, the known two-planet exosystems are known to cluster close to the Hill stability boundary (Barnes & Greenberg 2006). This motivates the hypothesis that all multiple planet systems may be dynamically "packed," in the sense that additional planets could not exist between the known planets on stable orbits (see Barnes et al. 2008 or Raymond et al. 2008b). The corollary of such a "packed planetary system" hypothesis (Barnes & Raymond 2004; Raymond & Barnes 2005; Raymond et al. 2006; see also Laskar 1997), of course, is that when planetary systems are observed not to be packed we should seek an additional "missing" (normally lower mass) planet in what appears to be an empty stable orbit.

It is easy to test numerically whether an observed multiple planet system is packed (e.g., Rivera & Lissauer 2000; Jones et al. 2001; Menou & Tabachnik 2003; Asghari et al. 2004; Raymond et al. 2008b; Kopparapu et al. 2009). In one interesting case, the HD 74156 system was dynamically mapped to reveal a narrow zone between the two known planets that was stable for Saturn-mass test planets but not for Jupiter-mass test planets (Raymond & Barnes 2005). Subsequently, Bean et al. (2008) presented evidence for a third planet in the system, in the stable zone and with a mass of 1.4 Saturn masses (Bean et al. 2008; Barnes et al. 2008).¹⁴ In a previous paper, we showed that pure planet–planet scattering naturally reproduces the observed distribution of dynamical configurations (Raymond et al. 2009b). Given the significant effects that planetesimal disks can have on planetary evolution, here we investigate whether this is still true at radii where the dynamical influence of disks is significant.

For two planets with masses M₁ and M₂ orbiting a star, it can be shown analytically that long-term dynamical stability requires

$$\begin{eqnarray} &&\frac{-2 (M_*+M_1+M_2)}{G^2 (M_1 M_2 + M_* M_1 + M_* M_2)^3} \, c^2 h \ge \nonumber \\ &&1 + 3^{4/3} \frac{M_1 M_2}{M_*^{2/3} (M_1+M_2)^{4/3}} - \frac{M_1 M_2 (11 M_1 + 7 M_2)}{3 M_* (M_1 +M_2)^2}, \end{eqnarray} \tag{ 11 }$$

where c and h represent the total orbital angular momentum and energy of the system, respectively (Marchal & Bozis 1982; Gladman 1993; Veras & Armitage 2004; note that this definition assumes that M₁ > M₂). Following Barnes & Greenberg (2006), we refer to the left side of Equation (11) as β and the right side as β_crit. The quantity β/β_crit therefore measures the proximity of a pair of orbits to the Hill stability limit of β/β_crit = 1. It is important to note that our β/β_crit analysis only applies for two-planet systems, because perturbations from additional companions can shift the stability boundary to values other than 1.

In Raymond et al. (2009b), we calculated β/β_crit by "observing" each system from a wide range of viewing angles, then using Equation (11) to calculate β/β_crit distributions. Although the β/β_crit value is sensitive to both the inclination between the planetary orbital plane and the line of sight, I, and the mutual inclination between the planets' orbits, Δi, we found that the effect of these angles was at the few percent level at most, i.e., negligible for our purposes. Additional tests have showed that a simple application of Equation (11) with no knowledge of viewing angle agrees remarkably well with the true and observed β/β_crit values as well as the cumulative distribution from a range of viewing angles. Thus, in this analysis we do not include the effects of viewing geometry.

Figure 16 compares the cumulative β/β_crit distributions for our unstable simulations with the observed two-planet extrasolar systems.¹⁵ As seen in Raymond al. (2009b), β/β_crit values from systems with equal-mass planets are far larger than those from systems with planetary mass ratios. This appears to be caused by the vastly larger number of encounters among equal-mass planets before a planet is destroyed.

Zoom In Zoom Out Reset image size

For the simulations with the most realistic mass distributions (i.e., the Mixed1 and Mixed2 cases), the effect of disks on the β/β_crit is small. Disks do not alter the prediction that scattering yields packed systems. Where there are differences they are in a surprising direction—after an instability, systems with planetesimal disks remain closer to the stability boundary (i.e., have lower β/β_crit values) than systems with no planetesimal disks. Since planetesimal disks damp eccentricities, they appear to leave the surviving planets in more dynamically compact configurations. This is true even for the low-mass systems such as N-N-N, for which planetesimal scattering almost always induces the migration of one planet out into disk. However, although the outer planet often migrates outward in N-N-N-type cases, the orbits of the inner planets are often compressed by scattering planetesimals outward. So, if we just consider the inner two planets in any system, they appear to end up closer together in the cases with planetesimal disks.

Without disks, the Mixed1, S-S-S and N-N-N simulations all provide good fits to the β/β_crit distribution—Table 3 contains p values from Kolmogorov–Smirnov (K–S) tests for each case. However, none of the sets provides a good match for β/β_crit ⩽ 1: the S-S-S and N-N-N cases are contaminated by systems which are themselves unstable on timescales longer than our 100 Myr integration time (these were removed by hand in Raymond et al. 2009b, but not here). With disks, the S-S-S simulations provide the best match to the known systems and the Mixed1 and J-J-J simulations also provide good fits. In addition, several sets match the distribution very well for β/β_crit ⩽ 1 with disks while they failed without disks. The cause of this increase in low-β/β_crit systems is the presence of MMRs. Systems in resonance tend to have β/β_crit < 1 (Barnes & Greenberg 2007), and these are preferentially set up in systems with planet masses comparable to M_J such as the Mixed1 and J-J-J simulations (see Section 4).

For MMR p + q:p, resonant arguments θ_i (also called "resonant angles") are of the form

$$\begin{equation} \theta _{1,2} = (p+q) \lambda _1 - p \lambda _2 -q \varpi _{1,2}, \end{equation} \tag{ 12 }$$

where λ are mean longitudes, ϖ are longitudes of pericenter, and subscripts 1 and 2 refer to the inner and outer planet, respectively (e.g., Murray & Dermott 2000). Resonant arguments measure the angle between the two planets at the conjunction point—if any argument librates rather than circulates, then the planets are in resonance. In fact, the bulk of resonant configurations are characterized by only one librating resonant argument (Michtchenko et al. 2008). Frequently, libration occurs around equilibrium angles of zero or 180°, but any angle can serve as the equilibrium. Different resonances have different numbers of resonant arguments, involving various permutations of the final terms in Equation (12). The order of a resonance is given simply by q. For example, the 3:1 MMR (q = 2, p = 1) is a second-order resonance that has three resonant arguments:

$$\begin{eqnarray} \theta _{1} &=& 3 \lambda _1 - \lambda _2 - 2 \varpi _{1}, \nonumber \\ \theta _{2} &=& 3 \lambda _1 - \lambda _2 - 2 \varpi _{2}, {\rm and } \nonumber \\ \theta _{3} &=& 3 \lambda _1 - \lambda _2 - (\varpi _{1}+\varpi _2). \end{eqnarray} \tag{ 13 }$$

In Raymond et al. (2008a), we found MMRs by targeting systems close to commensurabilities, then looking at individual resonant angles. In this paper, we take an automated approach. We devised a simple analysis script to analyze each of our simulations (with and without disks) and calculate all resonant angles for MMRs of up to fourth order. This includes the 2:1, 3:2, 3:1, 5:3, 5:2, and 4:1 MMRs. The script then determined which resonant angles were librating during the last 50 Myr of each simulation. This script is simple and reliable, but it misses the small fraction of cases that enter resonance between 50 and 100 Myr because it requires θ values to librate for the full interval (50–100 Myr). We impose a cutoff on the libration amplitude A of A ⩽ 150° to avoid false positives, which eliminates some additional systems that are just barely resonant.

In previous work, we found that a few to 10% of unstable systems could end in MMRs (Raymond et al. 2008a). The resonant fraction depends on the initial planetary mass distribution: MMRs were much more common in systems with radial mass gradients than in simulations that started with equal-mass planets. MMRs are populated by scattering simply because the density of resonant orbits is non-zero. In other words, the last close encounter leading to the destruction of one planet deposits the surviving planets onto orbits which have a chance of being resonant. This mechanism can populate high-order MMRs that are not populated by convergent migration in gaseous disks (Snellgrove et al. 2001; Lee & Peale 2002). The libration amplitudes of resonances populated by scattering are large, which also contrasts with the convergent migration scenario. Thus, as the population of resonant exoplanet systems increases, we may be able to differentiate MMRs caused by scattering versus convergent migration.

Figure 17 shows an example of an MMR caused by planet–planet scattering in a simulation that included a planetesimal disk. After 30,000 years, a close encounter with the outer (1.8 M_J) planet threw the lower-mass (0.4 M_J) planet outward into the planetesimal disk, where it remained on a stable but chaotic orbit for about 10 Myr. A final close encounter caused the lower-mass planet to be ejected and placed the surviving (1.9 and 1.8 M_J) planets in 5:2 MMR. In fact, the two surviving planets are in a special type of resonance called an apsidal corotation resonance in which all resonant angles librate and so does the apsidal alignment (Beaugé et al. 2003). This simulation was unusual because the entire planetesimal disk had been destroyed by the outer planet by the time of the last planetary scattering event such that the behavior was similar to simulations without disks.

Zoom In Zoom Out Reset image size

Most MMRs that occur in unstable systems with disks arise via a different mechanism. Usually, an instability occurs early in the simulation, removing one planet and placing the two surviving massive (M_p ∼ M_J) planets close to resonance. Subsequent interactions with the planetesimal disk act as a damping force that aligns the planetary orbits into MMR. In these cases, the outer surviving planet is massive enough to eject most of the planetesimals that it encounters, leading to some inward migration and a compression of the system. In fact, this mechanism has more in common with the generation of MMRs and MMR chains via planetesimal disk interactions (see Section 4.3) than MMRs from pure N-body scattering.

The frequency of MMRs from scattering is greatly reduced by the presence of the planetesimal disk (Table 4). In fact, only a small fraction of the MMRs with disks from Table 4 are even due to pure scattering—the majority are caused by disk damping as described above. The reason for the absence of MMRs in unstable systems with planetesimal disks is that the low-order MMRs populated by scattering occur preferentially in systems with high-mass ratios (for example, almost 10% of the unstable N-S-J simulations yielded 2:1 MMRs), meaning that one planet in the system was relatively low mass (usually ∼M_S) and subject to radial migration by planetesimal scattering after the instability placed the planets in resonance.

Thus, scattering appears to be a viable mechanism to populate resonances but its efficiency is reduced in regions which are strongly affected by planetesimal disks. Given that the observed exoplanets do not yet show an obvious signature of planetesimal scattering (see Section 5), some of the known systems could indeed have been populated by this mechanism (Raymond et al. 2008a). However, as we show below, interactions with planetesimal disks in systems that do not experience close encounters and scattering events are vastly more efficient at generating resonant systems.

We now consider the generation of MMRs in stable systems that do not experience close encounters between the planets. An empirical argument can be made that such systems are rare amongst the progenitors of currently observed systems, since the good agreement with the observed eccentricity distribution requires that a large fraction of systems are unstable. At larger radii, however, stable systems could be present and perhaps common. In stable systems with relatively massive planets (M_p ≳ M_J), the planetesimal disk can act as a damping force and induce planets to align into MMRs. This mechanism is not efficient for low-mass planets because the planetesimal disk has enough angular momentum to cause excessive radial migration and drive the planets apart. However, gentle migration for M_p ∼ M_J is very effective at forming resonant chains reminiscent of the Galilean satellites of Jupiter (Paper I—Raymond et al. 2009a).

Figure 18 (left panel) shows one example of a 4:2:1 resonant chain for a Mixed1 simulation with three ∼ Jupiter-mass planets. The θ₁ resonant argument is librating with a small amplitude for both the inner and outer pairs of planets, and for the inner pair θ₂ begins to librate after ∼20 Myr and its libration amplitude decreases over the next 10–20 Myr. The damping effect of the planetesimals is clearly seen as the system moves deeper into the resonance with time, i.e., the libration amplitude of the two resonant arguments decreases. The right panel of Figure 18 shows the same system with no planetesimal disk. It is clear that this system started with both pairs of planets close to the 2:1 MMR (see Figure 2) but without the disk the planets do not evolve into the resonance.

Zoom In Zoom Out Reset image size

Table 5 shows the number of resonances and resonant chains that were set up in each set of stable simulations with planetesimal disks, and that number as a fraction of the total number of stable simulations. For systems with planet masses M_S ≲ M_p ≲ 3 M_J, between roughly 50% and 75% of all stable systems ended up with at least one MMR. In addition, up to 40% of systems were in resonant chains involving all three planets. Recall that our definition of resonance requires at least one resonant argument to librate with an amplitude A ⩽ 150°. The second row for each set of simulations shows the number and fraction of systems that were deep in resonance, with A ⩽ 60°. Roughly half of the MMRs were deep in resonance, independent of their configuration in the system: both pairs of resonant planets were deep in resonance for about one quarter of the resonant chains.

The vast majority of resonant systems are in the 2:1 MMR, although for lower masses the 3:2 MMR is also common. This is simply because of the location of strong resonances in initial planetary separations versus mass (Figure 2). More (less) massive planets start closer to the 2:1 (3:2) MMR, and so given the limited amount of radial migration, tend to end up in that same resonance. The second-order MMRs were comprised almost exclusively of the 5:3 MMRs, although those were much rarer than first-order MMRs. Depending on the case, about half or slightly less of the resonances were deep, meaning that the libration amplitude was 60° or less. Systems that were extremely deep in the resonance, however, were rare: indeed only a single case with A < 10° was found in the entire ensemble of simulations. This is of interest since dynamical modeling of the GJ 876 system shows that one of the resonant angles librates with an amplitude that is, at most, θ = 7° ± 2° (Laughlin et al. 2005). Taken at face value, one concludes that capture into resonance due to gas disk migration is not only the favored explanation for the origin of the GJ 876 system itself, but that gas rather than planetesimal effects would also be needed to explain hypothetical analogs of GJ 876 that might be discovered very deep in resonance at large orbital radii. We caution, however, that we have not demonstrated that the distribution of libration amplitudes is independent of numerical details of our simulations, and in particular that it is independent of the mass of particles used to represent the disk. Further simulations would be needed to study the predicted distribution of A in more detail.

What causes the resonances in these systems? Figure 19 shows the semimajor axis shifts of the planets in the stable Mixed1 and Mixed2 systems: positive values indicate convergence and negative indicate divergence. There is a much more limited range of parameter space for convergence than divergence, since for massive planets the stability limit is located in the vicinity of the 3:2 MMR. Planet pairs were limited to converging ≲0.1 AU, but divergence was unlimited, at least for the outer pair of planets (too much divergence of the inner pair could impinge on the outer planet).

Zoom In Zoom Out Reset image size

Table 5 shows that MMRs generally occur at a higher frequency in systems that converge rather than diverge. However, for most giant planet configurations, systems for which all three planets converge are rare. The exception are the systems with positive mass gradients (S-J-3J and N-S-J), for which almost half the stable systems converged. This is because they contain a very massive outer planet capable of quickly ejecting planetesimals rather than scattering them inward, thereby migrating slightly inward. In contrast, systems with negative mass gradients (3J-J-S and J-S-N) favor divergence because the innermost planet is the ejector and the outer two preferentially scatter planetesimals inward and migrate outward. For massive equal-mass planets (J-J-J and 3J-3J-3J), it is often the case that the inner two planets converge while the outer two diverge. This is caused simply by the outer planet scattering planetesimals in to the middle planet. With higher relative velocities the middle planet can more easily eject the planetesimals, causing it to migrate inward, away from the outer planet but toward the inner one.

Among the simulations with radial mass gradients, MMRs are more common in systems with a more massive outer planet (those with positive mass gradients, S-J-3J and N-S-J). This is because the massive outer planet confines the planetary system, and the planetesimal scattering causes the system to converge, as a whole or for one pair of planets, much more frequently than for a lower-mass outer planet. Indeed, the frequency of MMRs as a whole correlates with the fraction of systems which converge rather than diverge (see Tables 5 and 6). Nonetheless, many MMRs are created in divergent systems, simply because planets were initially placed near strong MMRs (Figure 2).

For the Mixed1 simulations, we generated planet masses by sampling the observed planetary mass function randomly and independently, using mass limits of M_Sat < M < 3 M_J. Since the eccentricity distribution that results from these initial conditions matches the data closely (Figure 6)—irrespective of whether we include disks or not—it might appear as if there is no further problem to solve. This ignores the fact that the observed exoplanet eccentricity distribution is also mass dependent: lower-mass planets have lower eccentricities than higher-mass planets (Ribas & Miralda-Escudé 2007; Ford & Rasio 2008; Wright et al. 2009). The evidence for this correlation appears to be strong. If we divide the sample of known extrasolar planets into low-mass (M_psin i < M_J) and high-mass (M_psin i > M_J) bins, we find a K–S probability p of only 0.002 that they are drawn from the same distribution.

Assuming that this trend is real we have a clear problem: scattering does not predict that low-mass planets should be less eccentric than high-mass ones. Figure 20 shows eccentricity distributions for the Mixed1 and Mixed2 simulations divided into the same low-mass (M_p < M_J) and high-mass categories for cases with and without disks. The predicted trend is opposite to that observed. Indeed, once the sample has been split into mass bins it becomes clear that the good agreement of the overall distribution is somewhat fortuitous. It occurs because—for the Mixed1 simulations—the low-mass simulations match the observed high-mass curve and vice versa. Changing the lower-mass limit for the mass function sampling to 10 M_⊕ only makes matters worse. For the simulations with no disks, the more massive planets are even less eccentric than for the Mixed1 simulations simply because the typical mass ratio in a given scattering event is larger for the Mixed2 simulations. It is interesting to note that for the Mixed2 simulations with disks the high- and low-mass planets have virtually identical eccentricity distributions, meaning that the tendency for low-mass planets to acquire large eccentricities during scattering is almost exactly balanced by the eccentricity damping of the disk.

Zoom In Zoom Out Reset image size

The above problem raises a simple question, is it possible to simultaneously match the observed constraints on (1) the total eccentricity distribution, (2) the mass–eccentricity correlation, and (3) the mass function using a scattering model? Although our data are not able to provide a rigorous answer we contend that solutions are possible, either with or without appealing to the existence of planetesimal disks. In the absence of disks, the key assumption that we must drop is that planet masses in the same system are chosen independently from the overall mass function. If, instead, we allow for correlations between planet masses, the extra degree of freedom allows us to match what is seen observationally.

For a disk-free solution, we start with the equal-mass simulations and choose a weighted combination of the 3J-3J-3J, J-J-J, and S-S-S simulations that matches the observed mass distribution (i.e., that contains correspondingly more low-mass than high-mass systems). By construction this matches the mass distribution, and does a reasonable job at matching the total eccentricity distribution (the K–S p value is 0.09). However, this combination does not match the mass–eccentricity distribution; both the high-mass and low-mass planets have large eccentricities, and in fact, both provide a match to the M_p > M_J exoplanets. The high-mass planets have slightly higher eccentricities than the low-mass planets (see Figure 6) but the difference is small because the highest-eccentricity systems (3J-3J-3J) are very weakly weighted to match the observed mass distribution. Note that we count Jupiter-mass planets in both the high-mass and low-mass bins.

We now mix in systems for which low-mass planets have lower eccentricities than high-mass planets. In the absence of disks, this can be done if we assume that high-mass and low-mass planetary systems are fundamentally different. Specifically, we assume that high-mass systems form roughly equal-mass planets yielding large eccentricities, whereas low-mass systems form a much wider variety of systems and tend to have lower eccentricities. To show that this works, at least approximately, we assume that the J-S-N and N-S-J simulations yield an eccentricity distribution that is representative of the scattering among unequal-mass low-mass planets (with masses up to, but not beyond, M_J). We then add J-S-N and N-S-J simulations to the total and low-mass eccentricity distributions, but not to the high-mass distribution. We gave the J-S-N/N-S-J systems the same abundance as the S-S-S systems such that the low-mass systems are equally divided between equal- and unequal-mass cases (except for a small contribution from the J-J-J simulations).

The upshot of this cookery is shown in Figure 21. The blended distribution, which is roughly consistent with the mass function, provides a decent match to both the total eccentricity distribution and to the mass–eccentricity distributions. The K–S test p values are 0.25 (for the total distribution), 0.09 (for the low-mass planets), and 0.20 (for the high-mass planets). It is evident that matching the low-mass planets is by no means easy, and in fact our "solution" is close to being statistically unacceptable. Nonetheless, it is clearly possible to reverse the natural tendency of scattering to generate larger e for low-mass planets by strongly correlating the masses of planets. We further believe—though we do not have the data to demonstrate it—that the magnitude of the observed mass–eccentricity correlation lies within the range that could be accurately reproduced by fine tuning the mix of initial conditions used.

Zoom In Zoom Out Reset image size

An entirely different solution to the mass–eccentricity problem invokes the damping effects of planetesimal disks. Disks referentially damp the eccentricities of the smallest planets (Figure 8), but, as we have noted, the effect correlates most strongly with the total system mass rather than that of individual planets. We therefore consider the equal-mass cases, and (as before) blend the 3J-3J-3J, J-J-J, and S-S-S simulations such as to approximate a three-bin version of the mass function. We assume that a fraction F = 1/2 of the systems are affected by disks, while the rest are not. This combination provides reasonable matches for both the total eccentricity distribution (p = 0.08) as well as the low-mass (p = 0.02) and high-mass (p = 0.44) planets (red curves in Figure 21). The match for the low-mass planets can be improved by slightly increasing F, at the expense of the fit to the overall distribution.

Do the two weighted distributions match other observations? Using the same combination of different sets of simulations, we generated weighted distributions of β/β_crit to compare with observations, seen in the third panel of Figure 21. Remarkably, both combinations of simulations provide excellent fits to the distribution of separations in two-planet systems. The weighted case with disks was a better fit to the β/β_crit distribution (K–S p values of 0.20 with no disks and 0.72 with disks), largely because the cases with disks were able to fit the distribution of values at β/β_crit ⩽ 1 which are dominated by resonant systems. Thus, the β/β_crit distribution cannot differentiate between the two weighted combinations of simulations we have constructed. However, it is interesting that both weighted distributions can match all of the currently well-observed constraints.

What can we conclude from this proof-of-principle exercise? First, correlations between planet masses in the same system can qualitatively alter some of the predictions of scattering models. Although we have here allowed ourselves complete freedom to impose whatever correlations we need, our disk-free solution in fact involves physically conceivable variations in the outcome of planet formation in different disks. A massive disk of gas and solids would be expected to form several cores that go on to accrete substantial envelopes, such that the outcome is usually a system with several very massive planets. A less massive disk, on the other hand, may only form one relatively low-mass gas giant together with a number of low-mass core-dominated planets. Second, the effects of planet–planet scattering with mass correlations, and the effects of planetesimals, can be almost indistinguishable, at least for the distributions we have considered. This should be borne in mind, especially as observations probe larger orbital radii where it is more probable that planetesimal disks play a dynamical role.

Notwithstanding the formal degeneracy of our solutions, the disk-free solution is favored as the origin of currently observed eccentricities. The reason is simply that the current sample includes many planets at radii that are simply too small to have experienced damping from plausible planetesimal disks (note that our inner disk radius of 10 AU is already smaller than is often assumed). The match is not perfect but this is not surprising, since the comparison is not entirely fair. Both the observed and predicted eccentricities depend upon the Safronov number, defined as the ratio of the escape speed from a planet's surface to the escape speed from the planetary system at that location,

$$\begin{equation} \Theta ^2 = \left(\frac{M_p}{M_\star }\right) \left(\frac{a}{R}\right), \end{equation} \tag{ 14 }$$

where M_p and M_⋆ are the planetary and stellar mass, a is the orbital distance and R is the planetary radius (Ford & Rasio 2008). Our simulations probe the region at somewhat larger a than is currently observed, such that the Θ values are correspondingly higher. Only at higher planetary masses are the Θ values large enough that the variations in Θ between our simulations and the observed planets are too small to affect the eccentricities, and it is reassuring that our weighted fits match the high-mass planets quite well. Another point to bear in mind is that although the selection function of radial velocity surveys is approximately independent of eccentricity (Cumming et al. 2008), there are almost certainly systematic errors to add to the purely statistical errors that we have used (e.g., Shen & Turner 2008).

We do not expect the degeneracy between disks and mass correlations to persist to arbitrarily large orbital radii. Disk damping becomes increasingly more important as the semimajor axis increases, and as the planet mass decreases. For very low-mass planets, the resulting reduction in eccentricity is greater than could be explained by any plausible blend of simulations with different conditional mass functions. An ongoing trend to lower eccentricities for low-mass planets, as the orbital radius increases, would thus be a strong sign of dynamical influence from planetesimal disks.