TY - JOUR
AU - Stadel,, Joachim
AB - Abstract Terrestrial planets are thought to be the result of a vast number of gravitational interactions and collisions between smaller bodies. We use numerical simulations to show that practically identical initial conditions result in a wide array of final planetary configurations. This is a result of the chaotic evolution of trajectories which are highly sensitive to minuscule displacements. We determine that differences between systems evolved from virtually identical initial conditions can be larger than the differences between systems evolved from very different initial conditions. This implies that individual simulations lack predictive power. For example, there is not a reproducible mapping between the initial and final surface density profiles. However, some key global properties can still be extracted if the statistical spread across many simulations is considered. Based on these spreads, we explore the collisional growth and orbital properties of terrestrial planets, which assemble from different initial conditions (we vary the initial planetesimal distribution, planetesimal masses, and giant planet orbits.). Confirming past work, we find that the resulting planetary systems are sculpted by sweeping secular resonances. Configurations with giant planets on eccentric orbits produce fewer and more massive terrestrial planets on tighter orbits than those with giants on circular orbits. This is further enhanced if the initial mass distribution is biased to the inner regions. In all cases, the outer edge of the system is set by the final location of the ν6 resonance and we find that the mass distribution peaks at the ν5 resonance. Using existing observations, we find that extrasolar systems follow similar trends. Although differences between our numerical modelling and exoplanetary systems remain, we suggest that CoRoT-7, HD 20003 and HD 20781 may host undetected giant planets. chaos, methods: numerical, celestial mechanics, planets and satellites: dynamical evolution and stability, planets and satellites: formation 1 INTRODUCTION Numerical simulations suggest that our Solar system is inherently chaotic and that any small changes to the planets’ initial positions diverge exponentially with e-folding times on the order of 5 to 20 Myr (Sussman & Wisdom 1988, 1992; Laskar 1989, 1994; Quinn, Tremaine & Duncan 1991). However, the difference between chaotic and near-integrable orbits remains within the bounds of measurement uncertainties of the outer planets (Hayes 2008). Although the outer Solar system planets appear to be remarkably stable against developing crossing orbits on Gyr time-scales (Ito & Tanikawa 2002), this may not be the case for the inner Solar system planets (Laskar 2008; Laskar & Gastineau 2009). On such time-scales, chaos is mediated by overlapping resonances (Chirikov 1979; Wisdom 1980; Franklin, Lecar & Wiesel 1984; Laskar 1990; Laskar, Quinn & Tremaine 1992; Morbidelli & Moons 1993; Moons & Morbidelli 1995; Murray & Holman 1997; Nesvorný & Morbidelli 1998) to which minor bodies are particularly sensitive. Analytical and numerical work suggest that overlapping resonances account for the observed distribution of bodies in the asteroid belts (Wisdom 1982, 1983, 1985; Gladman et al. 1997; Moons, Morbidelli & Migliorini 1998), in the inner (Mikkola & Innanen 1995; Evans & Tabachnik 1999) and outer Solar system (Everhart 1973; Lecar & Franklin 1973; Franklin, Lecar & Soper 1989; Gladman & Duncan 1990; Duncan & Quinn 1993; Holman & Wisdom 1993; Holman 1995; Grazier et al. 1999a,b; Morbidelli et al. 2005), as well as within the Kuiper belt (Torbett 1989; Torbett & Smoluchowski 1990; Holman & Wisdom 1993; Levison & Duncan 1993, 1997; Duncan & Levison 1997).1 During the epoch of terrestrial planet formation, the Solar system environment was rather different than today. Set against a backdrop of migrating giant planets (Tsiganis et al. 2005; Morbidelli & Crida 2007; Morbidelli et al. 2007; Levison et al. 2011; Walsh et al. 2011) and being embedded in a dissipating gaseous disc (Haisch, Lada & Lada 2001; Mamajek 2009; Pfalzner, Steinhausen & Menten 2014), planetesimals are thought to grow collisionally and hierarchically into terrestrial planets (Safronov & Zvjagina 1969; Greenberg et al. 1978; Wetherill & Stewart 1989; Greenzweig & Lissauer 1990, 1992; Ida & Makino 1992a,b; Kokubo & Ida 1996, 1998, 2000; Weidenschilling et al. 1997). Planetesimals undergo perturbational encounters with each other (to within a few Hill radii) about once per orbit, whereas perturbations from resonant configurations with giant planets require hundreds of orbital periods to cause noticeable effects. Planetesimal disc dynamics resemble those of stellar systems, in which small orbital perturbations grow exponentially fast (Miller 1964; Kandrup & Smith 1991; Goodman, Heggie & Hut 1993; Valluri & Merritt 2000; Hut & Heggie 2002). This is the essence of stochasticity in planetesimal discs. Numerical simulations probe this regime by tracking the collisional evolution of planetesimals. Such simulations can address the formation and composition of the terrestrial planets (Chambers & Wetherill 1998; Chambers 2001; Raymond, Quinn & Lunine 2004, 2005b; Kokubo, Kominami & Ida 2006; O'Brien, Morbidelli & Levison 2006; Raymond, Quinn & Lunine 2006a; Raymond et al. 2009, hereafter R09; Morishima, Stadel & Moore 2010) and the extrasolar systems (Raymond, Quinn & Lunine 2005a; Raymond, Barnes & Kaib 2006b; Izidoro, Morbidelli & Raymond 2014b; Ogihara, Kobayashi & Inutsuka 2014).2 In these simulations, much of the dynamics of asteroids, planetesimals and embryos in the inner Solar system depends on the presence of giant planets, with interactions mediated most strongly by way of mean motion and secular resonances. The latter depend on the gravitational potential of the gas disc, so that their location sweeps inwards as the gas dissipates, which provides a credible mechanism to dynamically excite bodies in large parts of the inner Solar system. Attempting to reproduce the distribution of eccentricities and inclinations in the asteroid belt, Heppenheimer (1980), Ward (1981), Lecar & Franklin (1997), Nagasawa, Tanaka & Ida (2000) and Nagasawa, Ida & Tanaka (2001) demonstrated the mechanism of resonance sweeping on massless test particles. In the context of terrestrial planet formation, Chambers & Wetherill (1998) and Chambers (2001) numerically demonstrated the effect of secular interactions of embryos with Jupiter and Saturn (as well as among the embryos themselves), although their experiments did not consider a gaseous disc. Later, Levison & Agnor (2003) investigated the effect of various giant planet configurations on planetary embryos and described ‘secular conduction’, which allows embryos at secular resonances to excite embryos in other regions. Adding gas, Kominami & Ida (2004) numerically found that including of a decaying gas disc and giant planets leads to significantly shorter crossing times. Finally, Nagasawa, Lin & Thommes (2005) and Thommes, Nagasawa & Lin (2008) demonstrated that inward secular resonances sweep along embryos to the inner systems as gas dissipates. Extending preceding work, Morishima et al. (2010) demonstrated that the secular resonances also sweep planetesimals inward. In all cases, authors stress that the tendency of sweeping resonances to dynamically excite embryos and planetesimals is balanced by the dampening influence of hydrodynamic drag and dynamical friction (sometimes called gravitational drag). Irrespective of the detailed setup of simulations, their initial conditions (ICs) are usually generated by drawing realizations from some underlying solid mass distribution, possibly subject to stability constraints if planetary embryos are implanted directly (Chambers, Wetherill & Boss 1996; Yoshinaga, Kokubo & Makino 1999). As the system is inherently chaotic, we expect that different ICs drawn from the same underlying distribution will lead to different final systems, much like in simulations of stellar dynamics (Allison et al. 2010; Parker & Goodwin 2012). To date, few contributors have evolved multiple realizations of the same distribution as limited computational resources are typically focused on parameter studies.3 Those that did report distinctly different outcomes for different realizations of the same underlying distribution (Kokubo et al. 2006; R09; Walsh et al. 2011; Izidoro et al. 2014a). Therefore, just how reliable are simulations of the collisional growth of terrestrial planets – do they have predictive power? This question is at the heart of our paper and we tackle it in a twofold manner. First, we evolve identical realizations of a planetesimal disc multiple times. Due to differences in round-off errors, we will see that simulations terminate with different planetary configurations and we assess the statistical spread of several diagnostics. Secondly, we evolve the planetesimal disc in the absence and presence of Jupiter and Saturn (in two configurations). We also vary the initial mass and planetesimal distribution. Again, we compare diagnostics across runs and check whether trends in the diagnostics are visible or buried in the statistical spread. The paper is organized as follows. In Section 2, we outline numerical methods, ICs, the analytic gas model used and describe our simulations. In Section 3, we present results by covering some illustrative examples, describing how the planetesimal discs evolve, what kind of final systems they lead to, describing the driving dynamics, determining whether a simple mapping between initial and final surface density profile persists, and addressing variations in typical diagnostics used in terrestrial planet formation. In Section 4, we extend this by discussing caveats in our dynamical modelling and attempting to match our trends to observations. We conclude the paper in Section 5. 2 N-BODY METHODS, INITIAL CONDITIONS In this work, we use the Graphics Processing Unit (GPU) code Genga to follow the collisional evolution of planetesimal discs. We now describe the code, skirt the issue of program execution order, the ICs of the planetesimal disc, and the gas disc model. 2.1 Genga Genga (Grimm & Stadel 2014) is a hybrid symplectic integrator similar to Mercury (Chambers 1999), but running in parallel on GPUs. The integration scheme treats gravitational interactions between bodies as perturbations of their Keplerian orbits. Genga uses democratic coordinates (heliocentric positions, barycentric velocities) (Duncan, Levison & Lee 1998). This allows the code to separate close encounter pairs from the rest of the system, and integrate them separately with a direct N-body integrator up to machine precision. Outside of close encounters, the bodies are integrated with a symplectic integrator. The hybrid symplectic integrator has excellent energy conservation over a large number of orbits. Accelerations between bodies are computed directly. This requires |$\mathcal {O}(N^2)$| operations, which is more efficiently calculated on a GPU and is more accurate than a tree-based method. Genga supports an analytical gas disc model inherited from a patched version of Pkdgrav (Stadel 2001; Morishima et al. 2010) which we describe in Section 2.3. The code is available online.4 2.2 Forcing the order of program execution To exactly reproduce numerical results, the order of execution of steps within the program must be fixed.5 While this is easily controlled in single threaded applications, multi-threaded programs require additional logic. To ensure reproducibility of numerical experiments and probe the underlying mechanics of orbital divergence, we have implemented such logic in Genga. The most likely source of variations in the order of operations is the parallel sum operation. In Genga, this is implemented as a parallel reduction formula within one thread block and always operates in the same order. As such, all summation operations are excluded as the source of round-off error variations. The only remaining possible source is in the creation of the close encounter list. If a close encounter pair is found, a counting variable is increased through an _atomicAdd() operation. The order of this operation is not well defined across threads. A different order of the close encounter pair list leads to a different order of bodies in the direct N-body integrator. For multiple close encounter groups, this can lead to a different result. We prevent this behaviour through an additional sorting step, which reorders the close encounter list, but induces a performance penalty. The behaviour is controlled at compilation through the SERIAL_GROUPING flag. In this work, all simulations run with this flag disabled because we rely on variations in round-off errors to induce orbital divergence. Tests show that individual runs of Section 3 can be reproduced exactly if we enable SERIAL_GROUPING. We are presently preparing a companion paper which exploits this flag to determine the rate of orbital divergence. 2.3 Initial conditions, gas disc ICs are generated in the same way as in Morishima et al. (2010), where a number of samples (planetsimals of equal mass) are drawn from an underlying distribution of Keplerian elements. This generates a realization with a particular surface density profile and total mass. We generate realizations by drawing 2000 samples such that the radial surface density profile is \begin{equation} \Sigma _\mathrm{{Disc}}(r) = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\Sigma _{\mathrm{{Disc}},0} \left( \frac{r}{1 \, \mathrm{au}}\right)^{-p} &\qquad 0.5 \, \mathrm{au} < r < 4 \, \mathrm{au}, \\ 0 &\qquad \mathrm{otherwise}, \end{array}\right. \end{equation} (1) for total disc mass of \begin{equation} M_\mathrm{{Disc}} = 2 \pi \int _{0.5}^4 \Sigma (r) \, r \, \mathrm{d}r. \end{equation} (2) In this paper, we adopt p = 1 and 1.5 as well as MDisc = 5 MEarth and 10 MEarth (see Section 2.4 and Table 1). All planetesimals have mass M ∼ 0.04 MLunar (∼0.08 MLunar for 10 MEarth discs) and are on nearly circular (e < 0.02), low inclination (i < 0|$_{.}^{\circ}$|75) orbits. The planetesimals are embedded in a gas disc described by an analytical model. The gas surface density follows a power law and decays exponentially in time, i.e. \begin{equation} \Sigma _\mathrm{Gas}(r,t) = \Sigma _{\mathrm{Gas},0} \left( \frac{r}{1 \, \mathrm{au}} \right)^{-1} \exp \left( - \frac{t}{\tau } \right), \end{equation} (3) where τ is the decay time of the gas disc. For all simulations, τ = 1 Myr, ΣGas, 0 = 2000 g/cm2. After ∼4.6 Myr, only 1 per cent of the gas remains. Table 1. Orbital elements of Jupiter and Saturn as well as initial planetesimal disc conditions for our simulation sets. For the giant planets, the three angular arguments are initialized to zero. The CJSset corresponds to the ICs of the Nice model (Tsiganis et al. 2005). The EJS set corresponds to the present-day Solar system. In all cases, we start with 2000 planetesimals. Set . aJ (au)a . eJb . iJ (deg)c . aS (au)a . eSb . iS (deg)c . ΣDiscd . ΣDisc, 0 (g/cm2)e . MDisc (MEarth)f . NRunsg . NJS – – – – – – ∝r−1 6.1 5 12 EJS 5.2 0.048 1.30 9.55 0.056 2.49 CJS 5.45 0.0 0.0 8.18 0.0 0.5 EJS/Steep 5.2 0.048 1.30 9.55 0.056 2.49 ∝r−1.5 8.2 5 EJS/Heavy ∝r−1 12.4 10 CJS/Steep 5.45 0.0 0.0 8.18 0.0 0.5 ∝r−1.5 8.2 5 CJS/Heavy ∝r−1 12.4 10 Set . aJ (au)a . eJb . iJ (deg)c . aS (au)a . eSb . iS (deg)c . ΣDiscd . ΣDisc, 0 (g/cm2)e . MDisc (MEarth)f . NRunsg . NJS – – – – – – ∝r−1 6.1 5 12 EJS 5.2 0.048 1.30 9.55 0.056 2.49 CJS 5.45 0.0 0.0 8.18 0.0 0.5 EJS/Steep 5.2 0.048 1.30 9.55 0.056 2.49 ∝r−1.5 8.2 5 EJS/Heavy ∝r−1 12.4 10 CJS/Steep 5.45 0.0 0.0 8.18 0.0 0.5 ∝r−1.5 8.2 5 CJS/Heavy ∝r−1 12.4 10 aSemi-major axis. bEccentricity. cInclination. dSurface density profile. eSurface density at 1 au. fDisc mass. gNumber of independent runs. Open in new tab Table 1. Orbital elements of Jupiter and Saturn as well as initial planetesimal disc conditions for our simulation sets. For the giant planets, the three angular arguments are initialized to zero. The CJSset corresponds to the ICs of the Nice model (Tsiganis et al. 2005). The EJS set corresponds to the present-day Solar system. In all cases, we start with 2000 planetesimals. Set . aJ (au)a . eJb . iJ (deg)c . aS (au)a . eSb . iS (deg)c . ΣDiscd . ΣDisc, 0 (g/cm2)e . MDisc (MEarth)f . NRunsg . NJS – – – – – – ∝r−1 6.1 5 12 EJS 5.2 0.048 1.30 9.55 0.056 2.49 CJS 5.45 0.0 0.0 8.18 0.0 0.5 EJS/Steep 5.2 0.048 1.30 9.55 0.056 2.49 ∝r−1.5 8.2 5 EJS/Heavy ∝r−1 12.4 10 CJS/Steep 5.45 0.0 0.0 8.18 0.0 0.5 ∝r−1.5 8.2 5 CJS/Heavy ∝r−1 12.4 10 Set . aJ (au)a . eJb . iJ (deg)c . aS (au)a . eSb . iS (deg)c . ΣDiscd . ΣDisc, 0 (g/cm2)e . MDisc (MEarth)f . NRunsg . NJS – – – – – – ∝r−1 6.1 5 12 EJS 5.2 0.048 1.30 9.55 0.056 2.49 CJS 5.45 0.0 0.0 8.18 0.0 0.5 EJS/Steep 5.2 0.048 1.30 9.55 0.056 2.49 ∝r−1.5 8.2 5 EJS/Heavy ∝r−1 12.4 10 CJS/Steep 5.45 0.0 0.0 8.18 0.0 0.5 ∝r−1.5 8.2 5 CJS/Heavy ∝r−1 12.4 10 aSemi-major axis. bEccentricity. cInclination. dSurface density profile. eSurface density at 1 au. fDisc mass. gNumber of independent runs. Open in new tab Particles exchange angular momentum with the gas disc in three ways: (i) hydrodynamic drag due to differences in velocity, (ii) torques arising from spiral density waves launched by massive particles, and (iii) gravitational interactions between particles and the massive disc. Note that we artificially enhance hydrodynamic drag for particles with masses <0.01 MEarth to correct for the large initial planetesimal mass. For more details on the interaction between gas and particles, we refer to section 2.2 in Morishima et al. (2010), which we have implemented verbatim. We stress that these interactions are modelled analytically. Our simulations would benefit from a full hydrodynamic model,6 which may affect some of the results in this work. 2.4 Simulation setup, post-processing We generate planetesimal discs by drawing a single realization for three different distributions: (i) a reference disc, (ii) a disc with a steeper surface density profile, and (iii) a more massive disc. Note that even for steeper planetesimal disc profiles, we retain a gaseous disc with surface density profile ∼r−1, which is consistent with models of dust coagulation (Birnstiel, Klahr & Ercolano 2012). For each realization, we generate two sets of simulations with giant planets – one with Jupiter and Saturn on circular orbits (CJS) and one with Jupiter and Saturn on eccentric orbits (EJS). The EJS set corresponds to the present-day Solar system and the CJS set to the ICs of the original Nice model (Tsiganis et al. 2005). The Jovian planets are inserted at the start of the simulations. We also generate a set based on the reference disc without giant planets (NJS). Table 1 summarizes our ICs. Each set is evolved 12 times for a total of 84 runs covering nine billion steps (t ∼ 147.84 Myr, Δt = 6 d). Computing time per run is about a month on a NVIDIA GeForce GTX 590. We treat particle collisions as inelastic mergers. Particles are removed from the simulations if their heliocentric distance falls below 0.2 au or exceeds 20 au. Given our six day time-step, this makes for at least 5.4 steps per orbit. The relative energy error remains ΔE/E0 < 4.3 × 10−4 at all times. All the outputs are available online.7 In post-processing, we load the outputs from all 12 runs per set, calculate various diagnostics for each run, and plot/tabulate their median, as well as 10/90 and 25/75 percentile spreads across the 12 runs. In the text, we quote only median values and 10/90 percentile ranges. For simulations with giant planets we compute the location of the ν5, ν6, ν15, and ν16 secular resonances to overlay during the analysis. The calculation is implemented as in the appendix of Nagasawa et al. (2000) with the following caveats: (i) the locations are exact only for massless test particles with low eccentricity/inclination orbits, (ii) they are computed to first order and therefore only depend on the semi-major axis a (as well as the masses and semi-major axes of the giants), (iii) the gas disc is modelled as in equation (3), ignoring modifications of the potential by spiral density waves (whose effect is modelled analytically). We thus expect the location of the resonances in the simulation to slightly differ from the values derived in post-processing. Also note that the ν15 resonance does not appear in the region of interest, and that – as the gas dissipates – the ν16 resonance appears in two locations. See also fig. 4 in Nagasawa et al. (2000). For consistency, we apply the same colour-coding of simulation sets to all figures in this paper. Blue colours indicate NJS runs, red colours EJS runs, and green colours CJS runs. Depending on the context, the shade encodes different information. If only one run is shown (or multiple runs are aggregated in an otherwise discernable fashion), the shading correlates with the mass. If multiple runs are aggregated to describe statistics of simulation sets, the shading indicates the statistical spread, i.e. median values as well as relevant percentile offsets. 3 RESULTS We now present results from the runs described above. We begin by illustrating how two identical ICs diverge rapidly and lead to different systems as well as the principal mechanism that drives our planetesimal discs. Afterwards, we address the assembly of embryos from planetesimals, final system architectures, radial mass and mass–frequency distributions, the link between planetesimal distributions and system architecture, the long-term stability of our systems, how resonances sculpt the system, and the spread in frequently used diagnostics. 3.1 Some illustrative examples 3.1.1 Divergence of orbits and runs Fig. 1 shows a time sequence of the semi-major axis and eccentricity for two of the 12 simulations launched without Jupiter and Saturn. By 14.98 Myr, the simulations have clearly evolved different groups of planetary embryos as well as populations of remaining planetesimals. After 147.84 Myr, the simulations terminate with distinctly different terrestrial planets. Figure 1. Open in new tabDownload slide Semi-major axis and eccentricity for two identical ICs (columns) at four time slices (rows). Each marker represents a single body. Larger markers indicate more massive bodies. Figure 1. Open in new tabDownload slide Semi-major axis and eccentricity for two identical ICs (columns) at four time slices (rows). Each marker represents a single body. Larger markers indicate more massive bodies. Closer analysis reveals that initially identical orbits diverge exponentially fast with e-folding times of the order of a few to a few tens of years. Although the particle position and velocity vectors are initially identical, variations in round-off errors across simulation runs induce position variations at the level of floating point accuracy (one part in 1015, or about a millimetre for planets at 1 au) which grow exponentially due to the chaotic nature of the gravitational N-body problem (see e.g. Goodman et al. 1993). On time-scales of a few hundred years, two initially identical simulations diverge fully and the systems will undergo different collisional histories. The rate of divergence of initially nearby planetesimals depends on the mass resolution, i.e. the number of particles initialized to sample a given total disc mass. In general, increasing the number of particles accelerates divergence. This has potential implications for the resolution requirements. We cover this in a companion paper. 3.1.2 Sweeping secular and mean motion resonances Most of our runs host giant planets which primarily interact with the planetesimal population by way of mean motion and secular resonances. While mean motion resonances remain fixed in space (up to orbital variations in the giant planets), secular resonances can sweep through the regions populated by planetesimals as the protoplanetary gas disc dissipates. Doing so, they may shepherd planetesimals along to sculpt the final architecture of the terrestrial planets (Nagasawa et al. 2005; Thommes et al. 2008). Postponing discussion on the influence of sweeping secular resonances, we now recapitulate their origin and illustrate how they influence planetesimal dynamics. Due to their mutual gravitational interaction, the Jupiter–Saturn system has four eigenfrequencies labelled f1, f2, g1, and g2. Linear combinations of these drive secular variations in eccentricity and inclination (Murray & Dermott 1999). They depend on the planetary masses, orbital configuration, and gravitational potential of the gas disc in which they are embedded (Brouwer & Clemence 1961; Heppenheimer 1980; Ward 1981; Nagasawa et al. 2000). Given interactions with Jupiter and Saturn as well as the gravitational potential of the gas disc, planetesimals are subjected to perturbations with frequencies f and g. At particular locations in the disc, these frequencies match, driving resonances - labelled ν5 (g = g1), ν6 (g = g2), ν15 (f = f1), and ν16 (f = f2). They pump eccentricities (ν5, ν6) or inclinations (ν15, ν16). Fig. 2 shows the location of the secular resonances at different time slices. As time moves on, the gas dissipates and the secular resonances sweep inwards, pushing planetesimals in front of them. For the EJS configuration, ν5 and ν6 remain closer together than for the CJS case. In the EJS case, ν5 also settles at smaller semi-major axes a once the gas disappears. The net effect is that more material is delivered to smaller a (where growth time-scale is faster) in EJS configurations. As we will confirm below, this leads to faster assembly and growth of planetary embryos. Figure 2. Open in new tabDownload slide Semi-major axis and eccentricity for a single run with EJS (left) and CJS (right) conditions at six time slices (rows). Circles indicate planetesimals and planets with the size being proportional to the mass. Triangles indicate the location of (from left to right) 3:1, 5:2, 7:3, 2:1, and 3:2 mean motion resonances with Jupiter. Vertical bars indicate the location of the secular resonances ν5 and ν6, which sweep inwards as the gas disc dissipates. We do not show the ν16 resonance which settles ∼0.1 au beyond the ν6 resonance at ∼9.9 Myr. Animations are available at https://cheleb.net/astro/chaos15/media/. Figure 2. Open in new tabDownload slide Semi-major axis and eccentricity for a single run with EJS (left) and CJS (right) conditions at six time slices (rows). Circles indicate planetesimals and planets with the size being proportional to the mass. Triangles indicate the location of (from left to right) 3:1, 5:2, 7:3, 2:1, and 3:2 mean motion resonances with Jupiter. Vertical bars indicate the location of the secular resonances ν5 and ν6, which sweep inwards as the gas disc dissipates. We do not show the ν16 resonance which settles ∼0.1 au beyond the ν6 resonance at ∼9.9 Myr. Animations are available at https://cheleb.net/astro/chaos15/media/. The mean motion resonances occur at semi-major axes where planetesimals and giant planets periodically line up at fixed orbital phases. We indicate five of the lowest order mean motion resonances by triangles in Fig. 2. As the simulations evolve, their location remains approximately fixed. At early times, they are efficient at exciting planetesimals, although eccentricities are rapidly dampened by hydrodynamic drag. As time evolves, planetesimals are cleared from these resonant regions. Since there appears to be no difference in the early-time collision rate across simulations, we conclude that mean motion resonances do not drive the dynamics of collisional growth. They do, however, help to drive ejection of material on time-scales of 10 to 100 Myr. In CJS configurations, planetesimals still populate the regions covered by the 3:1, 5:2, and 7:3 mean motion resonances (2.5 ≲ a ≲ 3.2 au). Over time, angular momentum exchange with Jupiter excites surviving planetesimals in this region on to hyperbolic orbits. This ejects them from the system. 3.2 Disc evolution We now explore how the collisional evolution for a given set of simulations differs. Fig. 3 shows the total mass bound in (i) the disc, (ii) planetary embryos, and (iii) surviving planetesimals. As in Morishima et al. (2010), a planetary embryo is an object with mass M > MCut = 3.3 × 1026 g, which is about the mass of Mercury. Although embryos grow into planets, we do not label them separately. Objects with M < MCut are classified as planetesimals. For all simulation sets, we show the median mass per component. For the NJS runs, we also indicate the 10/90 percentile spread about the median. In all simulation sets, the spreads are comparable in magnitude to the NJS case, so we omit them for the other sets to remove visual clutter. Tables 2 and 3 tabulate the ranges indicated in Fig. 3. In Fig. 4, we additionally indicate the median amount of material that is ejected from the system or falls on to the star along with the 10/90 percentile range. Figure 3. Open in new tabDownload slide Remaining solid mass (total, planetesimals, embryos) over time as a fraction of the total initial mass. We show the mass as a fraction of the initial total solid disc mass. In general, we show the median across 12 runs. For runs without Jupiter and Saturn, we also show the 10/90 percentile spread about the median. The spread is similar for all runs, so we do not indicate it separately for EJS and CJS runs. Rows (top to bottom): total solid disc mass, planetesimal mass, and total mass in embryos/planets as a function of time (rows, top to bottom). Columns (left to right): masses for runs without Jupiter and Saturn (NJS), runs with Jupiter and Saturn on eccentric orbits (EJS), and runs with Jupiter and Saturn on circular orbits (CJS). For EJS and CJS runs, we show results for three different initial planetesimal discs. Figure 3. Open in new tabDownload slide Remaining solid mass (total, planetesimals, embryos) over time as a fraction of the total initial mass. We show the mass as a fraction of the initial total solid disc mass. In general, we show the median across 12 runs. For runs without Jupiter and Saturn, we also show the 10/90 percentile spread about the median. The spread is similar for all runs, so we do not indicate it separately for EJS and CJS runs. Rows (top to bottom): total solid disc mass, planetesimal mass, and total mass in embryos/planets as a function of time (rows, top to bottom). Columns (left to right): masses for runs without Jupiter and Saturn (NJS), runs with Jupiter and Saturn on eccentric orbits (EJS), and runs with Jupiter and Saturn on circular orbits (CJS). For EJS and CJS runs, we show results for three different initial planetesimal discs. Figure 4. Open in new tabDownload slide Fraction of initial solid disc mass that falls on to the star (heliocentric distance r < 0.2 au) or is ejected from the system (r > 20.0 au). We show the median over 12 runs for each set of simulations as the main bars. The small bars below and above indicate the 10/90 percentile spread of ejected and mass deposited on to the star, respectively. The 10/90 percentile ranges for NJS runs are unreliable due to file corruption. Figure 4. Open in new tabDownload slide Fraction of initial solid disc mass that falls on to the star (heliocentric distance r < 0.2 au) or is ejected from the system (r > 20.0 au). We show the median over 12 runs for each set of simulations as the main bars. The small bars below and above indicate the 10/90 percentile spread of ejected and mass deposited on to the star, respectively. The 10/90 percentile ranges for NJS runs are unreliable due to file corruption. Table 2. Total disc mass and total mass locked up in embryos. We show median mass, offsets from median to 25/75 percentile, and offsets from 25/75 to 10/90 percentile as sub/superscripts. Set . Time . MDisca,c . MEmbryosb,c . NJS 1 Myr |$4.87^{+0.04(+0.04)}_{-0.05(-0.09)}$| |$1.73^{+0.07(+0.04)}_{-0.09(-0.15)}$| End |$4.22^{+0.09(+0.04)}_{-0.15(-0.16)}$| |$3.30^{+0.08(+0.05)}_{-0.13(-0.14)}$| EJS 1 Myr |$4.91^{+0.07(+0.01)}_{-0.12(-0.05)}$| |$1.95^{+0.06(+0.03)}_{-0.12(-0.09)}$| End |$3.38^{+0.15(+0.06)}_{-0.27(-0.04)}$| |$3.38^{+0.13(+0.07)}_{-0.27(-0.04)}$| CJS 1 Myr |$4.92^{+0.06(+0.01)}_{-0.06(-0.02)}$| |$1.83^{+0.05(+0.08)}_{-0.07(-0.04)}$| End |$3.09^{+0.06(+0.26)}_{-0.15(-0.26)}$| |$2.99^{+0.08(+0.24)}_{-0.10(-0.28)}$| EJS/Steep 1 Myr |$4.26^{+0.08(+0.04)}_{-0.09(-0.04)}$| |$1.84^{+0.09(+0.07)}_{-0.08(-0.04)}$| End |$3.12^{+0.12(+0.10)}_{-0.21(-0.05)}$| |$3.11^{+0.12(+0.10)}_{-0.20(-0.05)}$| EJS/Heavy 1 Myr |$8.51^{+0.18(+0.18)}_{-0.19(-0.03)}$| |$3.24^{+0.26(+0.14)}_{-0.07(-0.30)}$| End |$4.15^{+0.35(+0.14)}_{-0.27(-0.33)}$| |$4.15^{+0.34(+0.15)}_{-0.28(-0.34)}$| CJS/Steep 1 Myr |$4.28^{+0.10(+0.06)}_{-0.09(-0.07)}$| |$1.70^{+0.13(+0.09)}_{-0.04(-0.07)}$| End |$2.83^{+0.13(+0.10)}_{-0.08(-0.11)}$| |$2.78^{+0.11(+0.11)}_{-0.14(-0.09)}$| CJS/Heavy 1 Myr |$8.95^{+0.09(+0.04)}_{-0.15(-0.07)}$| |$2.34^{+0.23(+0.16)}_{-0.09(-0.02)}$| End |$4.57^{+0.15(+0.23)}_{-0.24(-0.35)}$| |$4.54^{+0.16(+0.20)}_{-0.26(-0.38)}$| Set . Time . MDisca,c . MEmbryosb,c . NJS 1 Myr |$4.87^{+0.04(+0.04)}_{-0.05(-0.09)}$| |$1.73^{+0.07(+0.04)}_{-0.09(-0.15)}$| End |$4.22^{+0.09(+0.04)}_{-0.15(-0.16)}$| |$3.30^{+0.08(+0.05)}_{-0.13(-0.14)}$| EJS 1 Myr |$4.91^{+0.07(+0.01)}_{-0.12(-0.05)}$| |$1.95^{+0.06(+0.03)}_{-0.12(-0.09)}$| End |$3.38^{+0.15(+0.06)}_{-0.27(-0.04)}$| |$3.38^{+0.13(+0.07)}_{-0.27(-0.04)}$| CJS 1 Myr |$4.92^{+0.06(+0.01)}_{-0.06(-0.02)}$| |$1.83^{+0.05(+0.08)}_{-0.07(-0.04)}$| End |$3.09^{+0.06(+0.26)}_{-0.15(-0.26)}$| |$2.99^{+0.08(+0.24)}_{-0.10(-0.28)}$| EJS/Steep 1 Myr |$4.26^{+0.08(+0.04)}_{-0.09(-0.04)}$| |$1.84^{+0.09(+0.07)}_{-0.08(-0.04)}$| End |$3.12^{+0.12(+0.10)}_{-0.21(-0.05)}$| |$3.11^{+0.12(+0.10)}_{-0.20(-0.05)}$| EJS/Heavy 1 Myr |$8.51^{+0.18(+0.18)}_{-0.19(-0.03)}$| |$3.24^{+0.26(+0.14)}_{-0.07(-0.30)}$| End |$4.15^{+0.35(+0.14)}_{-0.27(-0.33)}$| |$4.15^{+0.34(+0.15)}_{-0.28(-0.34)}$| CJS/Steep 1 Myr |$4.28^{+0.10(+0.06)}_{-0.09(-0.07)}$| |$1.70^{+0.13(+0.09)}_{-0.04(-0.07)}$| End |$2.83^{+0.13(+0.10)}_{-0.08(-0.11)}$| |$2.78^{+0.11(+0.11)}_{-0.14(-0.09)}$| CJS/Heavy 1 Myr |$8.95^{+0.09(+0.04)}_{-0.15(-0.07)}$| |$2.34^{+0.23(+0.16)}_{-0.09(-0.02)}$| End |$4.57^{+0.15(+0.23)}_{-0.24(-0.35)}$| |$4.54^{+0.16(+0.20)}_{-0.26(-0.38)}$| aTotal mass in disc. bTotal mass in embryos. cIn Earth masses. Open in new tab Table 2. Total disc mass and total mass locked up in embryos. We show median mass, offsets from median to 25/75 percentile, and offsets from 25/75 to 10/90 percentile as sub/superscripts. Set . Time . MDisca,c . MEmbryosb,c . NJS 1 Myr |$4.87^{+0.04(+0.04)}_{-0.05(-0.09)}$| |$1.73^{+0.07(+0.04)}_{-0.09(-0.15)}$| End |$4.22^{+0.09(+0.04)}_{-0.15(-0.16)}$| |$3.30^{+0.08(+0.05)}_{-0.13(-0.14)}$| EJS 1 Myr |$4.91^{+0.07(+0.01)}_{-0.12(-0.05)}$| |$1.95^{+0.06(+0.03)}_{-0.12(-0.09)}$| End |$3.38^{+0.15(+0.06)}_{-0.27(-0.04)}$| |$3.38^{+0.13(+0.07)}_{-0.27(-0.04)}$| CJS 1 Myr |$4.92^{+0.06(+0.01)}_{-0.06(-0.02)}$| |$1.83^{+0.05(+0.08)}_{-0.07(-0.04)}$| End |$3.09^{+0.06(+0.26)}_{-0.15(-0.26)}$| |$2.99^{+0.08(+0.24)}_{-0.10(-0.28)}$| EJS/Steep 1 Myr |$4.26^{+0.08(+0.04)}_{-0.09(-0.04)}$| |$1.84^{+0.09(+0.07)}_{-0.08(-0.04)}$| End |$3.12^{+0.12(+0.10)}_{-0.21(-0.05)}$| |$3.11^{+0.12(+0.10)}_{-0.20(-0.05)}$| EJS/Heavy 1 Myr |$8.51^{+0.18(+0.18)}_{-0.19(-0.03)}$| |$3.24^{+0.26(+0.14)}_{-0.07(-0.30)}$| End |$4.15^{+0.35(+0.14)}_{-0.27(-0.33)}$| |$4.15^{+0.34(+0.15)}_{-0.28(-0.34)}$| CJS/Steep 1 Myr |$4.28^{+0.10(+0.06)}_{-0.09(-0.07)}$| |$1.70^{+0.13(+0.09)}_{-0.04(-0.07)}$| End |$2.83^{+0.13(+0.10)}_{-0.08(-0.11)}$| |$2.78^{+0.11(+0.11)}_{-0.14(-0.09)}$| CJS/Heavy 1 Myr |$8.95^{+0.09(+0.04)}_{-0.15(-0.07)}$| |$2.34^{+0.23(+0.16)}_{-0.09(-0.02)}$| End |$4.57^{+0.15(+0.23)}_{-0.24(-0.35)}$| |$4.54^{+0.16(+0.20)}_{-0.26(-0.38)}$| Set . Time . MDisca,c . MEmbryosb,c . NJS 1 Myr |$4.87^{+0.04(+0.04)}_{-0.05(-0.09)}$| |$1.73^{+0.07(+0.04)}_{-0.09(-0.15)}$| End |$4.22^{+0.09(+0.04)}_{-0.15(-0.16)}$| |$3.30^{+0.08(+0.05)}_{-0.13(-0.14)}$| EJS 1 Myr |$4.91^{+0.07(+0.01)}_{-0.12(-0.05)}$| |$1.95^{+0.06(+0.03)}_{-0.12(-0.09)}$| End |$3.38^{+0.15(+0.06)}_{-0.27(-0.04)}$| |$3.38^{+0.13(+0.07)}_{-0.27(-0.04)}$| CJS 1 Myr |$4.92^{+0.06(+0.01)}_{-0.06(-0.02)}$| |$1.83^{+0.05(+0.08)}_{-0.07(-0.04)}$| End |$3.09^{+0.06(+0.26)}_{-0.15(-0.26)}$| |$2.99^{+0.08(+0.24)}_{-0.10(-0.28)}$| EJS/Steep 1 Myr |$4.26^{+0.08(+0.04)}_{-0.09(-0.04)}$| |$1.84^{+0.09(+0.07)}_{-0.08(-0.04)}$| End |$3.12^{+0.12(+0.10)}_{-0.21(-0.05)}$| |$3.11^{+0.12(+0.10)}_{-0.20(-0.05)}$| EJS/Heavy 1 Myr |$8.51^{+0.18(+0.18)}_{-0.19(-0.03)}$| |$3.24^{+0.26(+0.14)}_{-0.07(-0.30)}$| End |$4.15^{+0.35(+0.14)}_{-0.27(-0.33)}$| |$4.15^{+0.34(+0.15)}_{-0.28(-0.34)}$| CJS/Steep 1 Myr |$4.28^{+0.10(+0.06)}_{-0.09(-0.07)}$| |$1.70^{+0.13(+0.09)}_{-0.04(-0.07)}$| End |$2.83^{+0.13(+0.10)}_{-0.08(-0.11)}$| |$2.78^{+0.11(+0.11)}_{-0.14(-0.09)}$| CJS/Heavy 1 Myr |$8.95^{+0.09(+0.04)}_{-0.15(-0.07)}$| |$2.34^{+0.23(+0.16)}_{-0.09(-0.02)}$| End |$4.57^{+0.15(+0.23)}_{-0.24(-0.35)}$| |$4.54^{+0.16(+0.20)}_{-0.26(-0.38)}$| aTotal mass in disc. bTotal mass in embryos. cIn Earth masses. Open in new tab Table 3. Percentile offsets from median: (i) total mass, (ii) total mass in embryos and (iii) total mass in planetesimals. The offsets are averaged (50 percentile) in time. Set . Percentile . MDisca,d . MEmbryosb,d . MPlanetesimalsc,d,e . xxh NJS 90 +0.13 +0.12 +0.08 75 +0.10 +0.10 +0.03 25 −0.15 −0.15 −0.04 10 −0.31 −0.27 −0.06 EJS 90 +0.27 +0.24 +0.18 75 +0.15 +0.13 +0.06 25 −0.12 −0.22 −0.05 10 −0.32 −0.32 −0.08 CJS 90 +0.22 +0.27 +0.15 75 +0.10 +0.08 +0.11 25 −0.13 −0.08 −0.07 10 −0.27 −0.22 −0.15 EJS/Steep 90 +0.12 +0.13 +0.13 75 +0.09 +0.09 +0.07 25 −0.19 −0.19 −0.03 10 −0.30 −0.32 −0.04 EJS/Heavy 90 +0.51 +0.50 +0.11 75 +0.35 +0.34 +0.05 25 −0.26 −0.26 −0.07 10 −0.60 −0.63 −0.14 CJS/Steep 90 +0.20 +0.16 +0.10 75 +0.14 +0.08 +0.04 25 −0.10 −0.15 −0.06 10 −0.20 −0.24 −0.12 CJS/Heavy 90 +0.34 +0.33 +0.37 75 +0.13 +0.23 +0.23 25 −0.27 −0.25 −0.23 10 −0.61 −0.53 −0.35 Set . Percentile . MDisca,d . MEmbryosb,d . MPlanetesimalsc,d,e . xxh NJS 90 +0.13 +0.12 +0.08 75 +0.10 +0.10 +0.03 25 −0.15 −0.15 −0.04 10 −0.31 −0.27 −0.06 EJS 90 +0.27 +0.24 +0.18 75 +0.15 +0.13 +0.06 25 −0.12 −0.22 −0.05 10 −0.32 −0.32 −0.08 CJS 90 +0.22 +0.27 +0.15 75 +0.10 +0.08 +0.11 25 −0.13 −0.08 −0.07 10 −0.27 −0.22 −0.15 EJS/Steep 90 +0.12 +0.13 +0.13 75 +0.09 +0.09 +0.07 25 −0.19 −0.19 −0.03 10 −0.30 −0.32 −0.04 EJS/Heavy 90 +0.51 +0.50 +0.11 75 +0.35 +0.34 +0.05 25 −0.26 −0.26 −0.07 10 −0.60 −0.63 −0.14 CJS/Steep 90 +0.20 +0.16 +0.10 75 +0.14 +0.08 +0.04 25 −0.10 −0.15 −0.06 10 −0.20 −0.24 −0.12 CJS/Heavy 90 +0.34 +0.33 +0.37 75 +0.13 +0.23 +0.23 25 −0.27 −0.25 −0.23 10 −0.61 −0.53 −0.35 aTotal mass in disc. bTotal mass in embryos. cTotal mass in planetesimals. dIn Earth masses. eOnly considers times <32.9 Myr. Open in new tab Table 3. Percentile offsets from median: (i) total mass, (ii) total mass in embryos and (iii) total mass in planetesimals. The offsets are averaged (50 percentile) in time. Set . Percentile . MDisca,d . MEmbryosb,d . MPlanetesimalsc,d,e . xxh NJS 90 +0.13 +0.12 +0.08 75 +0.10 +0.10 +0.03 25 −0.15 −0.15 −0.04 10 −0.31 −0.27 −0.06 EJS 90 +0.27 +0.24 +0.18 75 +0.15 +0.13 +0.06 25 −0.12 −0.22 −0.05 10 −0.32 −0.32 −0.08 CJS 90 +0.22 +0.27 +0.15 75 +0.10 +0.08 +0.11 25 −0.13 −0.08 −0.07 10 −0.27 −0.22 −0.15 EJS/Steep 90 +0.12 +0.13 +0.13 75 +0.09 +0.09 +0.07 25 −0.19 −0.19 −0.03 10 −0.30 −0.32 −0.04 EJS/Heavy 90 +0.51 +0.50 +0.11 75 +0.35 +0.34 +0.05 25 −0.26 −0.26 −0.07 10 −0.60 −0.63 −0.14 CJS/Steep 90 +0.20 +0.16 +0.10 75 +0.14 +0.08 +0.04 25 −0.10 −0.15 −0.06 10 −0.20 −0.24 −0.12 CJS/Heavy 90 +0.34 +0.33 +0.37 75 +0.13 +0.23 +0.23 25 −0.27 −0.25 −0.23 10 −0.61 −0.53 −0.35 Set . Percentile . MDisca,d . MEmbryosb,d . MPlanetesimalsc,d,e . xxh NJS 90 +0.13 +0.12 +0.08 75 +0.10 +0.10 +0.03 25 −0.15 −0.15 −0.04 10 −0.31 −0.27 −0.06 EJS 90 +0.27 +0.24 +0.18 75 +0.15 +0.13 +0.06 25 −0.12 −0.22 −0.05 10 −0.32 −0.32 −0.08 CJS 90 +0.22 +0.27 +0.15 75 +0.10 +0.08 +0.11 25 −0.13 −0.08 −0.07 10 −0.27 −0.22 −0.15 EJS/Steep 90 +0.12 +0.13 +0.13 75 +0.09 +0.09 +0.07 25 −0.19 −0.19 −0.03 10 −0.30 −0.32 −0.04 EJS/Heavy 90 +0.51 +0.50 +0.11 75 +0.35 +0.34 +0.05 25 −0.26 −0.26 −0.07 10 −0.60 −0.63 −0.14 CJS/Steep 90 +0.20 +0.16 +0.10 75 +0.14 +0.08 +0.04 25 −0.10 −0.15 −0.06 10 −0.20 −0.24 −0.12 CJS/Heavy 90 +0.34 +0.33 +0.37 75 +0.13 +0.23 +0.23 25 −0.27 −0.25 −0.23 10 −0.61 −0.53 −0.35 aTotal mass in disc. bTotal mass in embryos. cTotal mass in planetesimals. dIn Earth masses. eOnly considers times <32.9 Myr. Open in new tab Initially, all simulation sets evolve identically. After 105 years, the first planetary embryos become visible. It takes about four equal-mass collisions to reach the cutoff mass, so the observed early-time collision rate of ≈0.03 collisions per year per 2000 particles agrees with this build-up. After 1 Myr, all simulations have assembled between 1.5 and 2.0 MEarth into embryos, but have overall lost very little of their total solid disc mass. Less than half the systems have lost more than 0.1 MEarth, although isolated systems in configurations without giant planets have lost as much as 0.55 MEarth. By the 3 Myr mark, total mass loss across all runs has increased to between 0.54 and 1.14 MEarth. The mass loss is driven by hydrodynamic drag and type I migration, which by now have delivered particles to the inner edge of the disc, where we remove them from the simulation. In fact, the first particles already fall in around 2 × 105 years, but they do not carry significant mass. After 3 Myr, the gas disc is depleted by a factor of 20 and becomes dynamically irrelevant for migration (Tanaka, Takeuchi & Ward 2002) and secular resonance sweeping (Nagasawa et al. 2000), although it remains important for damping the eccentricities and inclinations of embryos (Kominami & Ida 2002; Tanaka & Ward 2004).8 At this stage, all simulations host between 5 and 18 oligarchic embryos at semi-major axes a ≲ 2 au. During the first 3 Myr, the differences between sets are small, largely because gas drag acts as an equalizer that dampens eccentricity excitements from resonant interaction between planetesimals and Jupiter. At later times, we observe three marked differences. First, simulations with giant planets are much more efficient at assembling embryos. At ∼50 Myr, EJS simulations have converted all low-mass particles. CJS simulations are slower, and retain ∼0.08 MEarth in the low-mass regime when the simulations terminate. The process is even slower in simulations without giant planets, which retain ∼1 MEarth in low-mass objects at termination. Secondly, simulations with giant planets continue to lose mass after 3 Myr, eventually bringing the remaining disc mass down to ∼3.38 MEarth (EJS) and ∼3.09 MEarth (CJS). Thirdly, there is a significant difference in how EJS and CJS simulations lose disc mass. For EJS, all lost mass is lost on to the host star, while CJS simulations eject about half the lost mass from the system. However, these processes are limited to late times. The spread of tracked mass ranges is never uncomfortably large, such that the 10 and 90 percentiles are never more than 15 per cent off the median. The exception to this is the total planetesimal mass, which suffers from small number statistics as planetesimals deplete at times ≳100 Myr. Overall, the tight spread suggests that the range of evolutionary paths available to individual runs in a set is limited. However, we do observe three key trends in statistical spread with simulation sets. First, the spread is consistently smallest in masses below the cutoff mass. Secondly, the spread across runs in total disc mass is larger in simulations that include giant planets. Thirdly, the probability of larger excursions is higher in the EJS and CJS simulations. For EJS sets, the 10/90 percentiles tend to be a factor of 2 further off the median than the 25/75 percentiles. In CJS sets, the difference grows and exceeds factors of 3 for the 90 percentile in mass above the cutoff. At small semi-major axis, relative velocities between planetesimals are higher, leading to shorter time-scales of collision and thus growth. Runs with steeper initial surface density profiles therefore start assembling embryos already after ∼2 × 104 years. As a consequence, more embryos are driven close to the star by type I migration, where they are removed from the simulation. This effectively stalls the conversion of planetesimals into embryos between 5 × 105 and 106 years. Although the conversion process picks up again, initially steep discs still tend to host ∼8 per cent less total mass in terrestrial planets than the reference discs. Apart from the short stall, the fraction of mass in embryos grows continuously until ∼10 Myr, whereafter it slows down significantly. This is irrespective of the orbit of the giant planets, although conversion of mass into embryos proceeds slightly faster if giant planets are on eccentric orbits. In massive discs, conversion of planetesimals into embryos begins earlier than in the reference run (but later than in runs with steep initial surface density profiles). This is likely a numerical artefact related to the doubling of the initial planetesimal mass.9 More robustly, we find that the fraction of mass in embryos either (i) stops increasing significantly after ∼1 Myr (eccentric giant planets) or (ii) proceeds slower overall (circular giant planets) for initially massive discs. However, in absolute terms, it does settle at a higher level (6 MEarth versus ∼3 MEarth for less massive discs). Despite the slower rate of embryo assembly and smaller total mass fraction of embryos, both steep and massive discs lose total mass earlier and more efficiently than the reference profiles. For the reference discs, mass loss only sets in after ∼1 Myr; almost 0.5 Myr later than when steep and massive discs begin losing material. Most extremely, systems hosting massive discs with eccentric giant planets begin losing mass almost immediately at the 2 × 104 yr mark. These systems must be dynamically more active than the reference profiles, such that more orbits can deliver mass on to the star, out of the system, or on to the giants. In fact, no single planetesimal or embryo collides with Saturn, and we record only three to four planetesimals (corresponding to 0.01 to 0.05 MEarth) impacting on Jupiter per run, but only in runs with massive planetesimal discs. There are occasional collisions with Jupiter in other runs, but these are isolated events and do not occur in all runs of a given configuration. This is not particularly surprising as Jupiter and Saturn have larger escape velocities than the host star (evaluated at their location). In other words, particles from the inner systems that have been thrown on to orbits crossing Jupiter or Saturn tend to be too fast to be caught by either (Ford & Rasio 2008; Raymond, Armitage & Gorelick 2010). Most of the mass loss is thus on to the star or through ejection from the system and the mechanics vary depending on the configuration (see Fig. 4). In systems without giant planets and with giant planets on eccentric orbits, almost all mass is lost on to the star, although massive discs in EJS configurations eject ∼10 per cent of the lost mass. For configurations of giants on circular orbits, between 20 and 50 per cent of the lost solid mass is ejected. For both giant planet configurations, initially massive discs appear to eject the largest fraction of their initial mass. Similarly to the reference runs, the spread of the masses remains well-bound. Except for outliers from small number statistics in the planetesimal mass at late times, the 10 and 90 percentiles of the masses remain within 15 per cent off the median. For eccentric giant planets, steep initial discs have a smaller and massive initial discs a larger spread than the reference case. For giants on circular orbits, steep profiles have comparable spreads to the initial profile. Massive planetesimal discs in CJS configurations show no definite trends either. For example, the 90 percentile mass in embryos in massive discs has a smaller offset from the median than the reference profile (7.2 versus 9.2 per cent) while the 75 per cent has a larger offset (5.17 versus 2.70 per cent). 3.3 Mass distributions and connection to initial conditions In terms of aggregate mass bound up in planetesimals, embryos, and the solid disc as a whole, runs within sets of initially identical conditions have variance below the 15 per cent level. But what about the radial mass distributions and mass–frequency distributions at the end of the simulations? Do the systems look architecturally similar? Is the mass distributed similarly? What is link between the initial and final mass distributions? Fig. 5 illustrates the diverse architectures of systems we have generated after 147.84 Myr. Considering the 12 runs within sets, we find the arrangement of final planets to vary substantially. Nevertheless, adding giant planets or adjusting their configuration has an even stronger impact that can easily be discerned visually over the different architectures generated by the same configuration. On the other hand, runs with different initial surface density profiles appear difficult to visually discern from the stochastic variations in architecture. For massive discs, this is easier because they tend to host massive terrestrial planets. Figure 5. Open in new tabDownload slide Distribution of planets/planetesimals with semi-major axis for all 7 × 12 initially identical runs after 147.84 Myr. Each line indicates a separate run. Darker colours and bigger circles indicate more massive particles. Columns (left to right): runs without Jupiter and Saturn, Jupiter and Saturn on their present-day eccentric orbits, Jupiter and Saturn on circular orbits. Rows (top to bottom): reference initial planetesimal disc, steep disc, more massive initial disc. Figure 5. Open in new tabDownload slide Distribution of planets/planetesimals with semi-major axis for all 7 × 12 initially identical runs after 147.84 Myr. Each line indicates a separate run. Darker colours and bigger circles indicate more massive particles. Columns (left to right): runs without Jupiter and Saturn, Jupiter and Saturn on their present-day eccentric orbits, Jupiter and Saturn on circular orbits. Rows (top to bottom): reference initial planetesimal disc, steep disc, more massive initial disc. To be more quantitative, we analyse the radial mass distributions (mass per semi-major axis) and mass–frequency distributions (particle counts per mass bin; also called the mass function in cosmology and galaxy formation) for various giant planet configurations and initial planetesimal distributions. In both cases, we use non-parametric kernel density estimates (KDE) to derive the distribution functions for all runs in a set. They are combined statistically in Figs 6 and 7. We also compute the final surface density profiles by grinding up the terrestrial planets, distributing their mass over annuli of suitable widths, and fitting power laws.10 In Fig. 8, we compare the range of power-law exponents for the planetary systems to the initial distribution of planetesimals. Figure 6. Open in new tabDownload slide Median (lines) and statistical spreads (shading) of the radial mass distribution each set of 7 × 12 runs shown in Fig. 5. The distributions are obtained as KDEs of the mass-weighted particle distribution dM/da along the semi-major axis a. Black lines indicate the median. The shadings indicate 25/75 and 10/90 percentile ranges. Triangles indicate (from left to right) the locations of the ν5, ν6, and ν16 secular resonances for runs with giant planets. Columns (left to right): runs without Jupiter and Saturn, Jupiter and Saturn on their present-day eccentric orbits, Jupiter and Saturn on circular orbits. Rows (top to bottom): reference initial planetesimal disc, steep disc, massive initial disc. Figure 6. Open in new tabDownload slide Median (lines) and statistical spreads (shading) of the radial mass distribution each set of 7 × 12 runs shown in Fig. 5. The distributions are obtained as KDEs of the mass-weighted particle distribution dM/da along the semi-major axis a. Black lines indicate the median. The shadings indicate 25/75 and 10/90 percentile ranges. Triangles indicate (from left to right) the locations of the ν5, ν6, and ν16 secular resonances for runs with giant planets. Columns (left to right): runs without Jupiter and Saturn, Jupiter and Saturn on their present-day eccentric orbits, Jupiter and Saturn on circular orbits. Rows (top to bottom): reference initial planetesimal disc, steep disc, massive initial disc. Figure 7. Open in new tabDownload slide Median (lines) and statistical spreads (shading) of the mass–frequency distribution (i.e. the number of particles per mass bin – dN/dM) for our simulation sets. The distributions are computed as KDEs. Columns (left to right): runs without Jupiter and Saturn, Jupiter and Saturn on their present-day eccentric orbits, Jupiter and Saturn on circular orbits. Rows (top to bottom): reference initial planetesimal disc, steep disc, more massive initial disc. Figure 7. Open in new tabDownload slide Median (lines) and statistical spreads (shading) of the mass–frequency distribution (i.e. the number of particles per mass bin – dN/dM) for our simulation sets. The distributions are computed as KDEs. Columns (left to right): runs without Jupiter and Saturn, Jupiter and Saturn on their present-day eccentric orbits, Jupiter and Saturn on circular orbits. Rows (top to bottom): reference initial planetesimal disc, steep disc, more massive initial disc. Figure 8. Open in new tabDownload slide Power-law slopes (β, where ΣPlanets ∝ aβ) for the projected surface densities ΣPlanets at semi-major axis a obtained by grinding up the final terrestrial planets. Smaller values of β indicate steeper profiles. Boxes indicate the 25/75 percentile range, their centreline the median, whiskers extend to 1.5 times the interquartile range, and crosses indicate outliers. We also show the slope of the initial planetesimal distribution with bold vertical lines. Figure 8. Open in new tabDownload slide Power-law slopes (β, where ΣPlanets ∝ aβ) for the projected surface densities ΣPlanets at semi-major axis a obtained by grinding up the final terrestrial planets. Smaller values of β indicate steeper profiles. Boxes indicate the 25/75 percentile range, their centreline the median, whiskers extend to 1.5 times the interquartile range, and crosses indicate outliers. We also show the slope of the initial planetesimal distribution with bold vertical lines. 3.3.1 System architectures and mass distributions Simulations without giant planets distribute mass most evenly with semi-major axis, while CJS and EJS simulations concentrate progressively more mass at smaller semi-major axes. The EJS systems tend to be truncated around 2 au. In CJS systems, terrestrial planets populate a region out to 3.5 au, but most mass is concentrated within 2.2 au. In simulations with giant planets, there appear to be regions that are preferably populated by terrestrial planets. Most notably, these are at 0.8, 1.2 (EJS) as well as at 0.9, 1.5, and 2.2 au (CJS), although the last peak is arguably weak. However, we observe significant variations in the distribution of mass and planets with semi-major axis within sets. These are smallest in systems without giant planets and largest in systems with giants on circular orbits. For example, the 90 percentile masses are (by median) 25, 50, and 69 per cent off the median mass (in order, NJS, EJS, CJS). The lack of giant planets severely stunts growth of terrestrial planets, and these systems retain larger numbers of planetesimals at their origin mass. In contrast, EJS configurations are the most efficient at assembling terrestrial planets with masses of the order of an Earth mass, while CJS runs populate a middle ground. Irrespective of presence and configuration of giant planets, there is significant variance in the mass–frequency distribution. For NJS runs, the scatter is most significant at intermediate mass ranges 0.03 to 0.3 MEarth, with the 25/75 percentiles off by a factor of about 2 from the median. The same holds for CJS configurations. For EJS configurations, the mass–frequency distributions are more difficult to interpret because almost no planetesimals remain. Focusing on planetary masses (0.5 to 1.5 au), we find a larger spread across EJS runs than for CJS runs (10 to 20 per cent versus factors of 2). Initializing the planetesimals to follow a steeper surface density profile has a surprisingly small effect. For the CJS case, we find no obvious discernible trends, although we do note two changes for EJS configurations. First, more mass is concentrated at smaller semi-major axis, which is a consequence of more mass being initially distributed here. This is only significant in EJS configurations because the semi-major axis range populated by terrestrial planets is tighter here. Secondly, the mass range covered by planets widens a bit and we find that more EJS runs retain some planetesimals in initially steep discs. This still holds for less than half of the runs, which is still more than in the reference profile. Increasing the initial disc mass has a more drastic effect on the distributions. Across the entire semi-major axis range populated by the terrestrial planets, both CJS and EJS configurations record a higher median mass as well as a larger scatter. Massive initial discs also generate a wider range of planet masses (at both the lower and upper mass end), which essentially flattens out the mass–frequency distribution, although the dip around 0.1 MEarth in CJS configurations remains and exhibits a large scatter. For EJS runs, we find more (but still not all) simulations to retain planetesimals, which puts power into the lower end of the mass–frequency distribution. Nevertheless, large spreads remain here as the number of surviving planetesimals per run is still in the single digits. Before closing, we wish to comment on the robustness of using KDEs to characterize radial mass and mass–frequency distributions. Fitting a KDE amounts to attempting to estimate the underlying distribution of planetesimals that is sampled by the simulation particles. For each sample (particle), the algorithm drops a Gaussian kernel of a width depending on the spacing of samples. The KDE is then sum of kernels. For NJS runs, this works well because the distribution per run is fit to between 150 and 200 samples. For CJS and EJS configurations, the number of available samples drops to between 10 and 20 (CJS), respectively 5 to 10 (EJS). The algorithm compensates by dropping much wider Gaussian kernels. Stacking runs and computing percentiles over these kernels lead to wide regions with jagged edges that are especially apparent in the mass–frequency distributions of EJS and EJS/Steep runs. We thus caution from attributing too much meaning to these regions because they are very sparsely sampled. 3.3.2 A memory of initial conditions? Frequently, attempts at reproducing the formation history of both the Solar system and exoplanetary systems are based around the concept of a Minimum Mass (Extra-)Solar Nebula. Here, the presently observed distribution of planets is ground up, distributed over their assumed feeding zones, and used to reconstruct a radial profile of the solid components that feed into the planets. By initializing simulations to follow such a profile, it is implicitly assumed that such a profile persists as planets form. For reference, Weidenschilling (1977) and Hayashi (1981) determine the Minimum Mass Solar Nebular to follow a surface density profile of ΣPlanets ∝ r β with β ∼ −1.5. For extrasolar systems, Kuchner (2004) and Chiang & Laughlin (2013) report β ∼ −1.6 to β ∼ −2.4. More recently, Raymond & Cossou (2014) report an even wider range of profiles that can be fit through observed systems (β ∼ −3.2 to β ∼ 0.5) and argue this indicates that there is no (or no tight range of) underlying initial profiles. We now show the ranges of final surface density profiles that can be obtained from a given initial distribution to determine if the initial profile indeed persists. We will find that – even when starting with identical ICs – there is a large scatter of generated profiles, such that the naive mapping of an initial to a final profile does not make sense. We now compute how the surface density of our terrestrial planets ΣPlanets varies with semi-major axis a at the end of each of our runs. We grind up each planet k over an annulus with boundaries determined by the geometric mean to the next planet on either side (Weidenschilling 1977; Kuchner 2004) to recover ΣPlanet,k and then fit a power law of the form ΣPlanets ∝ a β. We perform the least-squares fit in log-space and consider a coefficient of determination r2 > 0.6 to indicate a sufficient enough fit to allow such a power-law description. Fig. 8 summarizes the statistical spread of the obtained power-law exponents. In each of simulation set, at most three runs have r2 ≤ 0.6, which we exclude from the analysis. Note that we do not include NJS runs because they are dominated by a large number of small embryos or large planetesimals where the concept of ‘grinding up planets’ does not make sense. In ∼85 per cent of the final systems, the surface density profiles are steeper than the initial profiles, with rare occurrences of a shallower final profile existing in all initial planetesimal distributions and giant planet configurations. For the EJS runs, the final profiles are always steeper than the initial ones. The median steepening is largest for EJS/Heavy and CJS (142 and 94 per cent, respectively) and smallest for EJS/Steep and CJS/Steep runs (32 and 29 per cent, respectively). All in all, the results are sobering. Although a power law is good fit to over 75 per cent of the final distributions, there is a huge scatter in fitted exponents across runs within a set. Even in the most robust case (EJS), the difference between the shallowest and steepest profile exponents is about 1 with a 25 per cent scatter between the 10/90 percentiles and the median. For all other configurations, the average spread exceeds 40 per cent, but we note that the 10 percentiles tend to be within 30 per cent of the median, while the 90 percentile is usually 60 to 90 per cent off the median. The largest absolute difference between 10 and 90 percentile exponents is recorded for EJS/Heavy with a difference of 2.23. Clearly, generating a simple mapping from initial to a final distribution is out of the question. Interestingly enough, when Raymond et al. (2005b) characterized their difference between initial and final surface density profiles, they found little change in the slope of the final power law, except for profiles that had substantial amount of material in the vicinity a giant planet. However, they started from a distribution of embryos which did not grow from planetesimals in a gaseous disc. As such, their planets were relegated to forming locally based on the initial distribution. We therefore conclude that mapping initial to final profiles is only possible if the initial profile is considered at a stage when embryos have already formed, gas has dissipated, and there is no immediate influence of giant planets. Conversely, is it not possible to establish a link to the initial distribution of planetesimals in a gas disc by observing present-day systems. 3.4 Stability and orbital spacing On time-scales of a few 100 Myr, stability for two orbiting planets against close encounters requires their orbital separation Δa to exceed Δa > 10 RHill, where RHill = ((m1 + m2)/(3MSun))1/3 (a1 + a2)/2 is the mutual Hill radius of two neighbouring planets (Gladman 1993; Chambers et al. 1996). The threshold increases to Δa ≳ 20 for systems with 10 to 20 planets. By this criterion, we find about 90 per cent of all terrestrial planets to be on stable orbits although this number fluctuates with giant planet configuration and initial planetesimal disc profile. The configurations EJS/Steep and EJS/Heavy have no planet pairs with spacings Δa < 20 RHill. We therefore expect none of these runs to eject (or collide) any planets with a few hundred Myr after our simulations terminate. Of the remaining configurations, the CJS/Heavy simulations have the fewest pairs spaced Δa < 20 RHill, so that we expect at most one planet for every four runs to collide or be ejected. The systems in EJS, CJS, and CJS/Steep even admit a few extremely tightly planets with Δa < 3 RHill. Over 12 runs, we expect one planet to be ejected per every two to three runs. Finally, NJS configurations without giant planets host a large number of such tight separations. Here, we expect at least one to two planets per run to collide or be ejected within a few hundred Myr after our runs terminate. Yoshinaga et al. (1999) improve on this stability analysis by considering the influence of eccentric (and inclined) orbits. They use the scaled eccentricity |$\skew4\tilde{e} = e/h$| (with h = RHill/a) in systems of 10 protoplanets with mutual spacing 4 ≤ RHill ≤ 10 and |$0 \le \tilde{e} \le 4$| to find that eccentric orbits shorten the time it takes a set of protoplanets to become dynamically unstable by up to two orders of magnitude. Unfortunately, it is unclear to what degree their analysis can be applied to our simulations as they cover an entirely different dynamic regime. Except for NJS runs, we find systems that host four to eight planets with spacings 20 ≲ Δa/RHill ≲ 55 and eccentricities |$3 \lesssim \tilde{e} \lesssim 30$| (bounds are 10/90 percentiles). Attempting a simple extrapolation of their fits, we find a required spacing of 40 RHill for planets excited to |$\tilde{e} \sim 8$| (the median for our runs) to remain stable for a hundred Myr. For runs with giant planets, this is fulfilled for over 90 per cent of all cases which suggests stability to be in line with the analysis based on Chambers et al. (1996). For NJS runs, only 25 per cent of planets have such a large spacing, again emphasizing that we expect significant dynamical evolution in the systems after our simulations terminate. At the 90 percentile level, our simulations exhibit values between |$\tilde{e} \sim 14$| and |$\tilde{e} \sim 36$|⁠, which requires separations of ∼80 RHill and ∼160 RHill, which is not fulfilled by all our planets at such eccentricities. Overall, we conclude that the Yoshinaga et al. (1999) criteria are stricter than those of Chambers et al. (1996) and that their direct application would suggest a larger fraction of planets to not survive the next few hundred Myr in our simulations. However, the straightforward application of their criterion is questionable as our systems populate a different dynamical regime. It remains now to characterize what sets the spacing of terrestrial planets. In particular, the double peaked structures of the radial mass distribution for some of our simulations (EJS, CJS, and EJS/Steep) suggest preferred orbital spacings, which may be indicative of chains of mean motion resonances. However, out of a total of 35 554 planet pairs, we find that only 13 (0.036 per cent) are no further than 0.001 au away from the semi-major axis ratio required for any of the mean motion resonances considered (3:1, 5:2, 7:3, 2:1, and 3:2). Even if we are willing to relax the distance to the resonance to 0.01 au, only 117 (0.33 per cent) of planet pairs are in resonance. We conclude that resonances between the formed terrestrial planet are rare and do not drive the double peaked radial mass distribution seen in some configurations. 3.5 Final orbital parameters Numerically speaking, all particle distributions that are evolved from the same ICs are equally valid configurations of the planetary system. We thus argue that (at present) we have no means of determining the actual configuration of a simulated planetary system. This leaves us with two choices of analysis. We could (i) try to determine the most likely actual configuration or (ii) explore the range of permitted configurations. If we do the former, we may be tempted to declare the median of all configurations to be the most likely true configuration. This is by no means well-justified, and would require a large number of simulation runs to escape the confines of small number statistics. Given this complication, we settle for a survey of possible configurations. We determine the range of the valid configurations at the end of the simulations by stacking data from all runs in a given set. We restrict our analysis to particles M > MCut, but make no distinction between planets and planetary embryos. Fig. 9 shows the raw data as well as median and 10/90 percentile ranges of the semi-major axis a, mass M, orbital eccentricity e, and orbital inclination i. Results are also tabulated in Table 4 where we list the number of planets per system. We include the Solar system in Table 4 purely for reference and stress that none of our simulations was set up with the explicit intent of reproducing the Solar system. Figure 9. Open in new tabDownload slide Final distribution of (rows, from top to bottom rows) mass, eccentricity, and inclination with semi-major axis of all terrestrial planets in the simulations. Small markers indicate the formed terrestrial planets (we exclude remaining planetesimals). The large overlaid markers indicate the median and 10/90 percentile range. Note that there exist points outside the view of the figure, although their fraction of the total is at or below the 10 per cent level. We chose a zoomed view to focus on the differences in median and percentile spreads. Columns (left to right): reference initial planetesimal disc with different giant planet configurations, Jupiter and Saturn on present-day eccentric orbits for different planetesimal discs, Jupiter and Saturn on circular orbits for different planetesimal discs. Figure 9. Open in new tabDownload slide Final distribution of (rows, from top to bottom rows) mass, eccentricity, and inclination with semi-major axis of all terrestrial planets in the simulations. Small markers indicate the formed terrestrial planets (we exclude remaining planetesimals). The large overlaid markers indicate the median and 10/90 percentile range. Note that there exist points outside the view of the figure, although their fraction of the total is at or below the 10 per cent level. We chose a zoomed view to focus on the differences in median and percentile spreads. Columns (left to right): reference initial planetesimal disc with different giant planet configurations, Jupiter and Saturn on present-day eccentric orbits for different planetesimal discs, Jupiter and Saturn on circular orbits for different planetesimal discs. Table 4. Median terrestrial planet mass M, semi-major axis a, eccentricity e, inclination i across all 12 simulations of a set. N is the median number of planets per simulation in a set. Sub/superscripts are 10/90 percentiles. Set . M (MEarth) . a (au) . e . i (deg.) . N . NJS |$0.28_{\,0.04}^{\,0.61}$| |$2.40_{\,0.49}^{\,4.55}$| |$0.06_{\,0.02}^{\,0.16}$| |$2.89_{\,1.01}^{\,6.92}$| |$12_{\,10}^{\,12}$| EJS |$0.78_{\,0.40}^{\,1.27}$| |$0.86_{\,0.40}^{\,1.77}$| |$0.08_{\,0.03}^{\,0.24}$| |$5.18_{\,2.42}^{\,13.93}$| |$4_{\,3}^{\,5}$| EJS/Steep |$0.80_{\,0.32}^{\,1.33}$| |$0.79_{\,0.45}^{\,1.66}$| |$0.06_{\,0.03}^{\,0.14}$| |$3.31_{\,1.74}^{\,9.00}$| |$4_{\,3}^{\,4}$| EJS/Heavy |$1.14_{\,0.08}^{\,2.45}$| |$0.93_{\,0.49}^{\,2.02}$| |$0.08_{\,0.04}^{\,0.20}$| |$4.06_{\,1.68}^{\,16.48}$| |$3_{\,2}^{\,5}$| CJS |$0.52_{\,0.12}^{\,1.02}$| |$1.20_{\,0.44}^{\,2.29}$| |$0.08_{\,0.03}^{\,0.24}$| |$4.23_{\,1.76}^{\,11.03}$| |$5_{\,4}^{\,6}$| CJS/Steep |$0.49_{\,0.10}^{\,1.01}$| |$1.17_{\,0.41}^{\,2.33}$| |$0.08_{\,0.03}^{\,0.20}$| |$3.14_{\,1.48}^{\,9.37}$| |$5_{\,4}^{\,6}$| CJS/Heavy |$0.75_{\,0.20}^{\,1.74}$| |$1.03_{\,0.38}^{\,2.48}$| |$0.05_{\,0.01}^{\,0.15}$| |$2.84_{\,0.95}^{\,9.51}$| |$5_{\,3}^{\,6}$| SS |$0.46_{\,0.07}^{\,0.94}$| |$0.86_{\,0.49}^{\,1.37}$| |$0.06_{\,0.01}^{\,0.17}$| |$1.93_{\,1.61}^{\,5.10}$| |$4_{\,4}^{\,4}$| Set . M (MEarth) . a (au) . e . i (deg.) . N . NJS |$0.28_{\,0.04}^{\,0.61}$| |$2.40_{\,0.49}^{\,4.55}$| |$0.06_{\,0.02}^{\,0.16}$| |$2.89_{\,1.01}^{\,6.92}$| |$12_{\,10}^{\,12}$| EJS |$0.78_{\,0.40}^{\,1.27}$| |$0.86_{\,0.40}^{\,1.77}$| |$0.08_{\,0.03}^{\,0.24}$| |$5.18_{\,2.42}^{\,13.93}$| |$4_{\,3}^{\,5}$| EJS/Steep |$0.80_{\,0.32}^{\,1.33}$| |$0.79_{\,0.45}^{\,1.66}$| |$0.06_{\,0.03}^{\,0.14}$| |$3.31_{\,1.74}^{\,9.00}$| |$4_{\,3}^{\,4}$| EJS/Heavy |$1.14_{\,0.08}^{\,2.45}$| |$0.93_{\,0.49}^{\,2.02}$| |$0.08_{\,0.04}^{\,0.20}$| |$4.06_{\,1.68}^{\,16.48}$| |$3_{\,2}^{\,5}$| CJS |$0.52_{\,0.12}^{\,1.02}$| |$1.20_{\,0.44}^{\,2.29}$| |$0.08_{\,0.03}^{\,0.24}$| |$4.23_{\,1.76}^{\,11.03}$| |$5_{\,4}^{\,6}$| CJS/Steep |$0.49_{\,0.10}^{\,1.01}$| |$1.17_{\,0.41}^{\,2.33}$| |$0.08_{\,0.03}^{\,0.20}$| |$3.14_{\,1.48}^{\,9.37}$| |$5_{\,4}^{\,6}$| CJS/Heavy |$0.75_{\,0.20}^{\,1.74}$| |$1.03_{\,0.38}^{\,2.48}$| |$0.05_{\,0.01}^{\,0.15}$| |$2.84_{\,0.95}^{\,9.51}$| |$5_{\,3}^{\,6}$| SS |$0.46_{\,0.07}^{\,0.94}$| |$0.86_{\,0.49}^{\,1.37}$| |$0.06_{\,0.01}^{\,0.17}$| |$1.93_{\,1.61}^{\,5.10}$| |$4_{\,4}^{\,4}$| Open in new tab Table 4. Median terrestrial planet mass M, semi-major axis a, eccentricity e, inclination i across all 12 simulations of a set. N is the median number of planets per simulation in a set. Sub/superscripts are 10/90 percentiles. Set . M (MEarth) . a (au) . e . i (deg.) . N . NJS |$0.28_{\,0.04}^{\,0.61}$| |$2.40_{\,0.49}^{\,4.55}$| |$0.06_{\,0.02}^{\,0.16}$| |$2.89_{\,1.01}^{\,6.92}$| |$12_{\,10}^{\,12}$| EJS |$0.78_{\,0.40}^{\,1.27}$| |$0.86_{\,0.40}^{\,1.77}$| |$0.08_{\,0.03}^{\,0.24}$| |$5.18_{\,2.42}^{\,13.93}$| |$4_{\,3}^{\,5}$| EJS/Steep |$0.80_{\,0.32}^{\,1.33}$| |$0.79_{\,0.45}^{\,1.66}$| |$0.06_{\,0.03}^{\,0.14}$| |$3.31_{\,1.74}^{\,9.00}$| |$4_{\,3}^{\,4}$| EJS/Heavy |$1.14_{\,0.08}^{\,2.45}$| |$0.93_{\,0.49}^{\,2.02}$| |$0.08_{\,0.04}^{\,0.20}$| |$4.06_{\,1.68}^{\,16.48}$| |$3_{\,2}^{\,5}$| CJS |$0.52_{\,0.12}^{\,1.02}$| |$1.20_{\,0.44}^{\,2.29}$| |$0.08_{\,0.03}^{\,0.24}$| |$4.23_{\,1.76}^{\,11.03}$| |$5_{\,4}^{\,6}$| CJS/Steep |$0.49_{\,0.10}^{\,1.01}$| |$1.17_{\,0.41}^{\,2.33}$| |$0.08_{\,0.03}^{\,0.20}$| |$3.14_{\,1.48}^{\,9.37}$| |$5_{\,4}^{\,6}$| CJS/Heavy |$0.75_{\,0.20}^{\,1.74}$| |$1.03_{\,0.38}^{\,2.48}$| |$0.05_{\,0.01}^{\,0.15}$| |$2.84_{\,0.95}^{\,9.51}$| |$5_{\,3}^{\,6}$| SS |$0.46_{\,0.07}^{\,0.94}$| |$0.86_{\,0.49}^{\,1.37}$| |$0.06_{\,0.01}^{\,0.17}$| |$1.93_{\,1.61}^{\,5.10}$| |$4_{\,4}^{\,4}$| Set . M (MEarth) . a (au) . e . i (deg.) . N . NJS |$0.28_{\,0.04}^{\,0.61}$| |$2.40_{\,0.49}^{\,4.55}$| |$0.06_{\,0.02}^{\,0.16}$| |$2.89_{\,1.01}^{\,6.92}$| |$12_{\,10}^{\,12}$| EJS |$0.78_{\,0.40}^{\,1.27}$| |$0.86_{\,0.40}^{\,1.77}$| |$0.08_{\,0.03}^{\,0.24}$| |$5.18_{\,2.42}^{\,13.93}$| |$4_{\,3}^{\,5}$| EJS/Steep |$0.80_{\,0.32}^{\,1.33}$| |$0.79_{\,0.45}^{\,1.66}$| |$0.06_{\,0.03}^{\,0.14}$| |$3.31_{\,1.74}^{\,9.00}$| |$4_{\,3}^{\,4}$| EJS/Heavy |$1.14_{\,0.08}^{\,2.45}$| |$0.93_{\,0.49}^{\,2.02}$| |$0.08_{\,0.04}^{\,0.20}$| |$4.06_{\,1.68}^{\,16.48}$| |$3_{\,2}^{\,5}$| CJS |$0.52_{\,0.12}^{\,1.02}$| |$1.20_{\,0.44}^{\,2.29}$| |$0.08_{\,0.03}^{\,0.24}$| |$4.23_{\,1.76}^{\,11.03}$| |$5_{\,4}^{\,6}$| CJS/Steep |$0.49_{\,0.10}^{\,1.01}$| |$1.17_{\,0.41}^{\,2.33}$| |$0.08_{\,0.03}^{\,0.20}$| |$3.14_{\,1.48}^{\,9.37}$| |$5_{\,4}^{\,6}$| CJS/Heavy |$0.75_{\,0.20}^{\,1.74}$| |$1.03_{\,0.38}^{\,2.48}$| |$0.05_{\,0.01}^{\,0.15}$| |$2.84_{\,0.95}^{\,9.51}$| |$5_{\,3}^{\,6}$| SS |$0.46_{\,0.07}^{\,0.94}$| |$0.86_{\,0.49}^{\,1.37}$| |$0.06_{\,0.01}^{\,0.17}$| |$1.93_{\,1.61}^{\,5.10}$| |$4_{\,4}^{\,4}$| Open in new tab The presence and configuration of giant planets have the largest effect on the final configuration of terrestrial planets. In the absence of giant planets, we find 10 to 12 terrestrial planets per system. Of these, 90 per cent have mass M < 0.6 MEarth, and half have M < 0.3 MEarth. By number, they are spread evenly across a wide range of semi-major axes with more massive planets closer to the host star. If giant planets are present, the number of terrestrial planets per systems drops to 4 or 5. They are more massive with 50 per cent having M > 0.5 to 0.8 MEarth, depending on the giant planet configuration. Most (90 per cent) are below 1.2 MEarth, but we find isolated 1.6 MEarth planets. They also cover a much narrower range of semi-major axes and are closer to the host star. Systems with Jupiter and Saturn on eccentric orbits (EJS) tend to have (by median) 50 per cent more massive planets than systems with Jupiter and Saturn on circular orbits (CJS). They place terrestrial planets at smaller semi-maxis axis. Note that these are broad trends. Both configurations can produce terrestrial planets M > 1 MEarth, although EJS simulations are more likely to do so. Simulations with eccentric giants planets also tend to have one or two fewer planets per system. Except for a few outliers, all planets are in orbits with low inclination and eccentricity. Overall, 50 per cent of all planets are at inclinations i < 5° and eccentricities e < 0.1. More inclined and eccentric orbits are restricted to simulations with giant planets, where 10 per cent of planets can reach i > 15° or e > 0.25. These observations are consistent with those of Levison & Agnor (2003), who report systems hosting embryos on more eccentric orbits to result in fewer, more massive, and closer-in terrestrial planets. Excitation of planetesimal (and ultimately terrestrial planet) eccentricity increases in the presence of giant planets, especially when those are on eccentric orbits. Raymond et al. (2004) find similar correlations and also note that the total mass of terrestrial planets tends to be slightly lower in systems with more eccentric Jupiters. We observe the opposite (cf. Fig. 3), but point out their use of a different initial surface density profile. While the baseline simulations in Raymond, Scalo & Meadows (2007) do not host giant planets, they have performed some tests on their influence. Overall, they find the effect of adding a Jupiter-mass planet to be very small (the mean terrestrial planet mass increases by 10 per cent). However, their runs are initialized at a stage when planetary embryos have formed and the gaseous disc has dissipated. For a given giant planet configuration, changing the slope of the initial planetesimal mass distribution has a much smaller effect than adjusting the mass of planetesimal disc. In particular, changing the reference profile to a steeper profile changes the median mass and semi-major axis of the formed terrestrial planets by less than 10 per cent for both giant planet configurations. Additionally, the width of the range in mass and semi-major axis populated by 80 per cent of terrestrial planets also remains similar. Orbital inclinations decrease by 26 (CJS to CJS/Steep) to 36 (EJS to EJS/Steep) per cent. Eccentricities are only visibly stunted for EJS runs with the median eccentricity for EJS/Steep about by 25 per cent below that of EJS runs. Previously, Raymond et al. (2005b) also modelled the effect of varying the initial surface density profiles. They noted that steepening the profile increases the (mean) mass of the most massive planet, accelerates their formation, increases mean number of plants in a given system, and leads to more planets at smaller semi-major axes. They did not observe a conclusive trend for the eccentricity of the final planets. Some of these trends are in line with ours (more massive planets at smaller semi-major axes) while others are not (number of planets per system) with an additional dependency on our particular configuration of giant planets. However, besides similar initial surface density profiles, they initialize simulations in a distinctly different way. Their runs launch with fully formed embryos (no gas and no planetesimals) and have a single giant planet on a circular orbit. In our runs, sweeping resonances drive the final architecture of our systems. This mechanism is absent in their runs and their final architecture is driven solely by planet–planet scattering and interaction with the giant at mean motion resonances. Predictably, heavier planetesimal discs produce more massive terrestrial planets. The most massive planets are 3.71 MEarth (EJS) and 3.05 MEarth (CJS). At the 90 per cent level, terrestrial planets in initially massive discs exceed those in the reference by 93 (70) per cent in EJS (CJS) configurations. The situation is less clean-cut in semi-major axis, eccentricity, and inclination. For systems with eccentric giant planets, there is weak trend for terrestrial planets to end up at larger semi-major axis, although median only differs by 8 per cent. For systems with giants on circular orbits, the trend is reversed such that the terrestrial planets tend to form 15 per cent closer to the stars. Dynamically, planets tend to be less excited in initially massive discs. This holds especially in CJS configurations where eccentricities and inclinations decrease by 35 per cent. In EJS runs, the effect is weaker. Here, eccentricities are approximately the same. Although inclinations are reduced by 20 per cent for most planets, they in fact increase for at least a tenth of all planets. Although starting from entirely different ICs, both Kokubo et al. (2006) and Raymond et al. (2007) also find the planet mass to scale with the solid disc mass. However, Kokubo et al. (2006) have very different ICs. They start with a distribution of 16 to 32 embryos between 0.5 and 1.5 au and do no include giant planets. Qualitatively, this is close to the state of our EJS simulations after 10 Myr when the gas has dissipated and inward sweeping secular resonances have set the final architecture of the system. We do not reproduce their observation that the number of planets decreases with the solid disc mass, likely because the number of embryos is set by sweeping secular resonances for us, whereas Kokubo et al. (2006) treat this number as a free parameter. The simulations in Raymond et al. (2007) use between 75 and 190 embryos, but also spread them out over a wider range of semi-major axes. They are therefore slightly more comparable to ours, but again do not include the effect of giant planets or assembly of embryos from planetesimals in the presence of gas. 3.6 How sweeping resonances sculpt the disc In Sections 3.2, 3.3, and 3.5, we explored the key differences between our simulation runs. We found that the presence, absence, or orbital configuration of the giant planets has the strongest influence on the forming terrestrial planets. In runs with giant planets, terrestrial planets form on tighter orbits, and are less numerous, more massive, and dynamically hotter than in those without. In EJS configurations, orbits are tighter and planets more massive and with slightly higher eccentricities and inclinations than in CJS configurations. Changing the initial planetesimal distribution also had noticeable (but less drastic) effects. In particular, more massive initial planetesimals discs lead to more massive terrestrial planets. Steeper discs are more fickle, producing more massive planets on tighter orbits for EJS runs, but the opposite in CJS configurations. We now disseminate how the interactions between planetesimals, giant planets, and the massive gas disc drive these trends. 3.6.1 Trapping planetesimals in resonances The bulk of the architectural properties of the final planetary systems are determined during the first few Myr when planetesimals – embedded in a gaseous disc – first grow into embryos through mutual collisions. At this point, their evolution is driven by their mutual gravitational interaction, their interaction with the gaseous disc, and with their interaction with the giant planets. Once the gas dissipates and the embryos have cleared out their immediate spheres of influence (a few Hill radii), few changes occur in the system architectures. The most defining difference between simulations is the presence and orbital configuration of giant planets. Through mean motion and secular resonances, they can transfer significant angular momentum on to the planetesimals efficiently. This exchange also proceeds without resonances, but is less effective. Unless dampened, thus excited planetesimals can be launched into the inner Solar system or be ejected from the system. Planetesimals can also be repelled by or trapped in resonances. Those repelled are effectively blocked off from accessing regions of phase-space. Planetesimals which are locked in resonances find their eccentricities and inclinations perpetually excited unless dampened.11 Although these processes occur within the first few million years, they induce a signature of a dynamically hotter state which persists until the end of our simulations. Here, we find embryos and planetesimals in EJS and CJS runs at higher eccentricities and inclinations than those in NJS runs. As illustrated in Section 3.1.2, planetesimals that become trapped in secular resonances are swept along with the resonances, eventually settling near their final locations once the gas has dissipated (∼10 Myr). In EJS and CJS configurations, the inward sweeping ν5 resonance delivers large amounts of material to its final location (a ∼ 0.9 au). This significantly enhances the amount of material available for embryos to accrete from, thereby peaking the radial mass distribution. Planetesimals that manage to avoid being swept along by the ν5 resonance face the ν6 resonance soon after, which sweeps most of the remaining planetesimals inward, although some manage to avoid this fate. Although it does not drag planetesimals along, the ν16 resonance also sweeps outward (Thommes et al. 2008) and settles near ν6. Once the gas is gone, their location demarcates the outer boundary of the region populated by terrestrial planets. Those trapped in their vicinity find themselves excited on to eccentric orbits that either deliver them to the inner system to collide with embryos, fall on to the star, or be ejected from the system. After 147.84 Myr, most of the planetesimals beyond the ν6 and ν15 resonances have been removed, although stragglers remain. By being separated from the inner regions, they have avoided being accreted on to larger embryos, thus covering a wide range of masses. In systems without giant planets, on the other hand, no secular resonances exist, and planetesimals remain in region they are initially placed in (modulo inward drift due to hydrodynamic drag and gravitational scattering). The net effect of inward sweeping ν5 and ν6 resonances is then to deliver planetesimals from the disc and trap them in the orbital region between. The region between these resonances is narrower in EJS configurations than in CJS ones. Naturally, this leads to a higher mass concentration in EJS runs as well as a more complete conversion of planetesimals into planets, i.e. only a very weak bimodal signature in the mass–frequency distribution. Configurations with giant planets on circular orbits trap less mass per au in the inner regions, which retards conversion of planetesimals into planets. As such, they retain a bimodal mass–frequency distribution after 147.84 Myr. As the conversion of planetesimals into embryos is more complete and sourced from a tighter orbital region with dynamically more excited planetesimals (two resonances are close-by), it is unsurprising to find more massive terrestrial planets in EJS configurations. In systems with no giant planets, material is very much spread out, so that embryos only accrete slowly. This leads to low-mass terrestrial planets which are much more spread out in semi-major axis. In discs with initially steep density profiles, sweeping resonances appear to be more efficient in shepherding material inwards than in the reference case. In this configuration, the systems are initialized with fewer planetesimals at larger semi-major axis so that (i) the resonances have to sweep less material along and (ii) the radial mass distributions are already biased towards smaller semi-major axes. The net result then is for systems with steep initial density profiles to have radial mass distributions and system architectures to appear almost indistinguishable from the reference discs. On the other hand, the initial bias of material to the inner regions also means that fewer planetesimals are dynamically excited by the sweeping resonances, resulting in a dynamically colder state with smaller eccentricities and inclinations. If the initial planetesimal disc is massive, the dynamical evolution of planetesimals is more complicated. Although the initial excitement from resonances is similar to configurations with the reference and steep profiles, dampening by hydrodynamic drag is much less effective due to the larger planetesimal mass (see fig. 1 in Morishima et al. 2010). This has three consequences. First, sweeping secular resonances move fewer planetesimals inward.12 Secondly, the richly excited planetesimals actually remove angular momentum from Jupiter through dynamical friction [e.g. Kokubo & Ida (2012)], although this also happens to a lesser degree for discs following the reference and steep profiles. For example, the mean eccentricity of Jupiter during the first 312 kyr is eJ, Ref = 0.0678 and eJ, Heavy = 0.0738 for Run01 in the EJS and EJS/Heavy sets. After 10 Myr, the 312 kyr averaged eccentricities are eJ, Ref = 0.0474 and eJ, Heavy = 0.0385, corresponding to a decrease of 43 and 91 per cent, respectively. Thirdly, a small population of planetesimals actually manage to move outwards past the orbits of Jupiter and Saturn. Although most are promptly ejected (on time-scales of ∼105 years), a single planetesimal (across 12 runs) is caught by Jupiter as a moon!13 The final radial mass distributions and mass–frequency distributions in initially massive discs still resemble those of the reference configurations, which may be counterintuitive. As sweeping resonances move fewer planetesimals, we would expect proportionally less mass to end up in the inner system. However, this appears to be amply compensated for by the larger mass of the planetesimals. The only striking remaining difference is that a larger amount of mass remains in the disc simply by virtue of starting out with more mass. Dynamically, the effect of more massive initial planetesimals is also somewhat opaque. While more massive planetesimals are more difficult to excite by sweeping resonances, their excitations are also more difficult to dampen through hydrodynamic drag due to their size. Taking these considerations together, it appears that the more excited planetesimals are ejected from the system (or fall into the star) while the less excited ones are swept along into the inner system to build up terrestrial planets. 3.6.2 Statistical spread of runs The dynamical outcome of a planetesimal passing through a resonance depends sensitively on the initial orbit. For two initially almost identical orbits, the more resonances the particles pass through, the farther their orbits can diverge. Therefore, simulations where particles are exposed to a larger number of resonances have more evolutionary pathways available. The wider spread of EJS and CJS runs with respect to NJS is then hardly surprising. By the same argument, we naively expect EJS runs to have a larger spread than CJS runs because more particles cross the sweeping ν6 resonance. While this holds at the 25/75 percentile level, it does not always hold at the 10/90 percentile level. Closer inspection of Fig. 2 reveals that in CJS runs, the 3:1 and 5:2 mean motion resonances with Jupiter remains within the planetesimal disc, and the final location of the ν6 and ν16 (not shown) resonances is within 0.07 au of the 7:3 resonance. Over Myr time-scales, such close stacking provides pathways for particle orbits to diverge and promotes ejection of bodies from the disc. As the effects of secular resonances and mean motion resonance cannot be readily distinguished, we caution from attributing to them any differences across CJS and EJS runs. 3.7 The spread of typical diagnostics (AMD, RMC) Simulations of terrestrial planet formation typically invoke diagnostics such as the Angular Momentum Deficit (AMD) and Radial Mass Concentration (RMC) to quantify how close they are to reproducing some observed reference system (which is usually the Solar system). Due to their stochastic nature, we expect initially identical simulations to produce different values for the AMD and RMC. If values obtained from a given IC exhibit a large spread, they may overlap with values generated from a different IC. This makes it unclear whether the outcome arises from a different ICs or is merely a reflection of the stochastic nature of these simulations. This suggests that individual runs lack predictive power. We now characterize the spread in AMD and RMC obtained in our simulations. The RMC is given by \begin{equation} \mathrm{RMC} = \max \left( \frac{\sum _j m_j}{\sum _j m_j \left[ \log _{10} \left( a / a_j \right) \right]^2} \right), \end{equation} (4) where the sum is over the mass mj and semi-major axis aj of the terrestrial planets, and we search for the maximum in a semi-major axis range covering all terrestrial planets (Laskar 1997). The AMD is \begin{equation} \mathrm{AMD} = \frac{\sum _j m_j \sqrt{a_j} \left( 1 - \cos (i_j) \, \sqrt{1 - e_j^2} \right)}{\sum _j m_j \sqrt{a_j}}, \end{equation} (5) where mj, aj, ej, and ij are the mass, semi-major axis, eccentricity, and inclinations of the terrestrial planets (Chambers 2001). Fig. 10 shows the AMD and RMC from seven of our nine simulation sets as well as corresponding simulations from R09 who have used similar initial surface density profiles as well as giant planet orbits. Box plots are generated from (i) 12 runs per set (our runs), or (ii) four runs per set (R09). Figure 10. Open in new tabDownload slide Box plots of the AMD and RMC obtained from (i) seven (out of nine) of our sets (NJS, EJS, EJS/Steep, CJS, CJS/Steep) and (ii) corresponding simulations from R09 (R09 EJS1, R09 EJS15, R09 CJS1, R09 CJS15). For visual aid, we indicate similar giant planet configurations through the background colour. Boxes indicate the 25/75 range with the median marked in the box. Whiskers extend to 1.5 times the inter-quartile range and crosses indicate outliers. Figure 10. Open in new tabDownload slide Box plots of the AMD and RMC obtained from (i) seven (out of nine) of our sets (NJS, EJS, EJS/Steep, CJS, CJS/Steep) and (ii) corresponding simulations from R09 (R09 EJS1, R09 EJS15, R09 CJS1, R09 CJS15). For visual aid, we indicate similar giant planet configurations through the background colour. Boxes indicate the 25/75 range with the median marked in the box. Whiskers extend to 1.5 times the inter-quartile range and crosses indicate outliers. Deferring a comparison to the simulations of R09 to Section 4.2, we find that our simulations without giant planets have small spreads in the AMD and RMC values. Simulations including giant planets have larger spreads, especially in the AMD. Values derived from EJS and CJS sets also overlap, and only the overall statistics reveal a trend towards smaller RMC values in CJS runs. Spreads in runs with a steeper planetesimal distribution tend to be smaller except for RMC in EJS/Steep runs. Overall, there is significant spread in all diagnostics, which emphasizes the need for a statistical approach. Running only one simulation per set, pathological cases with reversed trends could arise, leading us to draw potentially wrongful conclusions. 4 DISCUSSION Based on the above results, we now address a serious caveat associated with the paradigm of sweeping secular resonances, demonstrate large statistical spreads in typical diagnostics used in terrestrial planet formation, and attempt a preliminary link of our numerical results to observations. 4.1 Sweeping resonances: summary and caveats Above, we have reinforced the notion of sweeping secular resonances to drive the dynamics of planetesimal discs. Depending on the configuration of giant planets, the ν5, ν6, ν15, and ν16 secular resonances with Jupiter and Saturn sweep inward as the gas dissipates until they settle in at their final positions (the inner ν16 resonances ‘falls into the star’, as it were). Based on their dynamics, the secular resonances end up truncating the planetesimal disc which constrains growth of embryos to particular regions. These dynamics are in line with previous work on secular resonances (as outlined in the Introduction), although we demonstrate that they extend to the dynamics of planetesimals instead of only those of embryos and asteroids (which are, respectively, more massive or massless as far as giant planets are concerned). However, a caveat exists in all numerical implementations above (including ours) – the feedback of the gaseous disc on to the giant planets is neglected (Raymond & Morbidelli 2014). As such, the eccentricity and inclinations of the giant planets will be dampened even beyond the current transfer of angular momentum on to the planetesimals and embryos. This modifies the time evolution and strength of the sweeping resonances which in turn affects the final configuration of terrestrial planets. For our simulations, this means that we may expect the EJS runs to become a little more like CJS runs. In other words, we most likely overestimate the mass and underestimate the semi-major axes of the final planets. Future work should clearly take into account a more sophisticated prescription of the gaseous protoplanetary disc that feeds back on to both the giant planets and the planetesimals and embryos. 4.2 Spread of RMC and AMD in solar system formation In Section 3.7, we presented the spreads of the RMC and AMD in our simulation runs. As alluded therein, we now extend this to compare our spreads to those reported in previous work. So far, we have deliberately omitted reference to the Solar system, but note that the AMD and RMC diagnostics are frequently used to assess how close a given model manages to reproduce the Solar system. As such, any comparison to previous work that reports on the AMD and RMC inevitably compares our simulations to those aimed at reproducing the Solar system. This warrants a short overview of available models. Although our ICs are based on simulations for exploring the formation of the Solar system, our goal was not to explicitly assess their viability for this. They are in fact not very well suited. For example, three key constraints in the Solar system are the comparatively low mass of Mars, the compactness of terrestrial planet spacing, and the small eccentricities and inclinations of the terrestrial planets. While CJS configurations generate low-eccentricity planets, they also spread out the terrestrial planets too much and produce too massive Mars-analogues. In EJS configurations, the terrestrial planets tend to be more concentrated and low-mass Mars-analogues are abundant, but they tend to be too eccentric. In terms of the RMC and AMD, CJS configurations roughly match constrains for the AMD but not for the RMC; and vice versa for EJS configurations. Barring unsuccessful attempts to sculpt extended planetesimals discs with extremely eccentric giant planets (R09) or excitement of terrestrial planets by a late migration of giant planets (Brasser et al. 2009; Agnor & Lin 2012; Brasser, Walsh & Nesvorný 2013), the most suitable ICs for generating terrestrial planets (up to factors of a few) fulfil AMD and RMC constraints of the Solar system require an initial planetesimal disc of extending only between 0.7 and 1.0 au (Morishima et al. 2008; Hansen 2009) instead of between 0.5 and 4.0 au. It remains unclear whether the outer edge of the annulus results from a truncation by migration of giant planets (Walsh et al. 2011) or is simply indicative of the initial mass distribution of solids in the protosolar nebula (Chambers & Cassen 2002; Jin et al. 2008; Izidoro et al. 2014a). The inner edge may be sculpted by material being trapped and piling up in a pressure bump (Youdin & Chiang 2004; Johansen, Youdin & Klahr 2009). For a more thorough exploration of ICs suitable to reproduce the Solar system terrestrial planets, we refer the interested readers to Izidoro et al. (2014a,b), Walsh et al. (2011), Morishima et al. (2010), R09, Kokubo et al. (2006), and Chambers (2001) as well as the review of Raymond et al. (2014). A number of authors report AMD and RMC for their Solar system formation simulations, but have only run a single instance for each set of simulation parameters (Chambers 2001; Morishima et al. 2010). To our knowledge, only Kokubo et al. (2006), R09, as well as Izidoro et al. (2014a) account for stochastic variations across runs by running similar ICs multiple times.14|$^,$|15 Of these, R09 consider and present simulations most similar to our own, and we include some of their results here for comparison (cf. Fig. 10). Although the ranges in AMD and RMC populated by our runs and those of R09 tend to be misaligned, they generally show overlap. More importantly, at least some trends resulting from changing the orbits of the giant planets and initial planetesimal mass distribution appear to be robust. For example, in CJS runs, moving to steeper initial planetesimal distribution decreases the AMD while the RMC remain approximately the same. For EJS runs, steeper initial profiles also decrease the AMD in both cases, although the RMC responds differently in the runs of R09 vis-à-vis our runs. Changing the orbits of the giant planets is much more robust and exhibits the same trends in both sets (when judging by the median). Overall, it appears that runs in R09 follow clear trends and are not as burdened by overly wide spreads and outliers as ours. Unfortunately, this clarity appears to be a consequences of simplifications made to conform to the computational resources available at the time. In particular, recall that depending on the giant planet configuration and initial planetesimal distribution, inward sweeping secular resonances generate a particular distribution of embryos after the gas disc has dissipated. It is only at this point that R09 initializes their simulations, by necessity with a more generic embryo distribution. This distribution lacks the imprint of the important dynamics that took place as the gas dissipated. As much of the shaping of the final system has already happened at this point, it should not be surprising that their runs produce smaller spreads in diagnostics.16 It appears that modelling of the initial phase of planetesimals embedded in a dissipating gas disc is essential. Omitting this stage neglects much of the dynamics that shape the final architecture of the system. 4.3 Preliminary comparison to observations At this point, we may wonder whether the key trends found in our simulations also hold in extrasolar planetary systems. In other words, are the observed terrestrial planets17 more massive, less numerous, and closer to the host star if the system also hosts giant planets? Deferring caveats of our admittedly naive comparison for now, Fig. 11 shows the mass and semi-major axis of all sub-Neptune mass planets found in multiple planet systems, grouped into whether the system also hosts giant planets or not. We also summarize their statistical properties in Table 5. Exoplanet observations are loaded from the Extrasolar Planet Encyclopaedia18 data base (Schneider et al. 2011) using a snapshot from 2015 August 02. We only consider systems hosting at least two terrestrial planets and require complete data on semi-major axis and mass for all planets in system. Raw data from a total of 1228 planetary systems are reduced to 62 systems, 43 of which host giant planets. Figure 11. Open in new tabDownload slide Mass and semi-major axis for all sub-Neptunes (M < 17.14 MEarth) in the reduced (see text) sample of observed exoplanets. Small points are observational data. Light and dark grey points are planets in systems with and without detected giant (M ≥ 50 MEarth) planets in the same system. Large circles and lines indicate the median as well as 10/90 percentiles for the two sets. Figure 11. Open in new tabDownload slide Mass and semi-major axis for all sub-Neptunes (M < 17.14 MEarth) in the reduced (see text) sample of observed exoplanets. Small points are observational data. Light and dark grey points are planets in systems with and without detected giant (M ≥ 50 MEarth) planets in the same system. Large circles and lines indicate the median as well as 10/90 percentiles for the two sets. Table 5. Median and 10/90 percentiles of mass M, semi-major axis a, and number of planets NP for the sub-Neptune (M < MEarth) population of two sets of exoplanet systems. The number of system fulfilling the filtering criteria is NS, and fS the corresponding fraction out of 1228 total systems. Set . M (MEarth) . a (au) . NP . NS . fS (%) . No giants |$6.61_{-3.84}^{+6.27}$| |$0.12_{-0.07}^{+0.23}$| |$2_{-1}^{+1}$| 43 3.5 Giants |$8.58_{-5.68}^{+4.52}$| |$0.09_{-0.05}^{+0.16}$| |$1_{-0}^{+2}$| 19 1.5 Set . M (MEarth) . a (au) . NP . NS . fS (%) . No giants |$6.61_{-3.84}^{+6.27}$| |$0.12_{-0.07}^{+0.23}$| |$2_{-1}^{+1}$| 43 3.5 Giants |$8.58_{-5.68}^{+4.52}$| |$0.09_{-0.05}^{+0.16}$| |$1_{-0}^{+2}$| 19 1.5 Open in new tab Table 5. Median and 10/90 percentiles of mass M, semi-major axis a, and number of planets NP for the sub-Neptune (M < MEarth) population of two sets of exoplanet systems. The number of system fulfilling the filtering criteria is NS, and fS the corresponding fraction out of 1228 total systems. Set . M (MEarth) . a (au) . NP . NS . fS (%) . No giants |$6.61_{-3.84}^{+6.27}$| |$0.12_{-0.07}^{+0.23}$| |$2_{-1}^{+1}$| 43 3.5 Giants |$8.58_{-5.68}^{+4.52}$| |$0.09_{-0.05}^{+0.16}$| |$1_{-0}^{+2}$| 19 1.5 Set . M (MEarth) . a (au) . NP . NS . fS (%) . No giants |$6.61_{-3.84}^{+6.27}$| |$0.12_{-0.07}^{+0.23}$| |$2_{-1}^{+1}$| 43 3.5 Giants |$8.58_{-5.68}^{+4.52}$| |$0.09_{-0.05}^{+0.16}$| |$1_{-0}^{+2}$| 19 1.5 Open in new tab Overall, the observations appear to corroborate the results of our simulations. In systems hosting giant planets, the median terrestrial planet mass is 30 per cent larger, the median semi-major axis 25 per cent smaller, and the terrestrial planets are generally less numerous. Although encouraging, these results clearly suffer from two major caveats. First, the observed planets are much more massive than those in our simulations. Secondly, the observed distribution of planets peaks around 0.1 au while our simulations concentrate planets at about 1.0 au. Superficially, the differences in planetary masses can be resolved by imposing that the observed planets formed in a more massive initial solid disc. In Section 3, we found that doubling the initial mass in planetesimals tends to increase the median mass of the final terrestrial planets by about 50 per cent. By this measure, we would require the initial protoplanetary discs to host about 80 MEarth of solid material (four doublings), although much larger required enhancements have also been suggested (Schlichting 2014). This does not account for the accretion of gas on to a rapidly formed core, which would account for some of the additional mass. A steeper initial surface density profile can further decrease the required total mass in solids throughout the disc by concentrating more material in the inner regions. However, surface density profiles steeper than ΣSolids ∝ r−1.5 are not supported by observations (Hughes et al. 2008; Andrews et al. 2009, 2010; Isella, Carpenter & Sargent 2009; Williams & Cieza 2011). A more massive initial disc and steeper surface density profile would significantly shorten the formation time-scale of terrestrial planets. This especially holds for those that may have formed in situ, i.e. at their detected orbital location around 0.1 au. Absent gas and giant planets, the time-scale for a planet to grow to isolation mass scales roughly as tIsolation ∝ a3/(ΣSolids)1/2 (Ida 1990; Kokubo & Ida 2012). In other words, while increasing the total solid disc mass by four doublings shortens the formation time-scale by only a factor of 4, in situ formation at 0.1 au (instead of at 1.0 au) accelerates accretion of planetary embryos by a factor 1000. In this case, terrestrial planets would have formed much faster than the giant planets, which implies that their inward sweeping secular resonances would have little effect on the already formed planets by virtue of their large mass. On the other hand, planetesimals accrete while being immersed in the gaseous disc. Depending on the disc structure, the presence or absence of giant planets, the gas dissipation rate, and the mass of the forming embryos, radial migration of planetesimals, embryos, and formed terrestrial planets is an almost inescapable consequence (Nagasawa, Lin & Thommes 2005; Masset, D'Angelo & Kley 2006; Terquem & Papaloizou 2007; Cresswell & Nelson 2008; Thommes et al. 2008; McNeil & Nelson 2010; Cossou et al. 2014). Planetesimals move inwards either from hydrodynamic drag, secular resonance sweeping, or shepherding while migration of embryos and planets is mostly driven by disc interactions. This scenario is conceptually much more in line with our simulations and provides sufficiently long formation time-scales for giant planets to affect the inner planets as they form. In either scenarios, embryos and terrestrial planets migrate inwards as they interact with the dissipating gas disc with the rate and direction of migration sensitively depending on the structure of the gas disc (Paardekooper et al. 2010; Paardekooper, Baruteau & Kley 2011; Bitsch et al. 2013, 2014a,b). Typically, migration terminates at the inner edge of the disc (∼0.1 au) where planets tend to pile up in resonant chains (Terquem & Papaloizou 2007; Raymond, Barnes & Mandell 2008; Ogihara & Ida 2009; Cossou et al. 2014), although this is inconsistent with observed period ratios (Lissauer et al. 2011). In our simulations, we remove planets approaching to within 0.2 au of the host star to meet integration accuracy requirements. This means that we potentially neglect a population of planets interior 0.2 au. To test for this, we have carried out additional test runs (at lower accuracy) where we only remove material that approaches the host star to within two solar radii (∼0.005 au). In these tests, we indeed observe a clustering of planets around the inner edge of the gas disc, which we have placed at 0.1 au.19 As such, we are conceptually able to reproduce the observed mass distribution of extrasolar terrestrial planets. In these trials, however, the formation of planets interior to 0.2 au appears to only be loosely coupled to the evolution of the giant planets. Initially, they grow solely by accretion of local material while migrating inwards until their eventual pile-up at inner edge of the gaseous disc. Compared to embryos beyond 0.2 au, they grow faster, have a smaller final mass,20 and source comparatively small amounts of material from the regions affected by inward sweeping resonances.21 Given the differences in mass and distribution of terrestrial planets between our simulations and observations, a direct comparison between the two may not be well-justified as yet. Nevertheless, it is encouraging that basic trends persist, such that we expect more suitable simulations to be able to match the basic properties of at least some of the exoplanetary systems discussed above. In such a scenario, we also expect to be able to predict as yet undetected giant planets (and their orbital configuration) by considering outliers in the system architectures of the terrestrial (sub-Neptune mass) planets. For example, if the terrestrial planets in a system without detected giant planets are significantly more massive, less numerous, and closer-in than the median values for this population, we may expect undetected giant planets to be present. A cursory exploration of this scenario suggests that CoRoT-7, HD 20003, and HD 20781 may host undetected giants. These systems should be the first targets in a future analysis. 5 SUMMARY AND CONCLUSIONS In this work, we have addressed the stochastic nature of terrestrial planet formation. Orbits which are initially identical up to one part in 1015 (∼double precision floating point accuracy) diverge exponentially until their separation becomes limited by the accessible phase-space. The divergence is driven by differences in the geometry of successive close encounters with other particles. We find that the rate of divergence increases with the number of simulation particles. After a few hundred years, orbits have diverged far enough for the collisional history of two simulations with initially nearby orbits to differ. They are now fully diverged, and produce distinctly different planetary systems within our simulation times of ∼147.84 Myr. If the position of a single planetesimal is changed by less than one millimetre, the positions, orbits and masses of the resulting planets are all different. We find that nearby orbits diverge exponentially with an e-folding time-scale of a few to a few tens of years, which decreasing as the number of particles increases. There is no reason to expect that this behaviour does not continue to much smaller scales. Perhaps if our early Solar system had contained one extra molecule, the Earth would not have formed at all. Even if simulations are initialized with identical ICs, variations in round-off errors quickly seed differences in the planetesimal orbits, which are amplified by sequences of close encounters. The variations are caused by the ill-defined behaviour of certain intrinsic functions available in parallel programming. They cause loss of control over the order of operations in the code. Although this can be mitigated, the extreme sensitivity of the order of operations implies the existence of multiple – numerically equally valid – final configurations for any given IC. Each of these configurations can give wildly different values for typical diagnostics. For example, we find variations in the final surface density profiles and note that the final profiles are not correlated with the initial ones. Individual simulations therefore lack predictive power, and we are relegated to considering distributions and stacks of multiple runs. Unfortunately at present, it remains unclear how many runs are required to adequately sample the systems. We advocate at least eight runs per IC. This would generate two samples per quartile given an underlying uniform distribution. Analysing our simulations, we find that varying the configuration of giant planets in systems has statistically robust results in line with previous work. Planetesimal dynamics are driven primarily by sweeping secular resonances and the extend of the terrestrial planet is set by the final locations of the resonances. Systems with giant planets tend to form fewer, more massive, and more eccentric terrestrial planets at smaller semi-major axes than those without. Although we probe different mass and orbital regimes, observations of sub-Neptune mass planets appear to support these trends. A proof-of-concept extraction of outliers in the observations suggests that CoRoT-7, HD 20003, and HD 20781 may host as yet undetected giant planets. Acknowledgments We are most grateful to the anonymous referee as well as Sean Raymond for a host of useful comments that improved the quality of this manuscript. We also wish to thank Ryuji Morishima, Rok Roškar, Thomas Peters, George Lake, Michael Rieder, and Joanna Dra̧żkowska for useful discussions, as well as Doug Potter for computational support. Simulations were run on the Tasna and Vesta GPU clusters at the Institute for Computational Science, University of Zürich. We also acknowledge use of Numpy, Scipy, IPython/Jupyter, Matplotlib, and Colorbrewer2 for post-processing. 1 " An extensive overview on resonances is given in Lecar et al. (2001). 2 " The literature on terrestrial planet formation can be overwhelming. Excellent reviews include that of Raymond et al. (2014) and Morbidelli et al. (2012). 3 " Instead of parameter sweeps, Richardson et al. (2000) tackled numerical issues and pushed the number of massive particles to 106. 4 " https://bitbucket.org/sigrimm/genga 5 " In computer arithmetic, a + b + c ≠ a + c + b because storage space for each number is finite, such that round-off errors will differ. 6 " In particular, spiral density waves launched from multiple massive bodies will locally modify the hydrodynamic drag, mutually interact, and affect the gravitational potential of the disc. Neither of these effects is captured in an analytic model. 7 " https://cheleb.net/astro/sp15 8 " Following equations (70) and (49) in Tanaka et al. (2002) and Tanaka & Ward (2004), respectively, yields – for a 1 MEarth planet at 1 au in gas disc depleted by a factor of 104 – a migration time-scale |$\dot{a}/a \approx 10^9$| years, but a damping time-scale of |$\dot{e}/e \approx 10^7$| years, which is well-within the bounds for relevance found by Kominami & Ida (2002). 9 " In massive discs, only three instead of four equal-mass collisions are required to reach the cutoff mass for classification as embryo. To avoid such numerical effects, future simulations that vary the total disc mass should keep a constant mass resolution. 10 " Profiles of surface density and radial mass distribution are related, but different methods of describing how mass is distributed. The radial mass distribution is a more general method relying on a non-parametric fit (generated by dropping Gaussians of adaptive width, summing them up, and then normalizing the integral) whereas the surface density requires more manual intervention. Although a more general method is generally preferred, much of the literature is based on surface densities. 11 " In fact, the balance between dynamical excitement from a passing secular resonance and damping by hydrodynamic and gravitational interactions determines whether planetesimals are swept along in the first place (Nagasawa et al. 2005). 12 " If planetesimals have too large eccentricities, the damping time-scale from interaction with the gas disc is longer than the time it takes for the secular resonance to sweep past the planetesimals. They cannot be trapped in the resonance (Nagasawa et al. 2005). 13 " For videos and time-sliced figures of the planetesimal dynamics of all simulations, we refer the reader to the supplementary website. 14 " In contrast to our simulations, Kokubo et al. (2006) and R09 do not evolve identical ICs, but rather redraw different ICs from the same underlying distribution. We are uncertain what is done in Izidoro et al. (2014a). 15 " Levison & Agnor (2003) and Raymond et al. (2004) also report stochastic variations across simulations, but do not report the RMC or AMD due to their focus on a different aspect. 16 " Of course, we have been wholly ignorant about the fact that giant planets may not yet be present (or massive enough) on their orbits as we start our runs. Future work will likely improve upon this and cheerfully point out this shortcoming. 17 " We adopt working definitions of terrestrial and giant planets as planets with masses M < MNeptune ∼ 17.15 MEarth and M > 50 MEarth. 18 " http://www.exoplanet.eu 19 " We consider the gaseous disc to extend inwards to 0.1 au, which corresponds to the dust sublimation radius (Dullemond et al. 2007). Of course, that does not necessarily mean that the gas disc extends to the same inner radius. To within a factor of 2, observations indicate the inner edge of the gaseous disc to coincide with the Solar corotation radius (Najita et al. 2007), which is at 0.17 au for the present-day Sun. 20 " Within the first Myr, the 90 percentile mass of embryos located within 0.2 au rapidly grow to ∼0.25 MEarth, where growth stalls out. For embryos beyond 0.2 au, the 90 percentile mass reaches this threshold between 4 and 5 Myr, but these embryos keep accreting afterwards. 21 " By ∼7.7 Myr, 95 per cent of material interior to 0.2 au originates from within 1.35 au, whereas embryos in the region 0.2 < a < 1.0 au source over 95 per cent of their material from beyond 1 au. REFERENCES Agnor C. B. , Lin D. N. C., 2012 , ApJ , 745 , 143 Crossref Search ADS Allison R. J. , Goodwin S. P., Parker R. J., Portegies Zwart S. F., de Grijs R., 2010 , MNRAS , 407 , 1098 Crossref Search ADS Andrews S. M. , Wilner D. J., Hughes A. M., Qi C., Dullemond C. P., 2009 , ApJ , 700 , 1502 Crossref Search ADS Andrews S. M. , Wilner D. J., Hughes A. M., Qi C., Dullemond C. P., 2010 , ApJ , 723 , 1241 Crossref Search ADS Birnstiel T. , Klahr H., Ercolano B., 2012 , A&A , 539 , A148 Crossref Search ADS Bitsch B. , Crida A., Morbidelli A., Kley W., Dobbs-Dixon I., 2013 , A&A , 549 , A124 Crossref Search ADS Bitsch B. , Morbidelli A., Lega E., Crida A., 2014a , A&A , 564 , A135 Crossref Search ADS Bitsch B. , Morbidelli A., Lega E., Kretke K., Crida A., 2014b , A&A , 570 , A75 Crossref Search ADS Brasser R. , Morbidelli A., Gomes R., Tsiganis K., Levison H. F., 2009 , A&A , 507 , 1053 Crossref Search ADS Brasser R. , Walsh K. J., Nesvorný D., 2013 , MNRAS , 433 , 3417 Crossref Search ADS Brouwer D. , Clemence G. M., 1961 , Methods of Celestial Mechanics . Academic Press , New York Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Chambers J. E. , 1999 , MNRAS , 304 , 793 Crossref Search ADS Chambers J. E. , 2001 , Icarus , 152 , 205 Crossref Search ADS Chambers J. E. , Cassen P., 2002 , Meteoritics Planet. Sci. , 37 , 1523 Crossref Search ADS Chambers J. E. , Wetherill G. W., 1998 , Icarus , 136 , 304 Crossref Search ADS Chambers J. E. , Wetherill G. W., Boss A. P., 1996 , Icarus , 119 , 261 Crossref Search ADS Chiang E. , Laughlin G., 2013 , MNRAS , 431 , 3444 Crossref Search ADS Chirikov B. V. , 1979 , Phys. Rep. , 52 , 263 Crossref Search ADS Cossou C. , Raymond S. N., Hersant F., Pierens A., 2014 , A&A , 569 , A56 Crossref Search ADS Cresswell P. , Nelson R. P., 2008 , A&A , 482 , 677 Crossref Search ADS Dullemond C. P. , Hollenbach D., Kamp I., D’Alessio P., 2007 , in Reipurth B., Jewitt D., Keil K., eds, Protostars and Planets V . Univ. Arizona Press , Tuscon , p. 555 Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Duncan M. J. , Levison H. F., 1997 , Science , 276 , 1670 Crossref Search ADS PubMed Duncan M. J. , Quinn T., 1993 , ARA&A , 31 , 265 Crossref Search ADS Duncan M. J. , Levison H. F., Lee M. H., 1998 , AJ , 116 , 2067 Crossref Search ADS Evans N. W. , Tabachnik S., 1999 , Nature , 399 , 41 Crossref Search ADS Everhart E. , 1973 , AJ , 78 , 329 Crossref Search ADS Ford E. B. , Rasio F. A., 2008 , ApJ , 686 , 621 Crossref Search ADS Franklin F. , Lecar M., Wiesel W., 1984 , in Greenberg R., Brahic A., eds, IAU Colloq. 75: Planetary Rings . Univ. Arizona Press , Tuscon , p. 562 Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Franklin F. , Lecar M., Soper P., 1989 , Icarus , 79 , 223 Crossref Search ADS Gladman B. , 1993 , Icarus , 106 , 247 Crossref Search ADS Gladman B. , Duncan M., 1990 , AJ , 100 , 1680 Crossref Search ADS Gladman B. J. et al. , 1997 , Science , 277 , 197 Crossref Search ADS Goodman J. , Heggie D. C., Hut P., 1993 , ApJ , 415 , 715 Crossref Search ADS Grazier K. R. , Newman W. I., Kaula W. M., Hyman J. M., 1999a , Icarus , 140 , 341 Crossref Search ADS Grazier K. R. , Newman W. I., Varadi F., Kaula W. M., Hyman J. M., 1999b , Icarus , 140 , 353 Crossref Search ADS Greenberg R. , Hartmann W. K., Chapman C. R., Wacker J. F., 1978 , Icarus , 35 , 1 Crossref Search ADS Greenzweig Y. , Lissauer J. J., 1990 , Icarus , 87 , 40 Crossref Search ADS Greenzweig Y. , Lissauer J. J., 1992 , Icarus , 100 , 440 Crossref Search ADS Grimm S. L. , Stadel J. G., 2014 , ApJ , 796 , 23 Crossref Search ADS Haisch K. E. Jr, Lada E. A., Lada C. J., 2001 , ApJ , 553 , L153 Crossref Search ADS Hansen B. M. S. , 2009 , ApJ , 703 , 1131 Crossref Search ADS Hayashi C. , 1981 , Progress Theor. Phys. Suppl. , 70 , 35 Crossref Search ADS Hayes W. B. , 2008 , MNRAS , 386 , 295 Crossref Search ADS Heppenheimer T. A. , 1980 , Icarus , 41 , 76 Crossref Search ADS Holman M. J. , 1995 , in Kinoshita H., Nakai H., eds, Proc. 27th Symposium on Celestial Mechanics , p. 116 Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Holman M. J. , Wisdom J., 1993 , AJ , 105 , 1987 Crossref Search ADS Hughes A. M. , Wilner D. J., Qi C., Hogerheijde M. R., 2008 , ApJ , 678 , 1119 Crossref Search ADS Hut P. , Heggie D. C., 2002 , J. Stat. Phys. , 109 , 1017 Crossref Search ADS Ida S. , 1990 , Icarus , 88 , 129 Crossref Search ADS Ida S. , Makino J., 1992a , Icarus , 96 , 107 Crossref Search ADS Ida S. , Makino J., 1992b , Icarus , 98 , 28 Crossref Search ADS Isella A. , Carpenter J. M., Sargent A. I., 2009 , ApJ , 701 , 260 Crossref Search ADS Ito T. , Tanikawa K., 2002 , MNRAS , 336 , 483 Crossref Search ADS Izidoro A. , Haghighipour N., Winter O. C., Tsuchida M., 2014a , ApJ , 782 , 31 Crossref Search ADS Izidoro A. , Morbidelli A., Raymond S. N., 2014b , ApJ , 794 , 11 Crossref Search ADS Jin L. , Arnett W. D., Sui N., Wang X., 2008 , ApJ , 674 , L105 Crossref Search ADS Johansen A. , Youdin A., Klahr H., 2009 , ApJ , 697 , 1269 Crossref Search ADS Kandrup H. E. , Smith H., Jr, 1991 , ApJ , 374 , 255 Crossref Search ADS Kokubo E. , Ida S., 1996 , Icarus , 123 , 180 Crossref Search ADS Kokubo E. , Ida S., 1998 , Icarus , 131 , 171 Crossref Search ADS Kokubo E. , Ida S., 2000 , Icarus , 143 , 15 Crossref Search ADS Kokubo E. , Ida S., 2012 , Prog. Theor. Exp. Phys. , 2012 , 010000 Crossref Search ADS Kokubo E. , Kominami J., Ida S., 2006 , ApJ , 642 , 1131 Crossref Search ADS Kominami J. , Ida S., 2002 , Icarus , 157 , 43 Crossref Search ADS Kominami J. , Ida S., 2004 , Icarus , 167 , 231 Crossref Search ADS Kuchner M. J. , 2004 , ApJ , 612 , 1147 Crossref Search ADS Laskar J. , 1989 , Nature , 338 , 237 Crossref Search ADS Laskar J. , 1990 , Icarus , 88 , 266 Crossref Search ADS Laskar J. , 1994 , A&A , 287 , L9 Laskar J. , 1997 , A&A , 317 , L75 Laskar J. , 2008 , Icarus , 196 , 1 Crossref Search ADS Laskar J. , Gastineau M., 2009 , Nature , 459 , 817 Crossref Search ADS PubMed Laskar J. , Quinn T., Tremaine S., 1992 , Icarus , 95 , 148 Crossref Search ADS Lecar M. , Franklin F. A., 1973 , Icarus , 20 , 422 Crossref Search ADS Lecar M. , Franklin F., 1997 , Icarus , 129 , 134 Crossref Search ADS Lecar M. , Franklin F. A., Holman M. J., Murray N. J., 2001 , ARA&A , 39 , 581 Crossref Search ADS Levison H. F. , Agnor C., 2003 , AJ , 125 , 2692 Crossref Search ADS Levison H. F. , Duncan M. J., 1993 , ApJ , 406 , L35 Crossref Search ADS Levison H. F. , Duncan M. J., 1997 , Icarus , 127 , 13 Crossref Search ADS Levison H. F. , Morbidelli A., Tsiganis K., Nesvorný D., Gomes R., 2011 , AJ , 142 , 152 Crossref Search ADS Lissauer J. J. et al. , 2011 , ApJS , 197 , 8 Crossref Search ADS McNeil D. S. , Nelson R. P., 2010 , MNRAS , 401 , 1691 Crossref Search ADS Mamajek E. E. , 2009 , in Usuda T., Tamura M., Ishii M., eds, AIP Conf. Ser. Vol. 1158, Exoplanets and Disks: Their Formation and Diversity . Am. Inst. Phys. , New York , p. 3 Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Masset F. S. , D’Angelo G., Kley W., 2006 , ApJ , 652 , 730 Crossref Search ADS Mikkola S. , Innanen K., 1995 , MNRAS , 277 , 497 Miller R. H. , 1964 , ApJ , 140 , 250 Crossref Search ADS Moons M. , Morbidelli A., 1995 , Icarus , 114 , 33 Crossref Search ADS Moons M. , Morbidelli A., Migliorini F., 1998 , Icarus , 135 , 458 Crossref Search ADS Morbidelli A. , Crida A., 2007 , Icarus , 191 , 158 Crossref Search ADS Morbidelli A. , Moons M., 1993 , Icarus , 102 , 316 Crossref Search ADS Morbidelli A. , Levison H. F., Tsiganis K., Gomes R., 2005 , Nature , 435 , 462 Crossref Search ADS PubMed Morbidelli A. , Tsiganis K., Crida A., Levison H. F., Gomes R., 2007 , AJ , 134 , 1790 Crossref Search ADS Morbidelli A. , Lunine J. I., O’Brien D. P., Raymond S. N., Walsh K. J., 2012 , Ann. Rev. Earth Planet. Sci. , 40 , 251 Crossref Search ADS Morishima R. , Schmidt M. W., Stadel J., Moore B., 2008 , ApJ , 685 , 1247 Crossref Search ADS Morishima R. , Stadel J., Moore B., 2010 , Icarus , 207 , 517 Crossref Search ADS Murray C. D. , Dermott S. F., 1999 , Solar System Dynamics . Cambridge Univ. Press , Cambridge Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Murray N. , Holman M., 1997 , AJ , 114 , 1246 Crossref Search ADS Nagasawa M. , Tanaka H., Ida S., 2000 , AJ , 119 , 1480 Crossref Search ADS Nagasawa M. , Ida S., Tanaka H., 2001 , Earth Planets Space , 53 , 1085 Crossref Search ADS Nagasawa M. , Lin D. N. C., Thommes E., 2005 , ApJ , 635 , 578 Crossref Search ADS Najita J. R. , Carr J. S., Glassgold A. E., Valenti J. A., 2007 , in Reipurth B., Jewitt D., Keil K., eds, Protostars and Planets V . Univ. Arizona Press , Tuscon , p. 507 Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Nesvorný D. , Morbidelli A., 1998 , AJ , 116 , 3029 Crossref Search ADS O’Brien D. P. , Morbidelli A., Levison H. F., 2006 , Icarus , 184 , 39 Crossref Search ADS Ogihara M. , Ida S., 2009 , ApJ , 699 , 824 Crossref Search ADS Ogihara M. , Kobayashi H., Inutsuka S.-i., 2014 , ApJ , 787 , 172 Crossref Search ADS Paardekooper S.-J. , Baruteau C., Crida A., Kley W., 2010 , MNRAS , 401 , 1950 Crossref Search ADS Paardekooper S.-J. , Baruteau C., Kley W., 2011 , MNRAS , 410 , 293 Crossref Search ADS Parker R. J. , Goodwin S. P., 2012 , MNRAS , 424 , 272 Crossref Search ADS Pfalzner S. , Steinhausen M., Menten K., 2014 , ApJ , 793 , L34 Crossref Search ADS Quinn T. R. , Tremaine S., Duncan M., 1991 , AJ , 101 , 2287 Crossref Search ADS Raymond S. N. , Cossou C., 2014 , MNRAS , 440 , L11 Crossref Search ADS Raymond S. N. , Morbidelli A., 2014 , in Knežević Z., Lemaitre A., eds, Proc. IAU Symp. 310, Complex Planetary Systems . Cambridge Univ. Press , Cambridge , p. 194 Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Raymond S. N. , Quinn T., Lunine J. I., 2004 , Icarus , 168 , 1 Crossref Search ADS Raymond S. N. , Quinn T., Lunine J. I., 2005a , Icarus , 177 , 256 Crossref Search ADS Raymond S. N. , Quinn T., Lunine J. I., 2005b , ApJ , 632 , 670 Crossref Search ADS Raymond S. N. , Quinn T., Lunine J. I., 2006a , Icarus , 183 , 265 Crossref Search ADS Raymond S. N. , Barnes R., Kaib N. A., 2006b , ApJ , 644 , 1223 Crossref Search ADS Raymond S. N. , Scalo J., Meadows V. S., 2007 , ApJ , 669 , 606 Crossref Search ADS Raymond S. N. , Barnes R., Mandell A. M., 2008 , MNRAS , 384 , 663 Crossref Search ADS Raymond S. N. , O’Brien D. P., Morbidelli A., Kaib N. A., 2009 , Icarus , 203 , 644 (R09) Crossref Search ADS Raymond S. N. , Armitage P. J., Gorelick N., 2010 , ApJ , 711 , 772 Crossref Search ADS Raymond S. N. , Kokubo E., Morbidelli A., Morishima R., Walsh K. J., 2014 , in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI . Univ. Arizona Press , Tuscon , p. 595 Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Richardson D. C. , Quinn T., Stadel J., Lake G., 2000 , Icarus , 143 , 45 Crossref Search ADS Safronov V. S. , Zvjagina E. V., 1969 , Icarus , 10 , 109 Crossref Search ADS Schlichting H. E. , 2014 , ApJ , 795 , L15 Crossref Search ADS Schneider J. , Dedieu C., Le Sidaner P., Savalle R., Zolotukhin I., 2011 , A&A , 532 , A79 Crossref Search ADS Stadel J. G. , 2001 , PhD thesis , Univ. Washington Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Sussman G. J. , Wisdom J., 1988 , Science , 241 , 433 Crossref Search ADS PubMed Sussman G. J. , Wisdom J., 1992 , Science , 257 , 56 Crossref Search ADS PubMed Tanaka H. , Ward W. R., 2004 , ApJ , 602 , 388 Crossref Search ADS Tanaka H. , Takeuchi T., Ward W. R., 2002 , ApJ , 565 , 1257 Crossref Search ADS Terquem C. , Papaloizou J. C. B., 2007 , ApJ , 654 , 1110 Crossref Search ADS Thommes E. , Nagasawa M., Lin D. N. C., 2008 , ApJ , 676 , 728 Crossref Search ADS Torbett M. V. , 1989 , AJ , 98 , 1477 Crossref Search ADS Torbett M. V. , Smoluchowski R., 1990 , Nature , 345 , 49 Crossref Search ADS Tsiganis K. , Gomes R., Morbidelli A., Levison H. F., 2005 , Nature , 435 , 459 Crossref Search ADS PubMed Valluri M. , Merritt D., 2000 , in Gurzadyan V. G., Ruffini R., eds, The Chaotic Universe . World Scientific Press , Singapore , p. 229 Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Walsh K. J. , Morbidelli A., Raymond S. N., O’Brien D. P., Mandell A. M., 2011 , Nature , 475 , 206 Crossref Search ADS PubMed Ward W. R. , 1981 , Icarus , 47 , 234 Crossref Search ADS Weidenschilling S. J. , 1977 , Ap&SS , 51 , 153 Crossref Search ADS Weidenschilling S. J. , Spaute D., Davis D. R., Marzari F., Ohtsuki K., 1997 , Icarus , 128 , 429 Crossref Search ADS Wetherill G. W. , Stewart G. R., 1989 , Icarus , 77 , 330 Crossref Search ADS Williams J. P. , Cieza L. A., 2011 , ARA&A , 49 , 67 Crossref Search ADS Wisdom J. , 1980 , AJ , 85 , 1122 Crossref Search ADS Wisdom J. , 1982 , AJ , 87 , 577 Crossref Search ADS Wisdom J. , 1983 , Icarus , 56 , 51 Crossref Search ADS Wisdom J. , 1985 , Icarus , 63 , 272 Crossref Search ADS Yoshinaga K. , Kokubo E., Makino J., 1999 , Icarus , 139 , 328 Crossref Search ADS Youdin A. N. , Chiang E. I., 2004 , ApJ , 601 , 1109 Crossref Search ADS © 2016 The Authors Published by Oxford University Press on behalf of the Royal Astronomical Society
TI - Stochasticity and predictability in terrestrial planet formation
JF - Monthly Notices of the Royal Astronomical Society
DO - 10.1093/mnras/stw2856
DA - 2017-02-21
UR - https://www.deepdyve.com/lp/oxford-university-press/stochasticity-and-predictability-in-terrestrial-planet-formation-SApmCEYoPc
SP - 2170
VL - 465
IS - 2
DP - DeepDyve
ER -