PLANET–PLANET SCATTERING IN PLANETESIMAL DISKS

Different dynamical mechanisms are commonly invoked to explain the architecture of the outer Solar System and extrasolar planetary systems. In the Solar System, scattering of small bodies ("planetesimals") by the ice giants (Fernandez & Ip 1984; Ida et al. 2000; Kirsh et al. 2009) is thought to drive outward planetary migration and concomitant capture of Pluto and other Kuiper Belt Objects into resonance (Malhotra 1995; Murray-Clay & Chiang 2005). The effects of planetesimal scattering on the gas giants are smaller but still significant, for example, in the "Nice model" (Tsiganis et al. 2005; Gomes et al. 2005), where a divergent resonance crossing between Jupiter and Saturn triggers the late heavy bombardment. The presence of small bodies around other stars can be inferred from observations of debris disks (e.g., Wyatt 2008), but as yet there is no evidence for a dynamical role of planetesimals in known extrasolar planetary systems. At radii where tidal effects are negligible (roughly a ≳ 0.1 AU) the eccentricity distribution of extrasolar planets matches relatively simple models of gravitational scattering among a system of two or more massive planets that typically include neither planetesimals nor residual gas (Chatterjee et al. 2008; Jurić & Tremaine 2008).

The success of pure planet–planet scattering models does not imply that other dynamical processes can be ignored. The observed distribution of semimajor axes of extrasolar planets at small orbital radii requires the existence of an additional dissipative process (Adams & Laughlin 2003), most probably gas disk migration (Lin & Papaloizou 1986), which will itself affect planetary eccentricity (Moorhead & Adams 2005). At larger orbital radii simple arguments suggest that a dynamically significant external reservoir of planetesimals ought to be a common feature of young planetary systems. The formation of giant planets becomes increasingly difficult at large radii (Pollack et al. 1996; Kokubo & Ida 2002), and hence it is probable that disks of leftover debris surround the zone of giant planet formation in most young systems. The typical masses of planetesimal disks are unknown, but values of 30–50 M_⊕ that are comparable to those inferred for the early outer Solar System are plausibly typical, since they are consistent with disk masses estimated from astronomical observations of the youngest stars (Andrews & Williams 2005). The dynamical effect of such disks on currently observed extrasolar planetary systems would be small, since radial velocity surveys preferentially detect planets that are either massive (and hence largely immune to influence from planetesimal disks) or orbit at very small radii where the mass of leftover debris is negligible.

Pooling knowledge from the Solar System and extrasolar planetary systems motivates consideration of a model in which planet formation typically yields a marginally unstable system of massive planets in dynamical contact with both a residual gas disk and an exterior planetesimal disk. In this Letter, we ignore the gas disk and study the subsequent evolution under the combined action of planet–planet and planetesimal scattering. We do not attempt to model the full distribution of extrasolar planetary properties (which would require the inclusion of hydrodynamic effects), but rather focus on how planetesimal disks affect the final eccentricity of extrasolar planets at moderately large orbital radii.

We assume that the gas-dominated epoch of planet formation is sufficiently distinct from the subsequent phase of planet–planet and planetesimal scattering that it makes sense to study the latter with pure N-body simulations. We focus on two large ensembles of runs. The highmass set comprises 1000 integrations of three planet systems in which the masses of the planets are chosen randomly in the range M_Sat < M_p < 3M_J, with a distribution,

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

which matches that observed (Marcy et al. 2008). The observed distribution is derived from an incomplete sample that represents (in the context of our model) the distribution after scattering, but these subtleties do not matter for our purposes. The lowmass set is identical except that we sample a wider swath of the mass function between 10M_⊕ and 3M_J. The planets are initially placed in a marginally unstable configuration defined by circular, nearly coplanar orbits with a separation of 4–5 r_h,m, where the mutual Hill radius,

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

Here, a₁ and a₂ are the planets' semimajor axes, M₁ and M₂ are their masses, and M_⋆ is the stellar mass. With this spacing the instability timescale is relatively long (Chambers et al. 1996; Chatterjee et al. 2008; the median timescale before the first planet–planet encounter was 0.3 Myr for the highmass integrations without disks). Our initial conditions are only a small subset of the architectures predicted from giant planet formation models (Thommes et al. 2008b; Mordasini et al. 2009), though broadly consistent with scenarios in which giant planets are captured into mean-motion resonances during the late stages of gas disk evolution (Thommes et al. 2008a) prior to being removed from resonance by turbulent perturbations (Adams et al. 2008). Each integration is repeated twice, once with just the three planets⁵ and once with an external planetesimal disk whose inner radius of a_in = 10 AU is 2 Hill radii beyond the orbit of the outermost planet.⁶ The inner edge of the disk lies within the radius where a test particle in the restricted 3-body problem would be stable, so the disk is in immediate dynamical contact with the outer planet. The planetesimal disk is represented by 1000 bodies distributed between 10 and 20 AU with a Σ_disk ∝ r⁻¹ surface density profile and a total mass of 50M_⊕.

We integrate these systems using the MERCURY code (Chambers 1999) for 100 Myr. The integrator uses the symplectic Wisdom–Holman mapping (Wisdom & Holman 1991) for well-separated bodies, and the Bulirsch–Stoer method when objects are within N mutual Hill radii, where N = 3 for our case. Planets were removed if their orbital distances were smaller than 0.1 AU ("hit Sun") or exceeded 100 AU ("ejection"). Collisions were treated as inelastic mergers conserving linear momentum.

A large ensemble of simulations includes some cases that are much harder to integrate accurately than the majority. To make the best use of our computational resources we adopted a default time step (20 days) that results in accurate integrations (as measured by the fractional orbital energy conservation dE/E) for the typical case. We then identified those runs (about 10%) in which energy was not adequately conserved and re-ran them with a smaller time step. For runs without disks we re-ran cases with dE/E > 10⁻⁴ with a time step of 5 days, while for the runs with disks we re-computed cases with dE/E > 5 × 10⁻⁴ with a time step of 10 days. A small number of the re-run simulations (typically 15–35) still did not meet our energy criterion and were discarded.

In the absence of planetesimal disks our model planetary systems are typically unstable on Myr timescales. There are also systems that are stable over the 100 Myr duration of our runs. In our initial analysis, we assume that the typical outcome of giant planet formation is a system that, in the absence of a disk, would be unstable. We therefore analyze the subset of disk-less simulations that are unstable, and compare the results to the matched set of simulations that include disks. This is not a perfect one-to-one comparison, since the chaotic nature of the evolution means that disk-less planetary systems can display different instability timescales in the presence of even negligible perturbations. Nonetheless we do observe statistical differences between the evolution of systems (with disks) that correlate with the stability of the disk-free systems, and hence it makes sense as a first approximation to separately consider the results for stable and unstable cases.

The results of our disk-less simulations agree with prior studies (Chatterjee et al. 2008; Jurić & Tremaine 2008; Ford & Rasio 2008). Scattering from initially unstable initial conditions frequently leads to the loss of one or more planets via ejection or collisions and sets up a broad eccentricity distribution (Rasio & Ford 1996; Weidenschilling & Marzari 1996; Lin & Ida 1997). Scattering among equal-mass planets produces larger eccentricities than scattering of unequal-mass planets (Ford et al. 2003; Raymond et al. 2008). Figure 1 compares the final eccentricities obtained from the unstable highmass simulations and the observed distribution (Butler et al. 2006; Schneider 2009) of extrasolar planets. They are in good quantitative agreement. The eccentricity distribution from the unstable lowmass simulations without disks is shifted toward lower values (Raymond et al. 2008; Ford & Rasio 2008), and fits the observed distribution of extrasolar planets with M_p < M_J (Wright et al. 2009). Our model therefore exhibits evolution that is consistent with current observations of extrasolar planetary systems, which as we noted previously are mostly of systems at such small radii that planetesimal disks are dynamically unimportant.

Zoom In Zoom Out Reset image size

At larger radii we expect that both disks and planet–planet scattering will play a dynamical role. A wide range of outcomes is then possible. Exchange of energy and angular momentum between the planets and the planetesimal disk leads to planetary migration (Fernandez & Ip 1984; Murray et al. 1998; Ida et al. 2000; Kirsh et al. 2009), which can be either stabilizing or destabilizing. A low-mass planet adjacent to the disk scatters planetesimals inward, resulting in divergent migration that is often stabilizing unless resonance crossing excites eccentricity to the point of triggering instability. Alternatively, an outer massive planet interacting with the disk directly ejects planetesimals and migrates inward, compressing the system and leading to instability or resonant capture. An equally important effect is that the disk can act to recircularize the orbits of scattered planets after dynamical instabilities (Thommes et al. 1999; Ford & Chiang 2007). To illustrate how significant recircularization can be we ran a small additional set of idealized experiments in which a single planet with mass M_p on an orbit with a = 10 AU and e = 0.5 begins to interact with our initial planetesimal disk. For 10⁴–10⁵ years (longer for smaller M_p), e is damped roughly exponentially with a damping timescale t_e, defined via,

$$\begin{equation} \frac{1}{e} \frac{{d}e}{{d}t} \equiv - \frac{1}{t_e} \end{equation} \tag{ 3 }$$

of t_e ≈ 5 × 10⁴ years, independent of M_p. The subsequent evolution was highly mass-dependent: for low-mass planets, e continued to decrease on much longer timescales (0.36, 0.63, and 4.6 Myr to reach e ≲ 0.1 for M_p = 10M_⊕, 30M_⊕, and M_S, respectively). Massive planets disrupted the disk, halting dynamical friction. The total decrease in e for M_p ≳ M_J was 0.15 or less, corresponding to an increase of 1.5 AU or less in perihelion distance.

Figure 2 illustrates the diversity of outcomes from our simulations that include planetesimal disks. We split our simulations into three mass bins (the Solar System's giant planets fall into the middle bin) and three stability categories. In "stable" systems there are no close encounters between planets and no large-scale change in system architecture (the ordering of the planets is preserved and all planets survive). "Moderately stable" systems experience substantial perturbations—which may be due to resonance crossing in high-mass systems or close encounters in low-mass systems—that are nonetheless insufficient to alter the architecture. "Unstable" systems undergo close encounters leading to architectural change.

Zoom In Zoom Out Reset image size

Subsets of our runs show dynamics analogous to that studied for the Solar System and for extrasolar planetary systems. At high masses planetesimal disks stabilize about 30% of cases but planet–planet scattering leading to the loss of one or more planets is still common. Quantitatively, the median eccentricity is reduced (Figure 1) but many highly eccentric systems remain. As planet masses decrease the dynamical importance of planetesimals grows. For low-mass systems, the masses of the planets and the planetesimal disk are comparable and planetesimal scattering inevitably leads to migration. Divergent crossing of mean-motion resonances (one example of which is shown in the center panel of Figure 2) can result in abrupt changes to planetary semimajor axis and eccentricity that qualitatively resemble those seen in the Nice model (Tsiganis et al. 2005; Gomes et al. 2005). We also see behavior that resembles an alternative Solar System model in which Uranus and Neptune formed in the Jupiter–Saturn region and were scattered outward (Thommes et al. 1999; top left panel). At the lowest masses even highly unstable systems rarely destroy any planets because recircularization of scattered planets is efficient (top center panel), though re-ordering of planets is common. In summary, dynamic characteristic of the outer Solar System is common among low- to medium-mass planetary systems.

The main prediction of our model is the statistical distribution of planetary eccentricity as a function of planet mass. Figure 3 shows the distribution of the eccentricity of the innermost surviving planet as a function of the total mass in surviving planets. In the absence of disks we observe similar behavior across all system masses—the shift to smaller eccentricities for the lowmass runs, seen in Figure 1, is not visually apparent. When disks are included, the eccentricity distribution divides into two distinct regimes: a low-mass regime in which planetesimal dynamics dominates to yield low eccentricities and a high-mass planet scattering dominated regime where planetesimals play a minor role. For our specific parameters (inner edge disk edge at 10 AU, and a disk mass of 50 M_⊕ assumed to be typical for a stellar metallicity Z = Z_☉) the transition between these regimes occurs for system masses M_crit ≈ 1 M_J. We predict that systems whose giant planets orbit between 5 and 10 AU, and which have a total mass below 1 M_J, will typically have low eccentricity orbits. This critical mass should scale roughly linearly with the stellar metallicity, as we expect the initial planetesimal disk mass to be proportional to Z. We expect the same qualitative behavior even if the zone of giant planet formation extends to larger radii (Goldreich et al. 2004), though in this case the low eccentricity regime would only be observable further out.

Zoom In Zoom Out Reset image size

Although the eccentricity distribution of extrasolar planets is consistent with the hypothesis that all newly formed multiple systems are unstable in the absence of disks, this conclusion may also be biased by selection effects. Most known extrasolar planets probably suffered significant gas disk migration prior to scattering (Lin et al. 1996; Trilling et al. 1998; Bodenheimer et al. 2000), so the high incidence of instability may be a consequence of migration rather than formation. With this in mind we have separately analyzed those (previously excluded) systems that were stable in the absence of disks to see what impact disks have on them. As expected, low eccentricity outcomes predominate. Resonances are also common: about 70% of all stable highmass simulations and 1/3 of stable lowmass simulations include at least one pair of planets in the 3:2 or, more often, the 2:1 mean motion resonance. This is a much higher probability of resonance capture than occurs for pure planet–planet scattering without disks (Raymond et al. 2008), and it also exceeds the fraction of resonant systems that are expected to survive the gaseous disk phase in the presence of turbulence (Adams et al. 2008). Most surprisingly within the highmass set a substantial fraction (about 1/3) of stable systems become locked into mean-motion resonances that involve all three of the planets—analogs of the Laplace resonance among Jupiter's Galilean satellites. Chains of resonances⁷ arise preferentially in higher-mass systems and in systems where planetesimal-driven migration causes compression rather than divergent migration. Detection of high-mass planets at the relevant radii (between 5 and 10 AU) should soon be possible via astrometric or direct imaging techniques, and observation of resonant chains would be consistent with our model. Intriguingly, the recently discovered triple planet system HR 8799 may be in a 4:2:1 resonant chain (Morois et al. 2008; Fabrycky & Murray-Clay 2009). Determining whether capture into resonance was initiated by a gas or planetesimal disk may be possible via detailed comparison of the outcome of resonant capture in planetesimal (Murray-Clay & Chiang 2005) versus gas disks (Lee & Peale 2002; Adams et al. 2008).

Circumstantial evidence suggests that many observed properties of the outer Solar System (Malhotra 1995) and of extrasolar planetary systems (Rasio & Ford 1996; Weidenschilling & Marzari 1996; Lin & Ida 1997) may be attributable to the dynamical effects of planetesimal scattering and planet–planet scattering. Here, we have studied the predicted architecture of planetary systems that results from the joint action of both mechanisms. We have argued that this regime will be relevant once lower-mass extrasolar planets are discovered at larger orbital radii than those currently known. Generically we predict that a transition to "Solar-System-like" architectures, characterized by near-circular orbits and relatively stable planetary separations, will be observed once surveys detect planets in the regime where planetesimal disks play a dynamical role. Our simulations suggest that the transition is a surprisingly sharp function of total planetary system mass, and that it occurs for system masses a factor of several larger than the initial planetesimal disk mass.

Our initial conditions do not sample anything approaching the full range of initial planetary system architectures. We believe that the existence of a transition between typically eccentric and near-circular orbits is a general feature of joint models of planet–planet and planetesimal scattering, but the transition mass and minimum orbital radii at which planetesimal effects become manifest is of course a function of the poorly known masses and radial extent of planetesimal disks. Our results suggest that the transition might be seen for sub-Jovian mass planets at orbital radii of 5–10 AU, but the transition would be pushed to greater orbital radii if giant planet formation consumes planetesimals across a wider extent of the disk. We also find that the final system architecture varies substantially depending on the initial separation of the planets. In particular, if planet formation yields a mixture of massive systems in initially stable orbits, interaction with planetesimals drives a large fraction of systems into resonance. Whether such systems exist should be testable in the near future.

We thank Google for the large amount of computer time needed for these simulations. S.N.R. acknowledges support from NASA's Astrobiology Institute through the Virtual Planetary Laboratory lead team, and from NASA's Origins of Solar Systems program (NNX09AB84G). P.J.A. acknowledges support from the NSF (AST-0807471), NASA's Origins of Solar Systems program (NNX09AB90G), and NASA's Astrophysics Theory program (NNX07AH08G).