A Statistics of Our Simulations

Growing the terrestrial planets from the gradual accumulation of sub-meter sized objects


Building the terrestrial planets has been a challenge for planet formation models. In particular, classical theories have been unable to reproduce the small mass of Mars and instead predict that a planet near 1.5 AU should roughly be the same mass as the Earth. Recently, a new model called Viscous Stirred Pebble Accretion (VSPA) has been developed that can explain the formation of the gas giants. This model envisions that the cores of the giant planets formed from 100 to 1000 km bodies that directly accreted a population of pebbles — sub-meter sized objects that slowly grew in the protoplanetary disk. Here we apply this model to the terrestrial planet region and find that it can reproduce the basic structure of the inner Solar System, including a small Mars and a low-mass asteroid belt. Our models show that for an initial population of planetesimals with sizes similar to those of the main belt asteroids, VSPA becomes inefficient beyond 1.5 AU. As a result, Mars’s growth is stunted and nothing large in the asteroid belt can accumulate.

Solar System Formation — planetary dynamics

www.pnas.org/cgi/doi/ \copyrightyear2015 \issuedateIssue Date \volumeVolume \issuenumberIssue Number


Submitted to Proceedings of the National Academy of Sciences of the United States of America


The fact that Mars is so much smaller than both the Earth and Venus has been a long standing puzzle of terrestrial planet formation. Here we show that a new mode of planet formation known as ‘Viscous Stirred Pebble Accretion’, which as recently been shown to produce the giant planets, also naturally explains the small size of Mars and the low mass of the asteroid belt. Thus there is a unified model that can be used to explain the all of the basic properties of the our Solar System.


VSPA, viscously stirred pebble accretion


Classical models of terrestrial planet formation have a problem, the same models that produce reasonable Earth and Venus analogs tend to produce Mars analogs that are far too large [1]. The only existing proposed explanations for the small mass of Mars based on classical modes of growth require a severe depletion solids beyond AU [2] involving either not-well-understood nebular processes [3] or a complicated and dramatic migration of the giant planets [4] to solve this problem. Recently, however, it has been shown that a new mode of planet formation known as Viscous Stirred Pebble Accretion (VSPA) can successfully explain the formation of the giant planets [5, 6]. Here it is our hypothesis that Mars’s mass may simply be another manifestation of VSPA. To understand how, we need to describe the process.

1 Review of Pebble Accretion

After the formation of the protoplanetary disk, dust particles, which are suspended in the gas, slowly collide and grow because of electrostatic forces. Once particles become large enough so that their Stokes numbers (, where is the stopping time due to aerodynamic drag and is the orbital frequency) are between and 1, depending on the model, these so-called ‘pebbles’ can be concentrated by aerodynamic processes [7, 8, 9, 10]. Under the appropriate physical conditions (which might not have been satisfied everywhere in the disk), these concentrations become dense enough that they become gravitationally unstable and thus collapse to form planetesimals [11] with radii between and km [10, 12]. This process can occur very quickly — on the order of the local orbital period.

Recent research shows that planetesimals embedded in a population of pebbles can grow rapidly via a newly discovered accretion mechanism that is aided by aerodynamic drag on the pebbles themselves [13, 5, 14]. In particular, if a pebble’s aerodynamic stopping time is less than or comparable to the time for it to encounter a growing body (hereafter known as an embryo) then it is decelerated with respect to the embryo and becomes gravitationally bound. After capture, the pebble spirals towards the embryo due to aerodynamic drag and is accreted. The accretional cross section for this situation is


where , is the mass of the embryo, and is the relative velocity between the pebble and embryo [13]. For the growing planets, can be orders of magnitude larger than the physical cross section alone. Full -body simulations [6] show that as long as pebbles form continuously over a long enough time period such that embryos have time to gravitationally stir each other, this process can form the observed gas giant planets before the gas disk dissipates.

Figure 1: The value of as a function of heliocentric distance () and embryo radius () in our fiducial protoplanetary disk, assuming that the embryos have circular (Keplerian) orbits and the pebbles are on orbits as determined by aerodynamic drag. The top and bottom panels are for and , respectively. As in general increases with heliocentric distance for an object of a given size, objects that can grow in the inner regions cannot grow farther from the Sun. For reference, the white horizontal line in the top panel corresponds to the radius of (1) Ceres.

Our hypothesis that this process can also explain Mars’s small size and the low mass of the asteroid belt is based on the term in Eq (1), which says that pebble accretion becomes exponentially less efficient for small embryos because the encounter times for these objects () becomes short compared to . As a result, aerodynamic drag does not have time to change the trajectory of the pebbles, so they are unlikely to be accreted. Eq (1) therefore predicts a sharp cutoff between small embryos, which cannot grow, and larger objects, which can. In addition, because is a function of location in the disk, this cutoff also varies with location. Fig 1 shows the value of in our fiducial disk (which is described in the Methods section) for two values of that are consistent with the requirements of the two competing models of planetesimal formation. In particular, the top panel employs pebbles that are required for the so-called streaming instability [8, 10, 15], while the bottom panel uses the pebbles needed by the turbulent concentration models of Ref. [7, 9]. As in general increases with heliocentric distance for an embryo of a given size, embryos that can grow in the inner regions cannot grow farther from the Sun. For example, if pebbles have , an object initially the size of Ceres will grow at AU (where ), but not at AU (where becomes large). This argument implies that, initially, all planetesimals could have been the size of currently observed main belt asteroids; those bodies in the asteroid belt did not grow appreciably, while those at AU did. Therefore, we postulate that Mars’s small mass and the lack of planets in the asteroid belt might be the result of this cutoff. It is important to note that this figure just shows and does not represent how the entire pebble accretion process will behave. In order to ascertain that, we must turn to numerical calculations. The remainder of this work presents such simulations.

2 Methods

To test this hypothesis, we performed a series of -body calculations of terrestrial planet formation starting with a population of planetesimals embedded in our fiducial gas disk, adopted from Ref [6]. We assume a flaring gas disk with a surface density of [16] and a scale height [17], where is heliocentric distance in AU. Here we set initially to 9000 g/cm, which is roughly 5 times that of the so-called minimum mass solar nebula [18]. In this work we set the density of the planetesimals and pebbles, , to 3 g cm. Additionally, our disk is assumed to be turbulent with [19] (although see below).

We allowed the gas surface density to decrease exponentially with a timescale of Myr, which is motivated by observations [20]. The disk has solar composition so that the solid-to-gas ratio is 0.005 in the terrestrial planet region [21]. We convert a fraction of the solids into planetesimals at the beginning of the simulations. We draw our planetesimals from a distribution of radii, , of the form such that is between and . The values of and are assumed to be independent of semi-major axis. We set our fiducial value of to the radius of Ceres, km, because we expect little growth in the asteroid belt and we need to produce this object. However, we do vary from 100 to km to test the sensitivity of our results to this value. Because we are interested in building the terrestrial planets and using the asteroid belt as a constraint, we study the growth of planetesimals spread from 0.7 to AU (the presumed location of the snow line; [18]). As is typical for this type of simulation [1], we do not treat planet formation in the Mercury forming region in order to save computer time.

Figure 2: Panels (a) to (c) show the formation of a system of terrestrial planets. Each panel shows the mass of the planetesimals and planetary embryos as a function of semi-major axis, while color indicates their eccentricity. In addition, the ‘error bars’, which are only shown for objects larger than in order to decrease clutter, indicate the range of heliocentric distance that an object travels. The yellow box shows the region populated by planetesimals at . Panel (b) shows the distribution of planetary embryos after pebble accretion but before the dynamical instability. Panel (c) panel shows the final system. The panel (d) shows the fiducial giant planet system constructed in Ref. [6].

The remaining solids (assumed to be dust) are slowly converted into pebbles with a fixed initial , which is a free parameter in our simulations. Following Ref [6], we utilize a simple prescription to convert dust into pebbles over time that assumes that pebbles form at a rate proportional to the instantaneous dust mass, correcting for dust lost as the gas disk evolves and as pebbles form. The functional form of pebble production can be found in Eq. 9 of the Methods section of Ref [6]. We scale the function such that the median production timescale is roughly . We assume that all of the pebbles are produced in 2 Myr. For simplicity, we assume that pebbles are randomly created throughout the disk according to the surface density.

We are justified in pursuing the long timescales for pebble formation for the following reasons. While models of dust coagulation (Ref [22], for example) predict that pebbles should grow on timescales on the order of 100-1000 orbital periods, this result is observationally problematic. Millimeter and even centimeter sized particles, which should have been lost rapidly, are observed in disks of a range of ages (e.g Ref.[23]). While it is possible that the drift of these pebbles could be slowed by variations in the disk structure [24], these trapping models need large, as of yet unobserved, variations in the disk structure. An simpler alternative is that pebbles are continuously formed. Indeed, models in which pebbles slowly form from dust and then are lost due to drift matches some features of observed disks [33]. Therefore, we will assume an initial planetesimal population along with pebbles that are steadily produced by the disk over its lifetime.

We also assume that the pebbles involved in terrestrial planet formation formed within AU. By having a cutoff at this location, we are assuming that material drifting in from the outer Solar System is unable to penetrate the snow-line, presumably due to sublimation. But no matter the mechanism, this assumption is required because solids from the outer Solar System are too carbon rich to have contributed more than a few percent of the mass of the terrestrial planets (see Ref. [26]). Also, carbon can not be removed without heating the material to above 500K [27], a temperature not reached in the midplane until well within 1 AU in reasonable disk models.

The values of present during VSPA are dependent on which planetesimal formation model we assume and range from to [7, 8, 9, 10]. Ideally, here we would prefer to study pebble sizes that cover the complete range required by both models. However, the CPU time required to perform our calculations increases drastically as decreases because of two effects. First, the timestep required by our code scales with the pebble’s aerodynamic stopping time. Thus, a smaller requires a smaller timestep. In addition, pebbles with small ’s have slower radial drift velocities than their larger siblings, and this spend more time in the calculation. As a result, at any time there are more objects present in the simulation that the code needs to deal with. This significantly increases the required CPU time per timestep. Therefore, to keep the calculations tractable, we require to be larger than . This issue will be addressed again in the Discussion Section below.

Each system is evolved for Myr using a Lagrangian code, LIPAD [28]. LIPAD is the first particle-based (i.e. Lagrangian) code that can follow the collisional/accretional/dynamical evolution of a large number of sub-kilometer objects through the entire growth process to become planets. It is built on top of the symplectic -body algorithm known as SyMBA [29]. In order to handle the very large number of sub-kilometer objects required by these simulations, we introduce the concept of a tracer particle. Each tracer represents a large number of small bodies with roughly the same orbit and size, and is characterized by three numbers: the physical radius, the bulk density, and the total mass of the disk particles represented by the tracer. LIPAD employs statistical algorithms for viscous stirring, dynamical friction, and collisional damping among the tracers. The tracers mainly dynamically interact with the larger planetary mass objects via the normal -body routines, which naturally follow changes in the trajectory of tracers due to the gravitational effects of the planets and vice versa. When a body is determined to have suffered an impact, it is assigned a new radius according to the probabilistic outcome of the collision based on a fragmentation law for basalt by Ref. [30]. In this way, the conglomeration of tracers and full -body objects represent the size distribution of the evolving population. LIPAD is therefore unique in its ability to accurately handle the mixing and redistribution of material due to gravitational encounters, including planetesimal-driven migration and resonant trapping, while also following the fragmentation and growth of bodies. An extensive suite of tests of LIPAD can be found in Ref. [28]. For the calculations described here, we will use a version of LIPAD that has been modified to handle the particular needs of pebble accretion [31].

The calculations are performed in three stages. For the first Myr, the terrestrial planet region is evolved in isolation as pebbles continually form and drift inward. At the end of this first stage, all pebbles have either been accreted or lost, and thus no more mass will be added to the system. Because we are interested in constructing systems similar to the Solar System, we only continue simulations with the appropriate amount of material (between and ) inside AU. If this criterion is met, the simulation is cloned 6 times by adding a random number between to AU to each component of the position vector of each body. For this second stage, we also add Jupiter and Saturn in orbits consistent with their pre-migrated configuration [32]. In particular, they are placed in their mutual 3:2 mean motion resonance with Jupiter at AU and Saturn at AU. The evolution is then followed until Myr. For the final stage, Jupiter and Saturn are moved to their current orbits and the system was integrated for an additional Myr.

3 Simulation Results

In total, we performed 28 simulations to at least Myr, varying between 0.004 and 0.01, and between 100 to km. For comparison, note that asteroid (1) Ceres has a radius of km. Our small values for were driven by the fact that we wanted to create asteroid belts as close to the observed system as possible assuming that VSPA is not effective there. Unfortunately, however, the smaller we force the more tracers that we needed to resolve the system and thus the more CPU time required by the simulation. Our compromise was to choose so that the initial mass between 2.2 and AU is roughly between 20 and 50 times the mass currently observed there. These small values of imply that there was at most of planetesimals between 0.7 and AU; most of the mass in the final planets comes from accreting pebbles.

We followed 9 of the 28 simulations to completion, making a total of 54 systems including the clones. See the Supporting Information for a description of the statistics of our runs.

The evolution of a system that produced a reasonable Solar System analog is shown in Fig 2. For this simulation, , is initially between between 200 and km, and . Growth occurs first and fastest closest to the Sun. This happens for two reasons. As described above, is function of semi-major axis. In addition, although we create pebbles randomly throughout our computational domain following the surface density of the gas disk, these objects quickly spiral toward the Sun due to the effects of aerodynamic drag. As a result, a planetesimal encounters pebbles that were created outward of its semi-major axis and thus the planetesimals closest to the Sun encounter more pebbles.

The era of pebble accretion ends when all the pebbles have been generated and have either impacted the Sun or been accreted by an embryo. At this time (Fig 2a) there is a series of 23 embryos with masses greater than on quasi-stable, nearly circular orbits. There is a direct correlation between mass and semi-major axis at this time, with no object larger than Mars beyond AU and none larger than the Moon beyond AU. The largest object in the system is . This correlation leads to very little mass beyond AU. As expected from Fig 1, there was little growth beyond AU.

This system remains stable until Myr, at which time the orbits of its embryos cross and they accrete each other. This dynamical instability leads to the formation of two roughly Earth-mass objects at 0.7 and AU. It is during this instability that a 2.7 Mars-mass object is gravitationally scattered to AU by its larger siblings and is stabilized by gravitational interactions with smaller objects found there. This leads to the system shown in the Fig. 2c that contains analogs of the Earth, Venus, and Mars with roughly the correct masses and orbits. The basic evolution seen here, where the system first develops a series of small planets on nearly circular orbits that suffer an instability at a few tens of millions of years, is a common outcome in our simulations. Thus, this model predicts that the Solar System had an initial generation of terrestrial planets — consisting of a large number of small planets — that is now lost. Mars is likely a remnant of this early system. The late timing of the impact that formed the moon [41] resulted, in part, from this instability.

Figure 3: The finial distribution of our planets. a) The black dots show a compilation of the planets constructed during our simulations. In particular, we plot the fraction of total system mass in each planet as a function of its semi-major axis. Here we plot the ‘fraction of total system mass’ because the total mass of our systems varies from run to run due to the stochastic variations in the efficiency of pebble accretion. The ‘error bars’ indicate the range of heliocentric distance that a planet travels and thus are a function of eccentricity. The red squares indicate the real terrestrial planets in the Solar System. b) The cumulative mass distribution of Mars analogs. Planets to the left of the red line are considered good Mars analogs according to Ref. [35].

Although, due to the chaotic nature of planet formation, not all of our simulations produce good Solar System analogs, as shown in Fig. 3a. While interpreting this plot it must be noted that we did not make any attempt to uniformly cover parameter space. Each calculation took many weeks to perform and thus the survey of parameter space was limited and ad hoc. Thus, the observed distribution shown needs to be viewed with caution. However, we believe that some of the trends are robust. For example, it is common to produce reasonable Earth and Venus analogs. In addition, planets near AU are systematically smaller than their siblings interior to AU. No planets grow beyond AU. The objects seen in this region were scattered out during a dynamical instability. The majority of these objects are still on orbits that cross their larger neighbors and thus will eventually be removed. However, a few are on stable orbits. The anomalous and highly processed main belt asteroid (4) Vesta might be an object that was captured in this manner.

The natural question is whether VSPA produces good Mars analogs more frequently than the standard planet formation picture. Ref. [35] finds that the standard model produces a good Mars analog in 4% of the simulations with Jupiter and Saturn on circular orbits (their likely state at the time of terrestrial planet formation). They define a ‘Mars analog’ as the largest planet with a semi-major axis both in the range 1.25 – AU and outside of the Earth analog’s orbit. They define an ‘Earth analog’ to be the largest planet between 0.75 and AU and if there is no planet within this range, the Earth analog is the closest planet to AU. A ‘good’ Mars analog is defined to be an object that is smaller than or 11% of the total system mass. Keeping the above caveats about our parameter space coverage in mind, we plot the cumulative mass distribution of Mars analogs from our simulations in Fig. 3b. We find that 30% of our systems produce good analogs. This represents a significant improvement over the standard model. Moreover, Ref. [35] finds that the probability that the standard model produces systems with both a good Mars analog and a low mass asteroid belt is roughly 0.6%. By design, our systems always have a low mass asteroid belt.

4 Discussion

The above experiments show that Mars’s small mass is a natural outcome of the process of VSPA. Having said this, there are other issues that need to be considered.

4.1 Sizes of Pebbles

There are currently two independent models in the literature for planetesimal formation directly from pebbles that make very different predictions concerning the size of pebbles: the streaming instability and turbulent concentration. The streaming instability takes advantage of the fact that a clump of pebbles is less affected by aerodynamic drag, and thus moves at a different speed, than free-floating pebbles [8]. As a result, individual pebbles catch up with the clump, and the clump grows to the point where it becomes gravitationally bound and collapse to form a planetesimal [10]. This mechanism requires [15]. On the other hand, turbulent concentration scenarios suggest that objects with migrate to the outer edges of turbulent eddies due to centrifugal forces caused by the eddie’s rotation [7, 9]. There the pebbles form gravitationally bound clumps that collapse to create planetesimals. While we are agnostic on which of these are correct and would prefer to cover the whole range of values, as we described above, we cannot perform simulations with ’s smaller than because the required amount of CPU time makes these calculations impractical.

Unfortunately, there have been recent developments in our understanding of the coagulation of pebbles that calls into question whether objects with can grow in disks. In traditional models the growth of objects is limited by particle fragmentation [33]. With this limitation, particles should grow to given the known strengths of compacted silicate grains [34]. However, some coagulation simulations have identified another barrier to growth, known as the “bouncing barrier,” that halts growth of rocky pebbles at much small sizes [36, 37]. If this barrier is proven robust then pebble growth may halt at .

The fact that we cannot perform direct LIPAD calculations with does not imply that our basic model will not work in this regime. For example, the simple calculation shown in Fig 1 shows that, regardless of pebble size, there is a radial dependence on the efficiency of pebble accretion. The only difference is that the transition between efficient and inefficient growth occurs at much smaller planetesimal size. This suggests that we would be able to create the same general features seen in Fig 2c (a low mass Mars and asteroid belt) for a range of assumed pebble sizes as long as a decrease in pebble size is commensurate with a decrease in the maximum initial planetesimal size.

This trend is already apparent with the range of pebble size we can investigate with LIPAD (see the table in the Supporting Information). For example, when was , reasonable terrestrial analogs were found only when the largest planetesimals initially in the system were bigger than km. However when , simulations that initially contained planetesimals this big grew planets in the asteroid belt. For ’s of this value, systems that look like our Solar System occurred only when that the maximum initial planetesimal size was km. Indeed, recent work indicates that the size-distribution of the larger asteroids in the main belt could be explained by planetesimals of roughly this size accreting pebbles[42]. Therefore, we expect that an initial combination of pebbles and planetesimals of a few tens of kilometers would result in systems like our own.

Having said this, even if the bouncing barrier model is correct, it has some uncertainties and free parameters that might allow larger pebbles to grow. One of the most important being the assumed level of turbulence in the disk. In environments with very little turbulence () it has been demonstrated that small pebbles and dust can combine to form larger aggregates with [38].

In our simulations we chose , a value motivated by disk evolution timescales [20] with the assumption that uniform turbulence is driving angular momentum transport in protoplanetary disks. However, while angular momentum must be being transferred in these accretion disks, this does not necessarily imply that there is turbulence at the midplane where the pebbles and planetesimals are located. Indeed, recent models of angular momentum transport in disks suggest that turbulence may be limited to the upper and lower regions of the disk [39] or that disks may not be turbulent at all, and angular momentum may instead by transferred by magnetocentrifulgal disk winds [40]. Therefore we investigated how pebble accretion would behave in more quiescent disks that may be more consistent with our assumed pebble size. We find that the masses of the growing planetary embryos are independent of the assumed level of turbulence for ranging from to our fiducial value of (see Fig. 4). As a result, although we did not directly perform simulations in a low- environment where pebbles can form even if the bouncing barrier is important, our results are directly applicable there.

Figure 4: The cumulative mass distribution of objects at the end of pebble accretion (i.e at Myr) in three simulations that are identical except for . These runs have , km, and (our smallest value). Three values of where studied; (our fiducial value; black curve), (cyan curve), and (red curve). Note that all three curves are basically identical.

4.2 Implications for the Asteroid Belt

This model also has profound implications for both the location and history of the asteroid belt. One of the primary reasons why increases with heliocentric distance in Fig 1 is because we assume the disk flares (i.e the ratio between the gas disk scale height and the heliocentric distance increases with heliocentric distance). While it is observed that most disks are flaring in the outer regions [43], there are not direct observations constraining the disk shape in the terrestrial planet region. Indeed, if a disk is viscously heated and one were to assume that the opacity of the disk is constant with heliocentric distance, then the disk will be flat, with a constant aspect ratio (e.g [44]). However, full radiative transfer models of the structure of viscously heated accretion disks (e.g [45]) show flaring interior to the snow line (which is why we made the assumption we did). This occurs because the dust opacity is expected to increase with decreasing temperatures as water freezes out. In viscously heated disk this occurs interior to the mid-plane snow-line because the temperture decreases with height. The change in opacity causes the disk to flare. Thus, a region of increasing is likely associated with snow-lines, and VSPA might lead us to expect that asteroid belts are features that are naturally associated with regions interior to snow-lines.

Our VSPA calculations also produce asteroid belts are profoundly different than the standard model. In the standard model, the region between 2 and AU originally contained several Earth-masses of planetesimals —  of which were lost either through the dynamical effects of Jupiter and Saturn or by collisional grinding [46](the asteroid belt currently contains ). Since is always less than 0.1 in our simulations and little growth occurred beyond AU, there was never very much mass in the asteroid belt. Indeed, our final planetary systems have asteroid belts with masses between 0.5 and 155 of that currently observed in this region. The model in Fig. 2 has an asteroid belt 3.5 times more massive than the current belt. This is a reasonable value because we expect the asteroid belt to lose between roughly and of its mass due to chaotic diffusion, giant planet migration, and collisions over the age of the Solar System [4, 47].

As previously noted by Ref. [4], an early low-mass asteroid belt might explain several of its most puzzling characteristics. For example, a low-mass main belt is consistent with a limited degree of collisional evolution, as predicted by modeling work [48]. Similarly, the 530 km diameter differentiated asteroid (4) Vesta only has two 400–km diameter basins on its surface, and one of them is only 1 Gy old [49]. These data are a good match to a primordial main belt population that was never very massive. Finally, if Vesta formed as part of a primordial asteroid belt that contained hundreds of times the mass of the current population, then it would be statistically likely have hundreds of Vesta-like bodies as well. Even a short interval of collisional evolution would produce more basaltic fragments than can be accommodated by the existing collisionally and dynamically evolved main belt [50]. All of these observables could be explained by asteroid belts that never contained much mass, like those constructed in our models.

Additionally, we note that by varying the parameters in this model, particularly the initial planetesimal size and the pebble Stokes number, we produce a wide range of different planetary architectures. While we present parameters that produce systems similar to our own, the natural variability of this process leads us to speculate that VSPA might be able to explain the variety in the observed exoplanetary systems [51].

4.3 Unified Picture of the Solar System Formation

Finally for completeness, Fig 2d shows a giant planet system constructed within the same disk and using VSPA [6]. It shows two gas giant planets plus 3 ice giants. There is no growth beyond AU because is again large.

We note that we get growth in the giant planet region [6] and not in the asteroid belt for at least three reasons: 1) The initial planetesimals were probably larger in the outer Solar System (to be consistent with the larger sizes of object in the Kuiper Belt as compared to the asteroid belt). In particular, in this work we use a maximum planetesimal size in the terrestrial region based on Ceres, while in Ref. [6] we used a maximum planetesimal size based on Pluto in the giant planet region. 2) The pebble sizes were likely also bigger because because the pebbles are icy and thus stickier [52]. 3) the giant planets have access to a much larger reservoir of pebbles than the inner Solar System due to sublimation of ices at the snow line. In particular as we explained above, here we assume that pebbles from the outer Solar System do not contribute to the growth of the terrestrial planets. As a result, the terrestrial planets only have access to solids out to AU. In Ref. [6], we assumed that pebbles formed out to AU and thus the giant planets had access to the substantially larger amount of solids that were between 2.7 and AU.

The lower two panels in the figure, therefore show a consistent planetary system generated by VSPA. This model reproduces the basic structure of our entire planetary system — two roughly Earth mass objects between 0.5 and AU, a small Mars, a low mass asteroid belt, two gas giants, ice giants, and a primordial Kuiper belt that contains objects the mass of (134340) Pluto and (136199) Eris but did not form planets.

This research was supported by the NASA Solar System Exploration Research Virtual Institute. We would like to thank M. Duncan, S. Jacobson, M. Lambrechts, A. Morbidelli, and D. Nesvorný for useful discussions.

Appendix A Statistics of Our Simulations

Tables SS1 & SS3 list the simulations that we have completed in our investigation of terrestrial planet formation. All of these employed the following setup. We used a total surface density distribution of , where is the heliocentric distance in AU, and . The disk flares with a scale height of [17]. We assume a solid-to-gas ratio of 0.005 and a bulk density of 3 g/cm. Pebbles are slowly generated and follow the same spatial distribution as above. In particular, since pebble generation follows the evolution of the gas disk, the median time for which a pebble is generated is Myr.

There are three important free parameters that we varied. These are listed in the first three columns of the tables:

  1. The fraction of solids in the disk that are initially converted into planetesimals, . For the fiducial simulation, .

  2. The range of initial radii of the planetesimals. We draw the initial planetesimals from a distribution of radii, , of the form [53]. For our fiducial simulation km.

  3. The initial Stokes number of a pebble, . Note that as pebbles spiral toward the Sun, their Stokes number changes (generally decreases) because both the gas density and the orbital period vary. For our fiducial simulation .

The last four columns in the table present the basic characteristics of our systems. In particular, is the total mass found within AU at Myr. Recall that we only continued to integrate systems for which . We chose this range because it is slightly larger the total of the terrestrial planets () and thus allow for some subsequent loss. We created 6 clones of any run that satisfies the above criterion. The clone ID is listed in Column (5). Column (6) lists the total mass in the asteroid belt (here defined to be between 2.1 and AU) in terms of the currently observed asteroid belt mass (). Finally Columns (7) and (8) list the total number of planets with mass between 0.5 Mars-mass and () and the number with mass greater than (), respectively.

Table S1: Our completed simulations.
(1) (2) (3) (4) (5) (6) (7) (8)
(km) () ()
0.004 50-100km 0.02 3.3
0.004 50-100km 0.025 2.4
Clone 1 13 1 2
Clone 2 13 0 2
Clone 3 12 0 2
Clone 4 37 1 1
Clone 5 7 0 2
Clone 6 11 0 2
0.004 50-100km 0.03 1.6
0.004 50-100km 0.04 1.3
0.004 50-100km 0.06 1.1
0.004 100-300km 0.03 3.6
0.004 100-300km 0.06 3.1
0.004 100-300km 0.08 2.4
Clone 1 1.2 2 1
Clone 2 11 1 3
Clone 3 82 2 2
Clone 4 112 1 4
Clone 5 31 1 3
Clone 6 37 0 2
0.004 100-300km 0.1 1.7
0.004 200-450km 0.03 3.9
0.004 200-450km 0.06 3.5
0.004 200-450km 0.08 2.8
Clone 1 2 0 2
Clone 2 7 2 2
Clone 3 5 1 2
Clone 4 10 0 3
Clone 5 17 0 3
Clone 6 50 2 2
0.004 200-450km 0.1 2.1
Clone 1 81 0 2
Clone 2 12 1 3
Clone 3 72 2 2
Clone 4 9 3 1
Clone 5 18 0 2
Clone 6 6 2 2
Table S2: Our completed simulations (cont).
(1) (2) (3) (4) (5) (6) (7) (8)
(km) () ()
0.008 100-300km 0.03 5.7
0.008 100-300km 0.06 3.1
0.008 100-300km 0.08 2.4
Clone 1 82 3 1
Clone 2 32 2 2
Clone 3 110 4 1
Clone 4 31 2 2
Clone 5 17 3 1
Clone 6 19 3 1
0.008 100-300km 0.1 2
0.008 200-450km 0.03 5.9
0.008 200-450km 0.06 3.8
0.008 200-450km 0.1 2.7
Clone 1 182 3 2
Clone 2 12 2 2
Clone 3 15
Clone 4 15 1 2
Clone 5 224 2 2
Clone 6 113 3 2
0.008 200-450km 0.3 1.1
0.008 200-600km 0.1 2.6
Clone 1 57 1 2
Clone 2 208 3 2
Clone 3 16 2 1
Clone 4 42 2 1
Clone 5 22 3 0
Clone 6 17 3 2
0.01 100-300km 0.06 3.6
0.01 100-300km 0.08 3.1
0.01 100-300km 0.1 2.6
Clone 1 0.5 0 2
Clone 2 12 1 2
Clone 3 102 1 2
Clone 4 155 2 3
Clone 5 19 0 2
Clone 6 28 1 3
0.01 200-450km 0.03 6.3
0.01 200-450km 0.1 2.7
Clone 1 7 2 2
Clone 2 52 2 3
Clone 3 54 0 2
Clone 4 19 1 2
Clone 5 38 1 2
Clone 6 16 1 2
0.01 200-450km 0.3 1.2
Table S3: Our completed simulations (cont).


  1. To whom correspondence should be addressed. Email: hal@boulder.swri.edu


  1. Raymond SN, O’Brien DP, Morbidelli A, Kaib NA (2009) Building the terrestrial planets: Constrained accretion in the inner Solar System. Icarus 203:644–662.
  2. Hansen, BMS (2009) Formation of the Terrestrial Planets from a Narrow Annulus. Astrophys. J. 703:1131–1140.
  3. Izidoro A, Haghighipour N, Winter OC, Tsuchida M (2014) Terrestrial Planet Formation in a Protoplanetary Disk with a Local Mass Depletion: A Successful Scenario for the Formation of Mars. Astrophys. J. 782:31.
  4. Walsh KJ, Morbidelli A, Raymond SN, O’Brien DP, Mandell AM (2011) A low mass for Mars from Jupiter’s early gas-driven migration. Nature 475:206–209.
  5. Lambrechts M, Johansen A (2012) Rapid growth of gas-giant cores by pebble accretion. Astron. Astrophys. 544:A32.
  6. Levison HF, Kretke KA, Duncan M (2015) Growing the gas-giant planets by the gradual accumulation of pebbles. Nature 524:32–324.
  7. Cuzzi JN, Hogan RC, Paque JM, Dobrovolskis AR (2001) Size-selective Concentration of Chondrules and Other Small Particles in Protoplanetary Nebula Turbulence. Astrophys. J. 546:496–508.
  8. Youdin AN, Goodman J (2005) Streaming Instabilities in Protoplanetary Disks. Astrophys. J. 620:459–469.
  9. Cuzzi JN, Hogan RC, Shariff K (2008) Toward Planetesimals: Dense Chondrule Clumps in the Protoplanetary Nebula. Astrophys. J. 687:1432–1447.
  10. Johansen A et al. (2007) Rapid planetesimal formation in turbulent circumstellar disks. Nature 448:1022–1025.
  11. Goldreich P, Ward WR (1973) The Formation of Planetesimals. Astrophys. J. 183:1051–1062.
  12. Youdin AN (2011) On the Formation of Planetesimals Via Secular Gravitational Instabilities with Turbulent Stirring. Astrophys. J. 731:99.
  13. Ormel CW, Klahr HH (2010) The effect of gas drag on the growth of protoplanets. Analytical expressions for the accretion of small bodies in laminar disks. Astron. Astrophys. 520:A43.
  14. Morbidelli A, Nesvorny D (2012) Dynamics of pebbles in the vicinity of a growing planetary embryo: hydro-dynamical simulations. Astron. Astrophys. 546:A18.
  15. Bai XN, Stone JM (2010) Dynamics of Solids in the Midplane of Protoplanetary Disks: Implications for Planetesimal Formation. Astrophys. J. 722:1437–1459.
  16. Andrews SM, Wilner DJ, Hughes AM, Qi C, Dullemond CP (2010) Protoplanetary Disk Structures in Ophiuchus. II. Extension to Fainter Sources. Astrophys. J. 723:1241–1254.
  17. Chiang EI, Goldreich P (1997) Spectral Energy Distributions of T Tauri Stars with Passive Circumstellar Disks. Astrophys. J. 490:368–376.
  18. Hayashi C (1981) Structure of the Solar Nebula, Growth and Decay of Magnetic Fields and Effects of Magnetic and Turbulent Viscosities on the Nebula. Progress of Theoretical Physics Supplement 70:35–53.
  19. Owen JE, Ercolano B, Clarke CJ (2011) Protoplanetary disc evolution and dispersal: the implications of X-ray photoevaporation. Mon. Not. R. Astron. Soc. 412:13–25.
  20. Haisch KE, Lada EA, Lada CJ (2001) Disk Frequencies and Lifetimes in Young Clusters. Astrophys. J. Let. 553:L153–L156.
  21. Lodders K (2003) Solar System Abundances and Condensation Temperatures of the Elements. Astrophys. J. 591:1220–1247.
  22. Brauer F, Dullemond CP, Henning T (2008) Coagulation, fragmentation and radial motion of solid particles in protoplanetary disks. Astron. Astrophys. 480: 859–877.
  23. Ricci L, Testi L, Natta A, Neri R, Cabrit S, Herczeg GJ (2010) Dust properties of protoplanetary disks in the Taurus-Auriga star forming region from millimeter wavelengths. Astron. Astrophys. 512:A15.
  24. Johansen A, Youdin A, Mac Low MM (2009) Particle Clumping and Planetesimal Formation Depend Strongly on Metallicity. Astrophys. J. Let. 704:L75–L79.
  25. Birnstiel T, Klahr H, Ercolano B (2009) A simple model for the evolution of the dust population in protoplanetary disks. Astron. Astrophys. 539:A148.
  26. Jessberger EK et al. (2001) Properties of Interplanetary Dust: Information from Collected Samples, eds. Grün E, Gustafson BAS, Dermott S, Fechtig H. p. 253.
  27. Draine BT (1979) On the chemisputtering of interstellar graphite grains. Astrophys. J. 230:106-115.
  28. Levison HF, Duncan MJ, Thommes E (2012) A Lagrangian Integrator for Planetary Accretion and Dynamics (LIPAD). Astron. J. 144:119–138.
  29. Duncan MJ, Levison HF, Lee MH (1998) A Multiple Time Step Symplectic Algorithm for Integrating Close Encounters. Astron. J. 116:2067–2077.
  30. Benz W, Asphaug E (1999) Catastrophic Disruptions Revisited. Icarus 142:5–20.
  31. Kretke KA, Levison HF (2014) Challenges in Forming the Solar System’s Giant Planet Cores via Pebble Accretion. Astron. J. 148:109.
  32. Morbidelli A, Tsiganis K, Crida A, Levison HF, Gomes R (2007) Dynamics of the Giant Planets of the Solar System in the Gaseous Protoplanetary Disk and Their Relationship to the Current Orbital Architecture. Astron. J. 134:1790–1798.
  33. Birnstiel T, Klahr H, Ercolano B (2012) A simple model for the evolution of the dust population in protoplanetary disks. Astron. Astrophys. 539:A148.
  34. Wurm G, Paraskov G, Krauss O (2005) Growth of planetesimals by impacts at 25 m/s. Icarus 178:253–263.
  35. Fischer RA, Ciesla, FJ (2014) Dynamics of the terrestrial planets from a large number of N-body simulations. Earth and Planetary Science Letters 392: 28–38.
  36. Zsom A, Ormel CW, Güttler C, Blum J, Dullemond CP (2010) The outcome of protoplanetary dust growth: pebbles, boulders, or planetesimals? II. Introducing the bouncing barrier. Astron. Astrophys. 513:A57.
  37. Güttler C, Blum J, Zsom A, Ormel CW, Dullemond CP (2010) The outcome of protoplanetary dust growth: pebbles, boulders, or planetesimals?. I. Mapping the zoo of laboratory collision experiments. Astron. Astrophys. 513:A56.
  38. Ormel CW, Cuzzi JN, Tielens AGGM (2008) Co-Accretion of Chondrules and Dust in the Solar Nebula. Astrophys. J. 679:1588–1610.
  39. Gammie CF (1996) Layered Accretion in T Tauri Disks. Astrophys. J. 457:355–+.
  40. Bai XN (2014) Hall-effect-Controlled Gas Dynamics in Protoplanetary Disks. I. Wind Solutions at the Inner Disk. Astrophys. J. 791:137.
  41. Touboul M, Kleine T, Bourdon B, Palme H, Wieler R (2007) Late formation and prolonged differentiation of the Moon inferred from W isotopes in lunar metals. Nature 450:1206–1209.
  42. Johansen A, Mac Low MM, Lacerda P, Bizzarro M (2015) Growth of asteroids, planetary embryos and Kuiper belt objects by chondrule accretion. ArXiv e-prints.
  43. Kenyon SJ, Hartmann L (1987) Spectral energy distributions of T Tauri stars - Disk flaring and limits on accretion. Astrophys. J. 323:714–733.
  44. Garaud P, Lin DNC (2007) The Effect of Internal Dissipation and Surface Irradiation on the Structure of Disks and the Location of the Snow Line around Sun-like Stars. Astrophys. J. 654:606–624.
  45. Bitsch B, Morbidelli A, Lega E, Crida A (2014) Stellar irradiated discs and implications on migration of embedded planets. II. Accreting-discs. Astron. Astrophys. 564:A135.
  46. Lecar M, Franklin FA (1973) On the original distribution of the asteriods. I. Icarus 20:422–436.
  47. Walsh KJ, Morbidelli A, Raymond SN, O’Brien DP, Mandell AM (2012) Populating the asteroid belt from two parent source regions due to the migration of giant planets—”The Grand Tack”. Meteoritics and Planetary Science 47:1941–1947.
  48. Bottke WF et al. (2005) Linking the collisional history of the main asteroid belt to its dynamical excitation and depletion. Icarus 179:63–94.
  49. Marchi S et al. (2012) The Violent Collisional History of Asteroid 4 Vesta. Science 336:690–.
  50. Bottke WF (2014) On the Origin and Evolution of Vesta and the V-Type Asteroids. LPI Contributions 1773:2024.
  51. Batalha NM et al. (2013) Planetary Candidates Observed by Kepler. III. Analysis of the First 16 Months of Data. Astrophys. J. Supp. 204:24.
  52. Gundlach B, Blum J (2015) The Stickiness of Micrometer-sized Water-ice Particles. Astrophys. J. 798:34.
  53. Morbidelli A, Bottke WF, Nesvorný D, Levison HF (2009) Asteroids were born big. Icarus 204:558–573.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minumum 40 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description