A pebble accretion model for the formation of the terrestrial planets in the Solar System

Authors: Anders Johansen, Thomas Ronnet, Martin Bizzarro, Martin Schiller, Michiel Lambrechts, Åke Nordlund, Helmut Lammer

Abstract

A new theory for formation of Earth by accretion of small pebbles explains how our planet acquired its water and carbon.

INTRODUCTION

The long formation time scale of gas giants and ice giants in the outer regions of protoplanetary discs by traditional planetesimal accretion (1, 2) instigated the development of the pebble accretion theory in which the pebbles drifting through the protoplanetary disc are accreted rapidly by the growing protoplanets (3, 4). While pebble accretion clearly aids the formation of gas giants (5), the formation of terrestrial planets has so far, with a few exceptions (6–9), mainly been explored in classical N-body simulations where terrestrial planets grow by successive impacts between increasingly massive protoplanets (10).

Observations of protoplanetary discs reveal that very young stars are orbited by several hundred Earth masses of pebbles (11), embedded in a gaseous disc of 100 times higher mass. This pebble population vanishes on a characteristic time scale of a few million years (12), likely due to a combination of radial drift and planetesimal formation outside gaps in the gas caused by the gravity of the growing planets (13–15). The orbital speed of the gas in a protoplanetary disc is set by a balance between the gravity from the central star and the outward-pointing pressure force. The resulting sub-Keplerian orbital motion of the gas acts as a headwind on the pebbles, draining them of their orbital angular momentum and causing them to drift radially toward the star (16). The inner regions of protoplanetary discs, where terrestrial planets form, therefore witness a flow of hundreds of Earth masses of pebbles throughout the lifetime of the protoplanetary disc.

The millimeter-sized chondrules found within primitive meteorites likely represent pebbles that formed around the young Sun; the leading mechanism for heating and melting these pebbles to form igneous chondrules is shock waves in the solar protoplanetary disc (17). Alternatively, chondrules may have formed from the molten debris of collisions between massive protoplanets (18). Chondrules formed over the first 3 million years (Ma) of the evolution of the protoplanetary disc (19, 20) and dominated the mass budget when the parent bodies of the ordinary chondrite meteorite class formed in the inner regions of the solar protoplanetary disc ∼3 Ma after the formation of the Sun. The large fraction of millimeter-sized chondrules in both ordinary chondrites and carbonaceous chondrites (CCs) is a strong indication that planetesimals in the solar protoplanetary disc accreted in the gravitational collapse of pebble clumps (possibly consisting of a mixture of individual chondrules and larger aggregates of chondrules and dust) that could have formed through the streaming instability. This hydrodynamical instability arises in the radial drift flow of pebbles and causes the pebbles to form dense filaments that collapse gravitationally into planetesimals with a characteristic size of 100 km (7, 21). The largest planetesimals then continue to accrete pebbles and grow to protoplanets with masses between the Moon and Mars.

The radial flux of pebbles toward the star was recently demonstrated to determine the outcome of planetary growth in the inner regions of the protoplanetary disc (22): High pebble fluxes from the outer protoplanetary disc lead to formation of migrating chains of super-Earths, while any reduction in the pebble flux, e.g., by the emergence of giant planets in the outer Solar System whose gravity on the gas disc acts as a pressure barrier for pebble drift, strands the protoplanets at Mars masses. This was followed by a phase of giant impacts akin to classical terrestrial planet formation.

Here, we explore the possibility that the radial drift of pebbles continues to drive the growth from protoplanets to masses comparable to Earth and Venus by pebble accretion. The role of giant impacts for terrestrial planet formation is here reduced to the Moon-forming impact between Earth and Theia, an additional terrestrial planet that collided with Earth after an instability in the system of primordial terrestrial planets that orbited our young Sun (23). We show that a pebble accretion scenario for terrestrial planet formation provides explanations for several properties of the terrestrial planets in the Solar System, including (i) the masses and orbits of Venus, Earth, and Mars, (ii) the isotopic composition of Earth and Mars, and (iii) the delivery of carbon and water to Earth in amounts that are comparable to the inferred reservoirs.

Our model, depicted as a sketch in Fig. 1, provides important constraints into the timing and flux of primitive outer disc material to the inner Solar System. In particular, we conclude that Jupiter did not halt inward mass transport of outer Solar System dust aggregates to the accretion region of terrestrial planets, in contrast to models that invoke Jupiter as an impermeable barrier that separated the inner and outer Solar System (24).

RESULTS

We develop and present here a model for terrestrial planet formation driven dominantly by pebble accretion, but we include also self-consistently the contribution from planetesimal accretion. We start by identifying the water ice line as the most likely location for the formation of a first generation of planetesimals that acted as seeds for pebble accretion in the solar protoplanetary disc. As the stellar luminosity decreases with time, the ice line moves interior to the region of terrestrial planet formation. This clearly has implications for the delivery of water to terrestrial planets that form by pebble accretion.

Location of the water ice line in the solar protoplanetary disc

Planetesimal formation by the streaming instability requires that the ratio of the surface density of pebbles relative to that of the gas is elevated above a threshold metallicity, higher than the nominal solar value, to trigger the formation of dense filaments that collapse to form planetesimals (25). Such an increased pebble density has been demonstrated to occur early in the evolution of the protoplanetary disc by pileup of dust released by icy pebbles interior of the water ice line (26) or by the deposition of water vapor transported by diffusion to exterior of the ice line (27, 28). We therefore assume in our model that an early generation of planetesimals formed at the water ice line and that the planetesimals that grew to form Venus, Earth, Theia, and Mars were among them.

The orbital distance of the water ice line depends on the luminosity of the young Sun as well as on the efficiency of viscous heating by the gas that is accreted. We present here both a model where viscous heating is provided at all distances from the star as well as a more realistic model where viscous heating is only provided when the magnetorotational instability is active above T = 800 K (29). This “dead zone” model is described in detail in Materials and Methods.

We assume that the protoplanetary disc column density evolves as a viscous α disc from an initial accretion rate of M·*0=10−6 M⊙ year−1 at t = 0 down to M·*1=10−6 M⊙ year−1 after 5 Ma. This decrease in the accretion rate implies that the heating becomes more dominated by the stellar irradiation with time. We show the calculated temperature of the solar protoplanetary disc in Fig. 2. Additional details of the structure of the protoplanetary disc model with a dead zone are shown in fig. S1.

The ice line in the more realistic dead zone model sits initially in the region between 1.2 and 2.0 astronomical units (AU), whereas the model with viscous heating everywhere has a primordial ice line at 7 AU. As the luminosity of the Sun decreases with the contraction of the protostar, the water ice line moves interior of 1 AU after around 1.5 Ma to 2 Ma in both models. Because the limiting case where viscous heating is applied everywhere is not physically realistic, as demonstrated in magnetohydrodynamical simulations (30), we adopt the dead zone model when integrating the growth tracks of planets. Overall, we infer from the dead zone model that the first planetesimals in the inner Solar Systems likely formed in a region devoid of magnetically driven turbulence between 1.2 and 2.0 AU.

Analytical growth tracks of terrestrial planets

Protoplanets in the protoplanetary disc will grow at a rate M· given by the accretion of pebbles and planetesimals and migrate toward the star at a rate r·. The combination of these two rates yields the growth track equation dM/dr=M·/r·. This equation can be integrated analytically (5, 31) for the case of growth by two-dimensional (2D) pebble accretion (relevant when the scale height of the pebbles is smaller than the pebble accretion radius) and standard, inward type I planetary migration to yield the analytical growth track expressionM(r)=Mmax[1−(rr0)1−ζ]3/4(1)

We express the growth track here in the limit where M ≫ M0, where M0 is the starting mass. The gas temperature of the protoplanetary disc is assumed to depend on the distance from the Sun as T ∝ r^−ζ^ and we take a constant value of ζ = 3/7, valid for a protoplanetary disc that intercepts stellar irradiation at a grazing angle (32). This power law approximation for the temperature holds in the dead zone of the protoplanetary disc where the terrestrial planets grow and migrate. Equation 1 has two free r0 is the starting position of the protoplanet, and Mmax is the mass of the protoplanet when it has migrated to r = 0. We fit now the growth track to Mars and Venus, assuming that their current masses reflect their primordial masses when the solar protoplanetary disc dissipated. Earth has suffered a giant impact that led to the formation of the Moon, so the current mass of Earth does not represent its preimpact mass. Matching the masses and locations of Mars (index 1) and Venus (index 2) in Eq. 1 allows us to divide out the Mmax parameter to yieldr0ζ−1[r11−ζ−(M1/M2)4/3r21−ζ]=1−(M1/M2)4/3(2)

We see how the starting position, r0, is dominantly set by the planet with the smallest mass (r0 → r1 for M1 ≪ M2). For Mars and Venus, we obtain r0 = 1.59 AU, slightly exterior of Mars’ current orbit at r = 1.53 AU. This location is broadly consistent with the likely location of the water ice line during the first million years of protoplanetary disc evolution when considering that the terrestrial planet zone was heated mainly by the radiation from the central star. The maximum mass is then set by applying Eq. 1 to Venus, giving Mmax = 1.74 ME.

Our successful fit of a single growth track to Mars and Venus is always mathematically possible when the outer planet (Mars) is less massive than the inner planet (Venus). How physically plausible this growth track is lies in the value of the maximum mass Mmax. This maximum mass, or migration mass, depends on a number of parameters of the protoplanetary disc (5). We focus here on the pebble Stokes number, St, a dimensionless number that characterizes the frictional stopping time of the pebbles, and the ratio of the radial flux of pebbles relative to the flux of the gas, ξ=Fp/M·⋆, where M·⋆ is the gas accretion rate onto the star. The maximum mass furthermore depends on starting position r0, which we fix here to 1.6 AU, and the temperature profile for which we take a power law with index ζ = 3/7 and temperature T1 = 140 K at 1 AU; we neglect any temporal dependence of the stellar luminosity. We illustrate the dependence of the maximum mass on these parameters in Fig. 3. The maximum mass that yields the best fit to the orbit and mass of Venus is obtained for a range of St from 0.001 to 0.1 and a range of ξ from 0.004 to 0.008. These ranges represent nominal values of the Stokes number of millimeter-sized pebbles and radial pebble fluxes (5). From this, we conclude that the masses and orbits of Venus and Mars are consistent with growth by accretion of millimeter-sized pebbles combined with standard type I planetary migration.

Numerical growth tracks of terrestrial planets

After the successful application of the analytical growth tracks to Venus and Mars, we now turn to numerical integration of the masses and orbits of protoplanets undergoing pebble accretion, planetesimal accretion, and inward type I migration, using the code developed and presented in (5, 2). The protoplanetary disc model is described in Materials and Methods. Planetesimals with fixed radii of 100 km are present from the beginning of the simulation as a Gaussian belt of width 0.05 AU centered at 1.6 AU. The total planetesimal mass is approximately 0.5 ME; this mass is broadly consistent with models of early planetesimal formation at the water ice line in protoplanetary discs of solar metallicity (28).

We start protoplanets representing Venus, Earth, Theia, and Mars after t = 0.66 Ma, 0.93 Ma, 1.50 Ma, and 2.67 Ma, respectively, after the formation of the protoplanetary disc. The protoplanets are all given initial masses of M0 = 1.0 × 10^−3^ ME. The starting times are chosen to be consistent with the growth time from an initial mass function of planetesimals peaking at M ∼ 10^−7^ − 10^−6^ ME to our starting mass of M0 = 10^−3^ ME by planetesimal collisions and pebble accretion (7); the difference in growth time scales for the different bodies is essentially set by their different levels of eccentricity, which controls the pebble accretion rate, until the bodies are circularized by the transition to rapid pebble accretion at M ∼ 10^−4^ ME (7), approximately the mass of the largest asteroid Ceres. We consider for simplicity pebble accretion only in the Hill regime, valid from a factor 10 times lower than our starting mass (4). We ignore the gravitational interaction between the planets as they grow; in the Supplementary Material, we present the results of a full N-body simulation showing that the single-planet growth tracks are reproduced when including the gravity between the growing planets.

In Fig. 4, we show growth tracks starting at r0 = 1.6 AU. The different starting times of the protoplanets make a hierarchy of final masses and orbits. We confirm from the numerical integrations that Venus and Mars belong to the same growth track, although Mars grows faster than estimated because of accreting a substantial mass fraction (approximately 60%) in planetesimals from its birth belt. The other planets accrete a similar contribution from planetesimals, but their total mass budget is dominated by pebble accretion and, hence, they follow the analytical growth track better than Mars. Earth obtains a final mass of 0.6 ME, but this mass is augmented by our assumption that Theia was a fifth terrestrial planet with an original orbit between Earth and Mars. We construct Theia to have a mass of 0.4 ME. This mass is broadly consistent with models for Moon formation that invoke a very massive impactor to vaporize the Earth-Theia collision debris (33).

Isotopic composition of the planets

The mass-independent isotopic compositions of various elements such as O, Ti, Cr, and Ca have been measured for Earth as well as for Mars (from martian meteorites) and the major meteorite classes (34) and show that variability exists between the various classes of meteorites and terrestrial planets. This nucleosynthetic variability likely reflects either a time-dependent supply of grains with different nucleosynthetic heritage to the solar protoplanetary disc or, alternatively, a temperature-dependent destruction of either presolar grains or grains condensed in the local interstellar medium (35–37), resulting in a composition dichotomy between the inner, hot regions and the outer, cold regions of the solar protoplanetary disc.

Broadly speaking, the isotopic abundances of meteorites fall in two distinct groups when plotting any two elements against each other. One group is represented by the noncarbonaceous (NC) meteorites, including ordinary chondrites, and the other group by the CC. Earth sits on a mixing line between the two components, together with the enstatite chondrites (38), establishing that a compositional gradient exists between the NC and CC reservoirs. Therefore, Earth could have formed either from material akin to enstatite chondrites to yield the measured isotopic fingerprint (39) or from a combination of material akin to ordinary chondrites mixed with material akin to CC (36, 38). The latter picture nevertheless opens the question of how CC, which likely formed beyond the ice line in the solar protoplanetary disc, could have been a major source of material for the terrestrial planets.

On the basis of the mean isotopic abundances of the NC chondrite group on the one hand and of the CC group on the other hand (38), it has been estimated that Earth and Mars accreted between 30 and 43% and between 15 and 30% mass from the latter group, respectively, when considering simultaneously the isotopes ^54^Cr and ^50^Ti. On the basis of the isotopes ^54^Cr and ^17^O, Earth and Mars should have accreted between 18 and 32% and between 2.6 and 18% from the CC group, respectively. Using instead the abundance of ^48^Ca in the ureilites and CI chondrites as end-members for the mixing, Earth requires a contribution of 42% from pebbles of CI composition and Mars requires 36% (36). We used ^48^Ca (with end-members ureilites and CI chondrites) to match the composition of Earth and Mars in Fig. 4, but choosing other end-members from the NC and the carbonaceous groups would give qualitatively similar results.

To calculate the contribution of material from the outer regions of the solar protoplanetary disc to the terrestrial planets formed in our simulations, we assume that the planetesimals in the birth belt as well as the pebbles (and chondrules) that drift through the disc for the first tCI (a variable time) have the isotopic composition of the ureilite meteorites (an end-member of the NC meteorite group understood to represent the initial inner disc composition (36)), followed by an influx of dust-aggregate pebbles of composition similar to CI chondrites for the last 5.0 Ma − tCI, before the protoplanetary disc finally dissipates after 5 Ma. The time tCI can be chosen to fit best the measured compositions of Earth and Mars. In Fig. 4C, we match (36) in getting approximately 42% contribution of pebbles with CI composition to Earth for a transition time tCI = 3.8 Ma. In fig. S2, we show that the CI material resided beyond 10 AU after 10^5^ years of evolution of the protoplanetary disc, having a total gas and dust mass of 0.027 M⊙ out of the total disc mass of 0.044 M⊙ at this stage, and was subsequently pushed outward by the viscous expansion of the disc before falling back through the terrestrial planet region at tCI = 3.8 Ma. At this time, the gas flux through the protoplanetary disc is M·*≈2×10−9 M⊙ year−1, and the radial pebble flux is Fp ≈ 2.4 ME Ma^−1^. Venus and Theia obtain very similar compositions to Earth because all three planets accrete from the same drifting material.

Mars would get a similar fraction as the other planets of NC relative to carbonaceous material from pebble accretion; however, Mars has an isotopic composition that lies closer to the ordinary chondrites than Earth. This could indicate that Mars terminated its growth earlier than Venus, Earth, and Theia and hence avoided the incorporation of pebbles of CI composition that drifted in later (36). However, that raises the question of how Mars’ eccentricity and inclination could have been excited during the gaseous disc phase. Instead, we assume here that all the four planets accrete planetesimals within their birth belt. This contribution of material with inner Solar System ureilite-like composition affects Mars the most because it accretes the largest fraction of its mass in the planetesimal birth belt. Hence, we can qualitatively explain the compositional difference between Earth and Mars as a natural consequence of Mars being less massive than Earth and thus having a larger planetesimal accretion contribution to its mass. The compositions of Earth and Mars are broadly reproduced by a pebble composition transition time tCI in the range between 3.5 and 4.0 Ma.

The isotopic signature of iron in the mantle of Earth is very similar to that of iron in CI chondrites (40). This agreement is consistent with our model. The early-accreted iron from the first generation of pebbles came from the NC reservoir and was hence largely reduced, with iron present as metallic iron nuggets as in the ordinary chondrites. This iron accreted together with ice until the envelope reached the ice sublimation temperature at a mass of 0.02 ME after approximately 2 Ma (see sections on water delivery and carbon delivery by pebble accretion below). The iron arriving together with ice would have oxidized in contact with liquid water in the protoplanet, melted by the heating by ^26^Al, but this early oxidized iron only constitutes a small fraction of the total oxidized iron found in the mantle of Earth today. Material accreted later contained too little ^26^Al to melt, and hence, the iron remained in solid form in a primitive mantle overlying the earlier differentiated layers. The bulk of the metal of NC composition separated from the silicates and entered the core as the protoplanet finally melted by the accretion heat. The second generation of pebbles of CI composition contained a large fraction of oxidized iron that remains in the mantle today. This is in agreement with the abundances of siderophile elements in Earth’s mantle, which require increasingly oxidizing conditions during silicate-metal separation with increasing mass of the growing Earth (41). The oxidized iron in Earth’s mantle constitutes around 16% of the total iron on Earth, but the outer core likely contains a substantial fraction of oxidized iron as well (42). Given our approximately 40% contribution from pebbles of CI composition to Earth, these pebbles should be approximately 60% oxidized to account for the iron budget in Earth’s mantle. The iron in the CI chondrites is nearly 100% oxidized, but the oxidation of primordial FeS and metallic iron could have occurred by liquid water flow in the parent body (43). The total mass of the oxidized iron Earth’s mantle implies furthermore that the iron in Earth’s core consists of 25% iron derived from CI-like solids and 75% iron derived from NC material.

Jupiter has been proposed to play an important role in the separation of the NC and CC reservoirs in the solar protoplanetary disc (24). However, hydrodynamical simulations show that the gap formed by the combined action of Jupiter and Saturn in the protoplanetary disc is permeable to small dust aggregates (44), particularly if the gap edge develops turbulence, e.g., through the Rossby wave instability (45). Larger pebbles will become trapped in the peak of the gas pressure exterior to the gap—but fragmenting pebble collisions continuously create dust that can pass through the gap (46). Hence, in our view, Jupiter did not present an efficient barrier to the supply of dust particles from the outer regions of the protoplanetary disc to the terrestrial planet zone. Jupiter could nevertheless have acted to reduce the flux of pebbles, consistent with the low pebble flux needed to explain the masses and orbits of Venus, Earth, Theia, and Mars in our model.

Free parameters in our model

It is remarkable that it is possible to choose a starting position and four starting times of the protoplanets and match the orbits, masses, and compositions of Venus, Earth, Theia, and Mars by pebble accretion, including a contribution from planetesimal accretion in the birth belt, using realistic parameters for the pebble sizes and the radial pebble flux. Our model nevertheless has a number of parameters that are to a varying degree free to choose. We briefly discuss the choice of those here.

We distinguish between three classes of (a) free parameters that are unconstrained from observations, (b) free parameters that are constrained from observations, and (c) fixed parameters that are constrained by observations and set according to the most reasonable value. The truly free parameters in class (a) are t0 (the four starting times of the protoplanets when they have the starting mass M0), r0 (the starting position of the protoplanets), tCI (time of infall of pebbles of CI composition), Mpla (the total planetesimal mass in the birth belt), and ξ (the pebble mass flux relative to the gas mass flux, which we fine-tune around the solar composition to fit the growth tracks). The only free parameter in class (b) is Rp (the pebble size, we set it here to be 1 mm, similar to chondrule sizes and to pebble sizes inferred from observations of protoplanetary discs). The fixed parameters in class (c) are ζ = 3/7 (the temperature gradient in the protoplanetary disc), M·*0=10−6 M⊙ year−1 and M·*1=10−9 M⊙ year−1 (the initial and final accretion rates of the protoplanetary disc), tdisc = 5 Ma (the lifetime of the protoplanetary disc), α = 10^−2^ (the turbulent viscosity in the protoplanetary disc), and δ = 10^−4^ (the diffusion coefficient of the gas, which governs the degree to which the pebbles sediment to the midplane). We chose reasonable values, based on observations of protoplanetary discs, for all parameters in class (c) and did not attempt to vary them to obtain better fits to the observed properties of the terrestrial planets.

The combined set of parameters describes a model that has as outcome the masses, orbits, and compositions of the planets Venus, Earth, Theia, and Mars. Here, Theia is a special case because its mass is only constrained to add up with the mass of proto-Earth to obtain the current mass of Earth. We therefore consider the real prediction of our model to be MV (the mass of Venus), MM (the mass of Mars), ME + MT (the total mass of proto-Earth and Theia), rV (the orbit of Venus), rM (the orbit of Mars), fCI, E (the composition of Earth), fCI, T (the composition of Theia, must be equal to that of Earth), and fCI, M (the composition of Mars). These eight model predictions are constructed by essentially choosing three starting times of the protoplanets (considering the starting times of Earth and Theia to be degenerate), the starting position of the protoplanets, the radial the pebble flux relative to the gas flux, the time of infall of pebbles of CI composition, and the total mass of planetesimals in the birth belt. This gives a total of seven truly free parameters that are used to fit the data. The fact that the number of parameters is smaller than the number of successful predictions implies that the model has a genuine explanation power for the properties of the terrestrial planets in the Solar System.

Monte Carlo populations

We now explore the effect of varying the birth positions and times of the protoplanets as well as the width of the planetesimal belt. We start 1000 protoplanets with masses M0 = 10^−3^ ME and give them random starting times between 0.5 Ma and 5.0 Ma. The positions are drawn from a Gaussian belt centered at r0 = 1.6 AU and with a width of either W = 0.05 AU or W = 0.5 AU. The growth tracks of Fig. 4 were integrated using the narrow belt. The wider belt has 10 times more mass in planetesimals than the narrow belt, so a total of approximately 5 ME.

In Fig. 5, we show the final positions and masses of the considered protoplanets. We demonstrate results for two values of the pebble ξ = 0.0036 as in Fig. 4 and a higher value of ξ = 0.01. The narrow belt case faithfully reproduces the characteristic masses and orbits of Venus and Mars. Widening the planetesimal belt to 0.5 AU increases the maximally reached planetary masses to between 2 ME and 5 ME (super-Earths). The mass hierarchy with increasing planetary mass with decreasing distance from the star is maintained as in the narrow belt case, but the terrestrial planets no longer occupy characteristic positions in mass versus semimajor axis.

The bottom panel of Fig. 5 shows the results when increasing the pebble-to-gas flux ratio to ξ = 0.01. This is closer to the nominal flux in the absence of drift barriers, such as planets, further out in the disc. The narrow belt case now fits Earth better, but Venus is well below the characteristic growth track and planets at the location of Mars now reach much higher masses before they start to migrate toward the star. Super-Earths of masses up to 10 ME pile up at the inner edge of the protoplanetary disc. The decisive role of the pebble flux in determining whether protoplanets in the inner regions of a protoplanetary disc grow to super-Earths or “only” to Mars masses was already found in (22); here, we demonstrate that the combination of a low pebble flux and a narrow planetesimal belt leads to the formation of analogs of the terrestrial planets in the Solar System in terms of both their masses and orbits.

Sublimation of volatiles in the planetary envelope

The accreted pebbles will carry molecules that are in solid form at the temperature level of the surrounding protoplanetary disc. Because the water ice line is interior of the accretion region of Earth and Venus during the main accretion phase of the terrestrial planets, this implies that the pebbles will be rich in volatiles. However, the growing planets attract gas envelopes that are heated by the high pebble accretion luminosity, and this heat will in turn process the pebbles thermally while they are falling toward the surface of the protoplanet. Our methods for calculating the temperature and density of the envelope in hydrostatic and energy equilibrium are described in Materials and Methods.

In Fig. 6, we show the temperature and density of the envelopes of our Earth and our Mars analogs at t = 5 Ma as a function of the distance from the center of protoplanet, for three different values of the opacity κ0 (see Materials and Methods). The density of the hydrogen/helium envelope rises exponentially inward from the Bondi radius (defined as RB=GM/cs2, where cs is the sound speed of the gas in the protoplanetary disc). For Earth, the envelope temperature lies above the 2600 K of silicate sublimation close to the planetary core. This high surface temperature will lead to extensive melting of the protoplanet and to the formation of a magma ocean of at least 1000 km in depth (47). Our Mars analog does not reach silicate melting temperatures at the surface, but trapped heat released by ^26^Al within the protoplanet could heat the interior above the melting temperature.

We indicate in Fig. 6 the location of the water ice line with blue pluses labeled with the entropy, relative to the entropy level of the surrounding protoplanetary disc, at the depth of the water ice line. Hydrodynamical simulations have demonstrated that flows from the protoplanetary disc reach down to a relative entropy level of s/s0 ≈ 0.2 (48). Hence, for the high-opacity cases, with κ0 = 0.1 m^2^ kg^−1^ and 1.0 m^2^ kg^−1^, the water vapor that is released by sublimating the ice layers from the pebbles will be recycled back to the protoplanetary disc. Such a high opacity is easily reached for a dust-to-gas ratio around the solar value or higher.

The recycling flows will first transport water vapor back across the ice line, where the vapor nucleates on tiny dust particles that carry the dominant surface area. Whether these ice-covered particles eventually reach the protoplanetary disc depends on the time scale to coagulate to large enough sizes to sediment back through the gas flow. The coagulation time scale of the ice particles can be calculated from the approximate expressionτc∼aiceρ•(Δv)ρice(3)Here aice is the radius of the particles, ρ~•~ is their material density, Δv is the particle collision speed, and ρice is the spatial density of the nucleated ice particles. We then assume that the collision speed is dominated by Brownian motion. That yieldsτc∼aiceρ•kBT/mρgϵ−1∼4.2 years (aiceμm)5/2(T170 K)−1/2(ρg10−6 kg m−3)−1(ϵ0.01)−1(4)We defined here ϵ = ρice/ρg as the ratio of the density of ice particles to the gas. The coagulation time scale should now be compared to the time scale for the recycling flow to transport the particles from the ice line at Rice out to the Hill radius with the characteristic speed vrec ∼ ΩRice. That givesτrec∼(RH/Rice)Ω−1∼10Ω−1∼2.0 years(5)This time scale is comparable to the coagulation time scale around the ice line, but the coagulation time scale increases steeply in the decreasing gas density exterior of the ice line. Hence, the freshly nucleated ice particles do not coagulate appreciably during their transport from the ice line to the Hill radius with the recycling flows.

Water delivery by pebble accretion

Figure 7 shows the water fraction of our growing Earth analog, assuming that water can only be delivered to the protoplanet in the form of ice accreted before the envelope reaches the temperature to sublimate water ice. We assume that the accreted planetesimals are dry [because they formed early enough to dry out due to the heat from ^26^Al; see (49)] and that the accreted pebbles carry either the nominal 35% ice by mass, which assumes that some oxygen is bound in the volatile CO molecule (50), or a much dryer 10%. The water fraction peaks at a protoplanet mass just below 0.02 ME (approximately the mass of the Moon) and then falls steeply for higher masses as water ice is sublimated in the hot envelope. The final water fraction of our Earth analog is 3000 parts per million (ppm) for the nominal ice fraction (35%) and 1000 ppm for the low ice fraction (10%). The nominal case lands thus at the high end of the 2000 to 3000 ppm estimated for Earth (51). Theia, with a final mass of 0.4 ME, acquires 4000 ppm of water and would increase the total water fraction of Earth to 3500 ppm after the giant impact, but this value could be reduced by loss of volatiles in energetic the impact.

All the water is delivered as ice in the early stages of protoplanet growth when the protoplanet is still solid. The protoplanet will melt later by accretion heat and radioactive decay, and hence, the magma ocean inherits a substantial fraction of water. Parts of this water will later be released from the magma during crystallization (52) and likely mixed with water formed from hydrogen ingassed from the protoplanetary disc (53). Any water that is outgassed during the growth of the protoplanet will be retained as a high-density primary atmosphere below the hydrogen/helium envelope, protected against loss to diffusive convection by the gradient in the mean molecular weight between the atmosphere and the envelope. Our model thus shows that water could have been delivered to Earth from “pebble snow” accreted during the early growth stages when the protoplanet was still solid and the envelope cold enough for water ice to survive the passage to the surface.

Carbon delivery by pebble accretion

In Fig. 6, we also indicate the temperatures where organic compounds are vaporized. We take the range between 325 and 425 K for the vaporization of 90% of the carbon by pyrolysis and sublimation (54). Pyrolysis and sublimation of organics convert these molecules to very volatile carbon-bearing molecules such as CH4, CO2, and CO. These species will not condense in the envelope and can diffuse freely to the recycling zone. We assume that protoplanetary disc flows penetrate into the Hill sphere with the characteristic speed ΩRB at the Bondi radius RB. This results in a mass loss rate of carbon-bearing molecules asM·C=4πRB2ρCΩRB(6)

We now equate this expression with the mass accretion rate of carbon, fCM·, and isolate the density of carbon-bearing molecules, ρC, at the Bondi radius to beρC=fCM·4πRB3Ω=1.1×10−11 kg m−3(fC0.05)(M·0.2 ME Ma−1)(M0.5 ME)−3×(cs7.0×102 m s−1)6(rAU)3/2(7)This density is much lower than the gas density in the protoplanetary disc, ρg ∼ 10^−8^ kg m^−3^, and hence, the equilibrium ratio of volatile carbon molecules relative to the gas is much lower than the ambient density of CO, the main carbon-bearing molecule in inner regions of the protoplanetary disc. The equilibrium mixing ratio at the Bondi radius, ρC/ρg, will be maintained in the entire envelope by turbulent diffusion with coefficient D. This diffusion equilibrium is maintained after a characteristic time scale ofτ=RB2D(8)

The characteristic value of the gas mass in the envelope of terrestrial planets that grow by pebble accretion is Mg ∼ 10^−4^ ME. The time scale to match this gas mass with carbon-bearing molecules released in the organics sublimation/pyrolysis region isτC=MgfCM·=1×104 years (Mg10−4 ME)(fC0.05)−1(M·0.2 ME Ma−1)(9)

The turbulent diffusion time scale must be shorter than this value to avoid building up a layer in the envelope dominated by carbon-bearing molecules, which would lead to an increase in the mean molecular weight of the gas and hence potentially protect the carbon-bearing molecules from turbulent diffusion. Setting the turbulent diffusion coefficient D = vturb × RB, we get a turbulent speed ofvturb=RBτ=1.3×10−3 m s−1(M0.5 ME)(cs7.0×102 m s−1)−2(τ1×104 years)−1(10)This is an extremely modest speed compared with the sound speed, and hence, we conclude that the carbon-bearing species released by vaporization of organics will become transported to the recycling flows of the protoplanetary disc before they can build up a substantial mass fraction in the envelope.

Only approximately 10% of the carbon provided by the pebbles will therefore make it below the region where organics are vaporized. The remaining pure carbon dust burns at 1100 K (54). This temperature is reached approximately 10,000 km above the surface of our Earth analog. Assuming that the carbon is only retained if it makes it to the protoplanet’s surface, we estimate the total carbon delivery to our Earth analog as follows. We assume that the carbon contents of the early generation of pebbles are similar to ordinary chondrites (approximately 0.3%) and similar to CI chondrites for the late generation (approximately 5%). The result is shown in Fig. 7. The fraction of carbon is initially constant at 3000 ppm (set by the carbon contents of ordinary chondrites) but falls steadily after carbon vaporization starts in the envelope. Our Earth analog ends up with 600 ppm of carbon.

The total amount of carbon in our planet is relatively poorly constrained. Earth’s bulk carbon may reside in the core, adding up to a total of (3 − 6) × 10^21^ kg or 500 to 1000 ppm (51). This is within range of our estimate for the carbon delivery integrated over the growth of Earth. Our Mars analog does not reach high enough temperatures to burn carbon dust in its envelope and therefore obtains a carbon fraction as high as 2000 ppm. The more massive planets (Venus, Earth, and Theia) accrete a substantial amount of carbon in their earlier growth phases, before entering a phase of burning carbon dust in the envelope after growing beyond 0.2 ME. This carbon-free accretion phase lowers the resulting carbon fraction by approximately a factor of three.

DISCUSSION

We built our model for terrestrial planet formation by pebble accretion around fundamental physical processes that have been explored within the context of modern planet formation theory. The water ice line has been demonstrated as a likely site for early planetesimal formation; the location of the primordial water ice line in the 1- to 2-AU region is a direct consequence of the elevated luminosity of the young Sun, compared to its current value, and fits broadly with the location of Mars. Taking Mars as a planet that underwent only modest growth by pebble accretion as one end-member and Venus as a planet that experienced substantial growth and migration as the other end-member, we demonstrated that the masses and orbits of these two planets are consistent with growth by accretion of millimeter-sized pebbles in a protoplanetary disc with a relatively low radial flux of solid material. The radial pebble flux needed to match the masses and orbits of the terrestrial planets corresponds to approximately 10 ME Ma^−1^ at a stellar accretion rate of 10^−8^ M⊙ year^−1^ (after 1 Ma of evolution in our model) and 1 ME Ma^−1^ when the protoplanetary disc dissipates at an accretion rate of 10^−9^ M⊙ year^−1^ after 5 Ma. These low pebble fluxes correspond well to the estimated 1 − 10 ME of solids left in protoplanetary discs at the 3-Ma to 5-Ma evolution stages (12). Although such evolved protoplanetary discs may not be able to form gas-giant planets any longer, they could presently be in their main phase of forming terrestrial planets by pebble accretion.

The mass, orbit, and composition of Mercury do not fit within this picture. Instead, we propose that Mercury formed by accretion of metallic pebbles and planetesimals outside the sublimation front of Fayalite (Fe2SiO4), which we demonstrate in Fig. 2 crosses Mercury’s current orbit after 0.5 Ma of disc evolution. This process is similar to how icy pebbles trigger the streaming instability exterior of the water ice line (27). Growth by pebble and planetesimal accretion is very rapid so close to the Sun, so the current mass of Mercury is reached within a few hundred thousand years. Mercury then ceased to accrete as the inner regions of the solar protoplanetary disc underwent depletion by disc winds (55), stranding Mercury with a dominant metallic core and a mantle whose silicates were reduced by accretion of pebbles rich in sulfur released at the troilite (FeS) sublimation front (56). Mercury could thus be the oldest of the terrestrial planets, having compiled its bulk mass within a million years of the formation of the Sun.

Within our model framework, Earth reaches 60% of its current mass at the dissipation of the protoplanetary disc after 5 Ma. This mass of the proto-Earth agrees well with constraints from matching the ^36^Ar/^8^Ar, ^20^Ne/^22^Ne, and ^36^Ar/^22^Ne ratios by drag from the hydrodynamic escape of the primordial hydrogen/helium envelope (47). We augment the mass of Earth with the introduction of an additional planet Theia (40% Earth mass) between the orbits of Earth and Mars, to later collide with Earth to form our Moon. The radial drift of pebbles naturally provides Earth and Theia with very similar compositions, even when the isotopic signature of the drifting pebbles changes with time. Hence, Earth and the Moon inherit the indistinguishable isotopic compositions of their parent planets.

Models of protoplanetary discs show that the water ice line inevitably passes interior to 1 AU as the stellar luminosity decreases with time (32). This implies that protoplanets growing in the terrestrial planet zone will accrete a substantial mass fraction of ice. However, ice is sublimated in the hot gas envelope once the protoplanet mass reaches 0.02 ME and the water then escapes back to the protoplanetary disc with the recycling flows that penetrate into the Hill sphere. Volatile delivery filtered through a hot envelope gives a good match to the water contents of Earth. The early-accreted water ice becomes incorporated into the magma ocean after Earth melts by the accretion heat, and the water will later degas as a dense vapor atmosphere upon crystallization of Earth (50), forming the first surface water masses on our planet.

The water on our Earth analog is accreted as ice in the first 2 Ma of the evolution of the protoplanetary disc, before the planetary envelope becomes hot enough to sublimate the ice from the accreted pebbles. Because water is delivered with pebbles from the early generation to Venus, Theia, and Mars as well, we predict that the primordial water delivered to Mars should have the same D/H ratio as water on Earth. The early generation of pebbles in our model is represented by the NC meteorite group, which includes the ordinary chondrites. Water in the ordinary chondrites is heavier than Earth’s water (57), although a recent study found that enstatite chondrites have a D/H ratio similar to water in Earth’s mantle (58). The D/H ratio of the ordinary chondrites could have been increased significantly by oxidation of iron on the parent body, a process that releases light hydrogen (59). The relatively dry parent bodies of the ordinary chondrites could have formed at the warm side of the water ice line after 3 Ma, when the ice line was interior of the current orbit of Venus. Some of these would later have been dynamically implanted into the asteroid belt by scattering, now seen as the S-type asteroids in the inner regions of the asteroid belt (60).

Carbon, in the form of organics and carbon dust, is more refractory than water and can reach the planetary surface for envelope temperatures up to 1100 K (where carbon dust burns). Carbon is thus delivered to the growing protoplanets at masses up to 0.2 ME, after which the envelope becomes hot enough to vaporize all the incoming carbon-bearing molecules. Our Earth analog accretes carbon at up to a time of 3.5 Ma. As with the water, the isotopic composition of the carbon will therefore reflect material in the inner regions of the solar protoplanetary disc. The overall carbon delivery filtered through the planetary envelopes gives a good match to the estimated carbon reservoirs on Earth, including a carbon component residing in the core.

Hence, both water and carbon—essential ingredients for life—may have been delivered to Earth in the form of “pebble snow” in the early phases of the growth of our planet. This implies that volatile delivery need not be a stochastic effect of a few giant impacts (51) but may instead follow predictable patterns that can be calculated from the evolution of the protoplanetary disc and the radial drift of pebbles. If millimeter-sized pebbles are the main providers of volatiles to planets in the inner regions of protoplanetary discs, then we predict that super-Earths (more massive counterparts of terrestrial planets found frequently around solar-like stars) would have the same overall volatile budgets as terrestrial planets, being unable to accrete any additional volatiles due to the intense heat in their gas envelopes. The accretion of dense atmospheres of hydrogen and helium poses an additional challenge to the habitability of super-Earths.

We believe that our pebble-driven model of terrestrial planet formation has several advantages over the classical models that use collisions between planetary embryos as the main driver of planetary growth. While the latter successfully produce planetary systems similar to the Solar System’s terrestrial planets in terms of masses and orbits, incorporating outer Solar System material into the growing planets by giant impacts is very challenging. Simulations of water delivery to terrestrial planets use collisions with wet protoplanets originating beyond 2.5 AU—the assumed ice line based on the dichotomy between the S-type asteroids, associated with the ordinary chondrites, that dominate the inner regions of the asteroid belt and the C-type asteroids, associated with the CC, that dominate the outer regions. This position of the ice line is in clear conflict with models of irradiated protoplanetary discs showing that the ice line was likely situated interior of 1 AU for most of the lifetime of the disc. Notwithstanding this discrepancy, the classical models provide the terrestrial planets with a water fraction between 0.1 and 1%. This fraction can then be considered the maximum contribution from volatile-rich outer Solar System material in the classical model—and hence, it stands in contrast to the 40% contribution needed to explain the isotopic composition of Earth based on ^48^Ca (36) and particularly to the match of the iron isotopic composition of Earth’s mantle with the CI chondrites (40). Models invoking the inward migration of Jupiter through the asteroid belt to explain the small mass of Mars yield similarly low fractions of outer Solar System volatile-rich material in the terrestrial planets (41).

As an alternative view, Earth could have formed from a local reservoir of material that carried the exact same elemental and isotopic composition as measured on Earth. The enstatite chondrites have been proposed as candidate material matching Earth’s isotopic composition, although the iron isotopes stand out as an exception (40). Such an approach nevertheless necessitates that all the terrestrial planets formed from local reservoirs of distinct composition, perhaps a remnant of the condensation sequence in the initially hot solar protoplanetary disc (61), although this view conflicts with the volatile-rich mantle of Mercury (62). However, explaining the isotopic differences between the terrestrial planets by local reservoirs highlights the problem of the similarity of Earth and the Moon, as Theia would thus be expected to have a distinct composition different from that of Earth. In the pebble-driven model, Earth and Theia naturally get very similar compositions because they accrete from the same radial flux of drifting pebbles.

MATERIALS AND METHODS

Protoplanetary disc model

We consider a protoplanetary disc model where the temperature is set by both irradiation from the host star and by viscous heating from the magnetorotational instability. The latter is active when the temperature is above 800 K (29). We therefore increase the turbulent viscosity that controls viscous heating from a background level of αv = 10^−4^ to αv = 10^−2^, relevant for fully active turbulence caused by the magnetorotational instability, over the temperature range from 800 to 1000 K. We calculate the temperature structure iteratively by solving the transcendental equation T = T(αv) (32), where αv is a function of T. We follow the database of Baraffe et al. (63) to calculate the luminosity evolution of the young Sun. We fit the Sun’s luminosity evolution as L(t) = 2L⊙(t/Ma)^−0.7^ for the times 0.5 Ma to 5 Ma, where L⊙ is the luminosity of the modern Sun. We then multiply this value by a factor exp { − 1/[t/(0.15 Ma)]} to represent the increase in stellar luminosity (including emission from accretion and from deuterium burning, but not including the luminosity of the disc itself) in the early stages of the stellar collapse phase (64); this effectively gives a peak in the protostellar luminosity of L ≈ 3L⊙ at t = 0.2 − 0.3 Ma, followed by a power law decline ending at L ≈ 0.5L⊙ at t = 5 Ma. The protoplanetary disc evolves via the α viscosity prescription, and the pebble mass flux through the disc is set to a fixed fraction ξ of the radial gas flux; this assumption of a constant flux ratio is justified by growth bottlenecks in the viscously expanding outer regions of the protoplanetary disc (5). The temporal evolution of the structure of the protoplanetary disc is shown in fig. S1.

The pebble column density Σp follows from setting the pebble flux, Fp = 2πrvrΣp, equal to ξM·* where M·⋆ is the gas accretion rate through the disc. Here vr is the radial drift speed of the pebbles, which is calculated self-consistently from the temperature and pressure gradient of the gas in the protoplanetary disc. The pebbles have a fixed size of 1 mm; this choice of pebble size is motivated by the characteristic size of chondrules, particle growth to the bouncing barrier (65), and observations of protoplanetary discs (66). We have checked that the fragmentation-limited particle size is larger than 1 mm at 1 AU for all the evolution stages of the adopted protoplanetary disc, for a turbulence strength of δ = 10^−4^ and a critical fragmentation speed of 1 m/s (67). The pebble column density for ξ = 0.0036 is indicated in fig. S1. The pebbles have Stokes numbers in the range 0.001 to 0.01 (see Fig. 3); for these small sizes, the drift speed is low enough that ξ approximately represents the ratio of the column density of pebbles relative to the gas. The pebble accretion radius is calculated based on the expressions provided in (68) and are integrated in either 3D (when the pebble scale-height is larger than the accretion radius) or 2D (when the pebble scale-height is smaller than the accretion radius). To keep the model simple, we do not decrease the pebble accretion rate when the ice layers are sublimated in the envelope at protoplanet masses above 0.02 ME.

Calculating the equilibrium structure of the envelope

We solve the structure of the planetary envelope by integrating the standard equations of hydrostatic balance and outward luminosity transport from the Hill sphere down to the center of the protoplanet. The boundary conditions at the Hill sphere are the temperature and density of the protoplanetary disc. We set the temperature gradient to the minimum of the radiative temperature gradient and the convective temperature gradient. Our consideration of the equilibrium structure of the envelope ignores the time needed to heat the material in the envelope to its equilibrium value. The accretion luminosity is generally high enough to validate this approach. The radiative temperature gradient is defined by the opacity; we assume here that the opacity follows a power law with the temperature κ = κ0(T/100 K)^0.5^, with κ0 = 0.01, 0.1, or 1.0 m^2^ kg^−1^ (69). The mean molecular weight of the gas is set to a weighted average of a solar mixture of hydrogen/helium and silicate vapor at its saturated vapor pressure.

The total luminosity generated by the accretion, the latent heat of sublimation of silicates, and the decay of ^26^Al areL=GMM·Rsil−QsilM·+L26(11)

Here, Rsil is the either the location of the silicate sublimation front (which is calculated self-consistently by equating the saturated vapor pressure of silicates with the ambient pressure) or the surface of the planet (if the temperature does not reach silicate sublimation). The latent heat of sublimation of silicate rock (Qsil) reduces the accretion energy when the temperature reaches silicate sublimation in the envelope. We nevertheless ignore the latent heat of water and silicates in the calculations because silicate sublimation only happens at planetary masses approaching 1 ME, while water sublimation is followed immediately by nucleation of the outward transport water vapor. The luminosity generated by the decay of ^26^Al is set toL26=r26E26M(12)

Here, r26 = 1.3 × 10^5^ kg^−1^ s^−1^ exp ( − t/τ) is the decay rate of ^26^Al per total silicate mass (70), t is the time, τ = 1.03 Ma is the decay constant, E26 = 3.12 MeV is the energy released per decay (71), and M is the mass of the silicate core. We assume that all the heat from ^26^Al is released in the core because the total mass of silicate vapor in the envelope is much lower than the core mass for the planetary masses that we consider.