Free Access
Issue
A&A
Volume 530, June 2011
Article Number A62
Number of page(s) 23
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/201116456
Published online 12 May 2011

© ESO, 2011

1. Introduction

Circumstellar disks of gas and dust are expected to produce three broad classes of planets in radially-segregated zones (Kokubo & Ida 2002). The inner disk forms terrestrial (rocky) planets because it contains too little solid mass to rapidly accrete giant planet cores, which are thought to form preferentially beyond the snow line where the surface density in solids is higher because of ice condensation and the isolation mass is larger (Lissauer 1987). Terrestrial planets form in 10−100 million years (Myr) via collisional agglomeration of Moon- to Mars-sized planetary embryos and a swarm of 1 − 103 km sized planetesimals (Chambers 2001; Kenyon & Bromley 2006; Raymond et al. 2006a, 2009c; O’Brien et al. 2006). From roughly a few to a few tens of AU, giant planet cores grow and accrete gaseous envelopes if the conditions are right (Pollack et al. 1996; Alibert et al. 2005). Despite their large masses, gas giants must form within the few million year lifetime of gaseous disks (Haisch et al. 2001; Pascucci et al. 2006; Kennedy & Kenyon 2009) and be present during the late phases of terrestrial planet growth. Resonant interactions, arising not just from giant planets but also from the changing surface density of the gas disk itself (Nagasawa et al. 2005), thus likely play a role in terrestrial planet formation. Finally, in the outer regions of planetary systems the growth time scale exceeds the lifetime of the gas disk, and the end point of accretion is a belt of Pluto-sized (and smaller) bodies (Kenyon & Luu 1998).

Our understanding of planetary formation in these different zones is constrained by a variety of observations. The initial conditions of inner disks can be probed by observations of hot dust. Such observations have shown evidence for grain growth in individual protoplanetary disks (Bouwman et al. 2008), and larger samples have shown that hot dust disappears from these disks within a few Myr and that cooler dust takes longer to disappear (Haisch et al. 2001; Mamajek et al. 2004; Meyer et al. 2008). Hundreds of close-in planets – the outcome of planet formation – have been detected by radial velocity and transit observations, with masses down to a few M (e.g., Butler et al. 2006; Udry et al. 2007; Léger et al. 2009; Batalha et al. 2011). The frequency of planets on short-period orbits has been found to be a function of the planet mass; less massive planets are significantly more common than high-mass “hot Jupiters” (Mayor et al. 2009; Howard et al. 2010). The frequency of 3−10 M planets has been measured at 12%, and by extrapolation the frequency of 0.5−2 M is 23% (for orbital periods less than 50 days; Howard et al. 2010).

Despite the existence of hot Jupiters around  ~1% of solar-type stars (Howard et al. 2010), the vast majority of giant planets is found beyond 1 AU (Butler et al. 2006; Udry & Santos 2007). The absolute frequency of giant planets is poorly constrained; estimates range from 10% (Cumming et al. 2008) to more than 50% (Gould et al. 2010). These planets are characterized by their broad eccentricity distribution that includes several planets with e ≥ 0.9 Butler et al. (2006). This distribution is quantitatively reproduced if giant planets form in systems with multiple planets, and dynamical instabilities occurred in 70 − 100% of all observed systems (Chatterjee et al. 2008; Jurić & Tremaine 2008; Raymond et al. 2010). The onset of instability may be caused by the changing planet-planet stability criterion as the gas disk dissipates (Iwasaki et al. 2001), resonant migration (Adams & Laughlin 2003), or chaotic dynamics (Chambers et al. 1996). Whatever the trigger, the instability leads to a phase of planet-planet scattering and in most cases to the eventual removal of one or more planets from the system by collision or hyperbolic ejection (Rasio & Ford 1996; Weidenschilling & Marzari 1996). It is the surviving planets that match the observed distribution. The properties of the outer Solar System are consistent with the giant planets having formed in a similarly unstable configuration (Thommes et al. 1999), though the low eccentricities of Jupiter and Saturn (eJ ≈ eS ≈ 0.05), and the late timing of the Late Heavy Bombardment (Strom et al. 2005), hint that the instability that occurred in our Solar System was weak and occurred too late to affect terrestrial accretion (Morbidelli et al. 2010).

Sub-millimeter observations of dust disks around young stars provide information on the initial conditions in outer planet-forming disks, and by extrapolation to entire disks. These observations suggest that the typical protoplanetary disk has a mass of 0.001 − 0.1 Solar masses (Andrews & Williams 2005, 2007a; Eisner et al. 2008) and a radial surface density profile of roughly r − (0.5−1) (Mundy et al. 2000; Andrews & Williams 2007b; Isella et al. 2010). There appears to be a roughly linear trend between the dust mass and the stellar mass such that the typical protoplanetary disk contains a few percent of the stellar mass (Andrews & Williams 2007a), which has important implications for the planet frequency as a function of stellar mass (Ida & Lin 2005; Raymond et al. 2007; Greaves 2010).

Debris disks – warm or cold dust observed around older stars, typically at infrared wavelengths (λ ~ 10−100   μm) – probe outer disks after planet formation has completed. Debris disks provide evidence for the existence of leftover planetesimals because the lifetime of observed dust particles under the effects of collisions and radiation forces is far shorter than the typical stellar age, implying a replenishment via collisional grinding of larger bodies. Spitzer observations show that about 15% of solar-type stars younger than 300 Myr have significant dust excesses at 24   μm but that this fraction decreases to about 3% for stars older than 1 Gyr (Meyer et al. 2008; Carpenter et al. 2009). At 70   μm the fraction of stars showing significant excesses is roughly constant in time (at  ~ 16%) but the upper envelope of the distribution decreases with the stellar age (Trilling et al. 2008; Gáspár et al. 2009; Carpenter et al. 2009). A stars have a significantly higher fraction of dust excesses at both 24 and 70   μm and the excesses themselves are brighter than for FGK stars but the dust lifetime is shorter (Su et al. 2006). Very few debris disks are known around M dwarfs and the dust brightness is necessarily far lower than for higher-mass stars (Gautier et al. 2007).

Here we perform a large ensemble of N-body simulations to model the interactions between the different radial components of planetary systems: formation and survival of terrestrial planets, dynamical evolution and scattering of giant planets, and dust production from collisional grinding. By matching the orbital distribution of the giant planets, we infer the characteristics of as-yet undetected terrestrial planets in those same systems. We post-process the simulations to compute the spectral energy distribution of dust by treating planetesimal particles as aggregates in collisional equilibrium (Dohnanyi 1969) to calculate the incident and re-emitted flux (Booth et al. 2009). We then correlate outcomes in the different radial zones and link to two key observational quantities: the orbital properties of giant planets and debris disks.

The paper is structured as follows. In Sect. 2 we explain our choice of initial conditions, the integration method, and our debris disk calculations. In Sect. 3 we demonstrate the detailed evolution of a single simulation in terms of its dynamics, the formation of terrestrial planets, and the dust evolution. In Sect. 4 we present results from our fiducial set of 152 simulations and explore the correlations between terrestrial planet formation efficiency, giant planet characteristics and debris disks. In Sect. 5 we discuss the implications of our results: we scale our simulations to the known eccentricity-semimajor axis distribution of giant planets, discuss the observed debris disks in known exoplanet systems, and evaluate the Solar System in the context of our results. We conclude in Sect. 6.

In a companion paper (Raymond et al. 2011; referred to in the text as Paper II) we test the effects of several system parameters on this process: the giant planet masses and mass distribution, the width of the outer planetesimal disk, the mass distribution of the outer planetesimal disk, and the presence of disk gas at the time of giant planet instabilities. We also calibrate our results to the statistics of known debris disks to produce an estimate of the fraction of solar-type stars that host terrestrial planets (ηEarth in the Drake equation).

2. Methods

We use numerical simulations to study the formation of terrestrial planets from massive embryos, the dynamical evolution of fully-formed gas and ice giants, and the long-term evolution of debris disks, starting from an assumed set of initial conditions that are specified at the epoch when the protoplanetary gas disk dissipates. The goal is to study how these processes, occurring in spatially distinct regions of the planetary system, are coupled, and to quantify the expected dispersion in the final outcome that arises from the chaotic nature of the evolution. In what follows, we describe in detail the specific initial conditions that we adopt, together with the numerical methods. It is worth emphasizing at the outset that our initial conditions are not unique. Observationally, little is known about the radial distribution of planetesimal formation within protoplanetary disks, and as a result the single most important initial condition for subsequent planet formation is not empirically well-constrained. Our fundamental assumption is that planetesimals form across a wide range of disk radii (extending out, in particular, to encompass a cold outer debris disk), with a smooth radial distribution. We further assume that cores are able to form quickly enough to yield planets at least as massive as ice giants in typical disks1. Given these assumptions, it is likely that the system-to-system variation in the masses of giant planets – arising ultimately from dispersion in the masses and radii of the protoplanetary disks – exceeds that of the terrestrial planets or that of the debris disk, because small variations in the formation time scale of cores result in large changes to the final envelope mass.

Our models are an extension of what might be described as the “classical” scenario for planet formation. Substantially different scenarios are also possible. In particular, planetesimal formation is poorly understood (for a review, see Chiang & Youdin 2010), and may be coupled to the level of intrinsic turbulence within the gaseous disk, which probably varies with radius (Gammie 1996; Armitage 2011). If planetesimal formation efficiency varies dramatically with radius, the initial conditions for subsequent planet formation would be radically different from what we have assumed, leading to qualitatively different conclusions. For example, it could be true that the zones of terrestrial and giant planet formation are typically dynamically well-separated, contrary to what is implied by our initial conditions. In such a model, the coupling between the giant and terrestrial planets would be much weaker than occurs given our assumptions, such that only exceptionally violent giant planet scattering – such as occurs when instabilities drive e → 1 – would affect the terrestrial zone. Similar caveats apply to the outer debris disk region.

2.1. Initial conditions

The initial conditions for planet formation are expected to vary with stellar mass, due both to the relatively well-understood variation of the location of the snow line with stellar type (Kennedy & Kenyon 2008), and as a result of any systematic variations in the star-to-disk mass ratio (Alexander & Armitage 2006). We consider Solar-mass stars, and start our simulations from highly idealized initial conditions that represent the predicted state of a planetary system at the time of the dissipation of the gaseous protoplanetary disk. There are three radial zones: an inner zone containing 9   M in planetary embryos and planetesimals from 0.5 to 4 AU, three giant planets on marginally stable orbits from Jupiter’s orbital distance of 5.2 AU out to  ~10 AU (depending on the masses), and an outer 10 AU-wide disk of planetesimals containing 50 M. In the majority of simulations the giant planets underwent a violent phase of scattering but in a significant fraction they did not. Figure 1 shows an example set of initial conditions.

thumbnail Fig. 1

An example of the initial conditions assumed to exist shortly after the epoch of gas disk dispersal. The terrestrial planet formation zone contains planetary embryos and lower mass planetesimals, while the outer planetesimal belt contains a population of relatively low-mass (compared to the giant planets) bodies that are represented numerically by N equal mass particles. The intermediate zone contains 3 fully-formed giant planets whose spacing is such that they are marginally unstable against dynamical instability.

Open with DEXTER

2.1.1. Terrestrial planet-forming region

The terrestrial planet forming region extends from 0.5 to 4 AU. As in the bulk of terrestrial planet formation simulations the inner boundary was chosen as a compromise between simulation run time and capturing the dynamics of the inner disk (e.g. Chambers 2001). The outer boundary corresponds to the approximate stability boundary with the innermost giant planet (at 5.2 AU), depending on that planet’s mass, such that objects beyond this limit would be immediately destabilized (Marchal & Bozis 1982; Gladman 1993). Within this zone, we adopt initial conditions that are motivated by “chaotic growth” models, of the type summarized by Goldreich et al. (2004) and calculated in detail by Kenyon & Bromley (2006). Specifically, the terrestrial zone contains a total of 9 M in 49 planetary embryos and 500 planetesimals, with equal mass in each component. The embryo mass of 0.09   M is a factor of a few to ten larger than that used in recent simulations of late-stage terrestrial accretion (Chambers 2001; Raymond et al. 2006a, 2009c; Morishima et al. 2010), but is within a factor of  ~2 of that obtained from semi-analytic models of oligarchic growth (Chambers 2006) over the radial range between 1 AU and 3 AU at 3 Myr. Of course, the outcome of late-stage accretion has been shown to be insensitive to the embryo mass distribution but rather to depend mainly on the large-scale mass distribution of the disk (Kokubo et al. 2006). The mass is distributed following a radial surface density profile Σ ∝ r-1. This profile is consistent with an inward extrapolation of that measured for DG Tau by Isella et al. (2010), who infer a slope for the disk surface density via interferometric observations of thermal dust emission at mm wavelengths. We do not consider the possibility of variations in the disk mass, due either to dispersion in the surface density of the protoplanetary disk at early times, or arising from dynamical effects such as the passage of massive planets migrating inward due to type II migration (Fogg & Nelson 2007; Raymond et al. 2006a; Mandell et al. 2007). Nor do we consider the effect of the changing gravitational potential of the dissipating gas disk, which in some cases could influence terrestrial accretion (Nagasawa et al. 2005; Morishima et al. 2010). Modeling these effects would require starting our calculations at an early epoch, when the gas disk was still present, and would introduce additional physical uncertainties. Our disk is comparable in mass to the “minimum-mass solar nebula” (MMSN) model (Weidenschilling 1977), with 2.6 [5.1] M inside 1.5 [2.5] AU, though we do not adopt the MMSN surface density slope.

The embryo mass increases only slightly with orbital distance, in agreement with models for embryo growth (Chambers 2006). In practice, embryos are spaced by Δ mutual Hill radii RH,m: (1)where a1 and a2 are orbital distances of two adjacent embryos, m1 and m2 are their masses, and M   is the stellar mass. We chose Δ to lie between 8 − 16 but decreasing systematically with orbital distance to avoid large variations in embryo mass in different regions of the disk. In practice, we used Δ = 8 + δ/a, where δ is randomly chosen between zero and 8 and a is the orbital distance in AU. The mass resolution of the terrestrial component of our model is within a factor of 2 − 4 of the best current simulations of terrestrial planet formation (Raymond et al. 2009c; Morishima et al. 2010). Embryos were given randomly-chosen initial eccentricities of up to 0.02 and initial inclinations of up to 0.1°. Planetesimals were given initial eccentricities of up to 0.02 and initial inclinations of up to 1°.

2.1.2. Giant planets

Our giant planet distribution is motivated by observations of extrasolar planetary systems, which show that the giant planet population detected at a < 5 AU exhibits a broad eccentricity distribution (Butler et al. 2006). The eccentricity distibution is not strongly affected by selection effects (Cumming et al. 2008), which do however dominate the radial distribution (which is essentially unknown beyond about 5 AU). There is no unique theoretical interpretation of this result, but the most common explanation is that the eccentricity represents the endpoint of a dynamical scattering phase. Models built upon this assumption (Jurić & Tremaine 2008; Chatterjee et al. 2008) can reproduce quantitatively the observed f(e) for extrasolar planets. When the effects of outer planetesimal disks are included, they are also consistent with near-circular orbits being typical for low mass giant planet systems at larger orbital radii (Raymond et al. 2010, as is found in the Solar System).

Here, we assume that marginally unstable initial conditions are also the rule for orbital separations modestly larger than those needed to explain the current observations of massive extrasolar planets. This may not be true, but it represents the simplest extrapolation of current observational results. We model systems of three giant planets, with the innermost giant planet given an orbital semimajor axis of 5.2 AU, the same as Jupiter’s. This is an arbitrary choice which allows for easy comparison with the Solar System (assuming that Jupiter formed at 5.2 AU). Two additional planets are spaced outward with randomly-chosen separations of (in the fiducial runs) 4 − 5   RH,m. This separation was chosen to select for systems that will likely become unstable on a timescale of 100 000 years or longer (Marzari & Weidenschilling 2002; Chatterjee et al. 2008). This value of 105 years is the expected timescale for the final dissipation of the gaseous protoplanetary disk (Wolk & Walter 1996; Ercolano et al. 2011), although substantial structural changes to the disk occur over longer timescales (Currie et al. 2009). Because damping within the disk tends to stabilize planetary systems (Iwasaki et al. 2002), we expect the natural instability timescale to be on the order of the gas dissipation timescale.

The outcome of scattering is dependent on both the initial mass distribution (which is related to the well-constrained observed mass distribution Butler et al. 2006), and on the distribution of mass ratios between planets in individual systems (Ford et al. 2003; Raymond et al. 2010). For the results shown in the paper we focus on a fiducial set of runs (called mixed in Paper II), in which the planet masses are chosen to follow the observed distribution of extra-solar planets, (2)where masses are chosen between one Saturn mass and 3 Jupiter masses. This range is chosen so as to include the majority of well-observed exoplanet systems (more massive planets are relatively uncommon, while selection effects become stronger for very low mass giants), but there is nothing particularly special about our choice. Similar models, that fit the observed constraints equally well, could almost certainly be constructed assuming different mass ranges, although it might be necessary to adopt different numbers of planets or to assume correlations between their masses (Raymond et al. 2010). In our runs, the masses of individual planets are chosen independently.

The radii of giant planets affect their early dynamics, by altering the ratio of physical collisions to scattering events. We adopt a constant bulk density, ρ = 1.3   g   cm-3, independent of mass, that matches that of the (present-day) Jupiter. Since young giant planets will assuredly have larger radii, our assumption underestimates the true number of physical collisions. We do not, however, consider this to be a major source of error. Planets in our mass range, formed via core accretion, are only modestly inflated at ages as short as a Myr (Marley et al. 2007), and collisions are suppressed at a > 5   AU because of the relatively modest orbital velocities at these radii (Ford et al. 2001).

As noted above, the initial conditions that we adopt for the giant planets are motivated by strictly empirical conditions: the observed mass function of extrasolar planets, and the requirement (within the context of a scattering model) that most multiple planet systems are eventually dynamically unstable. However, they are also broadly consistent with first principles models of giant planet formation (Bromley & Kenyon 2011), which suggest that radial migration and dynamical instability ought to be common, even prior to the dispersal of the gas disk. The specific separation of planets that we have assumed matches that expected if multiple giant planets interact with the gas disk such that they end up trapped in mean-motion resonances (Snellgrove et al. 2001; Lee & Peale 2002; Kley et al. 2004), provided that the majority of these resonances are broken before or shortly after the dispersal of the gas disk. Under some conditions fluctuating gravitational perturbations from disk turbulence may be able to break resonances (Adams et al. 2008; Lecoanet et al. 2009), though the ubiquity of this process remains unclear since the true strength of turbulence within disks cannot yet be reliably determined.

In Paper II we test the effect of the mass distribution of giant planets, including systems with much lower-mass giant planets that can undergo Nice model-like instabilities and systems with equal-mass giant planets that produce the strongest instabilities for a given planet mass (Raymond et al. 2010).

2.1.3. Outer planetesimal belt

In all cases, the outer planetesimal disk contains 50 M spread over a 10 AU-wide annulus. The inner edge of the annulus is chosen to be 4 linear Hill radii exterior to the outermost giant planet, such that the innermost planetesimals are “Hill stable” at the start of the simulation. This disk’s mass and mass distribution is comparable to that invoked by recent models of the evolution of the Solar System’s giant planets (the “Nice model”) to explain the late heavy bombardment (Tsiganis et al. 2005; Gomes et al. 2005). Note that, for numerical reasons and also to account for the much larger total mass, the mass of a “cometary” planetesimal particle in the outer system is roughly an order of magnitude larger than an “asteroidal” planetesimal in the inner system, although in both cases each planetesimal particle is assumed to represent an ensemble of smaller bodies (see below).

In Paper II we test the effect of varying certain properties of the outer planetesimal disk mass, including the width of the outer planetesimal disk and the presence of  ~ Earth-mass “seeds” in the disk.

2.2. Integration method

Our simulations were run using the hybrid integrator in the Mercury integration package (Chambers 1999), which combines the speed of a symplectic mapping scheme (Wisdom & Holman 1991) while particles are well-separated with the Bulirsch-Stoer method during closer encounters. We used a timestep of 6 days for each of our simulations. Particles were removed from a simulation when they attained perihelion distances of less than 0.2 AU, at which point they were assumed to collide with the star, or if they reached aphelion at more than 100 AU, when they were assumed to be ejected. Collisions were treated an inelastic mergers. We performed extensive numerical tests to validate our choice of timestep as well as the threshold integration error in orbital energy above which simulations were rejected as unreliable. These tests are described in Appendix A.

2.3. Debris disk modeling

To calculate the dust flux in a planetary system from a simulation we follow the procedure outlined in Sect. 2 of Booth et al. (2009) with only a few small modifications. This approach makes the assumption that planetesimal particles act as tracers of a collisional populations of small bodies and adopts a simple model for the evolution of the population. The properties of this collisional population are drawn from previous studies that fit models of the collisional evolution of planetesimal belts to the statistics for the evolution of debris disks (Dominik & Decin 2003; Krivov et al. 2005, 2006; Wyatt et al. 2007b; Löhne et al. 2008; Wyatt 2008; Krivov 2010; Kains et al. 2011). Although the parameters in this simple model have a precise physical meaning for the model as it is presented below, in practice they represent “effective” parameters that determine the mass evolution and are calibrated to match observations. For example, although we adopt the classic single power law size distribution for a collisional cascade (Dohnanyi 1969; Williams & Wetherill 1994), the true size distribution of a collisional population is certainly more complicated (O’Brien & Greenberg 2005; Kobayashi & Tanaka 2010). We also assume a constant value for , the impact energy required to catastrophically disrupt an object of size D that is not realistic because  is undoubtedly a function of D (e.g., Benz & Asphaug 1999). However, this type of simple model does a reasonably good job at matching more detailed calculations of debris disk evolution that include a size dependent  and (consequently) a multiphase size distribution (see Fig. 11 of Löhne et al. 2008; Kenyon & Bromley 2008, 2010; Krivov 2010). Furthermore the parameters of the model have been tuned to match the observed debris disk evolution (Wyatt et al. 2007b; Kains et al. 2011).

We divide our planetesimal populations into two components: asteroidal and cometary: asteroids are simply planetesimals interior to the giant planets’ initial orbits and comets are exterior. We assume that each planetesimal population is in a collision-dominated regime from the start of the simulation such that its differential size distribution is represented by n(D) ∝ D2 − 3qd, where D is the object diameter and qd = 11/6 for an infinite collisional cascade (Dohnanyi 1969; Williams & Wetherill 1994). This distribution spans from the smallest assumed particles Dbl = 2.2   μm up to the largest Dc = 2000 km. Smaller particles than Dbl are assumed to be blown out of the system by radiation pressure on short enough timescales to not contribute to the size distribution (see, e.g., Wyatt 2005). The largest bodies in the distribution were chosen to match the largest known objects in the Kuiper belt, assuming that the primitive Kuiper belt represents a proxy for the accretion that would have occurred in planetesimal disks of this mass. If Dc were significantly smaller then the collisional timescales would be shorter and the collisional evolution would proceed more quickly.

The surface area of a population of bodies with a collisional size distribution can be easily calculated. Using our chosen values for Dbl, Dc and qd, and assuming a physical density for all objects of 1 g   cm-3, the total cross-sectional area in a given radial bin is related to the total mass in that bin by: (3)where σ(R) is in AU2 and M(R) is in M.

For each population of bodies (asteroidal and cometary), we calculate the collisional timescale tc, which is a function of the particle size D as well as the orbital properties of the population (Wyatt et al. 1999, 2007a). This represents the mean timescale between collisions that are violent enough to destroy bodies of a given size at a mean orbital distance Rm with an annular width dr. (4)where (5)Here tc is in years, M   is the stellar mass in solar masses and σtot is the total surface area in AU2. The orbital characteristics of the planetesimal population are represented by their mean eccentricities e and inclinations I in radians. The factor fcc(D) represents the fraction of the total size distribution that could cause a catastrophic collision with a particle of size D, and for our assumptions can be expressed as (6)where Dcc(D) = XcD where XcD > Dbl and Dcc(D) = Dbl otherwise, and (7) is the impact energy required to catastrophically disrupt and disperse a body of size D. Here we assume that and that it does not vary with D. As discussed above, this assumption is not realistic in terms of the collisional dynamics (Benz & Asphaug 1999), but for our purposes represents an effective strength that determines the dust production and mass loss rate from the planetesimal belt, and using a constant value allows for a good fit to the statistics of debris disks around A stars (Wyatt et al. 2007b). In addition, it may not be reasonable to assume the same for asteroidal and cometary planetesimals given that their compositions are likely to be different. However, the typical collisional timescales for the largest (2000 km) bodies are very short, just a few  × 104 years, for the asteroidal planetesimal population (compared with 1 − 3  ×  108 years for the cometary planetesimal population). Thus, asteroidal dust is severely depleted within the first 0.1 − 1 Myr of each system’s evolution and our choice of has virtually no effect on the results.

For a given simulation timestep, we calculate the spectral energy distribution (SED) of the system dust as follows (again, as in Booth et al. 2009 with only small modifications):

  • Calculate tc for the largest bodies in both asteroidal and cometary planetesimal populations. We then artificially decrease the mass of the planetesimals in each population by a factor of [1 + t/tc(Dc)]-1 to account for collisional mass loss. This represents a very slow decrease for the comets but a significant mass loss among the asteroids2.

  • Divide the radial domain into Nbin radial bins spaced logarithmically between 0.2 and 100 AU, which are the inner and outer boundaries of our simulation. We tested a range of Nbin values and found good convergence with Nbin ≥ 30 − 50 and so we use Nbin = 100 to be conservative.

  • Calculate the total planetesimal mass in each radial bin for asteroids and comets by sampling the orbit of each planetesimal Nsamp times along its orbit at intervals that are equally spaced in time (i.e., in mean anomaly). We calculate Nsamp = 1 + e/elimit, where elimit represents the minimum statistical eccentricity needed to cross between radial bins: elimit ≈ 0.06 in our case but we use half of that value. Thus, we sample circular orbits sparsely but more eccentric orbits up to 30 times per orbit to allow the dust to be spread over a range of orbital distances.

  • For each radial bin and for asteroids and comets, calculate the blackbody temperature3(8)where L   is the stellar luminosity, R is the radial distance, and Tbb is in Kelvin. Then, assuming that all objects radiate as blackbodies, we calculate the flux density (in Janskys) Fν seen by an observer at distance d: (9)where Bν is the Planck function in units of Jy   sr-1 and Xλ is a factor derived by Wyatt et al. (2007b) to account for a decrease in emission beyond  ~200   μm needed to match observed sub-mm measurements of debris disks: Xλ = 1 for λ < 210   μm and Xλ = λ/210 for longer wavelengths.

  • Sum the contributions from each radial bin over the whole SED, and then sum the asteroidal and cometary contributions. We consider the wavelength range from 1 to 1000 μm.

Using this method, we calculate the SED of each simulation from the orbits of all planetesimal particles at each simulation timestep. To make the SED useful for comparison with observations, the most convenient quantity is the dust flux relative to the stellar flux. The stellar flux Fν   in Jy can be calculated as (10)where Bν is again the Planck function in units of Jy   sr-1 and T   is the stellar effective temperature, which we take to be 5777 K for all cases. Armed with the stellar flux as well as the dust flux, we can isolate the dust-to-stellar flux ratio at any wavelength of interest such as those observed with Spitzer or Herschel.

Our simulations ran for 100 − 200 Myr, but we want to compare fluxes with stars at a range of ages. Thus, we need to extrapolate the expected dust fluxes into the future. We do this by making the assumption that there is no more significant dynamical evolution of the system, i.e. that the orbital characteristics of the planetesimals are not going to change in the future. Clearly, this assumption does not account for punctual events like the late heavy bombardment (Gomes et al. 2005) or very late dynamical instabilities. For a given time t after the end of the simulation, we decrease the planetesimal mass by a simple factor related to the collisional timescale at the end of the simulation, i.e., as . In terms of the global analysis of our sample of simulations, we retain snapshots of each simulation at wavelengths of 1, 3, 5, 15, 20, 25, 50, 70, 100, 160, 250, 350, 500 and 850 μm at time intervals of 1, 3, 10, 30, 100, 300, 1000, and 3000 Myr.

Our dust flux calculations compare reasonably well with observations and other models. We perform this comparison using systems that are dynamically calm (referred to as “stable systems” later in the paper), in which the giant planets remain on stable orbits throughout the simulation such that the outer planetesimal disk remains largely intact. In these stable systems, our dust flux calculations yield typical values for the dust-to-stellar-flux ratio F/Fstar of 0.1 − 0.5 at 25   μm and 10 − 35 at 70   μm after 1 Gyr of simulated dynamical and calculated collisional evolution. These values are broadly consistent with the fluxes detected around solar-type stars with observed excesses at 24 and 70   μm (Habing et al. 2001; Beichman et al. 2006; Moór et al. 2006; Trilling et al. 2008; Hillenbrand et al. 2008; Carpenter et al. 2009; Gáspár et al. 2009), although our stable simulations yield very few systems with F/Fstar(70   μm) ≈ 1 − 10, probably because of the relatively large masses in our outer planetesimal disks (note that our unstable systems can yield those flux levels). Compared with the more sophisticated models of dust production during planetary accretion of Kenyon & Bromley (2008, 2010), our calculated fluxes are larger by a factor of a few, notably at 70   μm. Our fluxes are only a factor of 2 − 3 larger than those calculated by Kennedy & Wyatt (2010).

thumbnail Fig. 2

Evolution of a system with a relatively late (21 Myr) instability among the giant planets. Each panel shows a snapshot in time of orbital eccentricity vs. semi-major axis for all particles; vertical bars denote sin(i) for terrestrial bodies with Mp > 0.2   M and i > 10°. The particle size is proportional to the mass1/3, but giant planets (large black circles) are not on this scale. Colors denote water content, assuming a Solar System-like initial distribution (Raymond et al. 2004). The surviving terrestrial planet has a mass of 0.72 M, a stable orbit within the habitable zone (semimajor axis of 0.96 AU), and a high eccentricity and inclination (and large oscillations in these quantities). A movie of this simulation is available in the electronic edition of the Journal.

Open with DEXTER

3. An example simulation

In this section we explore the dynamics and dust properties of a single simulation. In following sections we will consider the outcome of an ensemble of many simulations, so this is an opportunity to inspect the detailed evolution of a particularly interesting system. The initial conditions for the chosen simulation are shown in Fig. 1: the giant planet masses were 1.46 MJ (at 5.2 AU), 0.75 MJ (7.9 AU) and 0.55 MJ (11.3 AU).

Figure 2 shows the dynamical evolution of the system, which became unstable after 21 million years. Before the instability, the system evolved in the expected quiescent fashion. In the inner disk, embryos accreted planetesimals and other embryos from the inside-out. Embryos’ eccentricities and inclinations were kept small by the planetesimals via dynamical friction and viscous stirring. After 21 million years there were several almost fully-grown terrestrial planets, including three planets more massive than 0.6 M at 0.61, 0.93 and 1.34 AU and a dozen other embryos of 0.06−0.3 M extending out to 2−3 AU. During this phase of quiescent accretion, the giant planets’ orbits remained almost perfectly circular, with eccentricities of less than 0.05 for each planet (less than 0.03 for the innermost gas giant). Planetesimal scattering caused a slow drift in the giant planets’ semimajor axes, of 0.06 AU for the innermost (and most massive) giant planet, and less than 0.01 AU for the other two giants. During this time, the outer planetesimal disk was slowly sculpted by the giant planets. The inner edge of the planetesimal disk was eroded, and strong resonances with the outermost giant planet (low-order mean motion resonances) were slowly cleared out. At the time of the instability, roughly three quarters of the initial outer planetesimal disk was intact.

The instability began with a rapid eccentricity increase in the eccentricity of the two outer giant planets and was followed by a close encounter between those two. That set off a series of close encounters between all three giant planets that lasted 84 000 years, involved 35 planet-planet scattering events, and culminated with the ejection of the outermost giant planet. The behavior of the simulation was characteristic in terms of dynamical instabilities in that the least massive giant planet was ejected and the surviving planets were segregated by mass, with the most massive one closest to the star (at 3.65 AU) and the less massive one farther out (at 36.6 AU, e.g. Chatterjee et al. 2008; Raymond et al. 2010).

During and immediately after the giant planet instability, the inner disk was almost entirely thrown into the star. This happened by a rapid increase in embryos’ and planetesimals’ eccentricities mainly by secular forcing from the scattered giant planets4. Among the 3.9 M of terrestrial material that was destroyed were 13 embryos including a 0.95 M at 0.6 AU, a 0.62 M planet in the habitable zone at 1.34 AU, and five other embryos larger than 0.15 M. All but one of the embryos – the sole survivor – collided with the star within a few hundred thousand years of the start of the instability. The rocky planetesimals were all destroyed within 1 Myr. The outer planetesimal disk was completely destabilized by the instability, as the two lower-mass giant planets scattered each other to tens of AU. All planetesimals were either ejected from the system (about 85%) or collided with the star (15%). The vast majority were destroyed within 1 Myr of the instability – all but three within 5 Myr – but the last planetesimal took almost 25 million years to be ejected.

Three planets – two gas giants and one terrestrial planet – survived to the end of the simulation, all on excited, eccentric and inclined orbits. The long-term dynamics of the two surviving giant planets is shown in Fig. 3. The eccentricities of the surviving giant planets are considerable as are their orbital inclinations with respect to the starting plane. The two giant planets have a mutual orbital inclination of  ~36°. Although this is less than the formal limit of 39.2° required for the Kozai effect to take place (Kozai 1962; Takeda & Rasio 2005), out of phase, Kozai-like oscillations are evident in the planets’ eccentricities and inclinations.

thumbnail Fig. 3

The eccentricity (top panel) and inclination (bottom panel) evolution of the two surviving giant planets in the example simulation. The three curves in the right panel show the inclination of each planet relative to the starting plane of the system – presumably the stellar equatorial plane – in black and light grey and their mutual inclination, which is larger than each of their individual inclinations, in dark grey.

Open with DEXTER

The terrestrial planet’s final mass is 0.72   M and its semimajor axis of 0.935 AU places it in the circumstellar habitable zone (Kasting et al. 1993; Selsis et al. 2007). However, Fig. 4 shows that the planet’s orbit is strongly perturbed. Its eccentricity oscillates between 0.4 and 0.7 and its inclination between almost zero and more than 60° on a  ~10 000 year timescale (although there other longer secular frequencies in the oscillations). On Myr timescale the orbit has a minimum and maximum eccentricity of 0.39 and 0.73 and a minimum and maximum inclination of 0.03° and 63.3°. Given its large eccentricity, one would expect this planet’s climate to be highly variable during the year (Williams & Pollard 2002; Dressing et al. 2010). In addition, the changes in both its eccentricity and inclination – equivalent to changes in obliquity for a fixed spin axis – would cause variations on the secular timescales (Spiegel et al. 2010). Given that the orbit-averaged stellar flux increases with the orbital eccentricity (as ) and that the planet’s closest approach to the star is only 0.25 AU, this planet may not be habitable during its high-eccentricity phases. However, more detailed modeling of such planets is beyond the scope of the paper and the reader is referred to recent climate modeling papers that include the effects of varying the eccentricity and obliquity (Williams & Pollard 2002, 2003; Spiegel et al. 2009, 2010; Dressing et al. 2010).

thumbnail Fig. 4

The eccentricity (top panel) and inclination (bottom panel) of the surviving terrestrial planet in the example simulation (mass of 0.72 M, semimajor axis of 0.93 AU). The two curves in the right panel show the inclination relative to the starting plane of the system in grey and the mutual inclination with respect to the innermost giant planet in black.

Open with DEXTER

The evolution of the dust brightness follows from the dynamical and collisional evolution of the planetesimals in the system. Figure 5 shows the resulting spectral energy distribution at different snapshots in time. The collisional timescale for the largest planetesimals (2000 km) in the terrestrial planet forming region is only tcoll ~ 104 years, so the planetesimals in the inner disk are ground into hot dust within just a few million years. Thus, the total dust brightness has dropped sharply at all wavelengths by the 20 million year snapshot, most dramatically at wavelengths shorter than roughly 50 μm. For the rest of the simulation, the primary source of dust is the outer planetesimal belt for which tcoll ≳ 108 years because the inner disk planetesimals have been ground away5.

When the instability occurs, the spectral energy distribution changes dramatically and quickly (Fig. 5). Comets from the outer planetesimal disk are scattered onto highly eccentric orbits. Those that enter the inner Solar System can cause bombardments on the terrestrial planets akin to or often far more intense than the Solar System’s late heavy bombardment (Gomes et al. 2005; Strom et al. 2005), although given the much earlier timing of most instabilities, these bombardments could act to seed the terrestrial planet-forming region rather than impact fully-formed planets. The decrease in the comets’ perihelia introduces a large amount of warm dust into the system which lasts for the duration of the bombardment and leads to a spike in brightness at near- to mid-infrared wavelengths, shown in Fig. 6 for λ = 5,10,25,70,160 and 500   μm (see also Booth et al. 2009; Nesvorný et al. 2010). At shorter wavelengths the flux is more jagged than at longer wavelengths because the short wavelength flux is entirely due to hot dust that is produced sporadically by individual planetesimals entering the inner planetary system. In contrast, the long wavelength flux combines the flux from a much larger range of dust temperatures and therefore includes a much larger number of particles. The spike in flux associated with the instability is strong for wavelengths shorter than  ~ 50   μm, but at longer wavelengths the spike is weak or absent. As objects are dynamically removed from the system the system’s brightness drops precipitously and out of the detectable range in the few million years following the instability. We note that a more realistic model of dust production by new comets suggests that the mid-infrared peak during the bombardment in our model may be underestimated by a factor of a few (Nesvorný et al. 2010).

thumbnail Fig. 5

Evolution of the spectral energy distribution (SED) of the simulation from Fig. 2. Each curve shows the SED during a given simulation snapshot. The instability occurred at 21 Myr, and the SED evolved dramatically in the immediate aftermath, as icy planetesimals were scattered onto high-eccentricity orbits – thereby producing transient hot dust – and then ejected from the system. The dashed line represents the stellar photosphere. The system is viewed at a distance of 10 pc.

Open with DEXTER

4. Results for the ensemble of simulations

In this section we study the outcome of our ensemble of simulations (called mixed in Paper II). This set is of particular interest because the surviving giant planets match the observed eccentricity distribution without any fine-tuning (Raymond et al. 2009a, 2010). We explore the formation process of terrestrial planets, the orbital characteristics of terrestrial planets that formed, giant planet dynamics, dust production from planetesimals, and correlations between these. We also compare the simulations with the observed extra-solar giant planets and debris disks. In Paper II we explore the effect of a variety of parameters on the process.

4.1. Giant planet instabilities

Of the 152 simulations that ran for at least 100 Myr and met our energy conservation cutoff (see Appendix A), the giant planets were unstable in 96 systems (63%). The instability times ranged from 247 years to 180 Myr after the start of the simulation with a median of 91 600 years. This is encouraging because it is close to the value expected based on our initial giant planet separations of 4 − 5   RH,m (Marzari & Weidenschilling 2002; Chatterjee et al. 2008). This means that our assumption of no disk gas at the start of the instability is marginally acceptable, because the typical timescale for the gas dissipation is thought to be  ~ 105 years (Wolk & Walter 1996; Ercolano et al. 2011). However, in many cases instabilities are likely to have occurred in the presence of a residual gas disk. To test this assumption, we ran an additional set of simulations that include a simple prescription for the effects of gas damping that are presented in Paper II (and show no significant changes from the gas-free simulations). We note that later instabilities can certainly occur (e.g., the late heavy bombardment) but we expect the number of instabilities to decrease in time, and the required computing time made it impractical to search for instabilities beyond 100 − 200 Myr.

Although 75% of instabilities occur within 1 Myr, there is a tail that extends to longer timescales. One in six instabilities (16 out of 96) occurred after 10 Myr, one in ten (10/96) after 30 Myr, and one in fifty (2/96) after 100 Myr. Given that our simulations only lasted 100 − 200 Myr, we undersampled the fraction of unstable systems in the 100 − 200 Myr timespan and we expect that even later instabilities should certainly occur in a fraction of systems due to long-term chaotic diffusion. The timing of the instability is important in terms of the sizes of objects in the inner disk; instabilities later than the typical terrestrial planet formation time of 10 − 100 Myr may destroy fully-formed Earth-sized planets – sometimes with appreciable water contents – rather than embryos and planetesimals.

As shown in Fig. 7, the surviving giants in the unstable simulations provide a quantitative match to the observed extra-solar giant planets (p = 0.49 from a Kolmogorov-Smirnov test, consistent with previous work with very similar initial conditions but much larger samples; Raymond et al. 2008a, 2009b,a, 2010). Thus, one might imagine that the unstable systems represent the appropriate subsample of simulations that we should use to represent the known exoplanet systems, and that our stable simulations are unrealistic in that they somehow lack a trigger to make them unstable (e.g., Malmberg et al. 2011). However, as we discuss in Sect. 5.1, the observed exoplanets do appear to require a contribution from dynamically stable systems. We note that additional constraints exist in the exoplanet sample (e.g., the mass-eccentricity correlation). Combinations of different sets of simulations can match all of the observed characteristics of the exoplanet distribution, and we perform this exercise in Paper II (see also Sect. 5 in Raymond et al. 2010).

thumbnail Fig. 6

The ratio of the dust-to-stellar flux as a function of time for six different wavelengths from 5 to 500 microns, including a zoom during the instability – note that each main plot is on a log scale and each zoom-in is on a linear scale. The rough observational limits of the MIPS instrument on NASA’s Spitzer Space Telescope are shown for 25 and 70 μm with the dashed line Trilling et al. (2008). All planetesimal particles were destroyed as of 45 Myr via either collision or ejection so there is no more dust in the system after that point.

Open with DEXTER

A prediction of the planet-planet scattering model is that most planetary systems should contain multiple giant planets, i.e., additional giant planets should exist exterior to most of the known ones (e.g., Rasio & Ford 1996; Marzari & Weidenschilling 2002). Likewise, the vast majority of unstable simulations in our sample (85/96 = 89%) contain multiple giant planets; the remaining 11% contain just a single surviving giant planet. The outer giant planet is typically 5 − 10 AU more distant than the inner giant planet (with a tail to  > 30 AU), so long-duration observations are needed to follow up the known giant exoplanets. The distribution of separations between planets in scattered systems may actually provide constraints on their initial mass distribution because the surviving planets tend to be more widely-separated if they are equal-mass than if their masses differ markedly (Raymond et al. 2009b).

4.2. Terrestrial planet formation

The number and spacing of terrestrial planets that form is governed by the eccentricities of the planetary embryos from which they form (Levison & Agnor 2003). These eccentricities are a result of gravitational forcing both from nearby embryos and the giant planets. There are two differences in terrestrial planet formation between stable and unstable systems: 1) perturbations during the instability can play a major role in shaping the embryo distribution; and 2) the dynamical state of the giant planets after an instability is generally more excited than before the instability. In some unlucky cases the surviving giant planets provide an accretion-friendly environment that is empty because the instability has already removed all rocky bodies. Past simulations of terrestrial planet formation with giant planets on excited orbits have all neglected the planet-planet scattering phase during which the giant planets actually acquired their eccentricities (Chambers & Cassen 2002; Levison & Agnor 2003; Raymond et al. 2004; Raymond 2006; Ji et al. 2011), which clearly plays a very important role in the dynamics.

thumbnail Fig. 7

Eccentricity distributions of the surviving giant (solid black line) and terrestrial (dashed black line) planets in the unstable systems. The grey curve represents the known extra-solar planets beyond 0.2 AU – this helps to exclude planets whose orbits have been tidally circularized as well as low-mass planets.

Open with DEXTER

Giant planet perturbations span a continuous range but the outcome is quantized into a discrete number of terrestrial planets during the accretion process. If perturbations are weak – if the giant planets collide rather than scattering (or are dynamically stable or low-mass) – then embryos’ eccentricities remain small and feeding zones narrow and several terrestrial planets form. For stronger giant planet perturbations, feeding zones widen and fewer terrestrial planets form, although the total mass in planets tends to decrease because stronger perturbations imply that the giant planets were scattered closer to the terrestrial planet region so more embryos end up on unstable orbits. In systems where embryos’ radial excursions are comparable to the radial extent of the surviving disk only one planet forms, usually on an excited orbit. In the simulation from Fig. 2, the lone terrestrial planet did not accrete from a disk of excited embryos but rather was the only planet to survive the instability. Perturbations during, not after, the instability determined the outcome.

The strength of the giant planet perturbations is directly related to the smallest giant planet perihelion distance qGP,min. Figure 8 shows that the efficiency of terrestrial planet formation is directly related to qGP,min. This is not surprising: planets with smaller qGP,min have higher eccentricities (since all giant planets start with the same range of semimajor axes) and have therefore undergone more violent scattering events. Every system with a giant planet scattered to qGP,min < 1.3 AU destroyed all terrestrial material in the system, as did some simulations out to 3 AU. This is consistent with the results of Veras & Armitage (2006), who found a similar scaling between qGP,min and the survival of test particles in the inner disk. As expected, systems that form a single terrestrial planet are those with giant planet perturbations that are almost strong enough to completely empty the terrestrial material from the system. In these cases the embryos’ eccentricities are large and, as we will see below, the single planet that survives maintains an excited orbit.

thumbnail Fig. 8

The total mass in surviving terrestrial planets as a function of the minimum perihelion distance of a giant planet during the simulation. Black dots represent systems in which the giant planets were dynamically unstable and grey dots are systems that were stable. Filled black dots are systems in which a single terrestrial planet formed.

Open with DEXTER

4.3. Characteristics of surviving terrestrial planets

All terrestrial material was destroyed in 41 out of 96 unstable simulations (43%; Fig. 9). Likewise, 22 unstable systems (23%) formed a single terrestrial planet, 16 formed two terrestrial planets (17%), and and the remaining 17 systems (18%) formed three or more terrestrial planets. Among the 56 stable simulations, only 7 (12.5%) formed two planets, and the remaining 87.5% of simulations formed three or more planets.

thumbnail Fig. 9

Distribution of the number of surviving terrestrial planets (defined as Mp > 0.05   M) in our fiducial mixed set of simulations. The grey shaded histogram shows the unstable simulations and the dashed line shows the stable simulations.

Open with DEXTER

As a population, the surviving terrestrial planets have smaller eccentricities than the giant planets (Fig. 7); the median eccentricity is 0.10 for terrestrials and 0.21 for giants. The terrestrial planets that formed in systems with unstable giant planets have only slightly more eccentric and inclined orbits than the terrestrial planets that formed in stable systems. The median e and i were 0.10 and 5.1° for the unstable systems and 0.08 and 2.5° for the stable systems, respectively. The single-terrestrial planet systems were the most excited, with median e of 0.14 and i of 8.7°. If we neglect the single-terrestrial planet systems, the eccentricities and inclinations of the terrestrial planets that formed in stable and unstable systems is very close6.

Although their time-averaged orbits are similar, terrestrial planets in unstable systems undergo much stronger orbital oscillations (Fig. 10). The median peak-to-peak eccentricity and inclination oscillation amplitudes are 0.11 and 5° for the unstable systems and 0.043 and 2.8° for the stable systems, respectively. Again, the single-terrestrial planets are the most excited and undergo the largest oscillations, with median e and i amplitudes of 0.21 and 11.9°. The climates of the single planets are likely to vary dramatically on the secular timescale of 103−106 years (Spiegel et al. 2010).

thumbnail Fig. 10

The peak-to-peak oscillation amplitudes in eccentricity e and inclination i for the surviving terrestrial planets in simulations with stable (grey) and unstable (black) giant planets. The single-terrestrial planet systems are shown with filled black circles, and Earth is the grey star.

Open with DEXTER

thumbnail Fig. 11

The total mass in surviving terrestrial planets as a function of the eccentricity of the innermost surviving giant planet. The Solar System is shown as the grey star.

Open with DEXTER

4.4. Correlations with observable quantities: giant planets and debris disks

Observations of massive exoplanets can only very roughly diagnose the outcome of terrestrial accretion. Figure 11 shows a strong negative correlation between the efficiency of terrestrial planet formation and the eccentricity of the innermost giant planet eg. For stable systems eg tends to be very small (median eg = 0.008) while for the unstable systems eg spans a very wide range, from  < 0.01 to 0.8 (median eg = 0.21). The correlation between the total terrestrial planet mass and eg is expected but the range in terrestrial planet mass for different systems with a similar eg shows the importance of the orbital history.

thumbnail Fig. 12

Distribution of eccentricities of the innermost giant planet in simulations which formed zero (blue line), one (green), and two or more (red) terrestrial planets. The sum of the three distributions (black dashed line) provides a quantitative match (p = 0.49 from a Kolmogorov-Smirnov test) to the observed exoplanets (thick grey line).

Open with DEXTER

The orbits of surviving giant planets retain a memory of the strength of the instability, or lack thereof. Figure 12 breaks down the giant planet eccentricity distribution by outcome in terms of the number of terrestrial planets that formed. Multiple terrestrial planets form preferentially in systems with low giant planet eccentricities, because these represent very weak instabilities. By the same argument, highly eccentric giant planets tend to destroy all terrestrial material in their systems. The intermediate regime is represented by the single-terrestrial planet systems; for these systems eg is typically in the range 0.1 − 0.3 (median of 0.17). These eccentricities are close to the median of the exoplanet distribution (Butler et al. 2006; Udry & Santos 2007) and there is considerable overlap from systems with zero or many terrestrial planets. It is therefore very difficult to diagnose a single-planet system from a measurement of just the eccentricity of a giant exoplanet.

thumbnail Fig. 13

The dust-to-stellar flux ratio at 70 μm after 1 Gyr of dynamical and collisional evolution, plotted as a function of the eccentricity of the innermost surviving giant planet. Unstable systems are plotted in black and stable systems in grey (systems with a single surviving planet are shown with solid circles). The pileups very close to the x axis represent systems with virtually zero 70 μm flux. The dashed line represents an approximate threshold above which excess emission was detectable using Spitzer data (Trilling et al. 2008). The star shows the estimated flux from the Kuiper Belt at 1 Gyr, as calculated previously (Booth et al. 2009) based on dynamical simulations of the outer Solar System (Gomes et al. 2005).

Open with DEXTER

thumbnail Fig. 14

The dust-to-stellar flux ratio at 70   μm after 1 Gyr of dynamical and collisional evolution, plotted as a function of the total mass in terrestrial planets. Unstable systems are plotted in black and stable systems in grey. The pileups close to the x and y axes represent systems with either no terrestrial planets or virtually zero 70   μm flux. The dashed line represents an approximate threshold above which excess emission was detectable using Spitzer data (Trilling et al. 2008). The star shows the estimated flux from the Kuiper Belt at 1 Gyr, as calculated previously (Booth et al. 2009) using models based on dynamical simulations of the outer Solar System (Gomes et al. 2005). The Solar System falls into an intermediate region of parameter space: the giant planets may have been unstable, but only weakly so (there is no clear evidence for ejections or planetary collisions), while the Kuiper Belt would have been bright for several hundred Myr prior to the Late Heavy Bombardment.

Open with DEXTER

The outer disk’s evolution is also governed by giant planet perturbations. Icy planetesimals – whose collisional erosion creates debris disks – survive in dynamically calm environments where the giant planets were either stable, low-mass, or underwent a relatively weak instability. Indeed, Fig. 13 shows a strong anti-correlation between the 70   μm dust flux and the giant planet eccentricity. This arises simply because eccentric giant planets have increased apocenter distances and impinge on the planetesimal disk, thereby dynamically removing planetesimals via ejection. In addition, if planetesimals survive on highly-eccentric orbits their collisional lifetimes may be shortened (Wyatt et al. 2010).

thumbnail Fig. 15

The dust-to-stellar flux ratio F/Fstar after 1 Gyr of dynamical and collisional evolution as a function of the total mass in terrestrial planets, for six wavelengths between 5 and 500 microns. The Spitzer observational limits are shown for 25 and 70 μm with the dashed line (Trilling et al. 2008). In each panel, grey circles represent stable simulations and black circles represent unstable simulations.

Open with DEXTER

Given that the efficiency of terrestrial planet formation anti-correlates with the giant planet eccentricity (Fig. 11) and the dust flux also anti-correlates with the giant planet eccentricity, Fig. 14 shows that the efficiency of terrestrial planet formation correlates with the dust flux. Stars with bright dust almost all contain terrestrial planets: the median system for which F/Fstar > 10 contains 3.6   M in terrestrial planets, and every single system contains at least 1.85   M in terrestrial planets. Systems that are extremely bright at long wavelengths should therefore be considered prime targets in the search for terrestrial planets.

thumbnail Fig. 16

Left: the fraction of systems that would be detectable after 1 Gyr as a function of the total mass in surviving terrestrial planets for six different wavelengths between 5 and 500   μm. Each curve represents a horizontal slice through Fig. 15 at a different vertical height. The chosen “detection limits” in F/Fstar were: 10-8 at 5   μm, 10-3 at 15   μm, 0.054 at 25   μm, 0.55 at 70   μm, 10 at 160   μm, and 100 at 500   μm. The error bars are 1-σ and were calculated using binomial statistics (Burgasser et al. 2003). Note that the different curves are offset by up to  ± 0.1   M for clarity. The bins themselves are at: zero, 0−1   M, 1−2   M, 2−3   M, and  > 3   M. Right: the fraction of systems with 0.5   M or more in terrestrial planets as a function of F/Fstar(70   μm) (1 Gyr). 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. 14.

Open with DEXTER

thumbnail Fig. 17

The dust-to-stellar flux ratio F/Fstar at 70   μm as a function of the total mass in terrestrial planets for eight different times from 1 Myr to 3 Gyr. The terrestrial planet mass refers to the final value such that simulations move vertically in time on the plot. The Spitzer observational limit is shown with the dashed line (Trilling et al. 2008). Grey circles represent stable simulations and black circles represent unstable simulations.

Open with DEXTER

thumbnail Fig. 18

Left: the fraction of systems that would be detectable after 1 Gyr as a function of the eccentricity of the innermost surviving giant planet eg for six different wavelengths between 5 and 500   μm. Each curve represents a horizontal slice through plots such as Fig. 13 at a different vertical height. The chosen “detection limits” in F/Fstar were the same as in Fig. 16: 10-8 at 5   μm, 10-3 at 15   μm, 0.054 at 25   μm, 0.55 at 70   μm, 10 at 160   μm, and 100 at 500   μm. The error bars are 1-σ and were calculated using binomial statistics. Note that the different curves are offset for clarity. The bins themselves are logarithmically spaced between 10-3 and 1. Right: the fraction of systems with eg > 0.1 as a function of F/Fstar(70   μm) (1 Gyr). 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. 13.

Open with DEXTER

The debris disk-terrestrial planet correlation seen in Fig. 14 exists because the inner and outer planetary system are subject to the same dynamical environment: the violent instabilities that abort terrestrial planet formation also tend to remove or erode their outer planetesimal disks. As was the case for the giant planet eccentricities, single-terrestrial planet systems populate the intermediate area of the correlation and overlap with systems with no planets as well as those with several. The terrestrial planet – debris disk correlation is not perfect as there exist “false positives” with bright dust emission and no terrestrial planets, corresponding to systems with asymmetric, inward-directed giant planet instabilities. Conversely, “false negatives” with terrestrial planets but little to no dust are systems that underwent asymmetric but outward-directed instabilities. We discuss these caveats in detail in Paper II when we compare our simulations with the observed debris disk statistics and use them to derive a crude estimate of ηEarth.

The debris disk-terrestrial planet correlation is a function of wavelength. Figure 15 shows the correlation at 1 Gyr for six wavelengths between 5 and 500 μm. At 5   μm there is no correlation but mainly a scatter plot. The only hint of a correlation is a lower envelope of dust fluxes larger than F/Fstar ≳ 10-20 for systems with more than  ~2   M in terrestrial planets. This non-correlation comes from the fact that the short-wavelength flux can be dominated by either a small amount of hot dust or a large amount of colder dust. At 15   μm the debris disk-terrestrial planet correlation clearly exists, although there are still some outliers with large fluxes and small terrestrial planet masses. At longer wavelengths these outliers slowly disappear because the signal from small amounts of transient hot dust is overwhelmed by the large amount of cold dust produced in quiescent outer planetesimal disks. The debris disk-terrestrial planet correlation is evident for all wavelengths longer than  ~25   μm.

Figure 16 (left panel) shows histograms that represent horizontal slices through Fig. 15. For each wavelength we chose a “detection limit” and tabulated the fraction of systems that were deemed detectable as a function of the total terrestrial planet mass. These detection limits were chosen for illustrative purposes and are in many cases unrealistic. For example, the current detection threshold at 5   μm is probably closer to  ~10-3, but none of our simulations would be detectable at that limit. The reader is referred to other work discussing current dust detection limits at short (Absil et al. 2006; Akeson et al. 2009) and long (Smith et al. 2009) wavelengths (see also Wyatt 2008). The fraction of systems with no terrestrial planets that is detectable varies significantly between the different wavelengths and simply reflects the location of the detection limit with respect to the general trends in Fig. 15. The detection limits at 25   μm and 70   μm, which correspond to the actual Spitzer limits (Trilling et al. 2008), are close to the extremes: at 25   μm not a single terrestrial planet-free system is detectable while at 70   μm more than 40% of all systems are detectable.

At each wavelength, the fraction of systems that is detectable increases with the terrestrial planet mass (Fig. 16). At all wavelengths the trend is relatively flat to a certain point where it increases significantly, by several standard deviations between adjacent bins. At all wavelengths this transition occurs between either the 1−2   M and 2−3   M bins or the 2−3   M and the  > 3   M bin. Even at 5   μm, which showed no obvious debris disk-terrestrial planet correlation in Fig. 15, there is marked jump in the fraction of detectable objects in the last bin. This jump is also seen at 15   μm, although these two shortest wavelengths are the only ones for which some systems with more than 3   M in planets are not detectable. Although our chosen detection limits at certain wavelengths may be overly optimistic or pessimistic (e.g., F/Fstar(5   μm) ≥ 10-8 is likely to be difficult to achieve in the near term), this shows that the debris disk-terrestrial planet correlation is present for all of these wavelengths.

The right panel of Fig. 16 shows the fraction of systems with significant mass in terrestrial planets as a function of F/Fstar at 70   μm. This plot represents a vertical slice through Fig. 14. At small dust fluxes, a minority of systems contain terrestrial planets – these are referred to in Sect. 5 as “false negatives”. The fraction of systems with terrestrial planets increases dramatically beyond the detection limit, and every single one of the 58 systems with F/Fstar > 10 (including 7 unstable systems) contains at least 1.9   M in terrestrial planets, with a median value of 3.5   M. Of the 53 unstable systems above the detection threshold at 70   μm, 35 () contain terrestrial planets. Thus, in our simulations debris disks appear to represent signposts for terrestrial planets with a confidence of at least 66% at 70   μm.

The anti-correlation between debris disks and eccentric giant planets in our simulations also holds across many wavelengths. Figure 18 (left panel) shows the fraction of stars that is detectable as a function of the giant planet eccentricity for a range of different wavelengths, similar to Fig. 16. Again, the trends are stronger for wavelengths of 25   μm and longer, but are still clear at shorter wavelengths (although the detection thresholds at these wavelengths are certainly much smaller than in reality). Figure 18 (right panel) shows that the fraction of systems containing an eccentric giant planet is anti-correlated with the dust emission at 70   μm, simply because in our simulations the calmest giant planets destroy the fewest number of outer planetesimals and therefore produce the brightest debris disks.

The debris disk-terrestrial planet correlation is also a function of time. Figure 17 shows F/Fstar at 70   μm vs. the final terrestrial planet mass for all simulations at eight snapshots between 1 Myr and 3 Gyr. At 1 Myr all systems are above the detection threshold, even those that have already undergone violent instabilities, as the timescale for clearing out planetesimals is generally closer to a few to 10 Myr from the time of the instability. In a given snapshot, systems that have not yet become unstable but that eventually will are those for which the flux remains as high as the cluster of stable systems but for which the total terrestrial planet mass is low, indicating a strong future depletion of rocky material. After an instability occurs the flux drops (though differently at different wavelengths; Fig. 6) but the total terrestrial planet mass, measured at the end of the simulation, does not change, such that a given system moves vertically between snapshots. In time, instabilities remove systems at low terrestrial planet mass and high F/Fstar; the last two instabilities occurred at 168 Myr (in which one 2.4   M terrestrial planet survived on a highly eccentric orbit) and 180 Myr (in which four fully-grown terrestrial planets were destroyed including a 1.3   M planet in the habitable zone). Beyond the end of our simulations one can imagine that a fraction of stable systems could actually become unstable and quickly lose a large fraction of their flux (and perhaps their terrestrial planets as well).

Stable systems remain clustered at high fluxes at all times, but decreasing in time due to collisional grinding7. The same collisional grinding affects the planetesimals that survive in unstable systems. In some cases the outer planetesimal disk is only moderately perturbed by the instability, although in all cases it is somewhat disturbed as shown by the lack of unstable systems with fluxes as high as the stable systems. After an instability, the mass in planetesimals decreases, although the planetesimals’ eccentricities and inclinations both tend to increase (not always in a simple correlated fashion). Mass loss causes the planetesimal population’s collisional evolution to slow down and eventually stop, as is thought to have occurred in the Solar System’s asteroid belt (Petit et al. 2001; Bottke et al. 2005). Thus, the dust flux decreases to a roughly asymptotic value. Once the collisional timescale for the small particles becomes longer than the timescale for Poynting-Robertson our calculation breaks down although this occurs at a low enough flux that it should not affect our results (Wyatt et al. 2007a).

Figure 17 shows that the debris disk-terrestrial planet correlation holds in time, especially after 10 − 100 Myr. However, if the figure were plotted including the total terrestrial mass at a given time the correlation would hold at even earlier times because individual systems would not be restricted to move between panels on vertical lines.

thumbnail Fig. 19

The p value from a Kolmogorov-Smirnov test comparing the eccentricities of our simulations vs. the observed exoplanet distribution as a function of the fraction of stable systems included in the sample. The p value drops below an acceptable level of 0.01 for more than a 17% contribution from stable systems.

Open with DEXTER

5. Discussion

In this section we first scale our simulations (Sect. 5.1) to match the exoplanet semimajor axis distribution and infer the orbital distribution of terrestrial planets in the current sample (mainly drawn from the radial velocity sample). In Sect. 5.2 we compare our results with the statistics of known debris disks (including cases with known giant planets). In Sect. 5.3 we discuss the Solar System in the context of our results. In Sect. 5.4 we discuss the limitations of our simulations.

5.1. Scaling to the observed a e exoplanet distribution

The surviving giant planets in our unstable simulations provide a match to the observed eccentricity distribution (Fig. 7). However, in constructing a sample of simulations that represents the observed exoplanet systems we do not think that is realistic not to include any stable systems for several reasons. First, current attempts to de-bias the observed eccentricity distribution infer a substantial fraction of systems – up to  ~30% – with circular or near-circular orbits (Shen & Turner 2008; Zakamska et al. 2010). Second, there exist individual systems that show no obvious signs of instability, for example with planets in resonance (e.g., GJ 876; Rivera et al. 2010) or with many relatively closely-packed giant planets (e.g., 55 Cancri; Fischer et al. 2008).

Figure 19 (left panel) shows the effect of the contribution from stable systems on the goodness of fit to the exoplanet eccentricity distribution. The distribution is best matched when we do not include any stable systems in the sample, and drops below a nominal statistically acceptable limit of p = 0.01 for a contribution larger than 17%. We therefore construct our sample with a 10% contribution of stable systems, which we think is a reasonable compromise between matching the exoplanet eccentricity distribution and allowing for stable systems.

The giant planets in our simulations are generally at larger orbital distances than the observed giant planets. The observed sample of giant planets displays a rapid rise between 0.5 − 1 AU, a plateau to 2 AU, and then a decrease at larger orbital distances (Butler et al. 2006; Udry & Santos 2007). The rise and plateau are real but the decrease beyond 2 AU is an observational bias due to the long orbital periods of these planets. The actual distribution of giant planets at Jupiter-Saturn distances is unknown, although observational constraints predict that at least 10%, and perhaps more than 50%, of stars have giant planets at these separations (Cumming et al. 2008; Gould et al. 2010).

Our simulations can be scaled to match the combined semimajor axis-eccentricity distribution of giant planets beyond 0.5 − 1 AU. At a smaller orbital distance, a giant planet’s scattering power decreases. This can be quantified by the quantity Θ2, the ratio of the escape speed from the planet’s surface to the escape speed from the planetary system at that location (Safronov 1969; Goldreich et al. 2004): (11)where Mp and M   are the planetary and stellar mass, respectively, Rp is the planetary radius and a is the orbital semimajor axis. As Θ2 is inversely proportional to the a, close in giant planets scatter less strongly than at larger distances. When Θ2 drops below 1, collisions become more important than scattering. For our simulations, only the lowest-mass giant planets would drop to Θ2 < 1 by scaling them to 0.5 − 1 AU. In addition, the planet-planet scattering mechanism has been proven at this range of distances (Marzari & Weidenschilling 2002; Jurić & Tremaine 2008). Thus, we conclude that it is dynamically appropriate to scale the giant planets inward to 0.5 − 1 AU. We cannot, however, scale to closer distances because the dynamical regime is different and giant planet-planet collisions are more likely than scattering events.

We assume that the underlying distribution that is being probed by current (mainly radial velocity) observations rises linearly from zero at 0.5 AU to 1 AU, where it flattens off and remains constant to 5 AU. We draw randomly from this distribution and scale the innermost giant planets from ten randomly chosen simulations – nine that were unstable and one that was stable according to relative contribution of stable vs. unstable systems in our sample – to match this value. We then re-scale the terrestrial planets in each system to the same size scale as the inner giant planet’s orbit. Our sample was created by repeating this 100 000 times with different random numbers.

Figure 20 shows the radial distribution of terrestrial and giant planets after this scaling. In essence, this figure shows the expected radial distribution of terrestrial planets in the known extra-solar giant planet systems. In these systems, our simulations predict a factor of a few higher abundance of terrestrial planets at a few tenths of an AU than at 1 AU because, given the typical giant-terrestrial planet spacing, the formation of planets at 1 AU requires distant giant planets that are hard to detect by current methods. Planets within  ~0.1 AU are sparsely populated because of the assumed inner edge of the embryo disk at 0.5 AU. The peak in the frequency of terrestrial planets at a few tenths of an AU is due to a combination of our inner disk edge and the giant planets’ radial distribution.

As noted above, the scaling we performed is dynamically permissible: it does not change the regime of accretion of the terrestrial planets nor the scattering regime of the giant planets. However, our simulations each contained a constant mass in embryos and planetesimals in the inner disk. By scaling the giant and terrestrial planets, we are effectively re-scaling the initial disk masses such that the same amount of mass would initially have been placed into an annulus whose position and width can vary. In addition, if we had included initially closer-in material in our simulations, the peak in Fig. 20 would have shifted inward. Despite these limitations, in the regime that we have considered we think that the shape of the curve is physical, and we predict that, at least within the known sample of extra-solar giant planets, terrestrial planets at  ~0.3 AU should be several times more abundant than at 1 AU.

5.2. Comparison between our simulations and observed debris disks

Observations (mainly with NASA’s Spitzer Space Telescope) have shown that roughly 15% of solar-type stars younger than 300 million years have measurable dust fluxes at 24   μm (Gáspár et al. 2009) but that this fraction decreases in time and flattens off at  ~3% (Carpenter et al. 2009). At 70   μm, 16% of stars observable dust and there is no measured decrease in this fraction with age (Trilling et al. 2008; Carpenter et al. 2009). Considering the current exoplanet/debris disk sample, 9%  ±  3% of planet-hosting stars were detected at 24 or 70   μm compared with 14%  ±  3% for stars without planets (Bryden et al. 2009). An update of that study found debris disks around the same fraction,  ~15%, of stars with and without known giant planets (Kóspál et al. 2009). In addition, the strong correlation between stellar metallicity and the fraction of stars with planets (Gonzalez 1997; Santos et al. 2001; Fischer & Valenti 2005) does not hold for the current sample of debris disks, whose presence is metallicity-independent (Beichman et al. 2006; Greaves et al. 2006; Bryden et al. 2006; Kóspál et al. 2009).

thumbnail Fig. 20

Semimajor axis distribution of simulated terrestrial planets (dashed line) from our set of simulations, derived by scaling the innermost surviving giant planet in each simulation to match an assumed underlying distribution for relevant exoplanets that increases linearly from zero at 0.5 AU and is constant from 1 − 5 AU.

Open with DEXTER

Figure 13 shows that, in our simulations, the dust flux is anti-correlated with the giant planet eccentricity. Almost all lower-eccentricity (e < 0.1−0.2) giant planets are in systems with debris disks, but at higher eccentricities the fraction of dusty systems decreases as does the dust brightness itself. Figure 18 shows that the fraction of systems that are detectable at all wavelengths from 5   μm to 500   μm decreases with increasing eccentricity of the innermost surviving giant planet. In addition, the fraction of planets with eg > 0.1 decreases with F/Fstar(70   μm), including a dramatic drop for F/Fstar > 10.

We therefore expect to see a correlation between the orbital properties of exoplanets and the detectability of cold dust. The dataset of Bryden et al. (2009) detected debris disks at 70   μm around 13 planet-hosting stars out of 146 for a detection frequency of . We arbitrarily divide the sample in two based on the eccentricity of the innermost giant planet in each system at a value of 0.2. Debris disks are detected in 8 of 76 systems () in the low-eccentricity subsample and 5 of 70 systems () in the high eccentricity subsample. Although the detection rate is somewhat higher in the lower-eccentricity subsample, the difference is not statistically significant and we must wait for better statistics from larger surveys.

Our simulations do produce some systems with high-eccentricity giant planets and bright dust emission (Fig. 13), in agreement with the detected debris disks in known exoplanet systems (Moro-Martín et al. 2010). In these cases the dynamical instability tends to be asymmetric and confined to the inner planetary system and these are therefore not generally good candidates for terrestrial planets, which also agrees with the observed systems (Moro-Martín et al. 2010). The outcome of a given system depends critically on the details of the instability, which is determined by the giant planet masses (Raymond et al. 2010).

Our approach is not unique; other approaches based on dust production by collisional cascades can also reproduce the debris disk observations (Krivov et al. 2005, 2006; Wyatt et al. 2007b; Wyatt 2008; Kenyon & Bromley 2008, 2010; Krivov 2010; Kennedy & Wyatt 2010). In addition, there is a qualitative difference between models that consider planetesimal disks to be “self-stirred” (i.e., eccentricities are excited by accreting bodies within planetesimal disks; Kenyon & Bromley 2008, 2010; Kennedy & Wyatt 2010) or those that consider external sources for planetesimal stirring (e.g. Wyatt et al. 2010). Mustill & Wyatt (2009) showed that self-stirred planetesimal disks tend to be fainter than disks stirred by giant planets. Debris disks may be explained by some combination of these ideas.

Nonetheless, we conclude that our simulations are consistent with the known sample of debris disks in exoplanet systems. However, our initial conditions are biased in that all systems that could potentially produce debris disks also contain giant planets, which is not consistent with the observation that the frequency of debris disk + giant exoplanet systems is about the same as debris disks with no detected giant exoplanets (Bryden et al. 2009; Kóspál et al. 2009). In Paper II we use multiple sets of simulations to construct a sample that adequately matches the entire debris disk and exoplanet samples, and use that sample to infer the properties of as-yet-undetected terrestrial exoplanets.

5.3. Our Solar System in context

Our results suggest that the Solar System is unusual at the  ~15−25% level. This corresponds to the fraction of simulations that form three or four terrestrial planets. (Figure 9; including a 10% contribution from stable systems. By the same weighting 38% of simulations destroy all their terrestrial material.)

The Solar System lies at the very edge of the debris disk correlations in Figs. 13 and 14 because of its combination of a rich terrestrial planet system, a low-eccentricity innermost giant planet, and a low dust flux. To a distant observer, the Solar System’s faint debris disk would suggest a dust-clearing instability in the system’s past. However, Jupiter’s low-eccentricity orbit would imply that the instability was weak and that the system may in fact be suitable for terrestrial planets8. This naive argument is remarkably consistent with our current picture of the LHB instability as an asymmetric, outward-directed instability that included a scattering event between Jupiter and an ice giant but not between Jupiter and Saturn (Morbidelli et al. 2010).

The Earth provides an interesting test case. On long timescales Earth’s eccentricity oscillates between 0.0002 and 0.058 and its inclination between zero and 4.3° (Quinn et al. 1991; Laskar et al. 1993). When compared with the stable simulations, the Earth’s time-averaged eccentricity is significantly smaller than the median and its inclination is also smaller. However, Earth’s oscillation amplitudes are more than 50% larger than the median values for the stable systems. This situation is the same for Venus’ orbit. This presents a confusing situation; perhaps Earth and Venus’ e and i are lower than these simulations because we are limited in the numerical resolution needed to adequately model dissipative processes. But if that were the case, we would expect Earth and Venus’ oscillation amplitudes to also be smaller than the simulated planets’. One explanation for the Earth and Venus’ orbits is that they formed in a dissipative environment but that were later dynamically perturbed during the instability that caused the late heavy bombardment (Gomes et al. 2005; Brasser et al. 2009; Morbidelli et al. 2010). The perturbation was not large enough to disrupt the system’s stability but sufficient to increase the amplitude of orbital oscillations of Earth and Venus.

The instability that caused the LHB is not captured in our simulations nor in current exoplanet observations. The degree to which our simulations interpret that the Solar System is unusual depends on how well the simulations characterize the instabilities in extra-solar planetary systems. As our simulations reproduce the observed exoplanet eccentricity distribution with no free parameters, it appears that our simulations do in fact capture the essence of instabilities among the known exoplanets. However, an instability such as the one proposed by the Nice model leaves little to no trace because the giant planets’ eccentricities remain very small (eJ and eS are only  ~0.05). It is plausible that Nice model-type instabilities are common in outer planetary systems, although if they systematically destroy their outer planetesimal disks then debris disk statistics constrains the fraction of stars that undergo such instabilities more than about 10 Myr after stellar formation to be less than about 10% (Booth et al. 2009). In Paper II we present an additional set of simulations with larger mass ratios between the giant planets in which Nice model-type instabilities can occur. Such instabilities do not change our conclusions; Nice model instabilities can effectively be lumped in with the stable systems as their impact on inner planetary systems is small9.

5.4. Limitations of our approach

As with any numerical study, our simulations are limited in several ways.

Our most important assumption is that the inner and outer regions of protoplanetary disks are connected such that observations of debris disks in the outer parts of planetary systems can tell us something about terrestrial planets in the inner parts of these systems. However, there exist substantial uncertainties. For example, several relevant processes – notably the formation of planetesimals as well as giant planets – are only modestly-well understood. In addition, it is unknown whether the efficiencies and timescales of these processes vary with distance from the star. If there exists a systematic bias to create an imbalance between the inner and outer disk mass, it could qualitatively affect our results. For example, if planetesimal formation is much more efficient in the outskirts of planetary systems then the typical system may contain an outer planetesimal belt but no inner planetesimals or embryos from which to form terrestrial planets. Alternately, one can imagine that outer planetesimal disks might be systematically destroyed in systems that do form terrestrial planets, as was the case for the Solar System. If one of these scenarios is true, then the initial conditions for inner and outer planetary systems may not be coupled as strongly as we have assumed, and their outcomes may not correlate as strongly as in our simulations.

Despite these uncertainties, we think that our approach, in particular the assumption that inner and outer disks are connected, is the simplest interpretation of current observations and theory. The frequency of close-in planets is 12% in the 3 − 10 M range and, by extrapolation, 23% in the 0.5−2 M range (Howard et al. 2010). The observed frequency of debris disks around FGK stars of 16% (Trilling et al. 2008) represents a lower limit for the frequency of outer planetesimal disks. Preliminary results from the Herschel DUNES survey (Eiroa et al. 2010) suggest that roughly 1/3 of stars have debris disks (C. Eiroa, pers. comm.), which means that that the frequency of inner planets and outer planetesimal belts are probably within a factor of a few or less. Although this certainly does not prove that our initial conditions are correct, it does provide circumstantial support for our basic assumption of a connection between the inner and outer disk although we note that there is as yet no observed correlation between the two.

Our initial conditions, though chosen to match models of earlier phases of planetary growth, are ad hoc. All planetary systems in our sample contain the same mass in terrestrial embryos and planetesimals (9   M), form their innermost giant planet at 5.2 AU, and contain the same mass in an outer disk of leftover icy planetesimals (50   M). In reality, there is a spread of several orders of magnitude in the disk mass (e.g., Andrews & Williams 2007a) that affects both the types of planets that form and the amount of leftover planetesimals (e.g., Greaves et al. 2007; Thommes et al. 2008). In addition, recent observations suggest that low-mass disks are more compact (Andrews et al. 2010) so there may be a correlation between the disk mass and the location of planet formation as well (see also Kennedy & Kenyon 2008). The disks that we modeled are comparable in mass to the minimum-mass solar nebula model and are probably more massive than the typical disk (Eisner & Carpenter 2003; Andrews & Williams 2007a).

Our simulations are confined to Solar-mass stars, which are a small minority of all stars (Bochanski et al. 2010). To expand on this work we should consider additional stellar types. Debris disks are currently very difficult to detect around low-mass stars (Gautier et al. 2007) but are extremely interesting as planet hosts because they dominate the stellar population of the Galaxy and are very long-lived. In contrast, debris disks are much easier to detect around A stars, but these are relatively few in number and their lifetimes are much shorter. There may be interesting differences in the evolution of planetary systems around other stellar types.

There are several physical processes not included in our simulations. In particular, we did not include the effects of giant planet migration because population synthesis models are currently unable to reproduce the observed exoplanet mass and orbital distribution (Howard et al. 2010). Including migration would have the benefit of providing a natural trigger for giant planet instabilities (Adams & Laughlin 2003; Moorhead & Adams 2005), although this depends on the details of the depletion of the gaseous disk (Crida et al. 2008). However, given that the giant planet observations can be matched with little to no giant planet migration, we chose not to include it.

6. Conclusions

Our main results are as follows:

  • Giant planet instabilities are destructive to terrestrial planetformation. The survival of terrestrial embryos and planetesimalsdepends on the strength of the instability as measured by theminimum giant planet perihelion distance (Fig. 8,see also Veras & Armitage 2006).In about 40% of our unstable simulations all rockymaterial was removed from the system, in largepart by being thrown into the host star. About 1/5 ofour unstable simulations produced a system containing a singleterrestrial planet (Fig. 9).

  • Terrestrial exoplanets on excited orbits should be common. The median eccentricity of surviving terrestrial planets in our simulations was about 0.1, but the distribution extends above 0.5 (Fig. 7). The most excited orbits belong to single-terrestrial planet systems. Compared with systems with many terrestrial planets, single-planet systems have only slightly higher eccentricities and inclinations but their oscillations in e and i are far larger (Fig. 10).

  • Debris disks are anti-correlated with strongly-scattered giant planets. Strong scattering events produce eccentric giant planets with large radial excursions that dynamically deplete the outer planetesimal disk by exciting their orbits until they cross a giant planet’s at which point they are quickly ejected from the system. Thus, we expect continued observations to show an anti-correlation between the fraction of systems with debris disks and the giant planet eccentricity.

  • Debris disks correlate with a high efficiency of terrestrial planet formation. Strong scattering events yield large giant planet eccentricities, and these eccentric giant planets tend to disturb both the inner and outer planetary system. Thus, in a strongly perturbed system the giant planets tend to destroy both terrestrial planetary embryos – aborting the growth of terrestrial planets – and outer planetesimals – preventing the creation of debris disks by long-term collisional evolution. In contrast, in a calm system the giant planets will not strongly impede on the inner or outer planetary system, allowing for the formation of terrestrial planets and long-lasting cold dust. The debris disk-terrestrial planet correlation holds for all wavelengths we tested but is clearer for λ ≳ 15   μm (Figs. 15 and 16). The correlation also holds for all times later than about 10 − 30 Myr (Fig. 17) and probably even earlier.

  • Within the known sample of extra-solar giant planets, terrestrial planets at a few tenths of an AU should be several times more abundant than either terrestrial or giant planets at 1 AU. In Sect. 5.1 we scaled the outcomes of our simulations to match the observed semimajor axis and eccentricity distributions of giant exoplanets. This scaling produced a radial distribution of terrestrial exoplanets that peaks at a few tenths of an AU and drops below the giant planet frequency at 1.3 AU (Fig. 20). However, we note that the distribution is incomplete in its inner regions due to our initial conditions, in particular the inner edge in our embryo distribution at 0.5 AU.

  • The Solar System lies at the outskirts of several correlations, probably because of the instability that caused the Late Heavy Bombardment. The Solar System has a rare combination of multiple terrestrial planets, low-eccentricity giant planets, and very low dust content (currently, F/Fstar ≈ 2  ×  10-2; Booth et al. 2009). In addition, Earth and Venus’ orbits are more circular and coplanar than terrestrial planets in typical stable planetary systems but they have significantly larger-amplitude oscillations in those quantities. This can be explained if the Solar System’s formation was quiescent but it underwent a later punctual event that was strong enough to remove most of the outer planetesimal disk and give the inner Solar System a small kick but did not destabilize the inner Solar System or impart a large eccentricity to Jupiter. This is in agreement with the current picture of the instability that caused the LHB (Morbidelli et al. 2010). This type of instability is poorly constrained by observations, is much weaker than the instabilities inferred from the exoplanet eccentricity distribution, and is not captured in the simulations presented here.

One implication of our results is that exo-zodiacal dust clouds around stars with terrestrial planets may often be brighter than the Solar System’s zodiacal cloud. A major obstacle to the direct imaging of terrestrial exoplanets is the amount of bright dust close to those planets, i.e., their exo-zodiacal dust clouds (Cash 2006; Defrère et al. 2010; Noecker & Kuchner 2010). The Solar System’s zodiacal dust has been shown to derive almost entirely from Kuiper Belt comets that were scattered by the giant planets into the inner Solar System, where they partially sublimated to produce warm dust before eventually being ejected (Nesvorný et al. 2010). Around other stars, cold debris disks should trace the same population of comets that produces exo-zodiacal dust: debris disks represent planetesimals on stable orbits in the outer system and exo-zodiacal dust is generated by the small fraction of bodies that has been destabilized and is in the process of being removed from the system. Our results suggest that systems with bright debris disks are excellent targets in the search for terrestrial exoplanets. These systems contain at least a few M (and often more than 20 M) in surviving cometary material, 1 − 2 orders of magnitude more than the current Kuiper Belt (Gladman et al. 2001). If the comet flux scales with the number of outer planetesimals then systems with bright debris disks should also harbor bright exo-zodiacal clouds close to the terrestrial planet zone. However, the dynamics of the outer planetary systems – in particular the architecture and masses of the giant planets – are key in determining the fluxes of new comets in these systems as well as their residence lifetimes in the inner planetary systems. In addition, there is almost certainly a significant population of systems with terrestrial planets without bright debris disks, i.e. what we have called “false negatives”. Those systems may be harder to find because they lack debris disks to signpost the presence of terrestrial planets, but they could prove easier to image because they may contain much fainter zodiacal clouds. We plan to test these ideas in future work.

In a companion paper (Raymond et al. 2011; referred to in the text as Paper II) we explore the effect of several other parameters on the results obtained here. In particular, we present results of several other sets of simulations that vary the giant planets’ mass distribution and total masses, the width of the outer planetesimal disk, the existence of icy embryos within the outer planetesimal disk, and the presence of disk gas at the time of giant planet instabilities. In that paper we confirm the large-scale results presented here but with several important clarifications and dependencies. We also carefully match our simulations to the observed statistics of giant exoplanets and debris disks to obtain an estimate for the fraction of stars that host terrestrial planets, ηEarth in the famous Drake equation.

In a second companion paper (Raymond et al. 2011b), we explore the fate of bodies – planetesimals, planetary embryos, and giant planets – that are ejected from unstable planetary systems and pollute the galaxy. By matching giant exoplanet and debris disk statistics we estimate the abundance of this population of free-floating bodies and the chances for their unambiguous detection either in interstellar space or entering the Solar System.

7. Appendix A: numerical tests

We performed simple numerical tests to choose an inner boundary appropriate for this timestep by placing a particle at 1 AU on an orbit that was highly inclined with respect to a distant, Jupiter-mass planet10. In this configuration, the Kozai effect (Kozai 1962) drives a dramatic increase in the particle’s eccentricity. By tracking the integrator error at ever-smaller perihelion distance, q = a(1 − e), we saw that the fractional error in the semimajor axis da/a increased to greater than 10-4 inside roughly 0.2 AU, as shown in Fig. 21. We therefore chose 0.2 AU as our inner particle boundary, inside which objects are removed from the simulation and assumed to have collided with the star. We chose 100 AU as our corresponding outer boundary, beyond which bodies are assumed to have been dynamically ejected from the system.

thumbnail Fig. 21

Fractional error in the semimajor axis a of a Mars-mass particle with initial a = 1 AU as a function of perihelion distance q. The particle was initially placed on a circular, nearly polar orbit (initial inclination of 89.9°) in the presence of an outer giant planet. The Kozai effect forced the particle to fall into the star, and we tracked the integration error in time.

Open with DEXTER

thumbnail Fig. 22

Fractional error in the system’s energy budget dE/E as a function of the smallest perihelion distance of a giant planet for all simulations.

Open with DEXTER

thumbnail Fig. 23

Eccentricity of the innermost giant planet (left panel) and dust-to-stellar flux ratio at 70 μm after 1 Gyr of evolution (right panel) as a function of the fractional error in the system’s energy budget dE/E. This plot is only for systems with minimum giant planet perihelia of less than 2 AU (see Fig. 22). Our energy error cutoff at 0.01 is shown with the dashed line. There is no evidence for contamination of our sample.

Open with DEXTER

With our choice of timestep and inner boundary, the majority of our simulations maintained excellent (dE/E < 10-4) energy and angular momentum conservation properties11. However, as shown in Fig. 22, a fraction of the simulations in which the minimum giant planet perihelion distance was small (q < 2 AU) exhibited substantially worse conservation properties. Given this behavior, we adopted an empirical energy conservation threshold of dE/E < 1  ×  10-2, and discarded from the final sample any runs that failed to meet this limit.

Does our cutoff in energy at dE/E < 1  ×  10-2 bias our results? The true outcome of our numerical experiment depends on the details of the giant planets’ orbital evolution, and it is only relatively extreme cases with q ≤ 2 AU for which significant integration error occurs. It is for these close perihelion passages that all terrestrial material is destroyed, called mode 2 accretion in the main text. In the case of the five simulations with dE/E > 1  ×  10-2, two contained a single surviving terrestrial planet and three had destroyed all of their terrestrial planets. By removing these cases we are therefore slightly biasing our sample at the 5% level away from accretion modes 2 and 3, and so we include this small extra contribution in our estimates in the paper of the fraction of systems for which the different accretion modes occur.

We do not see any clear signature of other bias introduced in our analysis from our energy cutoff. To assess this possibility, we restrict ourselves to systems for which the minimum giant planet perihelion was less than 2 AU because this is where large errors occur. Figure 23 shows the giant planet eccentricity and the dust-to-stellar flux ratio at 70 μm after 1 Gyr as a function of dE/E for this subsample. Both panels are scatter plots, with no clear trend or any indication that the computed dE/E value changes the outcome in any way. We therefore conclude that our chosen cutoff in integrator energy, while less stringent than some other studies, has no measurable impact on our results.


1

Recent calculations of planet formation via core accretion support such an assumption (Movshovitz et al. 2010), though it should be noted that these are subject to substantial uncertainties because the physics of type I migration remains imperfectly understood (Lubow & Ida 2010).

2

Our procedure of decreasing the particle mass according to the collisional timescale is not completely self-consistent as our planetesimal mass did not change during the simulation itself. Given the long collisional timescales in the outer disk, this assumption has little to no effect on the outer planetary systems. This assumption also has little effect on cases in which the giant planets were unstable quickly. However, for stable simulations or simulations with delayed instabilities, this means that, in the inner disk, our dynamics is not completely true to our calculated dust flux. In other words, the amount of damping provided to growing embryos by planetesimals (via viscous stirring and dynamical friction) should realistically have been lower if we had accounted for collisional evolution of planetesimals. However, given that we compare our simulations with observed debris disks that are generally far older than the 10 − 100 Myr timescale for terrestrial planet formation, this does not affect our results. For a careful treatment of dust production during terrestrial planet formation, see Kenyon & Bromley (2004).

3

Small grains with sizes of 2.2   μm do not actually radiate as blackbodies. Rather, the effective temperature of small grains is likely to be higher than the blackbody temperature. The consequence of this in the context of our model is discussed extensively in Bonsor & Wyatt (2010) and Kains et al. (2011). To summarize, the actual flux may be slightly higher or lower than that predicted with our assumptions, since on the one hand the higher temperature means that we underestimate the dust flux for a given dust mass, however this is compensated to some extent by the fact that we also overestimate the evolutionary timescale (since the evolutionary timescale for a disk of observed temperature is set observationally).

4

In the simulation, any body that came within 0.2 AU of the star is considered to have collided with it. It is conceivable that in some cases, star-grazing embryos could have their orbits shrunk and re-circularized by tidal effects, although the efficiency of this process is uncertain and probably quite small (Raymond et al. 2008b).

5

We do not consider regeneration of small bodies via giant impacts but we note that the debris from giant collisions could cause short-lived spikes in the dust brightness, in particular at mid- to near-infrared wavelengths (Stern 1994; Grogan et al. 2001; Kenyon & Bromley 2005; Lisse et al. 2009). These peaks in brightness from collisions can only occur within a few AU because a very massive collision is needed to cause substantial brightening.

6

We expect that numerical resolution should play a role here. Given that they contain only 500 planetesimal particles, our simulations cannot fully resolve dynamical friction at late times (e.g. O’Brien et al. 2006; Raymond et al. 2006b). We expect that eccentricities and inclinations are somewhat overestimated in the stable systems, but since the instabilities remove most of the planetesimals, not in unstable systems.

7

Note that one stable system ejected all of its planetesimals. In that system the interactions between the giant planets excited eccentricities of  ~ 0.1, causing an eventual complete depletion in the outer planetesimal disk and a decreased terrestrial planet formation efficiency compared with other stable simulations. Although the planets were all roughly one Jupiter mass, the two inner planets were actually also driven slowly across their 2:1 mutual mean motion resonance after 130 Myr (which actually decreased the eccentricities).

8

Of course, it would take about a decade of precise observations for these aliens to pin down Jupiter’s orbital eccentricity.

9

Brasser et al. (2009) showed that some Nice model instabilities can destabilize the orbits of the terrestrial planets due to sweeping secular resonances and potentially cause collisions between the terrestrial planets. However, this process does not remove all the terrestrial planets and so, in the context of our results, nothing has changed.

10

We thank Hal Levison for suggesting this numerical test.

11

We are aware that removing particles from the simulation, at any radius larger than the physical one (approximately the size of the star), can in principle result in unphysical behavior, even when energy and angular momentum are well-conserved. Unfortunately, it is currently infeasible to run long-duration terrestrial planet formation simulations with dramatically shorter timesteps and smaller inner boundary radii.

Acknowledgments

We thank the referee, Scott Kenyon, for a very thorough review 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). Some of the collaborations vital to this paper (in particular, between S.N.R, P.J.A., A.M.-M., M.B. and M.C.W.) started during the Isaac Newton Institutes Dynamics of Disks and Planets program in Cambridge, UK. S.N.R. acknowledges funding from the CNRS’s PNP and EPOV programs and NASA Astrobiology Institute’s Virtual Planetary Laboratory lead team. P.J.A. acknowledges funding from NASAs Origins of Solar Systems program (NNX09AB90G, NASAs Astrophysics Theory and Fundamental Physics program (NNX07AH08G), and the NSFs Division of Astronomical Sciences (0807471). This paper is dedicated to office B329 in the University of Washington’s Physics and Astronomy building, which from 2000 − 2003 housed S.N.R., J.C.A., and A.A.W. and was the site of many ridiculous pursuits.

References

Online material

Movie associated to Fig. 2 (Access here)

All Figures

thumbnail Fig. 1

An example of the initial conditions assumed to exist shortly after the epoch of gas disk dispersal. The terrestrial planet formation zone contains planetary embryos and lower mass planetesimals, while the outer planetesimal belt contains a population of relatively low-mass (compared to the giant planets) bodies that are represented numerically by N equal mass particles. The intermediate zone contains 3 fully-formed giant planets whose spacing is such that they are marginally unstable against dynamical instability.

Open with DEXTER
In the text
thumbnail Fig. 2

Evolution of a system with a relatively late (21 Myr) instability among the giant planets. Each panel shows a snapshot in time of orbital eccentricity vs. semi-major axis for all particles; vertical bars denote sin(i) for terrestrial bodies with Mp > 0.2   M and i > 10°. The particle size is proportional to the mass1/3, but giant planets (large black circles) are not on this scale. Colors denote water content, assuming a Solar System-like initial distribution (Raymond et al. 2004). The surviving terrestrial planet has a mass of 0.72 M, a stable orbit within the habitable zone (semimajor axis of 0.96 AU), and a high eccentricity and inclination (and large oscillations in these quantities). A movie of this simulation is available in the electronic edition of the Journal.

Open with DEXTER
In the text
thumbnail Fig. 3

The eccentricity (top panel) and inclination (bottom panel) evolution of the two surviving giant planets in the example simulation. The three curves in the right panel show the inclination of each planet relative to the starting plane of the system – presumably the stellar equatorial plane – in black and light grey and their mutual inclination, which is larger than each of their individual inclinations, in dark grey.

Open with DEXTER
In the text
thumbnail Fig. 4

The eccentricity (top panel) and inclination (bottom panel) of the surviving terrestrial planet in the example simulation (mass of 0.72 M, semimajor axis of 0.93 AU). The two curves in the right panel show the inclination relative to the starting plane of the system in grey and the mutual inclination with respect to the innermost giant planet in black.

Open with DEXTER
In the text
thumbnail Fig. 5

Evolution of the spectral energy distribution (SED) of the simulation from Fig. 2. Each curve shows the SED during a given simulation snapshot. The instability occurred at 21 Myr, and the SED evolved dramatically in the immediate aftermath, as icy planetesimals were scattered onto high-eccentricity orbits – thereby producing transient hot dust – and then ejected from the system. The dashed line represents the stellar photosphere. The system is viewed at a distance of 10 pc.

Open with DEXTER
In the text
thumbnail Fig. 6

The ratio of the dust-to-stellar flux as a function of time for six different wavelengths from 5 to 500 microns, including a zoom during the instability – note that each main plot is on a log scale and each zoom-in is on a linear scale. The rough observational limits of the MIPS instrument on NASA’s Spitzer Space Telescope are shown for 25 and 70 μm with the dashed line Trilling et al. (2008). All planetesimal particles were destroyed as of 45 Myr via either collision or ejection so there is no more dust in the system after that point.

Open with DEXTER
In the text
thumbnail Fig. 7

Eccentricity distributions of the surviving giant (solid black line) and terrestrial (dashed black line) planets in the unstable systems. The grey curve represents the known extra-solar planets beyond 0.2 AU – this helps to exclude planets whose orbits have been tidally circularized as well as low-mass planets.

Open with DEXTER
In the text
thumbnail Fig. 8

The total mass in surviving terrestrial planets as a function of the minimum perihelion distance of a giant planet during the simulation. Black dots represent systems in which the giant planets were dynamically unstable and grey dots are systems that were stable. Filled black dots are systems in which a single terrestrial planet formed.

Open with DEXTER
In the text
thumbnail Fig. 9

Distribution of the number of surviving terrestrial planets (defined as Mp > 0.05   M) in our fiducial mixed set of simulations. The grey shaded histogram shows the unstable simulations and the dashed line shows the stable simulations.

Open with DEXTER
In the text
thumbnail Fig. 10

The peak-to-peak oscillation amplitudes in eccentricity e and inclination i for the surviving terrestrial planets in simulations with stable (grey) and unstable (black) giant planets. The single-terrestrial planet systems are shown with filled black circles, and Earth is the grey star.

Open with DEXTER
In the text
thumbnail Fig. 11

The total mass in surviving terrestrial planets as a function of the eccentricity of the innermost surviving giant planet. The Solar System is shown as the grey star.

Open with DEXTER
In the text
thumbnail Fig. 12

Distribution of eccentricities of the innermost giant planet in simulations which formed zero (blue line), one (green), and two or more (red) terrestrial planets. The sum of the three distributions (black dashed line) provides a quantitative match (p = 0.49 from a Kolmogorov-Smirnov test) to the observed exoplanets (thick grey line).

Open with DEXTER
In the text
thumbnail Fig. 13

The dust-to-stellar flux ratio at 70 μm after 1 Gyr of dynamical and collisional evolution, plotted as a function of the eccentricity of the innermost surviving giant planet. Unstable systems are plotted in black and stable systems in grey (systems with a single surviving planet are shown with solid circles). The pileups very close to the x axis represent systems with virtually zero 70 μm flux. The dashed line represents an approximate threshold above which excess emission was detectable using Spitzer data (Trilling et al. 2008). The star shows the estimated flux from the Kuiper Belt at 1 Gyr, as calculated previously (Booth et al. 2009) based on dynamical simulations of the outer Solar System (Gomes et al. 2005).

Open with DEXTER
In the text
thumbnail Fig. 14

The dust-to-stellar flux ratio at 70   μm after 1 Gyr of dynamical and collisional evolution, plotted as a function of the total mass in terrestrial planets. Unstable systems are plotted in black and stable systems in grey. The pileups close to the x and y axes represent systems with either no terrestrial planets or virtually zero 70   μm flux. The dashed line represents an approximate threshold above which excess emission was detectable using Spitzer data (Trilling et al. 2008). The star shows the estimated flux from the Kuiper Belt at 1 Gyr, as calculated previously (Booth et al. 2009) using models based on dynamical simulations of the outer Solar System (Gomes et al. 2005). The Solar System falls into an intermediate region of parameter space: the giant planets may have been unstable, but only weakly so (there is no clear evidence for ejections or planetary collisions), while the Kuiper Belt would have been bright for several hundred Myr prior to the Late Heavy Bombardment.

Open with DEXTER
In the text
thumbnail Fig. 15

The dust-to-stellar flux ratio F/Fstar after 1 Gyr of dynamical and collisional evolution as a function of the total mass in terrestrial planets, for six wavelengths between 5 and 500 microns. The Spitzer observational limits are shown for 25 and 70 μm with the dashed line (Trilling et al. 2008). In each panel, grey circles represent stable simulations and black circles represent unstable simulations.

Open with DEXTER
In the text
thumbnail Fig. 16

Left: the fraction of systems that would be detectable after 1 Gyr as a function of the total mass in surviving terrestrial planets for six different wavelengths between 5 and 500   μm. Each curve represents a horizontal slice through Fig. 15 at a different vertical height. The chosen “detection limits” in F/Fstar were: 10-8 at 5   μm, 10-3 at 15   μm, 0.054 at 25   μm, 0.55 at 70   μm, 10 at 160   μm, and 100 at 500   μm. The error bars are 1-σ and were calculated using binomial statistics (Burgasser et al. 2003). Note that the different curves are offset by up to  ± 0.1   M for clarity. The bins themselves are at: zero, 0−1   M, 1−2   M, 2−3   M, and  > 3   M. Right: the fraction of systems with 0.5   M or more in terrestrial planets as a function of F/Fstar(70   μm) (1 Gyr). 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. 14.

Open with DEXTER
In the text
thumbnail Fig. 17

The dust-to-stellar flux ratio F/Fstar at 70   μm as a function of the total mass in terrestrial planets for eight different times from 1 Myr to 3 Gyr. The terrestrial planet mass refers to the final value such that simulations move vertically in time on the plot. The Spitzer observational limit is shown with the dashed line (Trilling et al. 2008). Grey circles represent stable simulations and black circles represent unstable simulations.

Open with DEXTER
In the text
thumbnail Fig. 18

Left: the fraction of systems that would be detectable after 1 Gyr as a function of the eccentricity of the innermost surviving giant planet eg for six different wavelengths between 5 and 500   μm. Each curve represents a horizontal slice through plots such as Fig. 13 at a different vertical height. The chosen “detection limits” in F/Fstar were the same as in Fig. 16: 10-8 at 5   μm, 10-3 at 15   μm, 0.054 at 25   μm, 0.55 at 70   μm, 10 at 160   μm, and 100 at 500   μm. The error bars are 1-σ and were calculated using binomial statistics. Note that the different curves are offset for clarity. The bins themselves are logarithmically spaced between 10-3 and 1. Right: the fraction of systems with eg > 0.1 as a function of F/Fstar(70   μm) (1 Gyr). 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. 13.

Open with DEXTER
In the text
thumbnail Fig. 19

The p value from a Kolmogorov-Smirnov test comparing the eccentricities of our simulations vs. the observed exoplanet distribution as a function of the fraction of stable systems included in the sample. The p value drops below an acceptable level of 0.01 for more than a 17% contribution from stable systems.

Open with DEXTER
In the text
thumbnail Fig. 20

Semimajor axis distribution of simulated terrestrial planets (dashed line) from our set of simulations, derived by scaling the innermost surviving giant planet in each simulation to match an assumed underlying distribution for relevant exoplanets that increases linearly from zero at 0.5 AU and is constant from 1 − 5 AU.

Open with DEXTER
In the text
thumbnail Fig. 21

Fractional error in the semimajor axis a of a Mars-mass particle with initial a = 1 AU as a function of perihelion distance q. The particle was initially placed on a circular, nearly polar orbit (initial inclination of 89.9°) in the presence of an outer giant planet. The Kozai effect forced the particle to fall into the star, and we tracked the integration error in time.

Open with DEXTER
In the text
thumbnail Fig. 22

Fractional error in the system’s energy budget dE/E as a function of the smallest perihelion distance of a giant planet for all simulations.

Open with DEXTER
In the text
thumbnail Fig. 23

Eccentricity of the innermost giant planet (left panel) and dust-to-stellar flux ratio at 70 μm after 1 Gyr of evolution (right panel) as a function of the fractional error in the system’s energy budget dE/E. This plot is only for systems with minimum giant planet perihelia of less than 2 AU (see Fig. 22). Our energy error cutoff at 0.01 is shown with the dashed line. There is no evidence for contamination of our sample.

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.