EDP Sciences
Free Access
Issue
A&A
Volume 541, May 2012
Article Number A11
Number of page(s) 25
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/201117049
Published online 19 April 2012

© ESO, 2012

1. Introduction

The solar system’s distinctive architecture, in which rocky terrestrial planets lie interior to gas and ice giants, with the Kuiper Belt of smaller bodies beyond, is not unexpected. The mass in protoplanetary feeding zones increases with orbital distance, but the resulting tendency toward the formation of larger planets further out is eventually frustrated both by the lengthening time scale for accretion (e.g., Lissauer 1993; Kokubo & Ida 2002), and by the increased ability of planetary cores to scatter planetesimals inward (Levison & Stewart 2001). The competition between these effects plausibly leads to the relatively slow (~100 Myr) assembly of a handful of terrestrial planets inside a few AU, the faster growth of several planetary cores with M ≳ 5  M inside  ~10 AU, and the persistence of a belt of unconsolidated debris further out. Simple arguments of this kind fail to establish how often planetary cores grow fast enough to admit the formation of fully fledged gas giants, but empirical estimates based on extrapolations of radial velocity and microlensing surveys suggest that gas giant formation is common (Cumming et al. 2008; Gould et al. 2010). Debris disks (Wyatt 2008) are also observed around a significant fraction of young stars – despite the existence of both dynamical and collisional processes that can destroy them on a time scale short compared to the main sequence lifetime of Solar-type stars – but the abundance of Earth-mass planets remains to be measured.

For the solar system, we have access to a unique array of observational constraints. Even with these advantages the exact nature of the interactions between the giant planets, the terrestrial planets, and the Kuiper belt remain under debate. In the inner solar system, at a minimum secular resonances with the giant planets would have influenced terrestrial planet formation (Nagasawa et al. 2005; Raymond et al. 2009c). There could, however, have been stronger effects. Gas driven migration of the giant planets could have brought them (temporarily) closer to the Sun (Walsh et al. 2011), directly reducing the supply of raw material in the Mars region and preventing the growth of a larger planet (Hansen 2009). In the outer solar system, early studies focused on the dynamics of Neptune, which is of low enough mass that inward scattering of planetesimals can drive substantial outward orbital migration (Fernandez & Ip 1984). The migration can deplete the mass in the Kuiper Belt and result in the resonant capture of Pluto and other bodies by Neptune (Malhotra 1993). Subsequent work introduced the idea of larger-scale dynamical instability among the solar system’s giant planets (Thommes et al. 1999). In the most-developed models, early outer solar system evolution is characterized by a combination of planetesimal migration, close encounters between planets, and resonant interactions (Tsiganis et al. 2005; Levison et al. 2011).

Theoretically, attempts to construct equally detailed models for extrasolar planetary system evolution are hampered by uncertainties in the distribution of the initial disk conditions, and by our poor knowledge of the evolution of gas disks (Armitage 2011) and formation mechanism for planetesimals (Chiang & Youdin 2010). Several observed properties of extrasolar planetary systems, however, including the existence of hot Jupiters (whose orbits are sometimes misaligned with respect to the stellar spin axis) and the prevalance of eccentric orbits, favor scenarios in which large-scale orbital evolution of giant planets is the norm (Winn et al. 2010; Triaud et al. 2010; Schlaufman 2010). It is therefore of interest to determine the diversity of outcomes that could arise given initial conditions and dynamical processes similar to those of the solar system, and to examine how those outcomes depend upon parameters such as the masses of the giant planets, and the properties of primordial planetesimal belts. Doing so is the goal of the present paper. We are particularly interested in studying how terrestrial planets form, and debris disks evolve, in the presence of dynamically active giant planet systems. In an earlier paper (“Paper I” Raymond et al. 2011) we showed that if giant planets form in or near dynamically unstable configurations, there are striking correlations between the nature of the terrestrial planets that form, and the properties of outer debris disks whose ongoing collisional evolution can be observed out to ages of several Gyr (e.g. Wyatt 2008; Krivov 2010). The dynamically calm conditions that favor the formation of massive terrestrial planet systems also result in long-lived outer debris disks, that remain bright in cold dust emission (e.g. at λ ~ 70   μm) to late times. In systems that suffer more dramatic dynamical evolution, we identified a channel for the formation of unusual terrestrial planet systems in which a single planet exhibits large oscillations in eccentricity and inclination due to secular coupling to a scattered giant. Here, we consider a broader range of models within the same qualitative class, and study how robust our earlier conclusions are to changes in the poorly-constrained model parameters.

Our paper is structured as follows. In Sect. 2 an outline of our methods is presented (the reader is referred to Paper I for more details and numerical tests). In subsequent sections we present results of different sets of simulations to test the effect of: the giant planet mass and mass distribution (Sect. 3), the width and mass distribution of the outer planetesimal disk (Sect. 4), and the presence of gas during giant planet instabilities (Sect. 5). In Sect. 6, we discuss the implications of our models for debris disks and terrestrial planet systems. We conclude in Sect. 7.

2. Methods

The initial conditions for our simulations assume that the location and mass of the rocky material in the terrestrial planet zone, and the mass of planetesimals in the outer disk, are fixed at values similar to those employed in solar system models. We assume that the masses of the giant planets, on the other hand, have a broad dispersion (extending up to masses above those realized in the solar system) that is uncorrelated with the mass in either the terrestrial planet region or in the outer planetesimal disk. All of our simulations contain three radially-segregated components orbiting a solar-mass star:

  • 1.

    The building blocks of terrestrial planets:9   M in 50 planetary embryos and 500 planetesimals from 0.5 to 4 AU, with equal mass in each component and a radial surface density profile Σ ~ r-1. The initial eccentricities were chosen at random from 0–0.02 and the initial inclinations from 0−0.5°.

  • 2.

    Three giant planets at Jupiter-Saturn distances: the innermost planet is placed at 5.2 AU and the two others are spaced outward by 4–5 mutual Hill radii. We adopt three-planet initial conditions because this is the simplest plausible configuration that evolves dynamically to match the measured eccentricity distribution of massive extrasolar planets (Chatterjee et al. 2008). The planets were placed on initially circular orbits with randomly-chosen inclinations of 0−1°.

  • 3.

    An outer disk of planetesimals thought to be analogous to the primitive Kuiper belt. This belt consists of 1000 planetesimals with a total mass of 50 M. The belt starts 4 Hill radii beyond the outermost giant planet and extends radially for 10 AU, and also follows an r-1 radial surface density profile. The initial eccentricities were chosen at random from 0–0.01 and the initial inclinations from 0−0.5°.

Adopting these initial conditions amounts to making implicit assumptions about the typical outcome of planet formation. First, our terrestrial, giant planet, and outer disk zones are located such that they are in immediate dynamical contact with each other. This is reasonable only if planetesimal formation results in a smooth, gap-less distribution of bodies in 0.5     AU < a ≲ 20     AU, and if giant planet migration is limited. Substantial giant planet migration, of the kind envisaged in models by Masset & Snellgrove (2001), Walsh et al. (2011) and Pierens & Raymond (2011), could create dynamical separation between the giant and terrestrial planets prior to the gas-less phase of evolution that we simulate. Second, we assume non-resonant initial conditions for the giant planets. Mean-motion resonances can be established readily if there is significant migration, due to either gas disk torques or planetesimal scattering, and a plausible alternate class of models could be constructed in which fully resonant initial conditions were the norm (Morbidelli et al. 2007). We do not consider this possibility further here.

Embryo and giant planet particles feel the gravitational attraction of all other bodies in the simulation. Planetesimal particles, both in the inner and outer disk, feel the gravity of embryos and giant planets but do not self-gravitate. This commonly-used approximation allows for an adequate treatment of collective particle effects and dramatically reduces the required computation time. Our methods are outlined in detail in Paper I. Here we summarize the key points. Each simulation was integrated for 100–200 million years using the hybrid version of the Mercury code (Chambers 1999) with a 6 day timestep. Collisions between particles were treated as inelastic mergers. Particles were removed from a simulation when they either came within 0.2 AU of the central (solar-type) star at which point they are assumed to have collided with the star, or ventured farther than 100 AU from the central star (1000 AU for the widedisk runs discussed in Sect. 4.2), at which point the particle is assumed to have been ejected from the system.

Table 1

Summary of the simulations.

In Paper I we presented the results from our fiducial mixed set of simulations. In these simulations the giant planet masses are drawn randomly from the observed exoplanet mass distribution (Butler et al. 2006; Udry & Santos 2007): (1)where masses were chosen between one Saturn mass and 3 Jupiter masses. In these runs, the masses of individual planets were chosen independently.

We also computed a number of alternate models in which either the range of the mass function, the assumption of independent masses, or the properties of the outer disk were altered (see Table 1):

  • The lowmass simulations represent systems with low-massgiant planets. For these cases the giant planet masses also followthe observed exoplanet distribution, but with masses between10 M and 1 MJup. The initial conditions for the giant planets in these lowmass simulations are the same as the “mixed2” simulations in Raymond et al. (2008a, 2009b,a, 2010). Abnormally low giant planet masses, relative to the mass in terrestrial planet-forming material, might occur physically in disks around stars where stronger than average photoevaporation limits the disk lifetime.

  • The equal simulations comprise four sets of simulations, each containing three giant planets with fixed masses of 30   M, 1   MSat, 1   MJ, or 3   MJ. For the equal simulations the planets were placed in a slightly more compact configuration (separated by 3.5–4 mutual Hill radii rather than 4–5) to ensure that they would become unstable. This variation mimics the reality that the conditions that favor the growth of one gas giant to high masses – for example early core formation or a long disk lifetime – probably apply also to other planets in the same system. Past work suggests that these simulations should produce the most violent instabilities (Ford et al. 2003; Raymond et al. 2010).

  • The widedisk simulations test the effect of planetesimal disks that are 20 AU wide rather than 10, and twice as massive (so that there is the same mass in the first 10 AU of the annulus as the fiducial case). In these simulations, the radius beyond which an object is considered ejected was 1000 AU rather than 100 AU. The giant planets’ initial orbits are identical to the mixed set.

  • The seeds simulations test the effect of the mass distribution within the planetesimal disk by including five or ten equally-spaced equal-mass fully self-gravitating seeds of either 2   M or 0.5   M respectively. The total mass and width of the planetesimal disk was held fixed at 10 AU.

  • The gas simulations test the effect of the presence of a gas disk during and after the giant planet instabilities. (The methodology is discussed in Sect. 5.)

Each simulation was post-processed to calculate the spectral energy distribution (SED) of dust in the system following the method of Booth et al. (2009) with a few small changes (see Sect. 2.3 in Paper I). To do this each planetesimal particle was assumed to represent a population of objects with sizes between 2.2   μm and 2000 km. This population was assumed to be in collisional equilibrium such that the differential size distribution can be written as n(D) ∝ D-3.5 (Dohnanyi 1969). The radial distribution of dust was calculated by a simple combination of the planetesimal orbital distribution (and by sampling eccentric orbits at multiple intervals along their orbit equally spaced in mean anomaly). The SED was calculated by assuming that the dust grains in each radial bin emit as blackbodies based on their effective temperature. At each simulation timestep the collisional timescale tc was calculated for the largest objects (D = 2000 km) for the population of both asteroidal and cometary planetesimals. In practice, tc represents the mean timescale between collisions energetic enough to disrupt an object of a given size at a given orbital radius. tc is a function of the mass, width, and orbital distribution of the planetesimal belt as well as the physical properties of planetesimals themselves, in particular their , the impact energy needed to catastrophically disrupt a planetesimal of size D (for details, see Wyatt et al. 1999, 2007b; Booth et al. 2009; Raymond et al. 2011; Kains et al. 2011). Once tc was calculated for a given timestep, the effective dust mass of each population was decreased by a factor of [1 + t/tc(Dc)]-1. This decrease in the dust mass is not self-consistent because the planetesimal mass in the simulations is constant. This effect can be important for the asteroidal planetesimals because their collisional timescale, tc ~ 104 years, is short compared to the interesting timescales for dynamical evolution. The opposite ordering typically applies for the outer, cometary planetesimals, whose collisional timescale is tc ≳ 108 years. The dust fluxes used in the analysis later in the paper are dominated by the cometary component so this inconsistency has little to no effect on our results. In addition, for the case of low-mass giant planets that migrate due to planetesimal scattering (i.e., the lowmass simulations), the timescale for dynamical mass loss from the outer planetesimal disk is roughly an order of magnitude shorter than the timescale for the calculated collisional mass loss in that same region. Thus, our assumption that the planet-planetesimal disk dynamics is not affected by the collisional cascade appears reasonable for the outer disk.

This simple model is based on previous studies that fit the statistics of debris disks using models for the collisional evolution of planetesimals (Dominik & Decin 2003; Krivov et al. 2005, 2006; Wyatt et al. 2007b; Wyatt 2008; Löhne et al. 2008; Kains et al. 2011). Our model agrees to within a factor of 2–3 at 24   μm and 70   μm with more detailed calculation of dust production during the collisional evolution (Kenyon & Bromley 2008, 2010) and also with dust fluxes observed around solar-type stars (Habing et al. 2001; Beichman et al. 2006; Moór et al. 2006; Trilling et al. 2008; Hillenbrand et al. 2008; Gáspár et al. 2009; Carpenter et al. 2009). However, due to our incomplete knowledge of the physical properties of planetesimals, there remains uncertainty in the dust fluxes of up to an order of magnitude for a given system (see Booth et al. 2009). Our model does not include outgassing from comets (i.e., outer disk planetesimals) when they enter the inner solar system. Thus, the dust fluxes that we calculate during planetesimal bombardments are significantly underestimated. Indeed, the 24   μm dust flux during the late heavy bombardment calculated by Booth et al. (2009) with our method reaches a peak that is roughly an order of magnitude lower than that calculated by Nesvorný et al. (2010), who accounted for cometary dust production. Cometary outgassing is of importance at mid-infrared wavelengths (10 ≲ λ ≲ 50   μm) during and shortly after bombardments. As our results focus on the steady-state production of cold dust in outer planetesimal disks rather than on bombardments, we are not strongly affected by this effect.

An important point is the fact that our sets of simulations systematically over-predict the frequency of debris disks by a factor of roughly 2–4. Given Spitzer’s detection limits, the observed frequency of debris disks at 70   μm around Solar-type stars older than 1 Gyr is 16.4% (Trilling et al. 2008; Carpenter et al. 2009). At 24   μm, the observed frequency is much lower, only 2–3% for stars older than 300 Myr (Carpenter et al. 2009; Gáspár et al. 2009). Various aspects of and potential solutions to this issue will be discussed in more detail throughout the paper.

thumbnail Fig. 1

Evolution of a reference mixed simulation in which the giant planets underwent a violent instability after 42.8 Myr of evolution. The initial giant planet masses were, in order of increasing orbital distance, 0.96, 0.46 and 0.64   MJ. Left: orbital eccentricity vs. semimajor axis of each body in the simulation. The size scales with the mass1/3 and the color corresponds to the water content, using initial values taken to solar system data (Raymond et al. 2004) and re-calculated during impacts by mass balance. The giant planets are shown as the large black bodies and are not on the same size scale. Top right: the spectral energy distribution of the dust during five simulation snapshots, as observed from a distance of 10 pc. The dashed line represents the stellar photosphere. Bottom right: the ratio of the dust-to-stellar flux F/Fstar at 25   μm as a function of time. The rough Spitzer observational limit is shown with the dashed line (Trilling et al. 2008). A movie of this simulation is available online or at http://www.obs.u-bordeaux1.fr/e3arths/raymond/scatterSED_199.mpg.

Open with DEXTER

2.1. An example mixed simulation

Here we briefly present a simulation from the mixed set to facilitate comparison against example simulations from other sets that are presented in later sections. The chosen system started with giant planets of 0.96, 0.46 and 0.64 MJ in order of increasing orbital distance, and an outer planetesimal disk that extended out to 22.8 AU.

Figure 1 shows the evolution of the simulation’s dynamics and calculated dust flux. The system remained stable for 42.8 Myr at which time it underwent a strong dynamical instability that started with a close encounter between the middle and outer giant planets that triggered a series of planet-planet scattering events over the next 400 000 years. The instability culminated in the ejection of the middle giant planet, and the surviving two giant planets swapped orbits (i.e., the innermost planet became the outermost and vice versa). At the end of the simulation both planets’ eccentricities are large: einner oscillates between 0.61 and 0.83 and eouter between 0.15 and 0.37. Despite their large eccentricities, both of the planets’ inclinations with respect to the initial orbital plane remain modest (at least in comparison with the planets’ eccentricities): iinner oscillates between 3.5° and 14° and iouter between zero and 3.1°. With a semimajor axis of 2.55 AU, the orbit of the inner giant planet is well-represented by the more eccentric of the known exoplanets, while the outer planet would probably not be currently detectable.

Before the instability the inner planetary system was undergoing standard terrestrial accretion. Embryos’ eccentricities remained damped by dynamical friction from planetesimals and embryos grew by frequent planetesimal impacts and occasional giant embryo impacts. At the time of the instability the accretion process was relatively mature, as only eight embryos remained inside 2 AU with masses between 0.08 and 0.91   M, as well as three smaller (0.07−0.2   M) embryos in the asteroid belt. In the immediate aftermath of the instability, ten of the eleven embryos collided with the central star. The mechanism that drove the embryos into the star was strong eccentricity pumping by a combination of close encounters and secular pumping by the (initially outermost) inward-scattered giant planet. The one embryo that was not driven into the star was the outermost one, which had a pre-instability semimajor axis of 2.6 AU and was ejected after a series of close encounters with the scattered giant planets.

The pre-instability outer planetesimal disk was for the most part dynamically calm. The inner edge of the disk was slowly eroded during this period as the giant planets cleared out planetesimals that were unstable on long timescales. A few other resonances (e.g., the 3:1 resonance with the outer giant planet at 21.4 AU) also acted to increase the eccentricities of long-term stable planetesimals. A slow trickle of planetesimals was also destabilized by certain strong mean motion resonances, notably the 2:1 resonance. When the instability started, the (previously middle) giant planet was scattered out into the planetesimal disk on an eccentric orbit. Its eccentricity was further pumped by a series of close encounters with the (initially innermost) planet such that for a period of several hundred thousand years (several thousand cometary orbits) the giant planet’s orbit completely crossed the initial planetesimal disk. The outer planetesimal disk was entirely destabilized by secular interactions and close encounters. These planetesimals’ eccentricities increased drastically until they either hit the sun (this occurred about 25% of the time), were scattered out beyond 100 AU and removed from the simulation by presumed hyperbolic ejection (~75%), or collided with a giant planet (~0.5%). More than 80% of the outer disk planetesimals were destroyed within 500 000 years and 97% within 5 Myr. Only 23 planetesimal particles survived more than 5 Myr after the instability on orbits that were unstable on 10 Myr timescales, typically with high inclinations and eccentricities. At the end of the simulation a single planetesimal survived, although it is almost certainly unstable on longer timescales as its orbit crosses the outer giant planet’s.

The system’s SED – shown at right in Fig. 1 – reflects its dynamical evolution. The asteroidal planetesimals are quickly ground to dust, as their collisional timescales are only 104−105 years. This causes a rapid decline in flux at short wavelengths (λ ≲ 20   μm). The erosion of the inner edge of the outer planetesimal disk causes a continued decrease in flux at shorter wavelengths. The flux at λ ≳ 50   μm is dominated by the mass in outer disk planetesimals and is only very weakly affected by the grinding of asteroids or the slow erosion of the inner edge of the planetesimal disk. When the instability occurs and the outer planetesimal disk is destabilized, a large number of planetesimals are temporarily placed on high-eccentricity (and therefore small-periastron) orbits which produce hot dust. This burst in hot dust changes the shape of the SED by increasing the flux at short wavelengths. This is the cause of the spike in flux seen at 25   μm in the lower right panel of Fig. 1. As noted already, the magnitude of the spike is underestimated because we do not account for the outgassing that occurs when an icy body is heated (see Nesvorný et al. 2010). However, as the outer planetesimal disk is removed the flux drops dramatically. The dust flux in the 30 Myr after the instability is maintained at a relatively high level by a single planetesimal that survived from 43 to 66.5 Myr between the giant planets with a perihelion that dropped periodically below 3 AU (on a retrograde orbit). This single close-in planetesimal produced enough dust to keep the system above the 25   μm detection threshold during this period. We note that our dust production scheme does not allow for the collisional evolution of a single planetesimal so its dust flux is certainly overestimated.

This example is relatively extreme in terms of the planets’ final orbital eccentricities and in that all the terrestrial and cometary particles were destroyed. However, as we will see below, this simulation allows for a convenient comparison with upcoming examples because the instability is delayed and so the evolution of the dust flux from the quiescent outer disk is unperturbed at early times.

thumbnail Fig. 2

Eccentricity vs. semimajor axis for the surviving giant planets in the mixed (left panel) and lowmass (right) simulations. Black circles represent unstable simulations – defined as simulations in which at least one giant planet-planet scattering event occurred – and grey dots are stable simulations. The size of each circle is proportional to the logarithm of the planet mass. The dashed vertical lines represent the median outer edge of the planetesimal disk for each set of simulations. Note that the outer edge varied from simulation to simulation depending on the giant planet masses, from 20.6 to 27 AU (with a median of 23.7 AU) for the mixed simulations and 17.3 to 23.3 AU for the lowmass simulations (median of 19.7 AU).

Open with DEXTER

3. Effect of the giant planet masses and mass distribution

We now analyze two sets of simulations that explore alternate giant planet mass distributions than the mixed set analyzed in Paper I. The mixed set included three planet systems with the masses of the planets being chosen randomly and independently in the range between a Saturn mass and 3 Jupiter masses. In the equal simulations, planet masses within a given system are the same. We consider masses between 30   M and 3   MJ for different systems. Assuming that the masses within individual systems are perfectly correlated has the effect of maximizing the strength of dynamical instabilities, compared to systems where there is a range in masses. In the lowmass simulations, planet masses are drawn randomly from the observed distribution but only in the range of 10   M to 1   MJ (in contrast to the range of MSat to 3   MJ for the mixed set). The inclusion of the lower mass planets increases the fraction of systems for which dynamical interactions between planets and the planetesimal disk are important. These interactions can take two forms. First, “planetesimal-driven migration” changes the orbital radius of a planet due to the back-reaction of planetesimals that are gravitationally scattered by the planet, and thus changes a planet’s orbital radius while maintaining a small eccentricity (Fernandez & Ip 1984; Malhotra 1993; Hahn & Malhotra 1999; Gomes et al. 2004; Kirsh et al. 2009; Levison et al. 2010). Second, the orbit of an eccentric planet can be re-circularized by “secular friction”, a process by which an eccentric planet excites the eccentricities of the outer disk planetesimals and causes a corresponding decrease in the planet’s eccentricity (Thommes et al. 1999; Levison et al. 2008). A low-mass planet can therefore be gravitationally scattered by another planet in the inner part of a planetary system and have its eccentricity decreased on a much wider orbit by secular friction with the outer planetesimal disk (e.g., Thommes et al. 1999; Raymond et al. 2010).

The low mass and equal sets of simulations are of particular interest because a combination of the two can produce an alternate sample that matches the observed exoplanet distribution. We explore the two sets of simulations independently (Sects. 3.1 and 3.2) and later combine them into a sample to compare with exoplanet statistics (called case B in Sect. 6).

3.1. Systems with low-mass giant planets (the lowmass simulations)

Figure 2 shows that, because of planetesimal-driven migration and secular friction, surviving low-mass giant planets populate different regions of parameter space than more massive giant planets. Given that all of our simulated giant planets started in the 5–15 AU region, massive giant planets are only able to alter their semimajor axes by energy exchange during close encounters to essentially follow curves with perihelia or aphelia at the encounter distance. In other words, high-mass planets at large a necessarily have large e. However, low-mass planets can have large a and small e by either 1) planetesimal-driven migration, which maintains planets’ low e out to large a (to the outer edge of the planetesimal disk), or 2) being scattered outward in a dynamical instability but having their eccentricities damped by secular friction with the outer planetesimal disk. Similarly, massive planets that are scattered interior to 5.2 AU necessarily have large e but low-mass inner giant planets can undergo inward changes in a and survive on low-e orbits. Indeed many low-mass planets do just that, ending up at a = 2.5−4 AU. Note that secular friction is only relevant in unstable systems in which a planet acquires a large orbital eccentricity. On the other hand, planetesimal-driven migration is mainly relevant for stable systems although periods of migration may occur in some cases after secular friction has already re-circularized the orbit of a scattered low-mass giant planet.

Thus, the surviving high-mass giant planets retain a memory of their initial conditions: the stable planets and many unstable planets are clustered at their original locations. However, given the ease and inevitability of planetesimal-driven migration and secular friction for low-mass giant planets, the initial conditions are erased.

thumbnail Fig. 3

Evolution of a simulation with low-mass giant planets. Left: orbital eccentricity vs. semimajor axis of each body in the simulation. The size scales with the mass1/3 and the color corresponds to the water content, using initial values taken to solar system data (Raymond et al. 2004) and re-calculated during impacts by mass balance. The giant planets are shown as the large black bodies and are not on the same size scale. Top right: the spectral energy distribution of the dust during five simulation snapshots, as observed from a distance of 10 pc. The dashed line represents the stellar photosphere. Bottom right: the ratio of the dust-to-stellar flux at 25 microns as a function of time. The rough observational limit of the MIPS instrument on NASA’s Spitzer Space Telescope is shown with the dashed line (Trilling et al. 2008). A movie of this simulation is available online or at http://www.obs.u-bordeaux1.fr/e3arths/raymond/scatter_lowmass_54.mpg.

Open with DEXTER

Figure 3 shows the evolution of a lowmass simulation with initial giant planet masses of 12.4 M (inner), 18.6 M (middle), and 35.9 M (outer). In this simulation the giant planets underwent an instability after 33 000 years, which threw the inner planet into the outer planetesimal disk. Once interacting with the planetesimal disk, the planet scattered planetesimals inward and migrated outward for roughly 20 Myr, then slowed when it came within  ~3.5 Hill radii of the outer edge of the planetesimal disk. As this represents the approximate boundary for dynamical stability (Marchal & Bozis 1982; Gladman 1993; Chambers et al. 1996), the number of planetesimals available to be scattered by the planet decreased and the planet’s migration slowed drastically. The planetesimals that were scattered inward from the outer disk were for the most part subsequently scattered outward by the inner giant planets, causing the two inner planets to migrate inward. However, some of the inward-scattered planetesimals were trapped on low-eccentricity orbits between the two inner giant planets and the outer giant as the outer planet migrated outward by continued planetesimal scattering. This is similar to the mechanism that may have been responsible for populating the solar system’s asteroid belt during the outward migration of Jupiter and Saturn (Walsh et al. 2011). The inner giant planets’ inward migration was fueled by the planetesimals that ended up with high-eccentricity, low-perihelion orbits after being scattered inward, although the inner giants also scattered some embryos and planetesimals from the inner disk. The inner giant planet migrated in to 2.69 AU but maintained an eccentricity lower than 0.1 throughout and less than 0.05 during the last phases.

Two terrestrial planets formed in the simulation shown in Fig. 3: a 1.25 M planet at 0.61 AU and a 0.69 M planet at 1.11 AU. The eccentricities of the inner and outer planet are 0.06 and 0.11, respectively, with peak to peak oscillation amplitudes of 0.11 and 0.20. The inner planet underwent its last giant (embryo) impact after 40.6 Myr but the outer planet did not undergo any giant impacts after 3.8 Myr. The inner planet is wet: it accreted a small amount of water from material that originated in the inner asteroid belt (an embryo and two planetesimals from  ~2 AU) but the bulk of its water came from a single cometary impact. The outer planet did not accrete any material from beyond 2 AU and so is considered to be dry (see Raymond et al. 2004, 2007). The closest giant planet’s semimajor axis is only 1.6 AU larger than the outer terrestrial planet’s but given the giant planet’s small mass (19.2   M) and modest eccentricity (0.035 with oscillations of 0.02 in full amplitude) the system is stable.

In the simulation from Fig. 3, only a relatively small fraction of the initial terrestrial mass was incorporated into the two surviving terrestrial planets. The majority of the initial terrestrial mass (57%) was ejected from the system and an additional 10% collided with the central star. This contrasts with the unstable systems with higher-mass giant planets, in which terrestrial material preferentially collides with the star (as in the simulation from Fig. 1). This difference is due to the fact that more massive giant planets pump the eccentricities of terrestrial bodies efficiently and can thus drive down their perihelion distances on short timescales. Low-mass giant planets require longer to excite high eccentricities and also migrate in reaction to scattering such bodies, inward in this scenario. Thus, in systems with unstable low-mass giant planets terrestrial material is more easily transported outward than inward, to be ejected after many encounters with the outer giant planets.

In addition, in the simulation from Fig. 3 two Mars-sized embryos were scattered out and survived on distant orbits, at 29.2 and 34.8 AU, and one embryo collided with the middle giant planet.

At the end of the simulation there are two surviving planetesimal belts: a low-eccentricity belt between the two inner and the outer giant planets and an outer disk of higher-eccentricity objects exterior to the outer giant planet. The outer belt is analogous to the solar system’s scattered disk (Luu et al. 1997; Duncan & Levison 1997), having been scattered by the outward-migrating giant planet. The scattered belt contains 3.7 M in 72 particles with a median eccentricity of 0.27 and a median inclination of 11.5°. This scattered disk also contains two embryos from the inner disk, including one that originated inside 2 AU. The inner belt of planetesimals – located between roughly 8 and 14 AU – contains 1.3 M in 27 planetesimals with a median eccentricity of 0.12 and a median inclination of 16.8°. The orbital distributions and surface densities of these two populations are quite different, and we suspect that a wide diversity of planetesimal belt structures must exist around other stars.

The evolution of the system’s dust brightness is shown in Fig. 3. The SED of the system decreases systematically as the system loses mass, but changes shape after roughly 80 Myr when four separate icy planetesimals entered the very inner planetary system and remained on orbits interior to the innermost giant planet (with perihelion distances as small as 0.3 AU) for several tens of Myr before being ejected1. At wavelengths longer than  ~50 μm, the dust brightness decreased monotonically in time. However, shorter wavelengths (such as 25   μm; Fig. 3) show the additional structure caused by the icy planetesimals entering the inner planetary system because they are sensitive to hot dust.

As a whole, the lowmass simulations were extremely efficient at forming terrestrial planets and also at creating long-lasting debris disks. Out of the 86 total simulations, 82 (95.3%) formed terrestrial planet systems containing a total of at least 0.5 M. Of the four remaining systems, three destroyed their terrestrial planets entirely and the fourth formed a single planet of 0.25   M. Of the 86 lowmass simulations, 76 (88.4%) were above the Spitzer detection threshold at 70   μm after 1 Gyr (73 remained bright after 3 Gyr; recall that this is far higher than the observed frequency of 16.4%). Of the 54 (62.7%) unstable simulations, only 4 did not yield at least 0.5 M in terrestrial planets and only 10 were not detectable at 70   μm after 1 Gyr. The few systems with destructive instabilities were those that by chance contained several massive planets, and so essentially overlapped with the mixed distribution. Figure 4 shows that all 32 stable simulations finished with two or more terrestrial planets and also with bright debris disks. Note that in this figure we use a low mass cutoff of just 0.05   M in our definition of a terrestrial planet. Any surviving planetary embryo is therefore considered a planet. This allows for a consistent comparison with other sets of simulations including more violent instabilities (like the equal simulations) in which surviving embryos are common. In the solar system, Mars is thought to be a surviving embryo (Dauphas & Pourmand 2011).

thumbnail Fig. 4

Distribution of the number of surviving terrestrial planets for the unstable (grey) and stable (dashed) simulations in the lowmass set of simulations.

Open with DEXTER

Among systems with planets less massive than 50–100 M (≈0.5−1MSat), there was little difference in the final outcome between systems that underwent planet-planet scattering and those that did not. This is because the planetesimal disk provides strong enough damping to quickly decrease the planets’ eccentricities back to near zero. Systems containing a single relatively massive giant planet (M ≳ MSat) also ended in a dynamically calm state because instabilities caused the lower-mass giant planets to be scattered and, again, their eccentricities and inclinations are quickly damped. The only situation that preserved large eccentricities was the relatively infrequent combination of multiple massive giant planets in the same system. In those systems the large eccentricity caused by strong scattering between giant planets could not be damped (the low-mass giant planets in such systems are usually scattered, sometimes to be ejected or sometimes re-circularized in the outer planetesimal disk). Thus, unlike massive giant planets, the eccentricities of low-mass planets do not retain a memory of the system’s dynamical history. In addition, multiple giant planets must exist in the same system to yield eccentric giant planets. The abundance of observed eccentric planets (e.g. Wright et al. 2009) thus points to the frequency of strong instabilities (although alternate models exist; see Ford & Rasio 2008, for a thorough review).

thumbnail Fig. 5

In the lowmass planetary systems, the fraction of stable zones between pairs of outer giant planets that contain at least one planetesimal on a stable orbit as a function of the mass of the outermost giant planet (left panel) and the width of the stable zone Δ in units of mutual Hill radii (right panel). For the left panel, bins were evenly logarithmically spaced from 10   M to 1   MJ. For the right panel, bins were evenly spaced from Δ = 3.5 to 30. The error bars were calculated using binomial statistics following Burgasser et al. (2003). The dashed line marks the two planet stability boundary (Marchal & Bozis 1982; Gladman 1993).

Open with DEXTER

Most of the lowmass systems have large gaps between giant planets, and many of these gaps contain planetesimals on stable orbits, such as the belt between 8–14 AU in the system from Fig. 3. The existence of isolated belts alters the radial distribution of dust and can be inferred from the SED. There are several systems known to contain belts of dust that appear to be radially confined, presumably by known or as-yet undetected planets (e.g., Beichman et al. 2005; Lisse et al. 2008; Su et al. 2009).

The majority of the lowmass systems (60/85 = 71%) have gaps between the outer two giant planets with separations Δ of 10 or more mutual Hill radii (in one system just a single giant planet survived and so is not counted). Given that two planets must be separated by at least mutual Hill radii for long-term dynamical stability (Marchal & Bozis 1982; Gladman 1993), the existence of a zone between two planets’ orbits that is stable for planetesimals requires at a minimum an interplanetary separation Δ ≳ 10. In practice, somewhat wider gaps are needed in realistic systems. The separation between Saturn and Uranus, and between Uranus and Neptune, amounts to 14 mutual Hill radii, and there is only a small region that is stable over long timescales between the latter two planets (Holman & Wisdom 1993). In our runs we frequently find gaps that are not only wider than their solar system counterparts, but which are also populated with primordial material despite the relatively coarse sampling of the outer planetesimal disk population. The lowmass simulations yield a range of separations from Δ = 3.9 to 30 and outer giant planet masses of 11   M to 0.97   MJ. From the entire sample, the stable zone between the outer two giant planets contained at least one particle in half of all lowmass simulations (38/85 = 45%, although 5 stable zones contained just a single planetesimal). The closest separation between two planets for which planetesimals existed on stable orbits between the two was Δ = 11.9, and the widest separation for which no planetesimals existed between two planets was Δ = 21.5. The total mass in these planetesimal belts ranged from the mass of one planetesimal particle (0.05   M) to 5.9   M, and in one case an embryo from the inner disk survived in such a belt.

The two most important factors that determine whether a stable zone between outer giant planets contains a planetesimal belt are first, the width of the stable zone and second, the mass of the outermost giant planet. Figure 5 shows that the probability that a stable zone contains at least one planetesimal particle on an orbit that is stable for long timescales (i.e., whose orbit does not come within 4 Hill radii of any giant planet’s orbit) as a function of the mass of the outermost giant planet and the width of the stable zones (as quantified by Δ). Stable zones are preferentially filled for larger Δ values simply because there is a larger region of parameter space into which planetesimals can be scattered and survive, and the fraction of stable zones that is filled increases dramatically for Δ ≳ 15. Stable zones are also preferentially filled in systems with outermost giant planets less massive than  ~50   M, and almost 100% of stable zones are filled when the outermost giant planet is less than about 20 M (Fig. 5). As they interact with the outer planetesimal disk, lower-mass planets scatter planetesimals onto lower eccentricities than do higher-mass planets, and these planetesimals are more likely to avoid encountering giant planets and to remain on stable orbits than if their obits are more eccentric.

Do other system parameters influence the probability that a stable zone between giant planets will contain planetesimals? We tested the importance of several other parameters using a suite of Kolmogorov-Smirnov (K-S) tests that compared different characteristics of stable zones with and without planetesimals. Using a cutoff of p < 0.01 for a statistically significant difference between the two populations, we found that five parameters influence whether a stable zone will contain planetesimals. In order of most important (lowest p value) to least important (highest p value), these are 1) the width of the stable zone (Δ), 2) the mass of the outermost giant planet, 3) the semimajor axis of the outermost giant planet, 4) the mass ratio of the outer two giant planets (systems with a large inner/outer mass ratio preferentially contain planetesimal belts), and 5) the total mass in the surviving giant planets (the probability of stable zones being empty increases for higher total mass). The effect of these five parameters is statistically significant, although several are correlated (e.g., Δ correlates with the semimajor axis of the outermost giant planet, although Δ is about three orders of magnitude more important in determining whether a stable zone will be filled). We tested five additional parameters whose effect turned out to be unimportant (p > 0.1 in each case): the mass and semimajor axis of the giant planet marking the inner boundary of the stable zone, the eccentricity of the outermost giant planet, and the mass-weighted eccentricities of all surviving giant planets.

3.2. Systems with equal-mass giant planets (the equal simulations)

While the lowmass simulations represent a calm environment conducive to the production of both terrestrial planets and bright debris disks, the equal simulations were destructive on both counts. This is simply because scattering among equal-mass giant planets is the most violent planet-planet instability (Raymond et al. 2010), and strong instabilities destroy small bodies in both the inner and outer disks, typically by driving a large fraction of inner bodies into the central star and ejecting the majority of outer bodies. Indeed, the eccentricity distribution of the surviving equal giant planets with masses of MSat or larger was skewed toward higher values than the mixed sample, reaching values as high as 0.89 and with a median of 0.35 (compared with 0.21 for the mixed giant planets). An exception to this rule are systems with equal-mass giant planets that are themselves low-mass. In systems containing three giant planets of 30 M, the outcome was similar to the lowmass systems because planet-disk interactions trumped planet-planet scattering. Indeed, the low-mass equal systems were dominated by planetesimal-driven migration of the outer giant planet and behaved very similarly to the lowmass simulations. Given this strong dichotomy, we now consider just the high-mass equal simulations, as the lower-mass cases are more appropriately included with the lowmass set.

thumbnail Fig. 6

Distribution of the number of surviving terrestrial planets for the equal set of simulations. Systems with giant planets of MSat or larger are shown in grey, and systems with giant planets of 30   M are shown by the black dashed line.

Open with DEXTER

Figure 6 shows that, of the 53 simulations with giant planet masses of MSat, MJ, or 3   MJ, 36 (68%) destroyed all terrestrial material. Five simulations (9%) formed a single terrestrial planet, although these were all small, roughly Mars-mass  ~0.1   M planets, and 4/5 of these were lone surviving embryos. The remaining 12 simulations with high-mass giant planets (21%) formed two or more terrestrial planets. In contrast, not a single simulations with a 30   M giant planet destroyed all of its terrestrial material. Rather, these systems usually formed several terrestrial planets – only 3/15 systems formed just one.

Of the 53 equal systems, only 17 (32%) contained detectable amounts of cold dust at 70   μm after 1 Gyr (only 14/53 = 26% after 3 Gyr, although again note that this is higher than the observed frequency of 16.4%). Of the 17 systems with detectable dust at 1 Gyr, 11 (65%) formed terrestrial planets. Of the 36 systems with no detectable dust, only 6 (17%) formed terrestrial planets. Thus, there exists a natural connection between debris disks and terrestrial planets that spans the domain of giant planet mass. This same debris disk-terrestrial planet correlation was found in Paper I for the mixed simulations. The details of this correlation depend on the system parameters and are of course confined to the context of our initial conditions, in particular to systems with relatively massive outer planetesimal disks (see discussion in Sect. 6).

The equal systems are violently unstable by construction and so the outcomes tend to be extreme. This is precisely the type of behavior that is required by exoplanet observations, in particular the trend for more massive planets to have more eccentric orbits than lower-mass planets (Jones et al. 2006; Ribas & Miralda-Escudé 2007; Wright et al. 2009; Raymond et al. 2010). If planet masses within individual planetary systems were random (as in the mixed simulations), then lower-mass giant planets should have higher eccentricities than higher-mass giants (Raymond et al. 2010); this is the opposite of what is observed. Thus, the equal simulations provide a key ingredient in constructing a sample of simulations that matches the observed giant exoplanets, as discussed in Sect. 6.

4. Effect of the properties of outer planetesimal disks

We now turn our attention to the effect of the properties of outer planetesimal disks. We first examine the seeds simulations – subdivided into the bigseed and smallseed sets – that contained a population of  ~Earth-mass embryos in their outer planetesimal disks. We then test the effects of doubling the width (and total mass) of the planetesimal disk in the widedisk simulations.

4.1. The mass distribution of the planetesimal disk (the seeds simulations)

In the seeds simulations a small number of fully-interacting massive bodies were included in the outer planetesimal disk (in contrast to planetesimal particles, which interact gravitationally with massive bodies but not with each other). The disk maintained the same total mass and numerical resolution (the masses of individual planetesimals were decreased to maintain a constant total mass). In 50 simulations comprising the bigseed set, five icy embryos of 2   M each were included at 11, 13, 15, 17 and 19 AU. In 50 additional smallseed simulations 10 seeds of 0.5   M each were spaced with 1 AU of separation from 10.5 to 19.5 AU. These seed masses are broadly consistent with calculations of accretion in outer planetesimal disks (Kenyon & Bromley 2008, 2010).

thumbnail Fig. 7

Evolution of a simulation that includes 5 seed embryos (shown in grey) in the outer planetesimal disk. Left: orbital eccentricity vs. semimajor axis of each body in the simulation. The size scales with the mass1/3 and the color corresponds to the water content, using initial values taken to solar system data (Raymond et al. 2004) and re-calculated during impacts by mass balance. The giant planets are shown as the large black bodies and are not on the same size scale. The two surviving terrestrial planets have a = 0.71 and 1.27 AU, m = 1.3 and 1.58 M, and Myr-averaged e = 0.09 and 0.17, respectively. Both terrestrials have substantial water contents accreted from hydrated asteroidal material. The two surviving giant planets have a = 4.1 and 24.4 AU, m = 1.27 and 0.45 MJ, and Myr-averaged e = 0.32 and 0.33, respectively. Top right: the spectral energy distribution of the dust during five simulation snapshots, as observed from a distance of 10 pc. The dashed line represents the stellar photosphere. Bottom right: the ratio of the dust-to-stellar flux F/Fstar at 25   μm as a function of time. The rough Spitzer observational limit is shown with the dashed line (Trilling et al. 2008). An animation of this simulation is available online or at http://www.obs.u-bordeaux1.fr/e3arths/raymond/scatterSED_seed13.mpg.

Open with DEXTER

Figure 7 shows the evolution of one bigseed simulation that became unstable after a long delay and therefore allows for a comparison with the evolution of both stable (prior to the instability) and unstable systems. The evolution is qualitatively similar to the mixed simulations in the inner disk such as the system from Fig. 1, as embryos maintain smaller eccentricities than planetesimals via dynamical friction and accrete as their feeding zones begin to overlap. The instability started after 55.35 Myr with a close encounter between the two outer giant planets. After a series of close encounters lasting 40 000 years, the middle giant planet was ejected. At the end of the simulation, the system contains two giant planets each with e ≈ 0.3 and two roughly Earth-sized terrestrial terrestrial planets (details in caption of Fig. 7).

The evolution of the outer disk differs significantly from simulations without seeds. The presence of the seeds introduces an effective viscosity into the outer planetesimal disk that is driven by scattering of planetesimals by close encounters with the icy embryos. This causes the disk to spread out radially. Within 1 Myr the outer edge of the low-eccentricity portion of the disk has expanded from 22 to 26 AU, and beyond 30 AU in the next few Myr (Fig. 7). This outward expansion is balanced by the loss of a large portion of the total disk mass that is scattered inward to encounter the giant planets and be ejected from the system. The planetesimal disk is further depleted on a 20–50 Myr timescale by instabilities in the system of icy embryos that increases their eccentricities – these instabilities are analogous to but far weaker than instabilities between the giant planets (and are weaker still in the smallseeds systems). When the giant planets go unstable the bulk of the planetesimal disk is rapidly ejected and the last icy planetesimal is removed just after 100 Myr.

The dust production in the seeds simulations is also different than in simulations with calmer planetesimal disks. The planetesimal population closest to the giant planets – the inner several AU of the outer planetesimal disk – is depleted in a few Myr as the disk “viscously” spreads. Thus, warm dust has a far shorter lifetime than in simulations with no seeds. This is reflected in the evolution of the 25   μm flux in Fig. 7, that drops below the detection threshold after  ~45 Myr while the system is still stable. In contrast, in the example mixed (Fig. 1) and lowmass (Fig. 3) simulations, the 25   μm flux was more than an order of magnitude higher after 40–50 Myr of evolution.

thumbnail Fig. 8

Correlations with the dust flux after 1 Gyr of dynamical and collisional evolution for the seeds simulations, as compared with the mixed simulations. The top panels show F/Fstar at 70   μm vs. the eccentricity of the innermost surviving giant planet (left) and the total mass in surviving terrestrial planets (right). The bottom panels show the same comparisons but at 25   μm. The mixed simulations are in black, the smallseed simulations in green, and the bigseed simulations in red. Filled circles represent stable simulations and open circles unstable simulations.

Open with DEXTER

thumbnail Fig. 9

Left: the fraction of systems that would be detectable with Spitzer (with F/Fstar(70   μm) ≥ 0.55 after 1 Gyr of collisional and dynamical evolution) as a function of the eccentricity of the innermost giant planet eg for the mixed and combined seeds simulations. The error bars are based on binomial statistics  (see Burgasser et al. 2003). This essentially represents a horizontal slice through the top left panel of Fig 8. Right: the fraction of systems with 0.5   M or more in surviving terrestrial planets as a function of F/Fstar(70   μm) ≥ 0.55 (1 Gyr) for the mixed and seeds simulations. Systems with F/Fstar < 10-2 are included in the bin at F/Fstar ≈ 10-2. The Spitzer detection limit is shown as the dashed line. This represents a vertical slice through the bottom right panel of Fig. 8.

Open with DEXTER

The shortened lifetime of warm dust is reflected in its changing SED (Fig. 7). Compared with simulations without seeds, there is a much faster decrease in flux at λ ≲ 100   μm at early times as the region just exterior to the outermost giant planet is cleared much more efficiently and to a larger radial separation. At longer wavelengths the flux also decreases more rapidly due to depletion by scattering from icy embryos, although at long wavelengths this is counteracted in part by the expansion of the dust disk to larger radii and therefore increased surface area (though lower temperature). There are spikes in the flux (seen at 25   μm) when objects enter the inner planetary system to encounter and be ejected by the innermost giant planet. These spikes are due to the fact that each particle represents a distribution of smaller bodies and, with a higher numerical resolution, these spikes would be less pronounced. There is a large-scale decrease in flux after the instability and all dust disappears when the last planetesimal is ejected after  ~100 Myr.

The giant and terrestrial planet evolution differed only slightly between the mixed and seeds simulations. (Recall that the giant planets in the seeds simulations are identical to the mixed simulations in terms of their total mass and mass distribution.) However, the giant planet instabilities in the seeds simulations appeared slightly weaker. The mean innermost giant planet eccentricity was 0.21 for the seeds simulations compared with 0.27 for the mixed simulations2. The cause of this difference appears to be late encounters between a giant planet and a seed embryo, after the end of giant planet-giant planet scattering. To test this, we sub-divide the 39 unstable bigseed simulations into the 22 simulations in which a giant planet underwent at least one planet-seed scattering event after the completion of planet-planet encounters and the 17 simulations that did not (i.e., the last planet-planet scattering event in these 17 simulations came after the last planet-seed scattering event). The simulations with late planet-seed encounters had systematically shorter instabilities (as measured by the time between the first and last planet-planet scattering events) and lower final giant planet eccentricities. We think that what happens is that in these systems a planet-seed scattering event can lead to a small readjustment in one giant planet’s orbit that separates it sufficiently from other giant planets to stabilize the system. However, we note that these samples are relatively small and we cannot rule out that the weaker seeds instabilities are a product of small number statistics.

Given the weaker giant planet instabilities in the seeds simulations, terrestrial planet formation was correspondingly more efficient: only about 1/4 of the unstable seeds simulations destroyed all of their terrestrial material compared with more than 40% for the unstable mixed simulations. However, the differences between the planetary evolution in the seeds and mixed simulations are small compared with those between some of the other sets of simulations such as the lowmass and equal runs.

Despite their influence on the outer planetesimal disk, seeds underwent little accretion. Among all fifty bigseed simulations, no seed accreted more than five planetesimals, and there was just a single seed-seed collision and 5 giant planet-seed collisions. There was slightly more accretion among the seeds in the smallseed simulations, which had a comparable rate of planetesimal-seed impacts but a higher rate of giant impacts, with 4 seed-seed collisions and 12 giant planet-seed collisions among the fifty simulations, although we note that all but 3 of the giant planet-seed collisions occurred very early and were probably caused by the seed being placed on an orbit that was initially very close to a giant planet.

Figure 8 shows the correlations between the dust-to-stellar flux ratio F/Fstar(70   μm) after 1 Gyr, and either the innermost giant planet eccentricity or the total terrestrial planet mass. We plot results for the smallseed, bigseed, and mixed simulations. The most important difference between the seeds and mixed simulations is that the seeds simulations produce less dust at late times, especially warm dust that is observable at wavelengths shorter than about 100   μm. Among just the unstable subset of simulations, 53 of 96 mixed simulations () had dust fluxes that would be detectable with Spitzer after 1 Gyr of evolution, i.e., with F/Fstar(70   μm) ≥ 0.55 (Trilling et al. 2008). In comparison, 14 of 29 () unstable smallseed and 16 of 39 (= ) bigseed simulations were detectable at 70   μm after 1 Gyr. Given the statistical error bars, the decreased detection rate compared with the mixed simulations is only a 1σ result for the smallseed simulations and 2σ for the bigseed simulations.

The effect of the seeds is more apparent when considering the stable simulations. The stable systems remain detectable at 98% or higher rate for each of the mixed, smallseed and bigseed sets of simulations. However, the actual dust brightness decreases dramatically for simulations with seeds. The median F/Fstar(70   μm) at 1 Gyr was 26.2 for the stable mixed simulations, 5.3 for smallseed, and 1.7 for bigseed. Given the scatter, the difference between the mixed and seeds simulations is significant at the 3σ level and the difference between the bigseed and smallseed is significant at 5σ.

At 25   μm the differences are even more striking. In unstable systems, 12/96 () mixed simulations, 2/29 () smallseed simulations, and 0/39 (0 + 4.4%) bigseed simulations were above the Spitzer detection threshold of F/Fstar(25   μm) ≥ 0.054 after 1 Gyr (Trilling et al. 2008). This constitutes a 1−2σ difference. Among the stable systems, 51/56 () mixed simulations were detectable at 25   μm after 1 Gyr but not a single stable bigseed or smallseed was detectable (0+7.3% for the combined seeds simulations).

With Spitzer’s detection limits, debris disks at 70   μm vastly outnumber those at 24   μm. Around stars older than 300 Myr the frequency of 24   μm dust excesses was estimated at by Carpenter et al. (2009) and at 1.9% ± 1.2% by Gáspár et al. (2009). However, these estimates are based on just a handful of detections (from more than 100 targets). Among the 24   μm detections are several systems such as η Corvi (Wyatt et al. 2005) and HD 69830 (Beichman et al. 2005; Bryden et al. 2006) that appear to contain dust at  ~1 AU; this dust has been interpreted as either being transient (Wyatt et al. 2007a) or due to a very peculiar outcome of planet formation (Wyatt et al. 2010). Thus, the frequency of systems with 24   μm excess due to collisional equilibrium processes is significantly smaller than the quoted values. The frequency of dust excesses at 70   μm is  (Trilling et al. 2008). This is at least 5–8 times higher than at 24   μm, and the removal of systems with potentially transient dust can only cause this ratio to increase.

The mixed simulations produce an overabundance of 25   μm dust excesses. After 1 Gyr of evolution, the frequency of detectable dust at 25   μm was and (1 − σ error bars) for stable and unstable systems, respectively. At 70   μm the frequency of detectable dust was for stable and for unstable systems. Thus, the ratio of the fraction of systems that were detectable after 1 Gyr at 70   μm to 25   μm is 1.08 for the stable mixed systems and 4.4 for unstable mixed systems. The higher ratio for the unstable systems is a simple consequence of the instability preferentially clearing out the inner portion of outer planetesimal disks and leaving behind the colder part of the disk that does not emit much flux at 25   μm. Nonetheless, no combination of these ratios for the mixed simulations can match the observed ratio of 5–10 after 1 Gyr of evolution3.

The seeds simulations may explain the dearth of 25   μm dust excesses because there is a natural suppression of the 25   μm flux for both stable and unstable systems. From the 44 stable and unstable smallseed simulations, only 2 were detectable at 25   μm. None of the bigseed simulations were detectable at 25   μm. The reason for the lack of 25   μm flux is that the viscous-like spreading out of the outer planetesimal disk acts to both deplete the inner part of the disk by inward scattering and to push the outer part of the disk to ever colder temperatures. The net effect is the near complete frustration of the 25   μm flux. However, the detection rates at 70   μm are higher than 50% for both the smallseed and bigseed simulations, including both the stable and unstable systems. Thus, the existence of  ~ Earth-mass seeds in outer planetesimal disks may provide a natural explanation for the very low frequency of 25   μm excesses compared with 70   μm excesses.

As seen in Fig. 8, the anti-correlation between giant planet eccentricity eg and debris disks still holds for the seeds simulations. For eccentricities larger than 0.1, the dust fluxes show a rapid decrease for the sets of simulations with and without seeds because the dynamics of unstable giant planets dominates the survival of planetesimals. However, for smaller eccentricities there is a clear segregation: the bigseed systems have the smallest dust fluxes, the mixed systems have the largest, and the smallseed are in the middle. In this realm the stirring up of outer planetesimal disks dominates the dust flux, and as we’ve seen before the seeds simulations create a lower-mass and colder planetesimal disk than the mixed systems, leading to lower dust fluxes in proportion to the seeds’ mass (not number). In fact, many of the unstable bigseed systems with modestly eccentric giant planets (eg ~ 0.1) have dust fluxes as high as the stable systems. However, despite their differences, Fig. 9 shows that the fraction of systems that is detectable at 70   μm as a function of eg is very similar for the seeds and mixed systems.

The debris disk-terrestrial planet correlation still holds for the seeds simulations. As with the eg, the correlation between the dust flux at 70   μm after 1 Gyr and the total surviving terrestrial planet mass is less evident because stable seeds systems that efficiently form terrestrial planets produce far less dust than their mixed counterparts. Again, for the stable systems there is a clear segregation between the sets of simulations, with the largest seed mass bigseed having the smallest flux. However, Fig. 9 shows that the fraction of systems that form at least 0.5   M in terrestrial planets increases for both the seeds and mixed simulations. However, the seeds curve is 1−2σ higher than the mixed curve close to the detection limit. This reflects the fact that the seeds simulations deplete their outer planetesimal disks far more than the mixed simulations: only very calm systems preserve enough planetesimals to produce dust. The systems that are observed at a given dust flux therefore represent more stable systems for the seeds simulations than the mixed ones. Thus, the seeds simulations predict an even stronger correlation between stars with observed debris disks and yet-to-be-discovered terrestrial planets.

4.2. The width of the outer planetesimal disk (the widedisk simulations)

thumbnail Fig. 10

Orbital evolution of a widedisk simulation that underwent a Nice model-like instability. The evolution of the semimajor axes, perihelion and aphelion distances of the three giant planets are shown in black, and the semimajor axes of the 19 surviving outer planetesimals in grey (note that one additional planetesimal survived with a = 59 AU). Two terrestrial planets formed in this system, at 0.71 AU (1.4   M) and 1.76 AU (0.13   M): accretion was perturbed by the perturbations moderately eccentric inner giant planets at 4.7 AU (0.78   MJ, e oscillates between 0.03 and 0.17) and 7.8 AU (2.4   MJ, e oscillates between 0.06 and 0.1).

Open with DEXTER

thumbnail Fig. 11

Left: the dust-to-stellar flux ratio at 1 Gyr at 70   μm vs. the innermost giant planet’s eccentricity for the widedisk (grey) and mixed (black) simulations. Filled circles represent stable simulations and open circles unstable ones. The solar system is shown with the grey star. Right: histogram of the fraction of systems that are detectable at 70   μm as a function of the innermost giant planet’s eccentricity; this is essentially a horizontal slice through the left panel.

Open with DEXTER

The widedisk simulations allow us to test the effects of a higher-mass, wider outer planetesimal disk. Compared with the mixed simulations, the outer planetesimal disk in the widedisk simulations was twice as wide, 20 AU rather than 10 AU, and contained twice the total mass in planetesimals, 100   M instead of 50   M. The inner 10 AU of the outer planetesimal disk is therefore the same for the widedisk and mixed simulations (although the numerical resolution is halved in the widedisk runs) but the widedisk systems contain an additional 10 AU of planetesimals. In addition, the outer boundary of each simulation – the limit beyond which particles are considered to be ejected – was 1000 AU rather than 100 AU.

The widedisk simulations behaved similarly to the mixed simulations in most respects. However, the more massive outer planetesimal disk did cause a few notable differences between the widedisk and mixed simulations. Given the much larger angular momentum reservoir contained in the outer planetesimal disk,  ~Saturn-mass giant planets on eccentric orbits can be captured in the outer system by planetesimal scattering. Figure 10 shows the evolution of one such system, in which the two inner planets (Minner = 0.78   MJ,Mouter = 2.4   MJ) started just interior to the 2:1 mean motion resonance. The giant planets’ clearing out of the inner portion of the outer planetesimal disk drove the two inner giant planets across the 2:1 MMR after  ~0.1 Myr. This caused a perturbative increase in their eccentricities and a corresponding increase in the eccentricity of the outer, Saturn-mass planet. The outer planet’s semimajor axis increased quickly and over the next 10 Myr its eccentricity slowly decreased by secular friction with the outer disk of planetesimals. Despite the perturbative evolution of the system, 19 planetesimals totaling 1.9   M survived in the outer planetesimal disk (their orbital evolution is shown in grey in Fig. 10). Most of these started the simulation in the outer parts of the planetesimal disk (one at 18 AU). As the outer giant planet migrated outward, it shepherded many of these planetesimals in its 3:2, 2:1, and even its 3:1 mean motion resonances, located at 30.8, 37.2 and 48.9 AU at the end of the simulation, respectively. This is analogous to the shepherding of Kuiper belt objects such as Pluto during Neptune’s planetesimal-driven migration (Levison et al. 2008). The surviving planetesimal disk in the simulation from Fig. 10 is massive enough to remain detectable at 70   μm for 3 Gyr but is cold enough not to be detectable at shorter wavelengths. Planetesimal-driven migration of a massive planet therefore represents another mechanism – in addition to the presence of seeds in outer planetesimal disks – to deplete outer planetesimal disks and to push them outward to colder temperatures. However, this mechanism operated efficiently in only a small fraction of widedisk simulations.

thumbnail Fig. 12

Left: the dust-to-stellar flux ratio at 1 Gyr at 70   μm vs. the total terrestrial planet mass for the widedisk (grey) and mixed (black) simulations. Filled circles represent stable simulations and open circles unstable ones. The solar system is shown with the grey star. Right: histogram of the fraction of systems that contain 0.5   M or more in terrestrial planets as a function of F/Fstar(70   μm) after 1 Gyr; this is essentially a vertical slice through the left panel.

Open with DEXTER

Despite having identical giant planet initial conditions, a significantly smaller fraction of widedisk simulations went unstable compared with the mixed simulations ( for widedisk vs. for mixed). This is because of the larger angular momentum reservoir in the outer planetesimal disks. For cases when the instability starts in the outer portion of the planetary system, as the outer planet’s orbit becomes eccentric the disk’s ability to transfer orbital angular momentum (as well as energy) to damp the outer planet’s eccentricity and also to increase its semimajor axis increases for a more massive planetesimal disk. Thus, if all systems of giant planets formed on orbits that would be unstable in the absence of planetesimal disks, their long-term orbital evolution should vary with the outer disk mass. The distribution of outer disk masses may therefore play a critical role in shaping the dynamics of inner planetary systems.

Figure 11 compares the debris disk-giant planet eccentricity anti-correlation for the widedisk and mixed simulations. Both distributions follow the same shape as the other sets of simulations: at eg ≲ 0.05 the vast majority of systems were dynamically stable and therefore have high fluxes, and for eg ≳ 0.05 unstable systems dominate and the probability of having a significant dust flux decreases strongly with increasing eg. For stable systems, the 70   μm fluxes are generally 2–3 times larger for the widedisk systems because their outer planetesimal disks are more massive and the dust flux is dominated by the outer regions of the disk where the collisional evolution is slow. However, the ratio of the median F/Fstar(70   μm) after 1 Gyr for the widedisk to the mixed simulations was 2.7, but at 25   μm the value was only 1.6. This shows that in the inner regions of the outer planetesimal disk, where collisional evolution is faster, the higher-mass widedisk disks approach the dust production level of the lower-mass mixed disks.

The distribution of fluxes for the unstable systems is the same for the widedisk and mixed systems, meaning that the rate of survival of planetesimal disks is dominated by the giant planet dynamics rather than the initial conditions, even the total mass and radial extent. For these unstable systems, the distributions of instability times for the mixed and widedisk simulations were indistinguishable. The distribution of the fraction of systems that is detectable with Spitzer is almost identical for the two sets of simulations: in both cases there is a plateau at  ~100% for eg ≲ 0.05 and a sharp decrease toward higher eccentricities.

Figure 12 compares the debris disk-terrestrial planet correlation for the widedisk and mixed simulations. Again, the two distributions have the same shape but the dust fluxes for the stable systems (i.e., for those with the highest surviving terrestrial planet mass) are 2–3 times higher for the widedisk systems. The distributions are almost indistinguishable for systems with less than about 2   M in surviving terrestrial planets. For both sets of simulations the fraction of systems with 0.5   M or more in terrestrial planets increases monotonically with F/Fstar (70 μm) at 1 Gyr. However, above the detection limit a slightly higher fraction of mixed systems contain terrestrial planets compared with widedisk. This is the opposite of the effect that we saw in Sect. 4.1 for the seeds systems. The widedisk systems form very bright dust disks and require somewhat stronger giant planet perturbations to decrease the flux below the detection limit. Thus, systems with bright debris disks at 70   μm are less sensitive to the presence of terrestrial planets than the mixed systems. Note that this difference comes from the slightly higher fraction of widedisk systems with eccentric giant planets that produce 70   μm excesses (see Fig. 11).

5. Effect of the gas disk during instabilities (the gas simulations)

In the gas set of simulations we added additional forces to the Mercury integrator (Chambers 1999) that acted on planetesimal and embryo particles to mimic the effects of the dissipating gaseous protoplanetary disk from which the planets formed. In these simulations we included two effects: 1) aerodynamic gas drag due to the headwind felt by bodies orbiting at the Keplerian speed while gas orbits slower due to pressure support. Aerodynamic drag acts most strongly on small objects, i.e. planetesimals, and leads to a rapid decay in eccentricity and inclination as well as a slower decay of the semimajor axis; and 2) tidal damping (also called “type 1 damping”) due to gravitational interactions between objects and the disk. Type 1 damping increases for more massive bodies because it is caused by waves excited in the disk, meaning that this was an important source of dissipation for embryos but not for planetesimals. Aerodynamic drag was calculated using standard models (Adachi et al. 1976) assuming planetesimals to be spheres with radii of 10 km. Type 1 damping was included based on linear calculations for planets embedded within isothermal disks (Tanaka & Ward 2004). We included additional terms for large eccentricities and inclinations that were derived by Cresswell & Nelson (2008). We included type 1 damping but not type 1 migration both to maintain a clearer comparison with the simulations without gas forces and because eccentricity damping from the disk is roughly 2 orders of magnitude faster than radial migration (Tanaka & Ward 2004; Cresswell & Nelson 2008).

We assumed the presence of an underlying gas disk that corresponds to roughly half the minimum-mass solar nebula (Weidenschilling 1977; Hayashi 1981), with surface density profile Σ(r) = 875   (r/1   AU)−3/2      g   cm-2 and vertical density distribution , where z0(r) = 0.0472   (r/1   AU)5/4      AU (see also Thommes et al. 2003; Raymond et al. 2006a; Mandell et al. 2007). We note that the distribution of solid mass in our initial conditions is distributed according to a different density profile, Σ ∝ r-1, in both the terrestrial and outer planetesimal zones. We used a steeper radial surface density profile for the gas in order to increase the gas density in the inner disk to maximize the effect of gas drag on the survival of terrestrial bodies. However, we note that our initial conditions for the inner and outer disks – which were chosen as represent approximate guesses for the solar system’s primordial disk – do not even follow the same global surface density profile because there is far too little mass in the terrestrial zone. The solution to this problem may lie with variations in the efficiency of planetesimal formation at different orbital radii within protoplanetary disks (e.g. Chambers 2010).

To model the final stage in the lifetime of the gaseous disk, the disk’s surface density was dissipated linearly and uniformly in 500 000 years. This is slightly longer than most estimates of the final dissipation phase (Simon & Prato 1995; Wolk & Walter 1996; Chiang & Murray-Clay 2007; Currie et al. 2009) and should thus maximize the importance of the gas disk phase. This situation is roughly consistent with models for dynamical instabilities among planets in the presence of gas disks (Chatterjee et al. 2008; Moeckel et al. 2008; Matsumura et al. 2010; Marzari et al. 2010; Moeckel & Armitage 2012), which predict that instabilities should preferentially occur late in the disk phase. The power-law gas density profile probably overestimates the amount of gas interior to the giant planets (Crida & Morbidelli 2007), so these simulations should provide an upper limit to the effects of gas on terrestrial bodies. The initial conditions for the gas simulations were drawn directly from the mixed set of simulations with one important change: the middle giant planet’s eccentricity was increased to make its orbit cross the orbit of the innermost planet. This was to ensure that the system would be immediately unstable so that we could test the effects of the relatively short-lived gaseous disk.

thumbnail Fig. 13

Total mass in surviving terrestrial planets as a function of the minimum perihelion distance of any giant planet during the simulation for the unstable mixed (black dots) and the gas (grey dots) simulations.

Open with DEXTER

The goal of the gas simulations is to test the effect of damping from the gas disk on the dynamics and survival of rocky and icy bodies in the inner and outer planetary system. To accomplish this, we want the giant planet instabilities to be the same as for the mixed set, to isolate the effects of the disk. Thus, we neglected type 1 and type 2 radial migration of giant planets (and embryos) in the gas simulations, although they would certainly occur in a realistic minimum-mass disk (e.g., Lin & Papaloizou 1986; Ward 1997). Indeed, in a more self-consistent setting the giant planets would likely be trapped in resonance at early times and the instability would be triggered by either eccentricity excitation (Marzari et al. 2010; Libert & Tsiganis 2011) or the dispersal of the gas disk (Moeckel et al. 2008; Chatterjee et al. 2008; Moeckel & Armitage 2012).

The giant planets behaved similarly in the gas and the unstable mixed simulations (note that in this section we compare with only the unstable portion of the mixed simulations because all of the gas simulations were unstable). The median eccentricity of the surviving giant planets in the gas simulations is slightly higher – 0.26 vs. 0.22. The reason for the stronger instabilities in the gas simulations is that they were initially placed on strongly unstable, crossing orbits and thus were unable to undergo weaker instabilities that can occur when they develop more slowly (and recall that no damping was felt by the giant planets in the gas simulations). The gas giants provided a comparable fit to the exoplanet distribution (p value from K-S test of 0.67 for gas, 0.49 for mixed). A notable difference between the distributions is a population of low-eccentricity (e ≲ 0.05) giant planets that is significantly more abundant in the unstable mixed than the gas simulations. This population was generated by weakly unstable systems and these systems are the most efficient at both forming terrestrial planets and producing debris disks. Although the instabilities in the gas simulations occurred very early by design, the distributions of the duration of instabilities were virtually identical for the gas and unstable mixed simulations, extending from 104 to 107 years with a median of slightly more than 300 000 years. This is due in part to the fact that we have not included appropriate drag forces acting on the giant planets as these are difficult to estimate without hydrodynamical simulations (see e.g., Moeckel et al. 2008; Marzari et al. 2010).

thumbnail Fig. 14

Distribution of the number of surviving terrestrial planets per system for the unstable mixed (grey) and gas (dashed line) simulations.

Open with DEXTER

The effect of the giant planets on terrestrial planet formation was similar for the gas and mixed simulations. Figure 13 shows that the sculpting of the terrestrial zone by the giant planets – as measured by the minimum giant planet perihelion distance during the simulation – is the same for the unstable mixed and the gas simulations. The only slight difference comes from two gas simulations in which a giant planet came closer than 1 AU to the star (in one case for a prolonged period of almost 1 Myr) but that succeeded in forming a terrestrial planet. In both of these cases the surviving planet was roughly a Mars mass (although in one case the planet accreted another embryo) and underwent large-scale oscillations in eccentricity and inclination. The survival of these planets may be due in part to gas drag from the disk, although they may simply represent a tail of the mixed distribution.

Figure 14 shows that the distributions of the number of surviving terrestrial planets in the gas and the unstable mixed systems are very close (1σ is 5–8% for most bins as calculated using binomial statistics; see Burgasser et al. 2003). The gas and unstable mixed simulations each formed a mean of 1.2 planets planets per unstable system, although there is a slight tendency for more zero- and one-planet systems in the gas simulations. We attribute this to the fact that in the gas simulations the giant planets started on orbits that were already strongly unstable and the instability involved the innermost giant planet with the strongest influence on the terrestrial planet zone. Thus, outward-directed instabilities and weakly unstable systems were less frequent in the gas simulations.

Although the number of surviving planets was similar, the surviving terrestrial planets in the unstable mixed simulations were significantly more massive than in the gas simulations. The median terrestrial planet mass was 0.43   M for gas and 0.73   M for mixed (counting only simulations that were integrated for  >100 Myr and planets  >0.1   M). In addition, 28 of the 90 unstable mixed terrestrial planets were more than 1   M () compared with 4 of 33 () for the gas terrestrial planets. The larger masses of the unstable mixed simulations come from the contribution from weakly unstable systems, i.e. those with minimum giant planet perihelion distances larger than about 4 AU in Fig. 13. For systems where a giant planet entered within 4 AU of the star, the two sets had the same median terrestrial planet mass.

thumbnail Fig. 15

The fraction of systems that contain  ≥ 0.5    M in surviving terrestrial planets as a function of F/Fstar(70   μm) after 1 Gyr for the unstable mixed (black) and gas (grey) sets of simulations. Error bars are 1σ values calculated with binomial statistics.

Open with DEXTER

A small fraction of unstable systems produced asteroid belts without terrestrial planets. In these systems a number of rocky planetesimals were the only survivors in the inner planetary system, as all terrestrial embryos had been destroyed. This occurred in 14 of 299 (4.7 ± 1.2%) of unstable simulations across all the sets of simulations (excluding the lowmass simulations). The gas simulations had a slightly higher rate of production of asteroid belt-only systems (), presumably because in a few cases gas drag was able to stabilize the orbits of planetesimals that were marginally unstable. These asteroid belts are low-mass, containing only 1–20 asteroid particles (each 5 × 10-3   M), albeit typically on excited orbits with moderate eccentricities (e ~ 0.2−0.3) and inclinations (i ~ 20−30°). Given their rapid collisional evolution, these belts probably become quickly dominate by a few relatively large objects and are probably not detectable with current instruments.

The surviving terrestrial planets in the gas and the unstable mixed simulations had similar median eccentricities e and inclinations i (median values of e ≈ 0.1 and i ≈ 5−6°). However, the gas simulations tend to have significantly higher oscillation amplitudes in e and i. Although the median oscillation amplitudes are relatively close (median peak-to-peak eosc = 0.13 and iosc = 7° for gas vs. 0.10 and 5° for mixed), planets in the gas simulations are shifted to higher values. Again, this difference is simply due to the lack of weakly unstable systems in the gas simulations; when a cut in the minimum giant planet perihelion of  <4 AU is applied the oscillation amplitudes are a match.

The anti-correlation between the giant planet eccentricity and the dust flux at 70   μm is very similar between the gas and unstable mixed simulations. The positive correlation between the mass in surviving terrestrial planets and dust flux is also preserved in the gas simulations. Figure 15 shows the fraction of systems that formed at least 0.5   M in terrestrial planets as a function of F/Fstar(70   μm) at 1 Gyr. The two distributions are almost identical. The only slight divergence is at small F/Fstar values, where the mixed simulations are about 1σ higher than the gas simulations. This is explained by the fact that the gas simulations are inward-directed by design, because the instability is triggered by a close encounter between the inner two giant planets (simply because our initial conditions put the middle giant planet on an orbit that crosses the inner one’s). In contrast, the mixed instabilities include both inward- and outward-directed instabilities, i.e., instabilities that can be triggered in, and largely confined to, either the inner or outer parts of the system. Inward-directed instabilities that perturb the outer planetesimal disk are necessarily very strong, somewhat stronger than equivalent outward-directed instabilities that perturb the outer disk. This only appears to be important (and then only at 1σ) for the lowest values of F/Fstar, where the disk is most strongly depleted (note that the lowest bin includes systems with F/Fstar < 0.01). Thus, the presence of disk gas at the time of giant planet instabilities does not appear to have a significant effect on the debris disk-terrestrial planet correlation.

We conclude that there are no strong systematic differences between the gas and the unstable mixed sets of simulations.

6. Implications of dynamically active giant planets

We now explore the implications of our results for expected correlations between extra-solar terrestrial planets, giant planets and debris disks. We emphasize that our conclusions are of course determined in part by our chosen initial conditions, which are poorly-constrained observationally. Several aspects of the initial conditions have considerable uncertainties: 1) the mass and mass distribution in the terrestrial planet zone, 2) the mass, mass distribution and extent of the outer planetesimal disk, 3) the number, masses and initial spacing of the giant planets, and 4) the relative masses and spacing of these three components. We have chosen what we consider to be reasonable values of the different components, but these values certainly vary from system to system and cover a far wider range than the subset included here. In addition, some systems with qualitatively different properties probably exist, such as disks with very widely-spaced giant planets or in which planetesimals only form in narrow regions. We discuss the impact of these assumptions on our results in Sect. 6.4 below.

With these limitations in mind, we now perform a simple experiment based on the results of our simulations. The goal of the experiment is to use the giant planets to match the observed exoplanet mass and eccentricity distributions and to then test different correlations within that framework. We construct two samples of systems that provide an adequate match to observations using simple, non-fine-tuned mixes of different sets of simulations. We then explore the implications of those samples.

6.1. Observational constraints

Any sample we construct is constrained by observations of giant exoplanets, debris disks, and correlations between those two. We now summarize these key characteristics.

First, we are constrained by the distribution of known giant exoplanets, particularly those beyond 0.2 AU that have presumably not been affected by tidal interactions with their host stars. The mass distribution can be fit with a simple power law: dN/dM ~ M-1.1 (Butler et al. 2006; Udry & Santos 2007). This is the mass distribution of the entire sample. It is likely similar to the ensemble-averaged mass distribution of planets prior to scattering (with some modification due to mass-dependent planetary ejections and collisions with the star) but, as discussed below, it need not be the mass distribution prior to scattering in any individual system (i.e. on a system-by-system level, planet masses may be correlated). The frequency of giant planets is very low (<1%) close-in, but increases sharply at 0.5–1 AU and appears to remain at a high level out to at least 3 AU (Cumming et al. 2008; Mayor et al. 2011). The median giant exoplanet eccentricity is  ~0.25 and the distribution extends to above 0.9 (Butler et al. 2006; Udry & Santos 2007). The eccentricity distribution is independent of orbital distance (for planets not affected by tides Ford & Rasio 2008). In addition, observations show that more massive giant planets (Mp ≳ MJ) have higher eccentricities than lower-mass giants (Mp ≲ MJ; Jones et al. 2006; Ribas & Miralda-Escudé 2007; Ford & Rasio 2008; Wright et al. 2009).

thumbnail Fig. 16

The dust-to-stellar flux ratio at 70   μm after 1 Gyr of evolution as a function of the eccentricity of the innermost giant planet eg,in for cases A and B. Systems in which the innermost giant planet is exterior to 8 AU have been excluded from this plot.

Open with DEXTER

Second, we are constrained by debris disk statistics, and we focus on observations at 70   μm made primarily with Spitzer. Observations show that % of solar-type stars have detectable dust emission at 70   μm (Trilling et al. 2008). There is no observed variation in this fraction with age, although the upper envelope of actual fluxes decreases for stars older than 1 Gyr or so (Hillenbrand et al. 2008; Carpenter et al. 2009).

Finally, we are constrained by any connection that might exist between the presence of giant exoplanets and debris disks. Debris disks have been detected around more than 20 stars with known exoplanets. However, there is currently no correlation between the presence of planets and debris disks: the incidence of debris disks is about 15% for both stars with and without planets (see Table 1 in Kóspál et al. 2009; Moro-Martín et al. 2007; Bryden et al. 2009). We also note that the strong observed correlation between the fraction of stars with currently-known exoplanets (hot Jupiters in particular) and the stellar metallicity (Gonzalez 1997; Santos et al. 2001; Fischer & Valenti 2005) is not apparent in the sample of known debris disks (Beichman et al. 2006; Greaves et al. 2006; Bryden et al. 2006; Kóspál et al. 2009).

6.2. Consistency with known giant planet properties

The masses of individual giant planets, as well as the mass ratio between planets in a given system, are the most important factor governing the outcome of planet-planet scattering (Raymond et al. 2010). Equal-mass giant planets provide the strongest instabilities for a given mass, and the most eccentric surviving planets (Ford et al. 2003; Raymond et al. 2010). Scattered equal-mass planets are also more widely-spaced than planets with mass ratios of a few (Raymond et al. 2009a, 2010). And among equal-mass unstable systems, more massive planets yield larger eccentricities but smaller inclinations (Raymond et al. 2010). The dynamics of scattering is only weakly dependent on the planet masses; Neptune-mass planets at a few AU require far more close encounters to eject one another than Jupiter-mass planets but their final orbital distributions are similar. The number of giant planets also plays a role; in general, more giant planets lead to more scattering events and higher final eccentricities (Jurić & Tremaine 2008).

We construct two mixtures of our simulations to reproduce the observations:

  • Case A is based on the mixed simulations. The eccentricitydistribution of surviving giant planets in the simulations(considering just the innermost planet as it provides the closestmatch to radial velocity observations) provides a quantitativematch to the observed distribution with a probability value p of 0.49 calculated from a K-S test. The best match is found by including only unstable systems, but the p value is still an acceptable 5–25% if the giant planet sample includes a 5–10% contribution from stable systems, with a higher p for smaller contributions of stable systems (see Fig. 19 in Paper I). Case A includes a 10% contribution from stable systems. Note that, since case A is built on the mixed simulations and that several other sets of simulations share the same giant planet characteristics (smallseed, bigseed, widedisk, and gas), variations on Case A can be constructed by substituting a different set of simulations for the mixed set.

  • Case B is constructed from a combination of the equal and lowmass simulations. The simplest scenario to explain the observation that higher-mass planets have higher eccentricities is for massive planets to form in systems with multiple, roughly equal-mass planets (see Sect. 5 of Raymond et al. 2010). At low planet masses the eccentricities are larger than observed, so to balance the sample a contribution of systems with lower-mass (M ≤ MSat − MJup) planets with significant mass ratios is needed – these are represented with our lowmass set. The exact number of low-mass systems needed is poorly constrained. In practice, we divide the lowmass in two based on the mass of the innermost surviving giant planet (dividing at 50   M); case B includes an equal number of systems at M ≤ MJup from the equal and unstable lowmass simulations.

Cases A and B each provide a marginally acceptable match to the observed exoplanet eccentricity distribution. Case A matches the giant exoplanet distribution with p = 0.08 from K-S tests, and if we allow a 5–10% increase in the number of planets with e = 0 (as suggested by Zakamska et al. 2011) case B also matches the distribution, with p ≈ 0.1. Given the uncertainties in orbital fitting of exoplanet eccentricities (Shen & Turner 2008; Zakamska et al. 2011) we do not attempt to fine-tune our samples to better fit the observations. Case A naturally matches the observed mass distribution (except for a bias due to the mass dependence of ejected planets) and as the outcomes for the equal simulations were mostly mass-independent, weighting of different outcomes with the equal contribution is not necessary, and case B can be also considered to provide a match to the mass distribution (see Sect. 5 in Raymond et al. 2010,for more details). However, the two cases have different implications for the nature of planetary systems. In case A, all planetary systems experience the same qualitative evolution because nearly all of them become dynamically unstable. In case B, the evolution of systems is divided according to the planetary masses: high-mass planetary systems undergo extremely violent instabilities, but the evolution of lower-mass systems is much calmer and many such systems are dynamically stable. Current observations favor case B over case A because it reproduces the higher eccentricities of more massive planets (Ribas & Miralda-Escudé 2007; Wright et al. 2009) while for case A higher-mass planets have lower eccentricities than low-mass planets (Raymond et al. 2010).

6.3. Debris disk correlations in cases A and B

The expected trends – an anti-correlation between giant planet eccentricity and debris disks and a positive correlation between terrestrial planets and debris disks – are clearly seen in both the case A and B systems.

Figure 16 shows the dust flux at 70   μm after 1 Gyr vs. the innermost giant planet’s eccentricity for all the simulations in each case without the weighting described above (e.g., case B still only contains the unstable lowmass simulations with inner planet masses greater than 50   M but there is no weighting between the equal and lowmass components). The anti-correlation is clear in the case A simulations but the large intrinsic scatter in the lowmass component of case B makes the trend less evident for case B. The scatter is caused by the fact that the case B simulations are susceptible to planetesimal-driven radial migration that allows for a wide range in the depletion of the outer planetesimal disk depending on the outer planet’s orbital history.

The correlation between eg and F/Fstar(70   μm) is clearer when the data are binned. Figure 17 (left panel) shows that cases A and B are very similar in terms of the fraction of systems that is detectable at 70   μm as a function of eg: both cases show the expected clear anti-correlation between debris disks and eccentric giant planets including a rapid decrease in the detectable fraction for eg > 0.03−0.1. When considering the fraction of systems with eg > 0.1 as a function of F/Fstar(70   μm) (right panel of Fig. 17) there is an offset of roughly 1 − σ between the two cases that is again due to the larger inherent scatter in the case B simulations, in particular the high-mass, unstable component of the lowmass simulations.

thumbnail Fig. 17

Left: the fraction of systems that would be detectable with Spitzer (i.e., with F/Fstar(70   μm) ≥ 0.55 after 1 Gyr of collisional and dynamical evolution) as a function of the eccentricity of the innermost giant planet eg, for cases A (black) and B (grey). The error bars are based on binomial statistics (see Burgasser et al. 2003). This essentially represents a horizontal slice through Fig. 16. Note that cases A and B include all relevant simulations; there is no weighting of stable/unstable or equal/lowmass simulations. Right: the fraction of systems with eg > 0.1 as a function of F/Fstar(70   μm) ≥ 0.55 (1 Gyr), again for cases A and B. Systems with F/Fstar < 10-2 are included in the bin at F/Fstar ≈ 10-2. The Spitzer detection limit is shown as the dashed line. This represents a vertical slice through Fig. 16.

Open with DEXTER

We now turn to the debris disk – terrestrial planet correlation shown in Paper I. This correlation arises naturally from the destructive role of dynamically active giant planets for both terrestrial planet formation and the survival of outer planetesimal disks. Indeed, the giant planets’ role in terrestrial planet formation is almost purely antagonistic: giant planets may quench or stifle terrestrial planet formation but have not been shown to help it along in any significant way (see Paper I and also Chambers & Cassen 2002; Levison & Agnor 2003; Raymond et al. 2004; Raymond 2006; Raymond et al. 2006a, 2009c), except in some circumstances to promote runaway growth of planetesimals (Kortenkamp et al. 2001). Nonetheless, if we assume that giant planets generally form outside the snow line, moderately high eccentricities are needed before the impact on terrestrial planet formation becomes deleterious.

Figure 18 shows the terrestrial planet-debris disk correlation for cases A and B. As before, the correlation between F/Fstar(70   μm) after 1 Gyr and the total mass in surviving terrestrial planets is clearer for case A. Again, this comes from the low mass simulations’ much larger scatter in F/Fstar and somewhat lower typical values for F/Fstar.

When the data are binned, the terrestrial planet-debris disk correlation is clearer for both cases. Figure 19 (left panel) shows the fraction of systems with F/Fstar(70   μm) above the Spitzer detection limit as a function of the total terrestrial planet mass (i.e., a horizontal slice through Fig. 18). The correlations are similar for cases A and B, with only small differences for small terrestrial planet masses, for which a higher fraction of systems remains detectable for case A. We attribute this difference to the more dynamic evolution of case B systems, which deplete or destroy their outer planetesimal disks at a much higher rate than case A systems. Cases A and B follow nearly identical trends in terms of the fraction of systems with at least 0.5 M in terrestrial planets as a function of F/Fstar(70   μm) (i.e., a vertical slice through Fig. 18). The only modest discrepancy that at very low fluxes (F/Fstar(70   μm) ≲ 0.05) there are more case A systems with terrestrial planets. Once again, we attribute this to the more dynamic evolution of case B systems: the typical case B system that destroys its outer planetesimal disk is more likely to also destroy its inner rocky material than a corresponding case A system because in this regime case B is dominated by the very violently unstable equal systems.

thumbnail Fig. 18

The dust-to-stellar flux ratio at 70   μm after 1 Gyr of evolution as a function of the total mass in surviving terrestrial planets for cases A and B. Again, the solar system is represented by the grey star.

Open with DEXTER

thumbnail Fig. 19

Left: the fraction of systems that would be detectable with Spitzer (F/Fstar(70   μm) ≥ 0.55 after 1 Gyr) as a function of the total mass in surviving terrestrial planets for cases A (black) and B (grey). This represents a horizontal slice through Fig 18. Right: the fraction of systems with 0.5   M or more in terrestrial planets as a function of F/Fstar(70   μm) ≥ 0.55 (1 Gyr) for cases A and B. Systems with F/Fstar < 10-2 are included in the bin at F/Fstar ≈ 10-2. The Spitzer detection limit is shown as the dashed line. This represents a vertical slice through Fig. 19.

Open with DEXTER

To summarize, we conclude from Figs. 1618 that the debris disk – eccentric giant planet anti-correlation and the debris disk – terrestrial planet correlation are clear in both cases A and B.

6.4. Discussion

Could our two predicted correlations be artifacts of our initial conditions? Is there any reasonable scenario that could remove these correlations?

The debris disk-eccentric giant planet anti-correlation exists because giant planet instabilities tend to clear out outer planetesimal disks, mainly by dynamical ejection. For eccentric giant planets not to be anti-correlated with debris disks, something fundamental about our proposed scenario must change. To start with, other mechanisms have been proposed to explain the large eccentricities of the observed exoplanets (for an exhaustive list, see Sect. 1 of Ford & Rasio 2008). However, to date the planet-planet scattering model is the only mechanism that has been shown to be physically viable and to fully reproduce the currently-observed characteristics of the giant exoplanet population (Ford & Rasio 2008; Chatterjee et al. 2008; Jurić & Tremaine 2008; Raymond et al. 2010).

The initial conditions of our simulations certainly play a role in producing the correlations we find; for example, if there exist wide radial gaps between the giant planets and outer planetesimal disks then the giant planets’ influence on debris disks would be weaker. The giant planets would appear to be basically irrelevant for the existence of debris disks. However, even a distant giant planet has a destructive influence on a planetesimal disk by increasing eccentricities sufficiently that inter-particle collisions become destructive (Mustill & Wyatt 2009). A simple fit to the observed distribution of debris disks using a self-stirred model like the one presented in Sect. 2 shows that outer planetesimal disks are typically located at 15–120 AU (Kennedy & Wyatt 2010). Thus, it may be possible for debris disks and giant planets at a few AU or less to appear uncorrelated (as is the case in the current sample Bryden et al. 2009; Kóspál et al. 2009) but the anti-correlation between eccentric giants and debris disks should be clear for giant planets beyond roughly 5–10 AU. In addition, there should be a minimal orbital radius (and a maximum temperature) for an outer planetesimal disk in a given system, depending on the giant planet architecture. The statistics of debris disks in known giant planet systems are currently too poor to provide constraints (see discussion in Sect. 5.2 of Paper I).

The debris disk-terrestrial planet correlation exists for three reasons. First, terrestrial planet formation is less efficient in the presence of eccentric giant planets. Second, the destruction rate of outer planetesimal disks by ejection increases under the influence of eccentric giant planets. Third, we have assumed that any system that can form terrestrial planets can also produce a debris disk.

The first reason is well-established from many dynamical studies (Chambers & Cassen 2002; Levison & Agnor 2003; Raymond et al. 2004; Raymond 2006; Raymond et al. 2006a, 2009c; Morishima et al. 2010). As discussed above, the dynamics behind the second reason is sound and the only reasonable way to negate or temper the eccentric giant planet – debris disk anti-correlation is for there to exist a large radial gap between giant planets and outer planetesimal disks. The third reason is hard to constrain. As a blanket statement this assumption is almost certainly incorrect, as the observed frequency of debris disks of  ~16% (Trilling et al. 2008; Carpenter et al. 2009) is smaller than the currently-estimated frequency of close-in super Earths of 20–50% (Howard et al. 2010, 2011; Mayor et al. 2011). Thus, there are probably many systems that can form terrestrial planets but without a massive outer planetesimal disk. The reason for the lower frequency of outer disks is unclear; it could be related to external perturbations from passing stars during the embedded cluster phase (Malmberg et al. 2007, 2011). Or perhaps outer planetesimal disks simply are not a common initial condition for planet formation, although that assertion appears to be at odds with the high frequency of disks around young stars (e.g. Hillenbrand et al. 2008). One can imagine that there could also exist systems with outer planetesimal disks but very little mass in the inner disk, for example if a giant planet migrated through and depleted the inner disk (although Raymond et al. 2006b; Mandell et al. 2007, shows that this depletion is much weaker than one would naively expect). Indeed, in our own solar system Jupiter’s inward-then-outward migration may have removed most of the mass from the asteroid belt (Walsh et al. 2011). However, there is no evidence in the exoplanet distribution for systematic depletion of inner disks by giant planet migration. Indeed, the radial distribution of giant exoplanets increases sharply beyond about 1 AU (Mayor et al. 2011) such that terrestrial planets closer than 0.5–0.7 AU to their stars are probably only weakly affected by giant planets, at least those with low to moderate eccentricities. In contrast, inner disks certainly have a wide mass range and high-mass disks may form super Earths or Neptune-like planets rather than Earth-like planets (Ikoma et al. 2001; Raymond et al. 2008b).

We think that the most likely interpretation of current observations and theory is roughly as follows. Protoplanetary disks start with a variety of masses and mass distributions, and are subsequently divided into different regions by any giant planets that form. Formation models suggest that giant planets may preferentially form at a few to ten AU (Kokubo & Ida 2002; Thommes et al. 2008; Levison et al. 2010), essentially dividing their disks into the distinct regions we have assumed: the inner terrestrial zone, the giant planet zone and the outer planetesimal disk. In some systems these zones are not cleanly separated, for example if the giant planets migrate very far outward or inward or if a relatively close stellar encounter disrupts the outer planetesimal disk. Indeed, observations suggest that there is probably an abundance of systems with inner terrestrial planet-forming disks but without the outer planetesimals to produce debris disks. However, we see no evidence or clear theory to contradict our assumption that all or at least most stars with debris disks also have inner disks of protoplanets. Thus, we think that debris disks can indeed act as signposts for systems that should have formed terrestrial planets. In contrast, in systems without debris disks eccentric giant planets may act as signposts of terrestrial planet destruction.

To summarize, we see no compelling argument against the assumptions we have made and we think that our two proposed correlations are reasonable.

7. Conclusions

In Paper I (Raymond et al. 2011) we showed that old solar-type stars with bright cold dust correlate strongly with dynamically calm environments that are conducive to efficient terrestrial accretion. The fact that both the inner and outer parts of planetary systems are sculpted by what lies in between – the giant planets – yields a natural connection between terrestrial planet and debris disks. We predicted two observational correlations: an anti-correlation between eccentric giant planets and debris disks, and a positive correlation between terrestrial planets and debris disks. We also showed that the solar system appears to be a somewhat unusual case in terms of having a rich system of terrestrial planets but a severely depleted Kuiper belt with little cold dust, which is probably a result of the outward-directed, relatively weak dynamical instability that caused the late heavy bombardment (Morbidelli et al. 2010).

In this paper we used six new sets of simulations to test the effect of several system parameters on the results from Paper I. None of the parameters qualitatively changed the eccentric giant planet-debris disk anti-correlation or the debris disk-terrestrial planet correlation but their effects were diverse and interesting:

  • Low-mass giant planets undergo planetesimal-driven migrationthat radially spreads the giant planets (the lowmass simulations).Systems with low-mass giant planet are very efficient at formingterrestrial planets and also at producing abundant observabledust simply because their weak giant planet perturbations do notstifle these processes. Systems with low-mass outer planets oftenproduce isolated belts of planetesimals orbiting between thegiant planets. The probability of a system containing such a beltincreases strongly with decreasing outer giant planet mass andincreasing separation between the giant planets.

  • In contrast, systems with equal-mass giant planets undergo the strongest instabilities. In more than 2/3 of the equal systems, all terrestrial material was destroyed, and all planetesimal disks were ejected in a similar fraction of cases (though not with a one to one correspondence). Given the tendency for more massive exoplanets to have more eccentric orbits (Wright et al. 2009), the equal systems may be representative of the high-mass (Mp ≳ MJ) exoplanet systems.

  • The presence of 0.5−2   M objects in outer planetesimal disks cause the disks to radially spread (the seeds simulations). The inward-spreading part of the disk is removed by interactions with giant planets (or, in systems without giant planets, presumably by accelerated collisional evolution) whereas the outer part of the disk spreads to larger orbital distances. Such cold disks are detectable at long but not short wavelengths. The presence of seeds may therefore explain the very low frequency of stars with 25   μm excess compared with the frequency of 70   μm excesses (Bryden et al. 2006; Trilling et al. 2008; Carpenter et al. 2009).

  • The widedisk simulations showed that a more massive and extended outer planetesimal disk stabilizes a significant fraction of systems because planet-planetesimal interactions allow higher-mass unstable giant planet systems to radially spread out and damp their eccentricities. Thus, if all systems of giant planet systems form in similar, near-unstable configurations, the outer planetesimal disk mass may be a key factor in the fraction of systems that end up being unstable. More massive outer disks produce larger quantities of dust, in particular at long wavelengths that are dominated by the outermost part of the planetesimal disks where the collisional timescales are the longest.

  • The presence of disk gas for the first 0.5 Myr of the system’s evolution does not have a strong effect on the outcome (the gas simulations). The debris disk-terrestrial planet correlation and the eccentric giant planet-debris disk anti-correlation are not affected.

We constructed two samples of simulations called cases A and B that matched the observed exoplanet mass and eccentricity distributions. Case A was built on the mixed simulations and case B on a combination of the equal and the high-mass unstable component of the lowmass simulations. These cases attempt to bracket the likely initial conditions for planet-planet scattering in the known systems (e.g., Chatterjee et al. 2008; Jurić & Tremaine 2008; Raymond et al. 2010). Each of these cases clearly show the eccentric giant planet-debris disk anti-correlation and the debris disk-terrestrial planet correlation.

The only plausible alternative to the eccentric giant planet-debris disk anti-correlation is for outer planetesimal disks to be located far away from the giant planets and thus to be extremely cold (see Sect. 6). Given that the known exoplanets are dominated by relatively close-in planets (typically within a few AU) and that the planetesimals responsible for the known debris disks are at 15–120 AU (Kennedy & Wyatt 2010), it is not surprising that the expected eccentric giant planet-debris disk anti-correlation has not yet been detected (Bryden et al. 2009, see also Sect. 5.2 of Paper I). However, we predict that it will be seen in upcoming datasets that include more distant planets.

The only plausible alternative to the debris disk-terrestrial planet correlation is for star with debris disks to be systematically depleted in their inner (terrestrial) regions and for outer planetesimal disks to be systematically cold enough to be uncorrelated with eccentric giant planets. Although it is easy to imagine scenarios in which this may occur – for example, a system in which an inward-migrating giant planet depletes the inner disk – there is no observational evidence or self-consistent theoretical model that suggests that such scenarios should be common or correlated with debris disks. Thus, we predict that upcoming datasets will find a strong correlation between debris disks and the presence of terrestrial planets.

The frequency of debris disks around Sun-like stars older than 1 Gyr is  ~16% (Trilling et al. 2008), which is lower than the observed frequency of close-in super Earth planets (Howard et al. 2010, 2011; Mayor et al. 2011). Thus, there is certainly a large population of systems that can form terrestrial planets but without the outer planetesimal disks needed to produce debris disks. In these systems, direct radial velocity observations or constraints on the giant planet architecture may offer the best insight for constraining the existence of terrestrial planets. But stars with debris disks most likely also had abundant material for building terrestrial planets. Thus, debris disks can indeed act as signposts of (past) terrestrial planet formation.

To sum up, our main result is a prediction that debris disks should be anti-correlated with systems containing eccentric giant planets and correlated with the presence of terrestrial planets. Solar-type stars with bright cold debris disks and no giant planets are excellent candidates to search for Earth-like planets. In contrast, systems without debris disks and with eccentric giant planets are probably not good candidates for terrestrial planets. Upcoming observations will test our predictions.


1

Note that the plateau in brightness seen in the 25   μm plot of Fig. 3 may be slightly overestimated because our method does not adequately account for collisional grinding of bodies on extremely close-in orbits when they are isolated particles. In this case, these four isolated particles were at close enough distances to dominate the flux at wavelengths shorter than  ~50 μm from about 80–120 Myr.

2

A somewhat higher fraction of the bigseed simulations were unstable compared with the mixed and smallseed simulations. This difference is only moderately statistically significant and might be due in part to a small glitch in our initial conditions for the seeds simulations: the inner edge of the outer planetesimal disk was always 4 Hill radii exterior to the outermost giant planet but the icy embryos were always initially between 10–20 AU. Thus, in many cases the innermost one to two seeds were in immediate dynamical contact with the outermost giant planet. This preferentially occurred when the outermost planet was very massive (and hence on a more distant initial orbit) and this glitch does not appear to have contaminated our results. In fact, the median instability time was later for the bigseed simulations than the mixed simulations, which is the opposite of what would be expected if the instabilities were systematically driven by seed-giant planet interactions.

3

We note that the unstable mixed systems with detectable 25   μm flux are all quite close to the detection limit (Fig. 8), and the fraction of unstable systems that is detectable at 25   μm drops drastically to 2/96 = after 3 Gyr of evolution, and the ratio between the detectable frequency at 70   μm to 25   μm increases to 25.5, well within the range allowed by observations. However, what is lacking in the mixed simulations is the ability to account for stable systems without 25   μm flux, as the fraction of stable systems with detectable 25   μm flux remains very high, at .

Acknowledgments

We thank the referee, Hal Levison, for a careful and thorough report that helped us improve the paper. Simulations were run at Weber State University and at Purdue University (supported in part by the NSF through TeraGrid resources). S.N.R. acknowledges the CNRS’s PNP and EPOV programs, the Conseil Regional d’Aquitaine, and NASA Astrobiology Institute’s Virtual Planetary Laboratory lead team. P.J.A. acknowledges funding from NASA’s Origins of solar systems program (NNX09AB90G), NASA’s Astrophysics Theory program (NNX11AE12G), and the NSF’s Division of Astronomical Sciences (0807471). F.S. acknowledges support from the European Research Council (ERC Grant 209622: E3ARTHs). This paper is dedicated to S.N.R.’s son Zachary Max Raymond, whose birth on June 14, 2011 caused a long but absolutely worthwhile delay of the publication of this paper.

References

Online material

Movie of Figure 1 (Access here)

Movie of Figure 3 (Access here)

Movie of Figure 7 (Access here)

All Tables

Table 1

Summary of the simulations.

All Figures

thumbnail Fig. 1

Evolution of a reference mixed simulation in which the giant planets underwent a violent instability after 42.8 Myr of evolution. The initial giant planet masses were, in order of increasing orbital distance, 0.96, 0.46 and 0.64   MJ. Left: orbital eccentricity vs. semimajor axis of each body in the simulation. The size scales with the mass1/3 and the color corresponds to the water content, using initial values taken to solar system data (Raymond et al. 2004) and re-calculated during impacts by mass balance. The giant planets are shown as the large black bodies and are not on the same size scale. Top right: the spectral energy distribution of the dust during five simulation snapshots, as observed from a distance of 10 pc. The dashed line represents the stellar photosphere. Bottom right: the ratio of the dust-to-stellar flux F/Fstar at 25   μm as a function of time. The rough Spitzer observational limit is shown with the dashed line (Trilling et al. 2008). A movie of this simulation is available online or at http://www.obs.u-bordeaux1.fr/e3arths/raymond/scatterSED_199.mpg.

Open with DEXTER
In the text
thumbnail Fig. 2

Eccentricity vs. semimajor axis for the surviving giant planets in the mixed (left panel) and lowmass (right) simulations. Black circles represent unstable simulations – defined as simulations in which at least one giant planet-planet scattering event occurred – and grey dots are stable simulations. The size of each circle is proportional to the logarithm of the planet mass. The dashed vertical lines represent the median outer edge of the planetesimal disk for each set of simulations. Note that the outer edge varied from simulation to simulation depending on the giant planet masses, from 20.6 to 27 AU (with a median of 23.7 AU) for the mixed simulations and 17.3 to 23.3 AU for the lowmass simulations (median of 19.7 AU).

Open with DEXTER
In the text
thumbnail Fig. 3

Evolution of a simulation with low-mass giant planets. Left: orbital eccentricity vs. semimajor axis of each body in the simulation. The size scales with the mass1/3 and the color corresponds to the water content, using initial values taken to solar system data (Raymond et al. 2004) and re-calculated during impacts by mass balance. The giant planets are shown as the large black bodies and are not on the same size scale. Top right: the spectral energy distribution of the dust during five simulation snapshots, as observed from a distance of 10 pc. The dashed line represents the stellar photosphere. Bottom right: the ratio of the dust-to-stellar flux at 25 microns as a function of time. The rough observational limit of the MIPS instrument on NASA’s Spitzer Space Telescope is shown with the dashed line (Trilling et al. 2008). A movie of this simulation is available online or at http://www.obs.u-bordeaux1.fr/e3arths/raymond/scatter_lowmass_54.mpg.

Open with DEXTER
In the text
thumbnail Fig. 4

Distribution of the number of surviving terrestrial planets for the unstable (grey) and stable (dashed) simulations in the lowmass set of simulations.

Open with DEXTER
In the text
thumbnail Fig. 5

In the lowmass planetary systems, the fraction of stable zones between pairs of outer giant planets that contain at least one planetesimal on a stable orbit as a function of the mass of the outermost giant planet (left panel) and the width of the stable zone Δ in units of mutual Hill radii (right panel). For the left panel, bins were evenly logarithmically spaced from 10   M to 1   MJ. For the right panel, bins were evenly spaced from Δ = 3.5 to 30. The error bars were calculated using binomial statistics following Burgasser et al. (2003). The dashed line marks the two planet stability boundary (Marchal & Bozis 1982; Gladman 1993).

Open with DEXTER
In the text
thumbnail Fig. 6

Distribution of the number of surviving terrestrial planets for the equal set of simulations. Systems with giant planets of MSat or larger are shown in grey, and systems with giant planets of 30   M are shown by the black dashed line.

Open with DEXTER
In the text
thumbnail Fig. 7

Evolution of a simulation that includes 5 seed embryos (shown in grey) in the outer planetesimal disk. Left: orbital eccentricity vs. semimajor axis of each body in the simulation. The size scales with the mass1/3 and the color corresponds to the water content, using initial values taken to solar system data (Raymond et al. 2004) and re-calculated during impacts by mass balance. The giant planets are shown as the large black bodies and are not on the same size scale. The two surviving terrestrial planets have a = 0.71 and 1.27 AU, m = 1.3 and 1.58 M, and Myr-averaged e = 0.09 and 0.17, respectively. Both terrestrials have substantial water contents accreted from hydrated asteroidal material. The two surviving giant planets have a = 4.1 and 24.4 AU, m = 1.27 and 0.45 MJ, and Myr-averaged e = 0.32 and 0.33, respectively. Top right: the spectral energy distribution of the dust during five simulation snapshots, as observed from a distance of 10 pc. The dashed line represents the stellar photosphere. Bottom right: the ratio of the dust-to-stellar flux F/Fstar at 25   μm as a function of time. The rough Spitzer observational limit is shown with the dashed line (Trilling et al. 2008). An animation of this simulation is available online or at http://www.obs.u-bordeaux1.fr/e3arths/raymond/scatterSED_seed13.mpg.

Open with DEXTER
In the text
thumbnail Fig. 8

Correlations with the dust flux after 1 Gyr of dynamical and collisional evolution for the seeds simulations, as compared with the mixed simulations. The top panels show F/Fstar at 70   μm vs. the eccentricity of the innermost surviving giant planet (left) and the total mass in surviving terrestrial planets (right). The bottom panels show the same comparisons but at 25   μm. The mixed simulations are in black, the smallseed simulations in green, and the bigseed simulations in red. Filled circles represent stable simulations and open circles unstable simulations.

Open with DEXTER
In the text
thumbnail Fig. 9

Left: the fraction of systems that would be detectable with Spitzer (with F/Fstar(70   μm) ≥ 0.55 after 1 Gyr of collisional and dynamical evolution) as a function of the eccentricity of the innermost giant planet eg for the mixed and combined seeds simulations. The error bars are based on binomial statistics  (see Burgasser et al. 2003). This essentially represents a horizontal slice through the top left panel of Fig 8. Right: the fraction of systems with 0.5   M or more in surviving terrestrial planets as a function of F/Fstar(70   μm) ≥ 0.55 (1 Gyr) for the mixed and seeds simulations. Systems with F/Fstar < 10-2 are included in the bin at F/Fstar ≈ 10-2. The Spitzer detection limit is shown as the dashed line. This represents a vertical slice through the bottom right panel of Fig. 8.

Open with DEXTER
In the text
thumbnail Fig. 10

Orbital evolution of a widedisk simulation that underwent a Nice model-like instability. The evolution of the semimajor axes, perihelion and aphelion distances of the three giant planets are shown in black, and the semimajor axes of the 19 surviving outer planetesimals in grey (note that one additional planetesimal survived with a = 59 AU). Two terrestrial planets formed in this system, at 0.71 AU (1.4   M) and 1.76 AU (0.13   M): accretion was perturbed by the perturbations moderately eccentric inner giant planets at 4.7 AU (0.78   MJ, e oscillates between 0.03 and 0.17) and 7.8 AU (2.4   MJ, e oscillates between 0.06 and 0.1).

Open with DEXTER
In the text
thumbnail Fig. 11

Left: the dust-to-stellar flux ratio at 1 Gyr at 70   μm vs. the innermost giant planet’s eccentricity for the widedisk (grey) and mixed (black) simulations. Filled circles represent stable simulations and open circles unstable ones. The solar system is shown with the grey star. Right: histogram of the fraction of systems that are detectable at 70   μm as a function of the innermost giant planet’s eccentricity; this is essentially a horizontal slice through the left panel.

Open with DEXTER
In the text
thumbnail Fig. 12

Left: the dust-to-stellar flux ratio at 1 Gyr at 70   μm vs. the total terrestrial planet mass for the widedisk (grey) and mixed (black) simulations. Filled circles represent stable simulations and open circles unstable ones. The solar system is shown with the grey star. Right: histogram of the fraction of systems that contain 0.5   M or more in terrestrial planets as a function of F/Fstar(70   μm) after 1 Gyr; this is essentially a vertical slice through the left panel.

Open with DEXTER
In the text
thumbnail Fig. 13

Total mass in surviving terrestrial planets as a function of the minimum perihelion distance of any giant planet during the simulation for the unstable mixed (black dots) and the gas (grey dots) simulations.

Open with DEXTER
In the text
thumbnail Fig. 14

Distribution of the number of surviving terrestrial planets per system for the unstable mixed (grey) and gas (dashed line) simulations.

Open with DEXTER
In the text
thumbnail Fig. 15

The fraction of systems that contain  ≥ 0.5    M in surviving terrestrial planets as a function of F/Fstar(70   μm) after 1 Gyr for the unstable mixed (black) and gas (grey) sets of simulations. Error bars are 1σ values calculated with binomial statistics.

Open with DEXTER
In the text
thumbnail Fig. 16

The dust-to-stellar flux ratio at 70   μm after 1 Gyr of evolution as a function of the eccentricity of the innermost giant planet eg,in for cases A and B. Systems in which the innermost giant planet is exterior to 8 AU have been excluded from this plot.

Open with DEXTER
In the text
thumbnail Fig. 17

Left: the fraction of systems that would be detectable with Spitzer (i.e., with F/Fstar(70   μm) ≥ 0.55 after 1 Gyr of collisional and dynamical evolution) as a function of the eccentricity of the innermost giant planet eg, for cases A (black) and B (grey). The error bars are based on binomial statistics (see Burgasser et al. 2003). This essentially represents a horizontal slice through Fig. 16. Note that cases A and B include all relevant simulations; there is no weighting of stable/unstable or equal/lowmass simulations. Right: the fraction of systems with eg > 0.1 as a function of F/Fstar(70   μm) ≥ 0.55 (1 Gyr), again for cases A and B. Systems with F/Fstar < 10-2 are included in the bin at F/Fstar ≈ 10-2. The Spitzer detection limit is shown as the dashed line. This represents a vertical slice through Fig. 16.

Open with DEXTER
In the text
thumbnail Fig. 18

The dust-to-stellar flux ratio at 70   μm after 1 Gyr of evolution as a function of the total mass in surviving terrestrial planets for cases A and B. Again, the solar system is represented by the grey star.

Open with DEXTER
In the text
thumbnail Fig. 19

Left: the fraction of systems that would be detectable with Spitzer (F/Fstar(70   μm) ≥ 0.55 after 1 Gyr) as a function of the total mass in surviving terrestrial planets for cases A (black) and B (grey). This represents a horizontal slice through Fig 18. Right: the fraction of systems with 0.5   M or more in terrestrial planets as a function of F/Fstar(70   μm) ≥ 0.55 (1 Gyr) for cases A and B. Systems with F/Fstar < 10-2 are included in the bin at F/Fstar ≈ 10-2. The Spitzer detection limit is shown as the dashed line. This represents a vertical slice through Fig. 19.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.