# Contemporary formation of early Solar System planetesimals at two distinct radial locations

A. Morbidelli, K. Baillié, K. Batygin, S. Charnoz, T. Guillot, D. Rubie, T.

Kleine

► To cite this version:

A. Morbidelli, K. Baillié, K. Batygin, S. Charnoz, T. Guillot, et al. Contemporary formation of early Solar System planetesimals at two distinct radial locations. Nature Astronomy, 2022, 6 (1), pp.72-79. 10.1038/s41550-021-01517-7. hal-03553375

HAL Id: hal-03553375

https://hal.science/hal-03553375v1

Submitted on 6 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Contemporary formation of early Solar System planetesimals at two distinct radial locations

A. Morbidelli ID<fnref n="1" />✉, K. Baillié ID², K. Batygin³, S. Charnoz ID⁴, T. Guillot ID<fnref n="1" />, D. C. Rubie⁵ and T. Kleine ID⁶,⁷

The formation of planetesimals is expected to occur via particle-gas instabilities that concentrate dust into self-gravitating clumps<fnref n="1" />⁻³. Triggering these instabilities requires the prior pile-up of dust in the protoplanetary disk⁴,⁵. This has been successfully modelled exclusively at the disk's snowline⁶⁻⁹, whereas rocky planetesimals in the inner disk were only obtained by assuming either unrealistically large particle sizes¹⁰,¹¹ or an enhanced global disk metallicity¹². However, planetesimal formation solely at the snowline is difficult to reconcile with the early and contemporaneous formation of iron meteorite parent bodies with distinct oxidation states¹³,¹⁴ and isotopic compositions¹⁵, indicating formation at different radial locations in the disk. Here, by modelling the evolution of a disk with ongoing accretion of material from the collapsing molecular cloud¹⁶⁻¹⁸, we show that planetesimal formation may have been triggered within the first 0.5 million years by dust pile-up at both the snowline (at ~5 AU) and the silicate sublimation line (at ~1 AU), provided turbulent diffusion was low. Particle concentration at ~1 AU is due to the early outward radial motion of gas¹⁹ and is assisted by the sublimation and recondensation of silicates²⁰,²¹. Our results indicate that, although the planetesimals at the two locations formed about contemporaneously, those at the snowline accreted a large fraction of their mass (~60%) from materials delivered to the disk in the first few tens of thousands of years, whereas this fraction is only 30% for the planetesimals formed at the silicate line. Thus, provided that the isotopic composition of the delivered material changed with time²², these two planetesimal populations should have distinct isotopic compositions, consistent with observations¹⁵.

The goal of this work is to identify the conditions that may lead to the contemporary formation of iron meteorite parent bodies at two distinct radial locations in the disk; one of these locations has to be characterized by a higher temperature than the snowline, so as to form ice-free planetesimals. Our model is similar to that in refs. ⁹,¹⁶⁻¹⁸,²³, but comprises a viscosity parameter $\alpha$ that, instead of being held fixed, is reduced from $1 \times 10^{-2}$ to $5 \times 10^{-4}$ as the accretion rate of mass onto the disk, the disk's local temperature and the propensity to undergo gravitational instabilities decrease (Methods). With this improvement, the early viscous, radially spreading disk evolves over time towards a low-viscosity state, consistent with observations of the dust distribution in protoplanetary disks²⁴,²⁵.

In our simulations, the Sun starts with half of its current mass ($M_\odot$), consistent with a class-0 protostar, and material is delivered to the Sun-disk system at a rate decaying as $e^{-(t/0.1\text{Myr})}$ where $t$ is time.

The time-integrated infall of material brings the Sun-disk system to $1 M_\odot$ in a few hundred thousand years, a timescale comparable with that in refs. ¹⁶⁻¹⁸ but notably shorter than in ref. ⁹. The Sun is assumed to accrete the material that falls directly within 0.05 AU or is transported by the disk to within this limit. Previous work⁹,¹⁶⁻¹⁸,²³ has assumed that the angular momentum of infalling material increases rapidly with time, but modern magneto-hydrodynamical simulations highlight the importance of magnetic breaking in removing angular momentum from the infalling material²⁶. Hence, we test different parametrizations of the time-evolution of the effective distance where material falls onto the disk, known as the centrifugal radius (Methods). We find that, as long as the inflow of infalling material is vigorous, the radial velocity of the gas is positive (that is, directed away from the star) beyond the centrifugal radius, whereas, when the inflow wanes, the disk rapidly becomes an accretional disk with a negative radial velocity in its inner part. Because a positive radial velocity of the gas can help in trapping dust particles¹⁹, we look for disks that have a protracted phase of radial expansion in their inner part. Assuming a centrifugal radius decreasing as $R_c = 0.35\,\text{AU}/(M_{\text{Sun}}(t))^{0.5}$, where $M_{\text{Sun}}$ is the mass of the Sun at time $t$ relatively to its current mass, we obtain a disk that expands radially beyond 0.4 AU during the first 0.3 Myr (Fig. 1). The time-evolution of the disk temperature is also shown in Fig. 1, whereas the evolution of the surface density and viscosity are depicted in Extended Data Figs. 1 and 2.

In our model, all elements heavier than hydrogen and helium are injected into the disk together with the gas. They are assumed to be in solid form (dust) when the local temperature of the disk is below their condensation temperature $T_{\text{cond}}$. For simplicity, we first consider only two broad species: rocks ($T_{\text{cond}} = 1,400\,\text{K}$) and water ice ($T_{\text{cond}} = 170\,\text{K}$). Initially, the dust is micrometres in grain size and is transported outwards during the radial expansion of the disk, while also growing on a timescale proportional to the local dust/gas mass ratio and orbital period¹⁰ (Methods). We cap the maximal dust grain size to be 10 cm beyond the snowline and 5 mm within the snowline, in agreement with earlier studies on dust coagulation, bouncing and fragmentation²⁷. When dust drifts inwards across the snowline, we assume that the ice sublimates and the remaining 70% of the solid mass is redistributed in 5 mm grains⁶,⁷ (Methods). The diffusion of water vapour and its recondensation enhance the solid/gas density ratio at the disk's midplane beyond the snowline (Fig. 2a), as has been found previously⁷,⁸,¹¹.

Inside of the snowline, the solid particles drift towards the Sun until their radial velocity becomes positive because their entrainment in the radially expanding gas dominates over the headwind drag¹⁹. A small pressure bump also appears, and so particles pile up

<footnote n="1">¹Laboratoire Lagrange, Université Cote d'Azur, CNRS, Observatoire de la Cote d'Azur, Nice, France. ²IMCCE, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Universités, UPMC Univ Paris 06, Univ. Lille, Paris, France. ³Division of Geological and Planetary Sciences, California Institute of Technology, Pasadena, CA, USA. ⁴Université de Paris, Institut de Physique du Globe de Paris, CNRS, Paris, France. ⁵Bayerisches Geoinstitut, Universität Bayreuth, Bayreuth, Germany. ⁶Institut für Planetologie, University of Münster, Münster, Germany. ⁷Max Planck Institute for Solar System Research, Justus-von-Liebig-Weg 3, Göttingen, Germany. ✉e-mail: morby@oca.eu</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 1** | The radial distribution of the disk's temperature at different times. We assume that $R_c = 0.35\text{AU}/[M_{\text{sun}}(t)]^{0.5}$. The thick part of each curve shows the region where the radial velocity of the gas is positive (outward), whereas the thin part depicts the accretion part of the disk (negative radial velocity), as also indicated by the black and orange arrows. The horizontal dashed lines mark the condensation temperature of water ($T_{\text{cond}} = 170$ K; black) and rocks ($T_{\text{cond}} = 1,400$ K; red). The intersection of these lines with the various coloured curves identifies the location of the condensation/sublimation fronts of these elements as a function of time.</figcaption></figure>

just outside of this location (Fig. 2a). However, turbulent diffusion, characterized by the coefficient $D=\nu/Sc$, where $\nu$ is the gas viscosity and $Sc$ is the Schmidt number, smooths the radial distribution and impedes efficient settling towards the midplane of such small particles, even for $Sc=10$. Thus, a solid/gas volume-density ratio of order unity, as required to trigger the streaming instability¹⁰, is never achieved (Fig. 2a). For this to occur $Sc=100$ is needed, as ref. ⁹ has assumed, but such a large value has never been observed in hydrodynamical simulations.

The situation changes if silicate sublimation is also taken into account. We now consider three broad species: refractories ($T_{\text{cond}}=1,400$ K), silicates ($T_{\text{cond}}=1,000$ K) and water ice (Methods). We assume that at $T=1,000$ K half of the rocky mass sublimates and the grains break into millimetre-size particles of more refractory material. This change in particle size makes the radial flow of solids converge at the silicate-sublimation front, up to 0.35 Myr (Fig. 2b). Moreover, diffusion and re-condensation of the silicate vapour increase the density of solids beyond the sublimation front. Altogether, this creates a local strong enhancement of the solid/gas ratio in the disk's midplane even for $Sc=10$. When this ratio becomes larger than 0.5, we convert part of the solid density excess into planetesimals¹⁰ at each timestep (Methods).

Figure 3 shows the radial mass distributions of the planetesimal populations produced at the snowline and silicate-sublimation line as a function of time. About 4.5 Earth masses ($M_{\oplus}$) of silicate-rich planetesimals form in a ring extending from 0.75 to 0.9 AU during a time period from 0.33 to 0.38 Myr, whereas ~32 $M_{\oplus}$ of ice-rich planetesimals form beyond the snowline, from ~3 to 5.5 AU, during 0.1–0.5 Myr. Such a large mass in icy planetesimals can explain the rapid formation of Jupiter's core at the snowline¹⁵, while the concentration of rocky planetesimals in a narrow ring is needed to explain the small masses of Mercury and Mars relative to Earth and Venus²⁹. Interestingly, if we further reduce the minimum value of the viscosity parameter $\alpha$ to $1 \times 10^{-4}$ (instead of $5 \times 10^{-4}$), the total mass

of planetesimals produced at the silicate sublimation line exceeds $40 M_{\oplus}$. This planetesimal mass, although too large for the Solar System, could readily explain the formation of rocky super-Earths, which are frequently observed around other stars³⁰ but are difficult to produce starting from a uniform distribution of planetesimals throughout the disk³¹. Our model predicts that the formation of rocky planets should always be accompanied by the formation of more distant icy planets (Extended Data Fig. 3).

We now compare our results with the constraints from the meteorite record. Iron meteorites are fragments of the metallic cores of some of the oldest planetesimals of the Solar System, which formed within 1 Myr after Solar System formation (as defined by the time of formation of its first solids, inclusions rich in Ca-Al or CAIs)¹⁵. The iron meteorites can be subdivided into two isotopically distinct groups, which are termed the carbonaceous (CC) and non-carbonaceous (NC)¹⁵ groups. Of note, the parent bodies of the CC irons tend to have smaller relative core sizes and are characterized by lower Fe/Ni ratios than those of the NC irons (Methods), suggesting that the former formed in more oxidizing environments than the latter. As such, our working hypothesis is to identify the planetesimals formed at the snowline as the parent bodies of CC iron meteorites, consistent with their formation in a more oxidizing environment, and those formed at the silicate-sublimation line as the parent bodies of NC irons, consistent with the observation that they accreted at higher temperature and were water-ice free. A larger water-ice fraction in CC iron parent bodies also leads to a more protracted timescale of core formation, due to the lower concentration of heat-producing $^{26}$Al (ref. ¹⁴). This is consistent with the observed later core formation time of CC compared with iron parent bodies at ~3 Myr and ~1 Myr, respectively¹⁴,¹⁵.

Another important difference between planetesimals formed at the silicate-sublimation line and at the snowline is that our model predicts the former to have ratios of silicates (olivine + pyroxene) to refractory elements 10–35% higher than the protosolar value.
<figure><img src="image_4.png" /></figure>

**Fig. 2** | The density ratio between dust and gas on the disk's midplane as a function of heliocentric distance and time. **a,b,** Simulation where only the snowline is accounted for (a) and simulation where silicate sublimation and recondensation at T=1,000 K are also taken into account (b). The red and green curves show the location of the snowline and silicate line as a function of time. The white line separates the (r, t) domain where the radial velocity of the dust particles is negative from that where it is positive. It is different in the two panels because in **b** particles are smaller inside the silicate line and therefore are more coupled with the gas. In **a**, outward particle motion is due to the positive radial motion of the gas up to 0.29 Myr; after, there is a narrow pressure bump due to the rapid increase of viscosity with temperature that causes a drop in surface density (Extended Data Fig. 1). In **b**, the particles are too small to feel the pressure bump and drift inward with the gas for t > 0.35 Myr. Convergent migration and dust pile-up tend to occur at the outer white line but are contrasted by particle diffusion. When the dust/gas ratio (colour scale) exceeds 0.5 (yellow colour) planetesimal formation is assumed to occur.

This property results from recondensation of the gas that sublimated off refractory grains at high temperature$^{20,21,32}$. The Ni/Ir ratio inferred for most bulk NC cores is indeed larger than solar$^{13}$ (Ni condenses together with silicates, whereas Ir is refractory), but this property is not unique to NC irons$^{33,34}$. So, these data do not provide clear evidence for the formation of NC parent bodies at the silicate-sublimation line. Instead, the enhanced ratio of silicates to refractory elements predicted by our model is consistent with the suprasolar Si/Al ratios of NC chondrites, which are not observed in any CC chondrites. Chondrites are later-formed planetesimals, and so modelling their formation goes beyond the scope of this study (see Supplementary Note, section S5 for a discussion). Nevertheless, early formed planetesimals may well be their precursors via the subsequent generation of chondrules as collisional debris$^{35}$. Consequently, the chemical composition of NC chondrites may still reflect that of the first planetesimals formed at the silicate-sublimation line, modelled in this work.

The most important constraint is that of the aforementioned isotopic dichotomy between NC and CC irons$^{15}$. To test the ability of our model to satisfy this constraint, we distinguish between material accreted to the disk before and after the first 20 Kyr (denoted 'early

material' and 'late material' hereafter; Fig. 3). The choice of this time is justified in the Supplementary Note, section S1.5, and the relationship between the condensates of the early material and CAIs is discussed in Supplementary Note, section S4. We find that planetesimals formed at the snowline incorporate a larger fraction of early material than planetesimals at ~1 AU (Fig. 3). This is because early material is efficiently transported to the outer disk during the radial expansion phase, while it is substituted by late-infalling material in the inner disk. By the time the inner planetesimals form, the drift of early material back into the inner disk again raises the early-to-late material ratio at ~1 AU, but this ratio nevertheless remains below that of the outer disk (Fig. 4). Assuming that the early and late materials are isotopically distinct$^{22}$, the two populations of planetesimals produced in our model at distinct radial locations have distinct isotopic compositions, as observed for NC and CC irons$^{15}$. Moreover, the mixing ratios between early and late materials in our model are in good agreement with those derived from the isotopic offset between the NC and CC reservoirs (Methods and Extended Data Fig. 4). Finally, we note that, although the presence of a barrier against dust drift is not needed to explain the isotopic dichotomy between the two populations of early formed planetesimals
<figure><img src="image_5.png" /><figcaption>**Fig. 3** | The radial mass distribution of the planetesimal populations formed in different time intervals. Planetesimal masses formed in the same radial bin at subsequent time intervals (colour-coded) are plotted on top of each other, so that the upper border of the histogram represents the total mass ever produced. The light green shaded area shows the fraction of this mass derived from early infalling material (the material accreted onto the disk before 20 Kyr). This fraction is 46-70% for the icy planetesimals between 3.1 and 5 AU and increasing from the inside out (the mass-weighted average is 60%), and 28-30% for the rocky planetesimals just inward of 1 AU.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 4** | The relative proportions of early infalling material among the disk's refractory elements, as a function of heliocentric distance at different times. By assumption, the late-infalling material starts at $t=20$ Kyr; thus, until this time 100% of the disk is made of early infalling material. The arrival of late-infalling material progressively chases the early infalling material from the inner disk until 0.1 Myr. Then, diffusion and the inward drift of dust again raises the early/late material ratio in the inner disk, while this ratio still decreases beyond ~4 AU.</figcaption></figure>

modelled in this study, it is nevertheless needed before the disk is completely homogenized, because otherwise the NC-CC dichotomy could not be preserved for the later-formed parent bodies of

chondrites. The formation of Jupiter from the population of ice-rich planetesimals (not included in our model) would be the most obvious cause of the appearance of such barrier¹⁵.
Our model highlights the fundamental processes and properties needed to act in concert to account for the meteoritic evidence for the contemporaneous formation of two isotopically different planetary populations at distinct radial locations: (1) a small centrifugal radius for the material falling onto the disk, which is necessary to sustain a protracted radial expansion of the gas and delay the inward drift of dust particles into the Sun; (2) sublimation and recondensation of water and silicates at the respective phase-transition lines, together with stepwise changes in the maximal sizes of solid particles at each line, to enhance the local solid/gas ratio; (3) a reduced turbulent diffusion, allowing for sufficient particle pile-up and sedimentation towards the mid-plane, together with a quite large disk temperature so that the silicate sublimation line is initially near 1 AU; (4) a rapid change in isotopic composition of the material accreted onto the disk over time, to account for the radial isotopic gradient in the disk that results in the NC-CC dichotomy when planetesimals form in non-contiguous regions. Importantly, within the context of the proposed model, any derogation from (1)-(4) would lead to results inconsistent with the meteorite record (Supplementary Note, sections S1 and S2).

## Methods

**Code description.** *Structure.* Our code uses a one-dimensional grid, similar to refs. 11,23,36, describing the radial distribution of gas and dust and their properties. The grid samples a user-defined radial range (from $r_{min}=0.05$ AU to $r_{max}=100$ AU for the simulation presented in the main text) in logarithmic bins. We used 100 bins for the presented simulation, although different numbers of bins have been used in convergence tests.

*Accretion of mass onto the disk.* For numerical reasons, the disk is initialized with an arbitrarily small surface density and a temperature $T_{irr}=115$ K (r/AU)$^{1-3/7}$, corresponding to a passively irradiated disk$^{37}$. The gas is supplied at a rate

$$ \dot{M}(t) = \frac{M_{\odot} - M_{\text{Sun}}(0)}{\tau} e^{-\left[\frac{t}{\tau}\right]} \quad (1) $$

where $M_{\text{Sun}}(0)$ is its initial mass in the simulation (here $1/2 M_{\odot}$). Previous work$^{16-18}$ has assumed that mass infall rate (equation (1)) is constant but truncated the infall abruptly when the mass of the star-disk system reached $1 M_{\odot}$. We think it is more realistic to expect that the accretion rate decays over time, with a low-rate tail and no artificial truncation. We adopt the value $\tau=10^5$ yr in equation (1), which is of the order of the duration of the infall in refs. 12,16,17 (170 Kyr). The gas falling in a radial bin is$^{16-18,23}$:

$$ \dot{M}(r) = \left[ \left( 1 - \sqrt{\frac{r_-}{R_c(t)}} \right)^{\frac{1}{2}} - \left( 1 - \sqrt{\frac{r_+}{R_c(t)}} \right)^{\frac{1}{2}} \right] \dot{M}(t) \quad (2) $$

where $r_+$ and $r_-$ are the upper and lower boundary of each bin. $R_c(t)$ is called centrifugal radius or injection radius and its parametrization is given in input. In previous studies$^{16-18}$, $R_c$ has been assumed to grow as $R_c(t)=10$ AU ($M_{\text{Sun}}(t)/M_{\odot}$)$^{3}$, but in this work we test different parametrizations. The nominal simulation presented in the main text is obtained with $R_c(t)=0.35$ AU/($M_{\text{Sun}}(t)$)$^{0.5}$. The effect of this $R_c$ prescription is discussed in Supplementary Note, section S1.1.

*Computation of the disk temperature.* The midplane temperature $T$ in each ring of the disk is computed taking into account several contributions. The first is the energy released by the infalling material shocking at the surface of the considered disk ring:

$$ Q_{\text{infall}} = \frac{1}{2} \frac{GM_{\text{Sun}}(t) \dot{M}(r)}{r} \quad (3) $$

where $G$ is the gravitational constant. We take the conservative assumption that only $1/2$ of the final potential energy of the infalling gas is injected in the disk (hence the factor $1/2$ in equation (3)), the rest being lost during the infalling phase. The second contribution is the energy released by viscous heating$^{38}$:

$$ Q_{\text{visc}} = 2\pi r \delta r \frac{9}{4} \Sigma_g \nu \Omega^2 \quad (4) $$

where $\Sigma_g$ is the surface density of the gas in the ring of width $\delta r=(r_+-r_-)$, $\nu$ is the viscosity and $\Omega$ is the Keplerian frequency. The third contribution is the ring's cooling due to black body irradiation at its surfaces:

$$ Q_{-} = 2 \times 2\pi r \delta r \sigma_B T_s^4 \quad (5) $$

where $\sigma_B$ is the Stephan-Boltzman constant and $T_s$ is the temperature at the surface of the disk, related to the midplane temperature $T$ by the relationship$^{39}$:

$$ T_s^4 = \frac{4}{3} \frac{2T^4}{\kappa \Sigma_g}, \quad (6) $$

which is valid where the disk is optically thick. The opacity $\kappa$ is a function of temperature$^{39}$. The last contribution is that of energy exchange between adjacent disks' rings. A ring gains or loses energy at a rate $\delta F=F_+-F_-$, where $F_+$ ($F_-$) is the flux of energy across the boundary with the external (internal) adjacent ring$^{10}$:

$$ F = (2\pi)^{3/2} \frac{16\lambda\sigma_B}{\kappa\rho_g} \frac{dT}{dr} T^3 rH \quad (7) $$

where $\rho_g = \sigma_g/((2\pi)^{1/2}H)$ is the volume density of the gas, $H=(R/\mu T r)^{1/2}$ is the pressure scale height of the disk (R being the gas constant and $\mu$ the gas mean molecular weight), $\lambda$ is the flux-limiter$^{11}$ and all quantities are taken at the boundary between adjacent rings. For $\mu$ we assume 2.3 g mol$^{-1}$ and we take the approximation to keep this number constant across condensation lines. The quantity $(Q_{\text{infall}}+Q_{\text{visc}}-Q_{-}+\delta F)\delta t$ describes the change of internal energy of a ring over an integration timestep $\delta t$, and the change in temperature $T$ is obtained by dividing this quantity by the heat capacity of the ring $c_v=R/[(\gamma-1)\mu]$ where $\gamma=1.4$ is the adiabatic index. If the temperature falls below that of a passively irradiated disk, $T_{irr}$ we reset $T=T_{irr}$. The temperature is further modified during the advection step, described below.

*Viscosity prescription.* The viscosity $\nu$ is as usual defined as $\nu=\alpha H^2 \Omega$. The viscosity parameter $\alpha$ had been set constant and equal to $1 \times 10^{-2}$ in previous works$^{16-18}$. In this work, we change $\alpha$ (between a minimum $\alpha_{\text{min}}$ and a maximum $\alpha_{\text{max}}$) over time and radial location. We set:

$$ \alpha = \alpha_{\text{min}} + (\alpha_{\text{max}} - \alpha_{\text{min}}) \frac{\dot{M}(t)}{M(0)} \quad (8) $$

the rationale being that the infall of material onto the disk generates Reynolds stresses that act as a viscosity$^{36}$, which become weaker as the infall wanes. It is also known that at high temperature, typically above the silicate sublimation value, the disk becomes prone to ionization and to the magneto-rotational instability, which raises the turbulent viscosity significantly. Thus, for the rings with temperature $T>1,500$ K we set $\alpha=\alpha_{\text{max}}$ and for rings with $1,000$ K $< T < 1,500$ K we set $\alpha$ to an intermediate, $T$-dependent value between equation (8) and $\alpha_{\text{max}}$ computed as:

$$ f = \sin \left[ \frac{T - 1,000}{1,000} \pi \right], \quad \alpha(T) = e^{\left[ (1-f) \log \alpha + f \log a_{\text{max}} \right]} $$

Similarly, it is known that when the disk is gravitationally unstable or close to instability, the disk develops clumps and waves that also generate an effective viscosity$^{12}$. Thus, for the rings where Toomre's Q parameter$^{12}$ is less than unity (a criterion for gravitational instability) we set $\alpha=3 \times 10^{-2}$ and for rings with $1 < Q < Q_{\text{lim}}$ we set $\alpha$ to a Q-dependent value intermediate between equation (8) and $3 \times 10^{-2}$, given by:

$$ f = \sin \left[ \frac{Q-1}{2(Q_{\text{lim}}-1)} \pi \right], \quad \alpha(T) = (1-f) 3 \times 10^{-2} + f\alpha. $$

For the nominal simulation presented in the main text, we set $\alpha_{\text{max}}=1 \times 10^{-2}$, $\alpha_{\text{min}}=5 \times 10^{-4}$ and $Q_{\text{lim}}=10$. We discuss in the Supplementary Note, section S1.2, how the results change when these parameters are varied. We find that radial energy exchange (equation (7)), which previous codes have not included$^{16-18}$, is essential to stabilize the disk when $\alpha$ is allowed to change over time at different radii as in our model.

*Computation of gas evolution.* We now discuss how the surface density of gas in the disk evolves. Besides receiving mass at a rate given by equation (2), a ring can exchange material with neighbouring rings. The radial velocity of the gas at the boundary between two rings is due to the mutual viscous torques that they arise on each other due to the differential rotation and results in:

$$ v_r^g = - \frac{3}{\Sigma_g \sqrt{r}} \frac{d}{dr} (\Sigma_g \nu \sqrt{r}) \quad (9) $$

where all quantities are evaluated at the boundary. This speed is then modified to account for the back-reaction of dust onto gas, as will be discussed below. A ring gains or loses mass at a rate $\delta F_M=F_{M-}-F_{M+}$, where $F_{M+}$ ($F_{M-}$) is the flux of mass across the boundary with the external (internal) adjacent ring:

$$ F_M = 2\pi r v_r^g \Sigma_g \quad (10) $$

where $\Sigma_g$ is here the surface density in the ring that is supplying mass to the other ring and $r$ is the radial distance of the boundary between the considered
rings, where equation (9) has been computed. From equation (10), the change in surface density of a ring over an integration timestep is readily computed. For code stability, we impose that $δt/v_g^2 < 0.2δr$. As anticipated, together with the advection of mass, we also compute an advection of thermal energy associated to the flux of gases in/out of rings with different temperatures. This further modifies the ring's temperature, on top of the prescription described above.

In doing these calculations, the boundary conditions play an important role, in particular the inner one (the outer boundary being very far from the region of interest). We use open boundary conditions, assuming that the gas flowing inward of $r_{min}$ is immediately accreted by the central star along magnetic channels (we increase the stellar mass accordingly) and that the gas flowing beyond $r_{max}$ is immediately photoevaporated or stripped away by passing stars.

**Dust species and particle growth.** When the gas is supplied to the disk following equation (2), we also assume that 1% of its mass is supplied in condensable materials. We consider three types of material: ice, with a condensation temperature of $T_{cond}$=170 K; silicates, with a condensation temperature of $T_{cond}$=1,000 K; and refractories, with a condensation temperature of $T_{cond}$=1,400 K. This choice, lower than those canonically assumed for the condensation of silicates and refractory elements, is discussed in Supplementary Note, section S3. For simplicity we neglect that the condensation temperature depends on the partial pressure of the considered material. Then, we introduce three surface density functions, $Σ_{ice}$, $Σ_{sil}$ and $Σ_{ref}$, for these three materials, respectively. At injection, we assume that 30% of the condensable material is in ice (consistent with comet composition), 35% is in silicates and 35% in more refractory materials. When the temperature is larger than the corresponding condensation temperature, the material is considered to be in vapour form, while for lower temperatures it is assumed to have condensed into dust.

The dust has an initial size (diameter) of 1 µm, but then grows with a timescale

$$τ_{\text{gr}} = \frac{1}{Z\Omega} \qquad (11)$$

where Z is the local solid/gas surface density ratio (see ref. ¹⁰ for a derivation). For simplicity, we consider only one dust size in each ring, instead of a size distribution. The reason is that, in dust growth models, most of the dust's mass is concentrated in particles near the maximal dust size¹⁰. Because of this simplification, when new dust is created in a bin that already hosts partially grown dust (for instance, due to the injection of fresh material from the molecular cloud) we take the total-mass-weighted mean size between the pre-existing size and the new injected one¹⁴. The maximum dust size is set by the drifting, bouncing and fragmentation barriers. Based on previous work on the effect of these barriers⁴,⁵, we limit the maximal size of dust in the ice regime ($T<170$ K) to 10 cm, that in the silicate regime ($1,000 K > T > 170$ K) to 5 mm and that in the refractory regime ($1,400 K > T > 1,000$ K) to 1 mm. Silicate and refractory particles are thus much smaller than in refs. ¹⁰,¹¹, which is consistent with a reduced fragmentation energy⁸ than that considered in those works and the existence of a bouncing barrier¹². These sizes are also consistent at the order of magnitude with those of silicate and refractory particles observed in meteorites, such as chondrules, chondrule clusters and CAIs. The large size contrast between icy and silicate particles is due to the fact that warm ice (near the snowline, where planetesimal formation will take place) is stickier than silicates, so aggregates are expected to grow bigger¹³. The size contrast between silicate and refractory particles could be justified by fragmentation during silicate sublimation. We discuss in Supplementary Note, section S1.6, how the results change with different size contrasts. Once the dust size is set, the dust's Stokes number $S_i$ is computed from the local density of gas¹⁰. Thus, a particle drifting towards the Sun with constant size has its Stokes number progressively reduced because the disk's gas density increases.

**Evolution of the dust surface densities $Σ_{ice}$, $Σ_{sil}$ and $Σ_{ref}$.** When the corresponding material is in vapour form, we assume that its radial velocity is equal to that of the disk's gas (equation (9)). When it is in dust form, its radial velocity $v_r^d$ is computed as described in the appendix of ref. ¹⁰, which includes a modification of the gas radial velocity (equation (9)) due to the back-reaction of dust on gas. Different from ref. ¹⁰, this modification of the gas velocity affects the evolution of the gas. Thus, our model is able, in principle, to capture the so-called self-induced trap phenomenon⁴⁵. For the record, we never observe this phenomenon in our nominal simulations, because the radial velocity of the gas is positive during the planetesimal-formation stage, but we do observe it in a classic, viscous accretion disk model if the particle size is 10 cm.

In addition to the advection process, analogue to that of the gas described above (equation (10) with $Σ_d$ and $v_r^d$ instead of $Σ_g$ and $v_r^g$, where $Σ_d$ stands generically for $Σ_{ice}$, $Σ_{sil}$ or $Σ_{ref}$), the evolution of $Σ_d$ is also affected by a diffusion equation¹⁰ that, from the ring perspective followed in this code description, generates a supplementary mass gain/loss $δF_D=F_{D-}F_{D+}$, where $F_{D-}$ ($F_{D+}$) is the flux of mass due to diffusion across the boundary with the external (internal) adjacent ring:

$$F_D = 2πrD\Sigma_g \frac{d}{dr} \left( \frac{\Sigma_d}{\Sigma_g} \right) \qquad (12)$$

Here D is the diffusion coefficient and r, $Σ_g$ and the gradient of $Σ_d/Σ_g$ are evaluated at the boundary between the adjacent rings. We set $D=ν/Sc$, Sc being the Schmidt number. For a passive tracer of the gas, Sc can be⁶⁶ as large as 10, which is our nominal choice (see Supplementary Note, section S1.2, for a discussion of the effects of this number). We assume the same Sc for vapour and solid particles because²⁵ the Stokes number of our particles is always smaller than 0.1 for $r<8$ AU and $t<0.5$ Myr.

The use of a unique density function $Σ_d$ to describe both the vapour and solid phases of the same material, which just differ in advection radial velocity, is a simplified but effective way to treat the sublimation/recondensation process. If one treats vapour and dust separately, each of their respective density functions has a discontinuity, dropping to zero at the condensation/sublimation line. If, instead of assuming instantaneous sublimation/condensation as we do here for simplicity, one considers a non-zero condensation or sublimation timescale dependent on partial pressures⁸,¹¹,⁷⁷, the discontinuity becomes a gradient, but such a gradient is nevertheless very steep. Consequently, in both cases equation (12) would give a strong mass flux at the boundary between the vapour-dominated and the dust-dominated regimes. But, in reality, most of the diffusion of vapour through the condensation line is counterbalanced by the diffusion of dust in the opposite direction. If, instead, one uses a single density function for both vapour and dust, as we do here, equation (12) automatically describes the net mass flux across the sublimation/condensation boundary, which is the one that really matters. Our procedure is mathematically exact if one makes the simplifying assumption of instantaneous and complete sublimation/condensation at the critical temperature. We provide in Supplementary Note, section S1.7, a test of the expected differences in the results using the two approaches.

**Planetesimal formation.** Associated with the surface densities of gas ($Σ_g$) and dust ($Σ_d$) are the volume densities on the midplane $ρ_g = Σ_g/((2π)^{1/2}H)$ and $ρ_d = Σ_d/((2π)^{1/2}H_d)$, with $H_d=H_e[α/Sc/(α/Sc+S_J)]$¹⁵. In the calculation of $ρ_d$, we sum up the contributions of all three species: ice, silicate and refractory. Whenever in a ring $ρ_d/ρ_g > 0.5$, we assume that planetesimal formation can take place via the streaming instability in that ring²⁷. A fraction of 0.01% of the solids is converted into planetesimals per ring's orbital period¹⁰. A more elaborate prescription¹² is also tested in Supplementary Note, section S1.3. The corresponding mass is subtracted from the dust density functions in proportion to the relative abundances of the three species of dust. If a species is in vapour form, its density remains untouched. When the production of planetesimals reduces the dust/gas mass ratio below 0.5, planetesimal formation is stopped. This regulates planetesimal production and dust accumulation, keeping the dust/gas ratio typically below Unity (Fig. 2). Without planetesimal formation the ratio would increase further. Thus, the phenomenon of planetesimal formation is not very sensitive to the adopted dust/gas threshold ratio, although the resulting total mass of planetesimals does increase (decrease) if the adopted threshold is decreased (increased).

**Isotopic composition.** To study the isotopic composition of planetesimals we compute the evolution of populations of dust tracers in the disk. We follow the idea proposed in ref. ²¹, according to which the isotopic dichotomy between NC and CC planetesimals is due to the injection into the disk of isotopically distinct materials at different times. Thus, we define a switch time $t_{dich}$ and we split $Σ_{ref}$ into two functions, $Σ_{ref}^{(1)}$ and $Σ_{ref}^{(2)}$, describing the surface density of refractory material injected at $t<t_{dich}$ or at $t>t_{dich}$ respectively. The two distributions are referred to as 'early' and 'late infalling material' in the main text. Both $Σ_{ref}^{(1)}$ and $Σ_{ref}^{(2)}$ undergo the diffusion process described by equation (12), which induces their mutual mixing, as illustrated in Fig. 4. When a planetesimal forms at time t, its composition in terms of early versus late material is given by the $Σ_{ref}^{(1)}/Σ_{ref}^{(2)}$ ratio at that time and location in the disk. Averaging the compositions of planetesimals formed in the same radial bin at different times (mass-weighted average) produces Fig. 3.

**Determining the sizes of cores of NC and CC iron meteorite parent bodies.** The core sizes of iron meteorite parent bodies can be estimated using the abundances of highly siderophile elements (HSEs) inferred for the bulk cores. Owing to their strong siderophile character, the HSEs quantitatively partitioned into the core. Thus, the core mass fraction can be calculated by dividing the HSE concentrations of the bulk body (assumed to be chondritic) by the HSE concentrations of the bulk core. To this end, the HSE concentrations of the bulk cores are inferred by modelling fractional crystallization. For the refractory HSEs rhenium, osmium, iridium, ruthenium and platinum the resulting relative ratios are typically broadly chondritic, such that for each element similar core mass fractions are calculated. The HSE concentration data used for calculating core mass fractions are summarized in Supplementary Table 1 and Supplementary Table 2 summarizes the mean core mass fractions for each iron meteorite parent body. For the CC irons, we assumed that before core formation the bulk body had the composition of a CI-chondrite. This is the most appropriate composition, given that these bodies formed at the snow line and, therefore, incorporated water ice. For the NC irons we used an average ordinary (OC) or enstatite chondrite (EC) composition. However, using the same starting compositions for both NC and CC irons does not change the resulting core mass fractions by much, except that it would lead to overlapping values for the IIC and IVA irons.
The calculations reveal overall larger core mass fractions in NC compared with CC iron meteorite parent bodies, where NC cores were typically ~20% of the mass of the parent body, and CC cores were <15% (Supplementary Table 2). The smaller core sizes of the CC parent bodies are consistent with larger water-ice fractions in these bodies, which results in oxidation and, hence, a smaller fraction of iron partitioned into the core. This is also consistent with the systematically lower Fe/ Ni ratios inferred for CC cores compared with NC cores¹³. These observations are difficult to reconcile with a model¹⁶ where the parent bodies of both NC and CC iron meteorites would have formed at the snowline at different times. By contrast, these observations are fully consistent with our model in which NC and CC bodies formed in rocky and icy environments, respectively.

The existence of water ice in CC iron meteorite parent bodies drastically changes the thermal evolution models usually used to infer the accretion time of planetesimals from their measured differentiation times. Most thermal modelling studies have assumed the same composition for NC and CC iron meteorite parent bodies, resulting in a monotonic relationship between accretion time and differentiation time, which in turn implied that CC iron meteorite parent bodies accreted later than their NC counterparts¹⁴. However, a more recent study has shown that water ice delays the onset of melting and core formation¹³ and that, therefore, the later core formation time of CC iron parent bodies does not imply later accretion. In fact, once the effect of different water ice fractions is taken into account, the inferred accretion times of CC and NC iron meteorite parent bodies are indistinguishable, where both groups of bodies are constrained to have formed within the first 1 Myr of the Solar System¹³, in line with our predictions.

**Determining the mixing ratios of distinct materials using the isotopic properties of CAIs, CC and NC meteorites.** For all elements that display the NC-CC isotopic dichotomy, the CC reservoir is always between the isotopic compositions of CAIs and NC meteorites¹⁹. This observation has led to the proposal that the NC-CC dichotomy reflects different mixing proportions of two isotopically distinct disk reservoirs, which were characterized by similar, broadly chondritic bulk chemical compositions²¹,⁴⁸. One of these reservoirs is characterized by a CAI-like isotopic composition (termed IC for inclusion-like chondritic reservoir¹⁸) and corresponds to the early infalling material. Note that although this material has a CAI-like isotopic composition, its chemical composition is distinct from CAIs and instead is assumed to be chondritic. This assumption stems from the observation that the NC-CC dichotomy exists for refractory (for example, molybdenum, titanium) and non-refractory elements (for example, chromium, nickel) and that for all these elements the CC reservoir is always isotopically intermediate between NC and CAIs. The other reservoir is characterized by an NC-like isotopic composition, but its exact isotopic composition, termed NCᵢ, is not known. Within this framework, the isotopic compositions of the CC and NC reservoirs can be expressed as simple binary mixtures between IC and NCᵢ material as follows:

$$CC = x \times IC + (1 - x) \times NC_i$$

$$NC = y \times IC + (1 - y) \times NC_i$$

where x and y denote the fractions of early material in the CC and NC reservoir, respectively. These two parameters cannot be calculated independently, because the isotopic composition of the late infall, NCᵢ, is not known. However, the two equations can be combined by eliminating NCᵢ, so that x can be calculated as a function of y as follows:

$$x = \frac{y(\text{CC} - \text{IC}) + \text{NC} - \text{CC}}{\text{NC} - \text{IC}}$$

The NC-CC dichotomy is best defined for titanium, chromium and molybdenum isotopes, which therefore are most suitable to calculate the dependence of x and y. As the $ε^{50}$Ti and $ε^{54}$Cr isotope anomalies among NC meteorites are correlated, and because this correlation points toward the composition of the CC reservoir, using $ε^{50}$Ti or $ε^{54}$Cr returns the same results. For titanium we used the following values: $ε^{50}$Ti$_{IC}$ = +9, which is the average titanium isotope anomaly of CAIs⁴⁹; $ε^{50}$Ti$_{NC}$ = -1, which is the average $ε^{50}$Ti of NC meteorites, or its most extreme negative value $ε^{50}$Ti$_{NC}$ = -2; the Ti isotopic composition of CI chondrites, $ε^{50}$Ti$_{CC}$ = +2, which best represents the composition of the outer disk⁴⁸. Using different values for $ε^{50}$Ti$_{NC}$ or $ε^{50}$Ti$_{CC}$ within the compositional range of the NC and CC reservoirs does not change the result significantly. For molybdenum we use the characteristic molybdenum isotopic difference between the CC and NC reservoirs, which can be expressed as $Δ^{95}$Mo (see ref. ⁵⁰), with the following values: $Δ^{95}$Mo$_{IC}$ = +125 (ref. ⁴⁹); $Δ^{95}$Mo$_{CC}$ = +26 (ref. ⁵⁰); and $Δ^{95}$Mo$_{NC}$ = -9 (ref. ⁵⁰).

The relationship between the fraction of early material in the CC and NC reservoirs x and y, respectively, calculated using the titanium and molybdenum isotope anomalies are shown in Extended Data Fig. 4 together with the proportions predicted by our model for $t_{dich}$ = 20 Kyr (see above). This comparison shows that our model can reproduce the early-to-late material ratios in the NC and CC reservoirs quite well.

## Data availability

The compiled code, the input file and the ascii output files of our reference simulation including silicate condensation/sublimation (one file per output timestep ($10^4$yr) for a total of 100 files) are provided at: lagrange.oca.gov/images/LAGRANGE/pages_perso/morby/forNature.tar.gz. A readme file describes the content of each file.

## Code availability

The code for the calculation of the disk evolution is available on request from the corresponding author.

Received: 2 March 2021; Accepted: 17 September 2021; Published online: 22 December 2021

## References

1. Youdin, A. N. & Goodman, J. Streaming instabilities in protoplanetary disks. *Astrophys. J.* **620**, 459–469 (2005).

2. Cuzzi, J. N., Hogan, R. C. & Shariff, K. Toward planetesimals: dense chondrule clumps in the protoplanetary nebula. *Astrophys. J.* **687**, 1432–1447 (2008).

3. Simon, J. B., Armitage, P. J., Li, R. & Youdin, A. N. The mass and size distribution of planetesimals formed by the streaming instability. I. The role of self-gravity. *Astrophys. J.* **822**, 55 (2016).

4. Drażkowska, J. & Dullemond, C. P. Can dust coagulation trigger streaming instability? *Astron. Astrophys.* **572**, A78 (2014).

5. Yang, C.-C., Johansen, A. & Carrera, D. Concentrating small particles in protoplanetary disks through the streaming instability. *Astron. Astrophys.* **606**, A80 (2017).

6. Ida, S. & Guillot, T. Formation of dust-rich planetesimals from sublimated pebbles inside of the snow line. *Astron. Astrophys.* **596**, L3 (2016).

7. Schoonenberg, D. & Ormel, C. W. Planetesimal formation near the snowline: in or out? *Astron. Astrophys.* **602**, A21 (2017).

8. Drażkowska, J. & Alibert, Y. Planetesimal formation starts at the snow line. *Astron. Astrophys.* **608**, A92 (2017).

9. Drażkowska, J. & Dullemond, C. P. Planetesimal formation during protoplanetary disk buildup. *Astron. Astrophys.* **614**, A62 (2018).

10. Drażkowska, J., Alibert, Y. & Moore, B. Close-in planetesimal formation by pile-up of drifting pebbles. *Astron. Astrophys.* **594**, A105 (2016).

11. Charnoz, S. et al. Planetesimal formation in an evolving protoplanetary disk with a dead zone. *Astron. Astrophys.* **627**, A50 (2019).

12. Schoonenberg, D., Ormel, C. W. & Krijt, S. A Lagrangian model for dust evolution in protoplanetary disks: formation of wet and dry planetesimals at different stellar masses. *Astron. Astrophys.* **620**, A134 (2018).

13. Chabot, N. L. Composition of metallic cores in the early solar system. *49th Lunar and Planetary Science Conference* 2083 (ed. LPI) (2018).

14. Spitzer, F., Burkhardt, C., Nimmo, F. & Kleine, T. Nucleosynthetic Pt isotope anomalies and the Hf-W chronology of core formation in inner and outer solar system planetesimals. *Earth Planet. Sci. Lett.* **576**, 117211 (2021).

15. Kruijer, T. S., Burkhardt, C., Budde, G. & Kleine, T. Age of Jupiter inferred from the distinct genetics and formation times of meteorites. *Proc. Natl Acad. Sci. USA* **114**, 6712–6716 (2017).

16. Hueso, R. & Guillot, T. Evolution of protoplanetary disks: constraints from DM Tauri and GM Aurigae. *Astron. Astrophys.* **442**, 703–725 (2005).

17. Pignatale, F. C., Charnoz, S., Chaussidon, M. & Jacquet, E. Making the planetary material diversity during the early assembling of the solar system. *Astrophys. J.* **867**, L23 (2018).

18. Baillié, K., Marques, J. & Piau, L. Building protoplanetary disks from the molecular cloud: redefining the disk timeline. *Astron. Astrophys.* **624**, A93 (2019).

19. Batygin, K. & Morbidelli, A. Formation of giant planet satellites. *Astrophys. J.* **894**, 143 (2020).

20. Aguichine, A., Mousis, O., Devouard, B. & Ronnet, T. Rocklines as cradles for refractory solids in the protosolar nebula. *Astrophys. J.* **901**, 97 (2020).

21. Miyazaki, Y. & Korenaga, J. Dynamic evolution of major element chemistry in protoplanetary disks and its implications for Earth-enstatite chondrite connection. *Icarus* **361**, 114368 (2021).

22. Nanne, J. A. M., Nimmo, F., Cuzzi, J. N. & Kleine, T. Origin of the non-carbonaceous-carbonaceous meteorite dichotomy. *Earth Planet. Sci. Lett.* **511**, 44–54 (2019).

23. Appelgren, J., Lambrechts, M. & Johansen, A. Dust clearing by radial drift in evolving protoplanetary discs. *Astron. Astrophys.* **638**, A156 (2020).

24. Pinte, C. et al. Dust and gas in the disk of HL Tauri: surface density, dust settling, and dust-to-gas ratio. *Astrophys. J.* **816**, 25 (2016).

25. Dullemond, C. P. et al. The Disk Substructures at High Angular Resolution Project (DSHARP). VI. Dust trapping in thin-ringed protoplanetary disks. *Astrophys. J.* **869**, L46 (2018).

26. Lee, Y. N., Charnoz, S. & Hennebelle, P. Protoplanetary disk formation from the collapse of a prestellar core. *Astron. Astrophys.* **648**, A101 (2021).
27. Güttler, C., Blum, J., Zsom, A., Ormel, C. W. & Dullemond, C. P. The outcome of protoplanetary dust growth: pebbles, boulders, or planetesimals? I. Mapping the zoo of laboratory collision experiments. *Astron. Astrophys.* **513**, https://doi.org/10.1051/0004-6361/200912852 (2010).

28. Gole, D. A., Simon, J. B., Li, R., Youdin, A. N. & Armitage, P. J. Turbulence regulates the rate of planetesimal formation via gravitational collapse. *Astrophys. J.* **904**, 132 (2020).

29. Hansen, B. M. S. Formation of the terrestrial planets from a narrow annulus. *Astrophys. J.* **703**, 1131–1140 (2009).

30. Owen, J. E. & Wu, Y. The evaporation valley in the Kepler planets. *Astrophys. J.* **847**, 29 (2017).

31. Izidoro, A. et al. Formation of planetary systems by pebble accretion and migration: hot super-Earth systems from breaking compact resonant chains. *Astron. Astrophys.* **650**, A 152 (2021).

32. Morbidelli, A., Libourel, G., Palme, H., Jacobson, S. A. & Rubie, D. C. Subsolar Al/Si and Mg/Si ratios of non-carbonaceous chondrites reveal planetesimal formation during early condensation in the protoplanetary disk. *Earth Planet. Sci. Lett.* **538**, 116220 (2020).

33. Tornabene, H. A., Ash, R. D. & Walker, R. J. New insights to the genetics, formation and crystallization history of group IC iron meteorites. *52nd Lunar Planetary Science Conference* 2548. LPI (2021).

34. Hilton, C. D., Ash, R. D. & Walker, R. J. Crystallization histories of the group IIF iron meteorites and Eagle Station pallasites. *Meteorit. Planet. Sci.* **55**, 2570–2586 (2020).

35. Stewart, S. T. et al. 2019. Collapsing impact vapor plume model for chondrule and chondrite formation. *50th Annual Lunar and Planetary Science Conference*. LPI. (2019).

36. Lichtenberg, T., Dąrżkowska, J., Schönbächer, M., Golabek, G. J. & Hand, T. O. Bifurcation of planetary building blocks during Solar System formation. *Science* **371**, 365–370 (2020).

37. Chiang, E. & Youdin, A. N. Forming planetesimals in solar and extrasolar nebulae. *Ann. Rev. Earth Planet. Sci.* **38**, 493–522 (2010).

38. Dullemond, C. *Theoretical Models of the Structure of Protoplanetary Disks*. Lecture notes. Heidelberg Univ. (2013).

39. Bitsch, B., Morbidelli, A., Lega, E. & Crida, A. Stellar irradiated discs and implications on migration of embedded planets. II. Accreting-discs. *Astron. Astrophys.* **564**, A135 (2014).

40. Kley, W., Bitsch, B. & Klahr, H. Planet migration in three-dimensional radiative discs. *Astron. Astrophys.* **506**, 971–987 (2009).

41. Kley, W. Radiation hydrodynamics of the boundary layer in accretion disks. I. Numerical methods. *Astron. Astrophys.* **208**, 98–110 (1989).

42. Rafikov, R. R. Viscosity prescription for gravitationally unstable accretion disks. *Astrophys. J.* **804**, 62 (2015).

43. Birnstiel, T., Fang, M. & Johansen, A. Dust evolution and the formation of planetesimals. *Space Sci. Rev.* **205**, 41–75 (2016).

44. Ormel, C. W. An atmospheric structure equation for grain growth. *Astrophys. J.* **789**, L18 (2014).

45. Gonzalez, J.-F., Laibe, G. & Maddison, S. T. Self-induced dust traps: overcoming planet formation barriers. *Mon. Not. R. Astron. Soc.* **467**, 1984–1996 (2017).

46. Carballido, A., Stone, J. M. & Pringle, J. E. Diffusion coefficient of a passive contaminant in a local MHD model of a turbulent accretion disc. *Mon. Not. R. Astron. Soc.* **358**, 1055–1060 (2005).

47. Hyodo, R., Ida, S. & Charnoz, S. Formation of rocky and icy planetesimals inside and outside the snow line: effects of diffusion, sublimation, and back-reaction. *Astron. Astrophys.* **629**, A90 (2019).

48. Burkhardt, C., Dauphas, N., Hans, U., Bourdon, B. & Kleine, T. Elemental and isotopic variability in solar system materials by mixing and processing of primordial disk reservoirs. *Geochim. Cosmochim. Acta* **261**, 145–170 (2019).

49. Brennecka, G. A. et al. Astronomical context of Solar System formation from molybdenum isotopes in meteorite inclusions. *Science* **370**, 837–840 (2020).

50. Budde, G., Burkhardt, C. & Kleine, T. Molybdenum isotopic evidence for the late accretion of outer Solar System material to Earth. *Nat. Astron.* **3**, 736–741 (2019).

## Acknowledgements

A.M. and S.C. acknowledge support from programme ANR-20-CE49-0006 (ANR DISKBUILD). The work presented here has been performed in preparation for the proposal HolyEarth by A.M. and T.K., which has been funded by the European Research Council (grant No. 101019380). The authors thank R. Deienno and C. Ormel for constructive and detailed comments.

## Author contributions

A.M. conceived the project, wrote the code, ran the simulations and led the writing of the manuscript. K. Baillié wrote an earlier version of the code. K. Baillié, S.C. and T.G. contributed with their experience on disk evolution. K. Batygin stressed the importance of the radial expansion of the disk. D.C.R. and T.K. provided their experience on the chemical and isotopic composition of meteorites, which allowed for testing the model against measured constraints. All authors contributed to writing the manuscript and discussing the significance of the results.

## Competing interests

The authors declare no competing interests.

## Additional information

*Extended data* is available for this paper at https://doi.org/10.1038/s41550-021-01517-7.

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41550-021-01517-7.

**Correspondence and requests for materials** should be addressed to A. Morbidelli.

**Peer review information** *Nature Astronomy* thanks Rogerio Dienno, Chris Ormel and the other, anonymous, reviewer(s) for their contribution to the peer review of this work.

**Reprints and permissions information** is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© The Author(s), under exclusive licence to Springer Nature Limited 2021
<figure><img src="image_7.png" /><figcaption>**Extended Data Fig. 1** | Surface density of the disk. Surface density of the disk as a function of heliocentric distance at different times.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Extended Data Fig. 2 | Turbulent parameter α. Turbulent parameter α, for the nominal simulation presented in the main text.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>mass in rocky planetesimals (Earth mass)</figcaption></figure>

**Extended Data Fig. 3** | total masses of rocky and icy planetesimals. total masses of rocky and icy planetesimals for 4 values of αmin and two values of Qlim.
<figure><img src="image_10.png" /></figure>

**Extended Data Fig. 4** | Fraction of early material in CC and NC according to isotopic constraints. Relation between fraction of early material in CC and NC as given by Ti and Mo isotope anomalies in meteorites. The thick solid line assumes the average value for NC meteorites $\varepsilon^{50}\text{Ti}_{\text{NC}} = -1$ while the thin line assumes $\varepsilon^{50}\text{Ti}_{\text{NC}} = -2$ (that is the extreme value observed in NC). Orange-shaded area indicates the predicted fractions of our model: 0.275-0.3 for NC planetesimals and 0.450-0.70 for CC planetesimals.