Bifurcation of planetary building blocks during Solar System formation
Tim Lichtenberg, Joanna Drazkowska, Maria Schönbächler, Gregor J. Golabek, Thomas O. Hands
BBifurcation of planetary building blocks during Solar System formation T IM L ICHTENBERG , J OANNA D RA ¸ ˙ ZKOWSKA , M ARIA S CH ¨ ONB ¨ ACHLER , G REGOR
J. G
OLABEK , AND T HOMAS
O. H
ANDS Atmospheric, Oceanic and Planetary Physics, Department of Physics, University of Oxford, United Kingdom University Observatory, Faculty of Physics, Ludwig-Maximilians-Universit¨at M¨unchen, Germany Institute for Geochemistry and Petrology, Department of Earth Sciences, ETH Zurich, Switzerland Bayerisches Geoinstitut, University of Bayreuth, Germany Institute for Computational Science, University of Zurich, Switzerland
Submitted 2 March 2020; Accepted 10 December 2020; Published 21 January 2021 (Science 371, 6527)
ABSTRACTGeochemical and astronomical evidence demonstrate that planet formation occurred in two spatially andtemporally separated reservoirs. The origin of this dichotomy is unknown. We use numerical models toinvestigate how the evolution of the solar protoplanetary disk influenced the timing of protoplanet formationand their internal evolution. Migration of the water snow line can generate two distinct bursts of plan-etesimal formation that sample different source regions. These reservoirs evolve in divergent geophys-ical modes and develop distinct volatile contents, consistent with constraints from accretion chronology,thermo-chemistry, and the mass divergence of inner and outer Solar System. Our simulations suggest thatthe compositional fractionation and isotopic dichotomy of the Solar System was initiated by the interplaybetween disk dynamics, heterogeneous accretion, and internal evolution of forming protoplanets.Planetary systems, including the Solar System, form byaccretion from a protoplanetary disk of gas and dust. As-tronomical observations of these disks provide evidencefor rapid dust coagulation (Segura-Cox et al. 2020), showringed substructure (Andrews et al. 2018), and indicate adecrease in total dust mass with disk age (Ansdell et al.2016), to below the total masses in fully-assembled exo-planetary systems (Najita & Kenyon 2014). This suggeststhat planet formation starts early. At the onset of plane-tary accretion, small dust grains coagulate to form largeraggregates (pebbles). These pebbles drift inward underaerodynamic drag and gravitationally collapse when theyreach sufficient local over-density (Johansen et al. 2015), which preferentially produces birth planetesimals of ≈ (cid:38) – yr,once they have grown to larger embryos through mutualcollisions (Liu et al. 2019).Meteorites record planet formation in our own Solar Sys-tem and constrain the astronomical timescales: radiomet-ric dating of meteorites suggests iron core formation (dif-ferentiation) in planetesimals by ≈ Corresponding author: Tim [email protected] of Ca,Al-rich inclusions (CAIs) (Kruijer et al. 2014), the old-est known solids that formed together with the proto-Sun(Connelly et al. 2012). Core formation in first-generation,birth planetesimals is driven by internal radiogenic heatingfrom Al (half-life of ∼ × yr). Therefore, the time in-terval between planetesimal formation and differentiation( (cid:38) – yr) increases with later formation time, point-ing to the formation of the earliest planetesimals in the in-ner Solar System (cid:46) structure, meteorite data indicate spatial heterogeneity inthe isotopic composition of planetary materials ( Supple-mentary Materials ): studies on Ti, Cr, Mo, and other iso-tope systems show variabilities across individual plane-tary bodies from different orbits (Leya et al. 2008; Trinquieret al. 2009; Warren 2011). Combined with temporal con-straints on the accretion process, these divide the outerand inner Solar System reservoirs into chronologically andspatially distinct populations, also known as the carbona-ceous chondrite (CC) and non-carbonaceous (NC) reser-voirs, after their representative meteorite classes. Theisotopic signatures record the heterogeneous distributionof presolar dust from multiple sites of stellar nucleosyn-thesis, and trace transport and mixing processes during a r X i v : . [ a s t r o - ph . E P ] J a n Lichtenberg et al. planet formation. In contrast, the elemental abundancesand internal structure of planetary bodies were modifiedthrough geophysical evolution during planetary formation(Lichtenberg et al. 2019a). Earth’s nucleosynthetic iso-tope composition is NC-like and disparate from outer SolarSystem materials (
Supplementary Materials ), but its de-pletion pattern of moderately volatile elements is close toCCs (Braukm ¨uller et al. 2019). This has been attributedto impact delivery (Sch ¨onb ¨achler et al. 2010) or additionof pebbles (Schiller et al. 2018). The inner Solar Sys-tem planets are depleted in highly volatile elements, suchas hydrogen, carbon, and nitrogen, which influenced theavailability of surface water and the composition of theiratmospheres. This has been linked with the accretion ofthe terrestrial planets inside the snow line, where waterice is not stable during the disk phase, and later deliveryof volatile-rich, outer Solar System materials (Peslier et al.2018; Raymond et al. 2020).The observed CC/NC dichotomy has been suggestedto result from an early formation of proto-Jupiter (cid:46)
Supplementary Mate-rials ).FORMATION OF TWO DISTINCT PLANETESIMALPOPULATIONSMultiple planetesimal populations can be formed duringthe build-up and evolution of the solar protoplanetary disk(Dra¸ ˙zkowska & Dullemond 2018). We test whether suchmulti-stage planetesimal formation is consistent with theaccretion chronology as inferred from meteorites. We usethe output of previous simulations (Dra¸ ˙zkowska & Dulle- mond 2018, see
Supplementary Materials ) to explore theCC/NC dichotomy. The simulations cover the initial infalland build-up of gas and dust, during which the disk itselfaccretes from the surrounding molecular cloud (Class I),and include the drift, coagulation and compositional evolu-tion of the dust, and viscous evolution of the gas during thelater evolutionary stage of the disk (Class II). In these sim-ulations, the initial material accretes onto the innermost re-gions and continuously feeds more distant orbits throughviscous expansion of the disk. The snow line moves out-ward during the Class I stage because a dense, compactdisk forms, which heats up viscously. With the onset ofthe Class II stage, the gas density decreases, the temper-ature structure becomes dominated by stellar irradiation,and the snow line moves inward again.
Orbital distance, ( ) T i m e , () Terrestrial planetmigration corridor . . . . C l a ss II C l a ss I Watersnow line
33 32 31 30 29 28 27 26 ln - - - - - - Planetesimal formation rate, / ln Figure 1. Formation of two distinct planetesimal popula-tions in the disk simulation.
The rate of planetesimal formation(change in surface density per timestep) is shown on logarithmicscale in cgs and natural units. The solid purple line indicates theorbital migration of the snow line, red and blue areas highlight theformation regions of the two planetesimal populations (ReservoirI and Reservoir II). The dashed gray lines (labeled ’Terrestrialplanet migration corridor’, cf. Fig. S10) bound the range of possi-ble planet migration scenarios consistent with the present-day or-bits of the terrestrial planets. While smaller embryos move slowlyinward, more rapidly accreting (and thus more massive) embryoswould migrate faster toward the proto-Sun. The light gray arrowson the left side indicate the transition from the Class I to Class IIdisk stage.
The snow line location affects the redistribution of dustgrains. These are initially well-mixed with the disk gas,but radially drift due to coagulation and aerodynamic dragfrom the ambient gas. Inward-drifting, icy pebbles un-dergo rapid dehydration and size reduction at the snowline, which reduces their drift velocity and causes a pile-up of solid material. Outward diffusion of water vapor addi-tionally leads to its recondensation onto icy grains beyondthe snow line (the cold-finger effect), locally increasing thedensity of solids. When the conditions for the formation ofdense, gravitationally unstable dust filaments (the stream-ing instability) are met in the simulation, planetesimals areformed. Planetesimal formation in the simulation prefer-entially occurs around the snow line. At 5 Myr into thesimulation, which is the approximate upper lifetime of thesolar nebula (
Supplementary Materials ), we assume thatthe gas disk dissipates and planetesimal formation halts.In wind-driven disks with low levels of midplane turbu-lence the global angular momentum transport is domi-nated by near-surface layers (
Supplementary Materials ). ifurcation of planetary building blocks In this scenario, the outward-inward tacking snow linegenerates two distinct episodes of planetesimal forma-tion in different orbital regions and time intervals (Fig. 1).The first, early-formed planetesimal population (hereafterReservoir I) is triggered by the cold-finger effect between ≈ ≈ ≈ ≈ Supplementary Materials ). The build-upof both reservoirs depends on the local evolution of thepebble flux, but in turn influences the coagulation and driftof pebbles toward the inner disk (Figs. S1–S3). We there-fore assess the feedback between reservoir formation andcontinuing planetary accretion.HETEROGENEOUS ACCRETION AND RESERVOIRSEPARATIONWe aim to determine the dominant mode of solid masstransfer in this disk build-up scenario, and to evaluate itsconsistency with the timescales of planetary accretion de-rived from radiometric ages (Rudge et al. 2010; Dauphas& Pourmand 2011; Kruijer et al. 2014). To relate the nu-merical simulation to geochemical chronology, we equatetime zero in the disk model with the time of CAI forma-tion (Connelly et al. 2012,
Supplementary Materials ). Forfixed orbits at 2 and 15 au, approximately representa-tive of Reservoir I and II in the inner and outer disk, weevaluate the anticipated growth timescales for two sce- narios of protoplanet accretion (
Supplementary Materials ).We first consider collisional growth due to planetesimal-planetesimal interactions, during which the largest plan-etesimals grow by accreting smaller bodies. Second, weinvestigate the efficacy of pebble accretion, which is drivenby the abundance of dust grains that cross the proto-planet orbit. Consistent with dynamical evidence from theasteroid-belt (Delbo et al. 2017) and numerical simula-tions of the streaming instability (Johansen et al. 2015),we assume a maximum planetesimal radius of 300 km,and that collisional growth is dominated by planetesimalsof 50 km in radius. The orbit- and time-dependent evolu-tion of the pebble flux (i.e., the transfer of solid mass in thedisk), pebble size, and the gas disk characteristics are self-consistently derived from the disk and coagulation simu-
Figure 2. Pebble flux and planetesimal growth timescalesin the disk simulation. ( A ) Pebble flux at 2 (red) and 15 (blue)au over time. During the disk build-up stage, the pebble flux isdominated by outward-moving dust (dotted lines), after which thepebbles start drifting inward. Reservoir I and II in the simulationprogressively diverge in pebble flux by more than one order ofmagnitude due to accretion of Reservoir II (black arrows, labeled’Reservoir separation’). ( B ) Comparison of growth timescalesfor birth-sized planetesimals of 300 km radius by either pebbleor planetesimal (collisional) accretion ( Supplementary Materials )within the approximate lifetime of the solar protoplanetary disk(gray horizontal band, labeled ’ ≈ disk lifetime’). The red shadedarea in the bottom indicates the time interval during which peb-ble accretion is more effective than collisional growth (labeled’Pebble-aided growth’, cf. Fig. S11). lation described above (Dra¸ ˙zkowska & Dullemond 2018, Supplementary Materials ).Under these conditions, the pebble flux is dominated byoutward-moving dust grains (Fig. 2A) between ≈ ≈ Lichtenberg et al.
Planetesimal radius, (km) F o r m a t i on t i m e a ft e r C A I f o r m a t i on ( M y r) A
151 436 720 1005 1290 1574
Highest mean temperature, < > (K) H y d r o t he r m a l a c t i v i t y C o r e f o r m a t i on Rain-outC1 C2 B H2H1Percolation Hydrous rockdecompositionWater icemelting T decomp T hydr φ rain T perc Figure 3. Simulated thermochemical evolution of planetesimals. ( A ) Highest mean internal temperature (color bar), interpolatedfrom simulated planetesimal radii and formation times (dots). Reservoir I and II in the disk model are indicated by red and blue dots,respectively (cf. Fig. 1). Green and yellow lines (labeled H1–C2) indicate bodies for which 50 vol% (dashed lines) or 90 vol% (solidlines) of their interiors undergo hydrothermal activity (H1, H2) and iron core formation (C1, C2). ( B ) Schematic illustration of the sce-narios indicated by the lines in panel A at peak heating and the threshold parameters (see main text and Supplementary Materials ) todistinguish between regimes. H1, H2: Dark and light blue represent water ice and liquid water, respectively. Silicates are colored yellowto orange, indicating increasing temperature. Dashed black lines indicate the thresholds for water ice melting ( T hydr ) and hydrous rockdecomposition ( T decomp ). The light blue arrows indicate hydrothermal circulation (H1) or degassing of volatiles (H2). C1, C2: Dark grayrepresents Fe,Ni–S metals. Silicates are colored yellow to red, indicating increasing temperature (yellow to orange) and progressivemelting (red). Dashed black lines indicate the thresholds for metal percolation ( T perc ) and rain-out ( ϕ rain ). Dark gray arrows indicate coreformation by percolation (C1) or rain-out (C2). Dark red dashed arrows (C2) indicate convection in the internal magma ocean. the disk transitions into the Class II stage (Fig. 1). Sub-sequently, pebbles substantially contribute to planetesimalgrowth for a time interval of ≈ ≈ Supplemen- tary Materials ). This is because a substantial fraction ofsolids, which drift from the outer toward the inner disk, aretrapped in Reservoir II, where the pebbles gravitationallycollapse to form planetesimals starting from ≈ DIVERGENT GEOPHYSICAL EVOLUTION OFPLANETESIMAL POPULATIONSThese simulation results suggest that both planetesi-mals and pebbles contribute to planetary accretion. Rockyprotoplanets that form from the planetesimal populationsin both reservoirs inherit geochemical and geophysical properties from their precursor bodies (Lichtenberg et al.2019a). The two planetesimal reservoirs differ substan-tially in accretion time and thus in radiogenic, internal heat-ing from the decay of short-lived Al. Meteorite datahave shown that Al heating in early-formed Solar Sys-tem planetesimals drove their internal evolution, and thusaltered their structure and bulk composition (Sutton et al.2017; Alexander et al. 2018,
Supplementary Materials ).Therefore, we determine the internal and compositionalevolution of the Reservoir I and II planetesimal populationsformed in the disk simulation. We constrain the timing ofiron core formation and hydrothermal activity (aqueous al-teration and degassing) using geodynamic simulations ofthe thermochemical evolution of planetesimals. Our nu- ifurcation of planetary building blocks merical setup assumes isolated planetesimals that forminstantaneously and evolve according to a fluid mechani-cal model (Lichtenberg et al. 2019a, 2018, SupplementaryMaterials ).We perform 700 single-planetesimal simulations, span-ning the parameter range of planetesimal formation time, t form ∈ [0.1, 3.0] Myr, and planetesimal radius, R P ∈ [1,300] km (Fig. 3A). To quantify their compositional evolu-tion, we employ four thermochemical criteria (Fig. 3B) todetermine the onset and cessation of core formation andhydrothermal activity ( Supplementary Materials ). Core for-mation in planetesimals initiates with the percolation ofFe,Ni–S liquids before silicate melting occurs, and termi-nates once the interior is sufficiently melted to form an in-ternal magma ocean. At this stage metal droplets rain outfrom the convecting magma flow, producing complete dif-ferentiation between the metal core and liquid silicate man-tle (Lichtenberg et al. 2018). Quantitatively, we assumethat core formation occurs in regions with temperatures T higher than the threshold for metal-sulfide percolation, T perc ≡ K, and completes for silicate melt fractions ϕ above the rheological transition, ϕ rain ≡ . , where rocksstart to behave like a liquid rather than a solid. As temper-ature limits on hydrothermal (water-rock) activity, we as-sume water ice melts above T hydr ≡ K, and hydrousrock phases fully decompose above T decomp ≡ K.Averaged over the planetesimal volume, Fig. 3A indi-cates that larger and earlier-formed planetesimals expe-rience higher temperatures during their evolution, leadingto a greater degree of internal processing. Compared toReservoir II, Reservoir I planetesimals reach systemati-cally higher temperatures and undergo body-wide metal-silicate segregation. Reservoir I bodies larger than about30 km in radius also reach near-complete dehydration.Reservoir II planetesimals only experience limited degreesof heating and therefore less core formation and dehydra-tion. The peak temperatures (Fig. 3) only reveal approx-imate trends in the simulations, as heating and composi- tional evolution are highly variable on timescales of – yr (Figs. S4–S7).To compare the simulations with the meteorite record,we explore the time-dependent evolution of the planetesi-mal populations in both reservoirs. We assume the plan-etesimal number distributions follow that of the streaminginstability ( Supplementary Materials ) in the interval R P ∈ [1, 300] km, and take into account the generation of newplanetesimals in the simulation. We then evolve the birthplanetesimal populations in Reservoir I and II and derivethe fraction of each reservoir that falls within the thresholdcriteria, normalized to the final reservoir population after 5Myr in the disk simulation. A Metal-silicate segregation times in meteoritesCarbonaceous reservoir (CC)Non-carbonaceous reservoir (NC)
IC IIAB IIIABIVAIIIE IABIIIFIIDIIF IICIVB
Time after CAI formation, (Myr) F r a c t i on o f f i na l p l ane t e s i m a l popu l a t i on ( v o l % ) B Simulation
Reservoir Ipercolation Reservoir Irain-out Res IIpercolationRes IIrain-out
Differentiated materialTotal, Rate, / Figure 4. Comparison of metal-silicate separation timesin simulated planetesimal populations with the meteoriterecord. ( A ) Metal-silicate separation times in NC and CC me-teorite classes (Tab. S1). ( B ) Timing of core formation in Reser-voir I (red) or Reservoir II (blue). Solid lines with colored shadingrepresent the fraction of material undergoing metal-silicate sep-aration (cf. Fig. 3 and Supplementary Materials ). Dashed linesshow the total volumetric fraction for each scenario over time.Light and dark red/blue indicate percolation or rain-out, respec-tively.
Fig. 4 compares the radiometric ages of metal-silicateseparation in meteoritic measurements with the thermo-chemical evolution in our simulations. Reservoir I plan-etesimals undergo core formation between ≈ ≈ ≈ by multiple stages of percolation and rain-out, are brack-eted by these time intervals. Reservoir II planetesimalsundergo core formation at later stages, starting at ≈ ≈ ≈ ≈ Al. The timing of core formation in both simu-lated populations agree with the meteorite data within theuncertainties of the individual ages (Kruijer et al. 2014,2017; Hunt et al. 2018). The spread in early-formed NCs isreproduced within the uncertainties, with the exception of
Lichtenberg et al.
Time after CAI formation, (Myr) F r a c t i on o f f i na l p l ane t e s i m a l popu l a t i on ( v o l % ) Reservoir I Reservoir II B Simulation
Water phase in planetesimalsHydrous rock, < Water ice, < A Aqueous alteration times in meteoritesCarbonaceous reservoir (CC)Non-carbonaceous reservoir (NC) OC Tagish Lake
CMCVCR CI CO
Figure 5. Comparison of hydrothermal activity in simulatedplanetesimal populations with aqueous alteration times inthe meteorite record. ( A ) Aqueous alteration times in NC andCC meteorite classes (Tab. S1). ( B ) Fraction of planetary ma-terials that retain primordial water ice (dashed lines) or undergohydrothermal activity and retain water in hydrous silicate phases(solid lines with colored shading) in the simulations. the IAB meteorite group (Hunt et al. 2018). The fractionalvolumes in both simulated reservoirs differ, with ≈ (cid:38)
90 vol%.Fig. 5 compares aqueous alteration ages from mete-orites with our simulations. Simulated planetesimals inReservoir I experience a brief phase of hydrothermal ac-tivity between ≈ ≈ ≈
10 vol% of rock contains hydrous silicate phases, and (cid:46) Al, experience protracted hydrothermal activitylasting for several Myr with a peak at ≈ ≈ ≈ ≈
15 vol% water iceand ≈
75 vol% hydrous rock. The peak for hydrothermalactivity in Reservoir II reproduces the clustering of aque-ous alteration in the CC meteorite record. The single avail-able age for the NC population (an ordinary chondrite, OC)does not coincide with the peak in the Reservoir I popula-tion, which we discuss below.
COMPOSITIONAL DICHOTOMY BETWEEN INNERAND OUTER SOLAR SYSTEMOur simulations indicate that the formation of spatiallyand temporally distinct planetesimal populations resultedin differing evolutionary pathways for planet formation inthe inner and outer Solar System. The initial seed-population of the inner Solar System planets formed dur-ing the infall stage, while the outer Solar System planetes-imal population started to form later, at the beginning ofthe Class II disk stage. Fig. 6 shows our interpretationof the isotopic evolution during early Solar System forma-tion. Measurements of CAIs indicate that the earliest SolarSystem infall material was initially dominated by dust en-riched in supernovae-derived isotopes, then transitionedto a more depleted isotopic composition (Desch et al.2018; Nanne et al. 2019; Jacquet et al. 2019). The dis-tribution of isotopes evolved with time because of viscousdisk spreading and increasing specific angular momentumof infalling matter. As a result, the outer disk became pro-gressively dominated by early-infalling, enriched material,whereas the inner disk became dominated by late-infalling,depleted material (Desch et al. 2018; Nanne et al. 2019;Jacquet et al. 2019). We suggest that the distinct NCand CC reservoirs measured in the meteorite record re-sult from the combined influence of early planetesimal for-mation in the inner Solar System and the subsequent dustpile-up at the snow line. The latter effect suppressed therapid inward drift of enriched material and prevented iso-topic homogenization of the disk. It also led to the onsetof planetesimal formation and planetary accretion in theouter Solar System, which further limited mixing by dustdrift from the outer to the inner reservoir.Fig. 6 also illustrates our proposed timeline of com-positional evolution and accretion, based on our simula-tions. The first planetesimal population formed due to thecold-finger effect. These bodies incorporated a substan-tial amount of Al and thus dehydrated and differentiatedrapidly. Initial growth in the inner Solar System proceeded from mutual collisions among planetesimals to an interme-diate phase of pebble accretion. During the final phase be-fore disk dissipation, the pebble flux in the inner Solar Sys-tem diminished and pebble growth may have stalled (seebelow). The outer Solar System population began to formlater, during the Class II disk stage, and incorporated moretotal dust mass into its planetesimals. The initial phase ofcollisional growth was then succeeded by pebble accre-tion. This series of growth stages proceeded faster than inthe inner Solar System due to the higher local pebble flux.This chronology is consistent with meteoritic and astro-nomical evidence and makes several potentially testablepredictions. The rapid decrease in the pebble flux in theinner Solar System (Reservoir I), to less than one Earth ifurcation of planetary building blocks Enrichedinfall
EnrichedDepleted
Supernovae-derived isotopevariability
Isotopic evolution
CC NC volatile-richvolatile-poorWatersnow line H O-poor gasCollisionalgrowth Pebbleaccretion + Migration Dust dehydrationDisk spreadingProto-Sun
Time
Radiogenicdifferentiation+ dehydrationDepletedinfall
Growth chronology
Dust inward drift H O-rich gas Icy planetesimalformation
Figure 6. Schematic illustration of our proposed chronology of early Solar System accretion.
Nucleosynthetic isotope variability(left) across the disk due to varying composition of infall material is retained by the pile-up of inward-drifting dust grains at the snow line.The formation of two distinct planetesimal populations initiates divergent evolutionary pathways of inner and outer Solar System (right)due to the secular variation of local material composition, internal radiogenic heating, and dominant mode of planetary growth. mass per Myr, reduces the growth by pebble accretion dur-ing later disk stages. Parameter space exploration (
Sup-plementary Materials ) and previous work (Levison et al.2015) indicate that under such conditions pebble accre-tion onto inner Solar System planetary embryos is inef-ficient. This implies that the intermediate, pebble-driven accretion phase stalled, and subsequent growth of theterrestrial protoplanets was driven by collisional interac-tions during the final disk stages and afterwards. Suchprotracted accretion timescales for the terrestrial planetsare supported by radiometric evidence for both proto-Earthand proto-Mars (Dauphas & Pourmand 2011; Rudge et al.2010), and the absence of rocky planets larger than Earth(super-Earths) in the inner Solar System. The formation ofsuper-Earth exoplanets has been attributed to a high peb-ble flux and rapid inward migration with increasing planetmass (Raymond et al. 2020; Bitsch et al. 2019b; Liu et al.2019,
Supplementary Materials ). The secular transitionof an early, low pebble flux to a brief period of high peb-ble flux could amplify preceding mass differences betweenaccreting protoplanets and thus explain why Earth is the inner Solar System body with the closest composition toCI-chondrites (Trinquier et al. 2009; Sch ¨onb ¨achler et al.2010; Schiller et al. 2018; Braukm ¨uller et al. 2019).The secular transition between growth regimes, the in-ternal evolution of planetesimals, and the temporal varia-tion of the local pebble composition would lead to a hetero- geneous delivery of water and other highly volatile com-pounds to the inner Solar System (
Supplementary Ma-terials ): the birth planetesimals in Reservoir I form aninitial seed-population of water-rich bodies that subse-quently dehydrate from Al heating and accrete via col-lisions. With progressing time and planetary growth, theprotoplanets in the inner Solar System experience an in-flux of dry pebbles between ≈ Lichtenberg et al. ≈ (cid:38) ≈ et al. 2008) and abundance of presolar interstellar mate-rials (Davidson et al. 2014) in carbonaceous chondrites,and the scatter in cometary D/H ratios (Alexander et al.2018).In summary, the chronology of Solar System forma-tion (Fig. 6) we infer from our simulations links severalcharacteristics found by geochemical laboratory analysesand astronomical observations. We interpret the chemical(volatile-poor versus volatile-rich) and isotopic (NC versusCC) dichotomy as causally linked by the build-up of dis-tinct planetesimal populations in the inner and outer So-lar System. Mixing by dust drift between them is limitedto the earliest disk phases and declining with time due tothe progressive accretion of the outer Solar System plan- etary population. The temporal and spatial variation in themain volatile reservoirs and heterogeneous planet growthare intrinsically coupled to the disk structure, the redistri-bution and fractionation of volatile ices, the abundance ofshort-lived radionuclides, and the geophysical evolution ofprotoplanets during planetary formation. Acknowledgements — We thank C. P. Dullemond, A. C.Hunt, I. Pascucci, S. Ida, J. Wade, S.-J. Paardekooper,T. Birnstiel, S. M. Stammler, J. J. Barnes, A. Morbidelli, W.Kley, and members of the ERC EXOCONDENSE projectat Oxford for discussions; B. Liu ( 刘 倍 贝 ), J. F. J. Bryson,M. Ek, R. D. Alexander, and S. Charnoz for comments onearlier draft versions; T. V. Gerya for usage of the I2ELVIS code family; and C. P. Dullemond for usage of the diskevolution code. F
UNDING : T.L. received funding from theSimons Foundation (SCOL award
UTHOR CONTRIBUTIONS :Conceptualization, T.L., J.D., M.S., G.J.G.; Methodology,T.L., J.D., T.O.H., G.J.G.; Software, T.L., J.D., T.O.H.,G.J.G.; Validation, all authors; Formal Analysis, T.L., J.D.,T.O.H.; Investigation, T.L., J.D., M.S.; Resources, T.L.,J.D., M.S.; Writing, T.L., J.D., T.O.H., M.S.; Review & Edit-ing, all authors; Visualization, T.L.; C
OMPETING INTER - ESTS : We declare no competing interests. D
ATA AND MA - TERIALS AVAILABILITY : Simulation codes and output dataare available at osf.io/e2kfv. ifurcation of planetary building blocks S1 SUPPLEMENTARY MATERIALS
MATERIALS & METHODS
DISK EVOLUTION & PLANETESIMAL FORMATION
Planetesimals, the gravitationally bound seeds of the ac-cretion process and building blocks of planets, are thoughtto form in a multi-stage process with dust first growingto pebble sizes, which can then be concentrated by thestreaming instability (Johansen et al. 2014). The stream-ing instability is a two-fluid mechanism driven by the rela-tive flow of gas and dust which leads to the formation ofdense dust filaments in circumstellar disks. Under certainconditions, those filaments may become massive enoughto collapse under their own gravity, leading to rapid plan-etesimal formation (Youdin & Goodman 2005; Johansenet al. 2007). Planetesimals may be formed throughout theprotoplanetary disk as long as both gas and pebbles are present to create localized overdensities that breach the physical conditions for gravitational collapse. To calculatethe timing and integrated mass of planetesimals formedvia the streaming instability, we consider a physical modelof the nebular gas disk, dust evolution including growthand fragmentation (Birnstiel et al. 2016), and a prescrip-tion for turning pebbles into planetesimals that parameter-izes the effects of gravitational collapse and planetesimalformation as obtained from fluid dynamical models.We make use of previously-published simulations(Dra¸ ˙zkowska & Dullemond 2018) that include a one di-mensional model in which the protoplanetary disk isformed from the parent molecular cloud and undergoesviscous evolution. These simulations follow the rotatinginfalling cloud model (Hueso & Guillot 2005; Dullemondet al. 2006a,b). In this model there are essentially threecomponents: the parent molecular cloud, which rotates asa solid body, the central star, and the disk, both of whichare growing, fed by the infalling cloud. We use this sim-ulation to follow the gas and dust evolution from the infall(Class I) to the Class II disk stage, for which there is obser- vational evidence for declining dust mass with time (Ans-dell et al. 2016; Pascucci et al. 2016; Andrews et al. 2018;Tychoniec et al. 2020). The disk simulation assumes aninitial mass of 1 solar mass (M (cid:12) ), an isotropic tempera-ture of 10 K, and a rotation rate of × − s − for theparent molecular cloud of the Solar System. After about × yr this cloud forms a single central star surroundedby a circumstellar disk with a peak total mass of 0.2 M (cid:12) (Fig. S1). The temperature of the disk is calculated takinginto account heating due to viscosity and irradiation by thecentral star. The evolution of the gas surface density Σ g (Fig. S2) is computed using ∂ Σ g ∂t + 1 r ∂∂r ( r Σ g v r ) = S g , (S1) M a ss , () JupiterEarthMarsMoon
Proto-Sun H/He gasH O gas DustTotalplanetesimalsReservoir IReservoir II
Figure S1. Mass evolution for different reservoirs in the disksimulation (Dra¸ ˙zkowska & Dullemond 2018). The mass of theproto-Sun (yellow line), H/He gas (green), H O gas (blue), anddust (red) grow from the surrounding molecular cloud until theend of the Class I disk stage at ≈ where r denotes the radial distance to the central star, t the time, v r is the radial velocity of the gas, and S g rep-resents a source term due to matter infalling onto the diskfrom the molecular cloud. The viscous evolution of gas fol-lows the standard α -formalism (Shakura & Sunyaev 1973),where the gas viscosity is defined as ν = α v c s H g , (S2)where c s is the sound speed, H g is the gas scale height,and α v is the disk viscosity, parameterizing the phys-ical mechanism of angular momentum transport of thedisk. The viscosity is fixed to α v = 10 − . We refer to (Dra¸ ˙zkowska & Dullemond 2018) for the effects of a vary-ing α v on dust coagulation and planetesimal formation.The viscosity determines the radial gas velocity, which iscalculated from v r = − g √ r ∂∂r (cid:0) Σ g ν √ r (cid:1) . (S3)The infall of gas is complemented by the delivery of dust(Fig. S3), with a dust-to-gas ratio of 1%, as found in the in-terstellar medium (Ferri `ere 2001). We follow the evolutionof dust surface density Σ d using the advection-diffusionequation, ∂ Σ d ∂t + 1 r ∂∂r (cid:20) r (cid:18) Σ d ¯ v − ν Σ g ∂∂r (cid:20) Σ d Σ g (cid:21)(cid:19)(cid:21) = S d , (S4) Lichtenberg et al. ( g c m ) A Surface density of disk gas ( g c m ) B Surface density of dust , ( g c m ) C Surface density of H O in gas Time (Myr)0.20.3 0.72.5 5.0 , ( g c m ) D Surface density of H O in dust ( g c m ) E Surface density of planetesimals
Orbital distance, (au) , ( K ) F Midplane gas temperature
Figure S2. Radial distribution and time evolution of the composition in the disk simulation (Dra¸ ˙zkowska & Dullemond 2018).Shown are the surface densities of H/He gas (A), dust (B), H O gas (C), H O ice (D), planetesimals (E), and midplane gas temperatureat 0.2, 0.3, 0.7, 2.5, and 5.0 Myr (light to dark blue). ifurcation of planetary building blocks S3 where S d is the source term for the µ m-sized dust infallingfrom the molecular cloud and ¯ v is the mass-weighted av-eraged radial speed of dust grains, which depends on dustsize. Dust and gas are assumed to be initially well-mixed.The dust grains start to decouple from the gas as theygrow into larger (pebble) sizes. Dust grain growth is com-puted adopting an algorithm (Birnstiel et al. 2012) that re-lies on only tracking two populations of dust grains: thesmall and the large ones. Dust size is regulated by theinitial growth phase and halted either by fragmentation orradial drift.Alongside dust growth, ice sublimation and water vaporre-condensation are incorporated (Ciesla & Cuzzi 2006).The infalling dust is assumed to consist of 50% refrac-tory material and 50% water ice. When the ice subli-mates, its mass is added to the water vapor reservoir. Forthe gas disk content, the surface density of water vapourand hydrogen-helium is tracked separately. When the wa-ter vapor exceeds the equilibrium vapor pressure value,it condenses onto the existing grains, adding to their wa-ter ice content. In contrast, when the water vapor pres-sure is lower than the equilibrium pressure, sublimationtakes place. This way, the position of the snow line can bemeasured as the location where the ice content of grainschanges rapidly. The composition of grains alters theirsticking properties such that ice-rich grains are more stickycompared to refractory, ice-free grains (Wada et al. 2009,2011; Aumatell & Wurm 2014). Thus, the fragmentationthreshold velocity for dry aggregates inside of the snowline is set to v frag , in = 1 m s − and for icy aggregatesoutside of the snow line to v frag , out = 10 m s − .Planetesimal formation is included by assuming thatplanetesimals form through the streaming instability whenthe midplane dust-to-gas ratio of pebbles, which are char-acterized by the Stokes number St ≥ − , exceeds unity.Whenever this criterion is fulfilled, part of the surface den-sity of pebbles is transferred into planetesimals. For thisto take place, a dense midplane layer of sufficiently large pebbles is required. In the disk model, these large peb-bles only grow outside of the water snow line, where thedust is more sticky. Enhancements of the vertically inte-grated dust-to-gas ratio to trigger planetesimal formationhave been shown to build-up in the inner part of the pro-toplanetary disk during long-term dust evolution and massinflux from the outer parts of the disk (Dra¸ ˙zkowska et al.2016). The adopted disk model does not include a phys-ical mechanism for gas disk dispersal, such as internalor external photoevaporation (Ercolano & Pascucci 2017),which can affect the total mass of planetesimals that areformed during the final stages of disk evolution. How-ever, geochemical (Bollard et al. 2017) and paleomagnetic(Wang et al. 2017) proxies indicate a dispersal of the solar nebula on the order of 4–5 Myr, which is when we assumethe disk to disperse and planetesimal formation to cease.We focus on the physical processes taking place aroundthe snow line, which facilitates dust-to-gas overdensi-ties as the dry dust inside of the snow line reachessmaller sizes and thus drifts more slowly than the icydust outside, which leads to a pile-up effect (Ida & Guillot2016; Dra¸ ˙zkowska & Alibert 2017; Schoonenberg & Ormel2017). However, the build-up of this enhancement simi-larly requires mass transfer from the outer to the inner diskon a timescale of (cid:38) yr, depending on the assumeddisk parameters. Thus, planetesimal formation triggeredby this mechanism typically takes place in the Class IIdisk stage (Dra¸ ˙zkowska & Dullemond 2018). An additionalmechanism that operates on shorter timescales, but alsoleads to less pronounced enhancement of dust density, isthe cold-finger effect (Stevenson & Lunine 1988; Cuzzi & Zahnle 2004). This relies on the outward diffusion of wa- ter vapor from inside and its subsequent re-condensationonto the grains outside of the snow line. This mecha-nism produces a moderate enhancement of the verticallyintegrated dust-to-gas ratio during the protoplanetary diskbuild-up stage, and depends on the level of midplane gasturbulence.In the adopted inside-out infall model, the snow line firstmoves outward during the infall phase, reaching its fur-thest location when the disk obtains its peak mass (atabout 0.7 Myr) resulting from an increase in stellar lumi-nosity and gas density during disk build-up (Fig. S2A). Thetemperature is calculated taking into account both the irra-diation from the central star and dissipation of viscous en-ergy by turbulence. Stellar irradiation sets a shallow tem-perature profile in the whole disk, while viscous heatingincreases the temperature dominantly in the inner parts,which leads to a steeper temperature profile close to thestar compared to a profile with only irradiation (Fig. S2F).In the Class II phase, snow line movement is driven bothby the evolution of the slowly cooling disk, and by dust evolution. The position of the snow line is self-consistentlycomputed as the location at which the dust grains containless than 1% ice component, thus the snow line evolutionis an outcome rather than an assumption of the simulation.The pebble pile-up outside of the snow line (Figs. S2B andD) increases the time it takes the solid ice component ofpebbles to evaporate while they drift inward. As a result,the snow line position moves inward faster than it wouldjust from the temperature evolution.The position of the snow line is related to the temper-ature structure, which in the inner disk is dominated bythe dissipation of viscous energy by turbulence. From as-tronomical observations of other protoplanetary disks, thelevel of turbulence, and thus the major physical mecha- Lichtenberg et al. ( M y r) A Total dust flux ( M y r) B Inward dust flux Time (Myr)0.20.30.72.55.0 ( M y r) C Outward dust flux ( c m ) D Dust grain size
Orbital distance, (au) S t ( non - d i m . ) E Stokes number
Figure S3. Radial distribution and time evolution of additional dust parameters in the disk simulation (Dra¸ ˙zkowska & Dullemond2018). Shown are total dust flux (A), inward dust flux (B), outward dust flux (C), dust grain size (D), and the Stokes number (E) at 0.2,0.3, 0.7, 2.5, and 5.0 Myr (light to dark green). ifurcation of planetary building blocks S5 nism of angular momentum transport, remains unclear.Some observational attempts to measure turbulence havebeen inconclusive (Flaherty et al. 2015; Teague et al.2016). In a standard model of protoplanetary disk evo-lution (Balbus & Hawley 1991), the angular momentumtransfer is driven by isotropic turbulence, such that the ra-tio of midplane turbulence ( α t ) equals the disk viscosity,hence α v = α t . However, large regions of the protoplan-etary disk, mostly surrounding its midplane, are expectedto be free from turbulence (Gammie 1996). This is sup-ported by observational evidence for turbulence levels thatare lower than anticipated (Flaherty et al. 2017; Teagueet al. 2019; Flaherty et al. 2020). The disk model we em-ploy (Dra¸ ˙zkowska & Dullemond 2018) mimics a disk struc-ture in which the midplane turbulence is lower than that ex-pected from the global angular momentum transport rateby decoupling the ratio of midplane turbulence, α t , fromthe global viscous α -parameter, α v . This is a simple mim-icking of a potentially much more complex structure of anon-ideal magneto-hydrodynamic disk, in which large re-gions of the midplane are free from turbulence and the an-gular momentum can be removed vertically by magneticwinds or radially by laminar torques (Lesur et al. 2014; Baiet al. 2016; Bai 2017). Similar approaches were used inother recent models (Carrera et al. 2017; Ercolano et al.2017). A parameter study of α v and α t space has beenperformed for this model (Dra¸ ˙zkowska & Dullemond 2018),which concluded that early planetesimal formation in theinfall stage is only possible if α t (cid:28) α v because planetes-imal formation from the cold-finger effect is favoured bynear-laminar midplane conditions. We focus on the sce-nario with α t = 10 − and α v = 10 − (Dra¸ ˙zkowska &Dullemond 2018). However, our results would be quali-tatively similar for any model forming planetesimals bothduring the infall and the disk stage.Decoupling the α v and α t values aims to mimic the lay-ered accretion scenario in which the gas flow takes placein active layers above a laminar midplane. However, the disk model is one dimensional and vertically integrated.Thus, α v describes the density-averaged gas flow and isused to calculate the disk temperature in the viscous heat-ing regime. This leads to intermittent high temperaturesin the inner disk and pushes the water ice line to radii be-yond 10 au for a short time interval. Calculating the mid-plane temperature is generally a complex radiative transferproblem (Turner et al. 2014) and numerical models vary intheir conclusions on the midplane temperature in the innerparts of the disk (Flock et al. 2013; Schobert et al. 2019),associated dust redistribution and composition, and henceplanetesimal formation rates (Sato et al. 2016; Charnozet al. 2019; Ida et al. 2019). THERMOCHEMICAL EVOLUTION OF PLANETESIMALS
Once planetesimals are formed in the disk model, wefollow their internal evolution using published methodology(Lichtenberg et al. 2018, 2019a). We consider single plan-etesimals of radius R P and instantaneous formation time t form . The formation time defines the planetesimal internalheat budget assuming a disk-wide initial homogeneous ra-tio of ( Al/ Al)
CAI ≡ . × − (Kita et al. 2013) at thetime of CAI formation (Connelly et al. 2012), which definesour time zero.Our numerical model of the thermochemical evolutionof planetesimals builds on the fluid mechanical framework I2ELVIS (Gerya & Yuen 2003, 2007; Gerya 2019). We nu-merically solve the fluid dynamic conservation equationsusing the extended Boussinesq approximation (Gerya2019) in a two-dimensional infinite cylinder geometry ona Cartesian grid, namely the continuity, Stokes, Poisson,and energy conservation equations, ∂ρ∂t + ∇ ρ v = 0 , (S5) ∇ σ (cid:48) − ∇ P + ρ g = 0 , (S6) ∇ Φ = 4 πGρ, (S7) ρc P (cid:18) ∂T∂t + v i · ∇ T (cid:19) = − ∂q i ∂x i + H r + H s + H L , (S8)with density ρ , flow velocity v , deviatoric stress tensor σ (cid:48) ,pressure P , directional gravity g , gravitational potential Φ ,Newton’s constant G , heat capacity c P , temperature T ,heat flux q i = − k∂T /∂x i , thermal conductivity k , and ra-diogenic ( H r ), shear ( H s ) and latent ( H L ) heat productionterms. Employing a Lagrangian marker-in-cell techniqueto integrate the energy equation, we minimise numericaldiffusion and trace advection of non-diffusive flow proper-ties during material deformation (Gerya 2019). The com-putational stencil is formulated using a fully staggered-gridfinite-differences method, which captures sharp variations of stresses and thermal gradients in models with stronglyvariable viscosity and thermal conductivity.The thermal disk gas boundary conditions of the plan-etesimals are assumed to be constant in time, T disk ≡ K. Radiogenic heating from Al decays with time, H r ( t ) = f Al · (cid:18) Al Al (cid:19) CAI · E Al τ Al · e − ∆ t CAI /τ , (S9)with the chondritic abundance of aluminum f Al (Lodders2003), the decay energy E Al = 3 . MeV (Castillo-Rogez et al. 2009), the time relative to CAI formation ∆ t CAI , and the Al mean lifetime τ Al = 1 . Myr (Kitaet al. 2013). The thermal consequences of magma oceanstages beyond the rock disaggregation threshold (Costaet al. 2009) are parameterized using a scaled thermal con- Lichtenberg et al.
Planetesimal radius, (km) F o r m a t i on t i m e a ft e r C A I f o r m a t i on ( M y r) A Reservoir IIplanetesimals not formed
Reservoir Iplanetesimals
Retained hydrous phases, , (vol%) Time after CAI formation, (Myr) R e t a i ned h y d r ou s pha s e s , ( v o l % ) B Planetesimal radius , formation time (R1) 10 km, 0.3 Myr(R2) 30 km, 0.3 Myr(R3) 300 km, 0.3 Myr (T3) 100 km, 1.43 Myr(T2) 100 km, 0.72 Myr(T1) 100 km, 0.3 Myr
Figure S4. Parameter space and time evolution of retained hydrous rock phases in planetesimal simulations . ( A ) Analogous toFig. 3A, but for the retained volumetric fraction of hydrous rock phases in simulated planetesimals after 5 Myr (thermochemical criterion T < T decomp ). Gray dots (labeled ’not formed’) are used for the color interpolation, but do not contribute to Reservoir I or II. ( B )Time evolution for distinct combinations of formation time and planetesimal radius. Lines R1–R3 and T1–T3 indicate the evolution ofplanetesimal simulations labeled in panel A . ductivity, k eff = ( q conv / . / · α liq gc p − Si / (∆ T ρ s η num ) , (S10)with the convective heat flux q conv , the temperature differ-ence across nodes ∆ T , silicate density ρ s , thermal expan-sivity of molten silicates α liq , silicate heat capacity c p − Si ,local gravity g ( x, y ) , and lower cut-off viscosity η num . Nu-merical rock properties and the I2ELVIS treatment of rockmelting have been described previously (Lichtenberg et al.2016).To estimate the compositional evolution of planetesimalsthat form from an initial mixture of silicate rock, iron metal,and water ice, we follow an updated version of the pro-cedure described in (Lichtenberg et al. 2019a; Monteux et al. 2018) and derive the fraction of each planetesimalin our simulation grid that exceed certain thermochemicalcriteria. We employ four criteria, two related to water icemelting and water-rock reactions, and two related to metalcore formation. During radiogenic heat-up of a water-ice-rich planetesimal, the primordial ice melts and may reactwith the ambient rock to create hydrous silicates and oxi-dize available iron phases, resulting in loss of hydrogen tospace (Castillo-Rogez & Young 2017; Sutton et al. 2017).While in liquid phase, the aqueous fluid may undergopore water convection along a down-temperature gradient(Grimm & McSween 1993; Young et al. 1999), but a frac-tion of water can be trapped in hydrous silicate phases.Therefore, we derive approximate upper and lower ther- mal bounds of the thermochemical interior evolution bycalculating the fraction of the planetesimal body that ex-ceeds certain threshold temperatures for ice melting/rockhydration at T ≥ T hydr ≡
273 K, and dehydration and de-composition of the most stable hydrous silicate phases at T ≥ T decomp ≡ the host body and subsequent gravitational drainage (per-colation) toward the center, which we fix to the Fe-FeSeutectic temperature at 1 bar, T ≥ T perc ≡ ϕ (cid:38) . (Costa et al. 2009),which is equivalent to ≈ ϕ ≥ ϕ rain ≡ . , we assume closureof core formation due to efficient iron droplet rain-out fromthe internal magma ocean (Lichtenberg et al. 2018; Elkins-Tanton 2012; Patoˇcka et al. 2020).For each timestep in the disk evolution model we thencompute the fraction f P , crit of computational markers inthe planetesimal body that breach the aforementioned ifurcation of planetary building blocks S7 Planetesimal radius, (km) F o r m a t i on t i m e a ft e r C A I f o r m a t i on ( M y r) A Reservoir IIplanetesimals not formed
Reservoir Iplanetesimals
Retained water ice, , (vol%) Time after CAI formation, (Myr) R e t a i ned w a t e r i c e , ( v o l % ) B Planetesimal radius , formation time (R1) 3 km, 0.3 Myr(R2) 10 km, 0.3 Myr(R3) 300 km, 0.3 Myr (T3) 100 km, 2.5 Myr(T2) 100 km, 0.72 Myr(T1) 100 km, 0.3 Myr
Figure S5. Parameter space and time evolution of retained water ice in planetesimal simulations . Same as Fig. S4, but for theretained volumetric fraction of water ice (thermochemical criterion
T < T hydr ). Planetesimal radius, (km) F o r m a t i on t i m e a ft e r C A I f o r m a t i on ( M y r) A Reservoir IIplanetesimals not formed
Reservoir Iplanetesimals
Potential percolative core formation, , (vol%) Time after CAI formation, (Myr) P e r c o l a t i v e c o r e f o r m a t i on , ( v o l % ) B Planetesimal radius , formation time (R1) 10 km, 0.3 Myr(R2) 30 km, 0.3 Myr(R3) 300 km, 0.3 Myr (T3) 100 km, 1.3 Myr(T2) 100 km, 0.72 Myr(T1) 100 km, 0.3 Myr
Figure S6. Parameter space and time evolution of core formation via metal percolation in planetesimal simulations . Same asFig. S4, but for the volumetric fraction eligible for percolative core formation (thermochemical criterion T ≥ T perc ). thresholds, f P , crit ( t ) ≡ N markers , crit /N markers , total , (S11)with N markers , crit the number of markers in a planetesimalthat exceed the threshold, and N markers , total the total num-ber of markers in the planetesimal body. Fig. 3A showsthe peak temperatures for various combinations of plan-etesimal formation time and radius (as a volumetric meanper body); we show additional parameter ranges in Figs. S4–S7. In general, larger planetesimal radius and ear-lier formation time (due to more live Al) lead to highertemperatures in the interior. Peak temperatures are al-ways reached at the center of the planetesimal, which ismost insulated from the cold exterior. Typically, tempera-tures decrease from the inside out. However, because thepenetration depth of inward conduction at the planetesi-mal surface is very limited (see below and Castillo-Rogez& Young (2017); Lichtenberg et al. (2016, 2019b)), plan- Lichtenberg et al.
Planetesimal radius, (km) F o r m a t i on t i m e a ft e r C A I f o r m a t i on ( M y r) A Reservoir IIplanetesimals not formed
Reservoir Iplanetesimals
Metal rain-out core formation, , (vol%) Time after CAI formation, (Myr) M e t a l r a i n - ou t c o r e f o r m a t i on , ( v o l % ) B Planetesimal radius , formation time (R1) 13 km, 0.3 Myr(R2) 30 km, 0.3 Myr(R3) 300 km, 0.3 Myr (T3) 100 km, 1.0 Myr(T2) 100 km, 0.72 Myr(T1) 100 km, 0.3 Myr
Figure S7. Parameter space and time evolution of core formation via metal rain-out in planetesimal simulations . Same asFig. S4, but for the volumetric fraction eligible for rain-out core formation (thermochemical criterion ϕ ≥ ϕ rain ). etesimal interiors above ≈
50 km in size are approximatelyisothermal in their deeper interiors after an initial heat-upphase. Once high silicate melt fractions beyond the disag-gregation threshold are reached, the maximum tempera-tures are effectively buffered by convective heat transport.The computational models cover the evolution of singleplanetesimals. To process the time-dependent evolution ofthe planetesimals formed in the disk simulation, we followa procedure (Lichtenberg et al. 2018) to analyze the ther-mochemical evolution of a planetesimal population usinga Monte Carlo approach that parameterizes the birth plan-etesimal population from the streaming instability mecha-nism, as indicated by fluid dynamical models (Johansenet al. 2015; Simon et al. 2017; Abod et al. 2019; Nesvorn´yet al. 2019). For each timestep in the planetesimal forma-tion model (see above) we randomly generate a family of newly formed planetesimals that follow a radius power law ∂N P ∂R P ∝ R − q P , (S12)with the number of bodies in a specific radius bin N P , andpower law index q = 2 . (Johansen et al. 2015; Simonet al. 2017), from which we generate integer radii R min ≤ R P ≤ R max according to R P = || R min (1 − x rand ) − / ( q − || , (S13)with pseudo-random number x rand ∈ [0,1], approximatingthe total newly generated planetesimal mass per timestepand reservoir using ∂M P , res ∂t = (cid:90) m max m min d m P , res ( R P , t ) , (S14) where the planetesimal mass m P and the minimum andmaximum bounds, m min and m max , are related to theplanetesimal radius by a fixed initial mean density of ρ P =2750 kg m − . We choose the minimum planetesimal ra-dius in our parameter space to be R min = 1 km, and setthe upper bound of the birth planetesimal population to R max = 300 km, consistent with a tapered power law to-ward larger planetesimal radii (Johansen et al. 2015; Abodet al. 2019) and planetesimal mean birth sizes on the orderof ≈
100 km (Morbidelli et al. 2009; Delbo et al. 2017).To compute the compositional evolution of the wholeplanetesimal population per timestep and thermochemicalcriterion, we bin the planetesimal mass formed in the diskmodel into distinct radius bins according to the primordialplanetesimal size-frequency distribution (Eq. S12). Wethen calculate the fraction of the planetesimal population per reservoir in a given radius and time bin that fulfills aspecific thermochemical criterion as f pop , crit ( t ) = 1 N P , tot ( t ) (S15) · (cid:90) R max R min N P , bin · f bin , crit · d V bin ( t ) , with N P , tot ( R P , t ) the total number of planetesimals perreservoir, N P , bin the number of planetesimals per bin, and f bin , crit ( R P , t ) · d V bin ( t ) the volumetric fraction, linearly in-terpolated from the single planetesimal simulation grid de-scribed above. The volume change per timestep is nor-malized to the total planetesimal volume per reservoir andtime after 5 Myr in the disk simulation. The volumetric frac-tion increases when more planetesimals form or when a ifurcation of planetary building blocks S9 Time after CAI formation, (Myr) P l ane t e s i m a l popu l a t i on t e m pe r a t u r e , <> ( K ) Water ice melting = 273 KAmphibolite decomposition, = = ( ) = Serpentine breakdown
Reservoir I Reservoir II
Figure S8. Temperature evolution of the simulated planetesimal populations.
Shown are the evolution of the volumetric mean (solidlines) and maximum (dashed lines) temperatures over the whole planetesimal population. Reservoir I is indicated in red, Reservoir IIin blue. Overplotted are the employed thermochemical criteria for the stability of water ice ( T hydr ), hydrous rock ( T decomp ), percolativecore formation ( T perc ), and rain-out core formation ( ϕ rain ) with dark green dashed lines. Additionally we indicate the breakdown formore common serpentine phyllosilicate phases (gray dotted lines, labeled ’Serpentine breakdown’) in the interval of T ∈ [573–673] K(Nakamura 2006; Nakato et al. 2008). planetesimal sub-volume starts to satisfy a specific crite-rion (e.g., hot enough for percolative core formation), anddecreases when a planetesimal sub-volume does not sat-isfy a given criterion anymore (e.g., too hot for water icestability). In Fig. S8 we display the combined temper-ature evolution of the planetesimal populations formed inthe disk model in Reservoir I and Reservoir II together withthe thermochemical criteria. The effects of minimum andmaximum radii and choice of power law have been dis- cussed elsewhere (Lichtenberg et al. 2018). EXTERIOR HEATING EFFECTS ON PLANETESIMALSDUE TO DISK EVOLUTION
Our planetesimal models employ a constant tempera-ture boundary condition, whereas in the disk simulationsthe local gas temperature can vary over time (Fig. S2). Ingeneral, the exterior disk temperature has minor effects onplanetesimal evolution (Golabek et al. 2014). In this sec-tion we seek to determine whether it is negligible in oursimulations. To estimate the effect of the time-dependentdisk temperature on the thermal evolution of planetesi-mals, we developed a thermal diffusion code that solvesthe radial 1-D heat diffusion equation (Hevey & Sanders ∂T∂t = kρc P r ∂∂r (cid:18) r ∂T∂r (cid:19) , (S16)where T is the temperature, t is the time, k is the ther-mal conductivity, ρ is the density, c P is the heat capacityand r is the radial distance from the center. We solvethe heat diffusion equation using the implicit finite differ-ence method (Gerya 2019). For consistency with the plan-etesimal calculations above, we assume that thermal con- ductivity, density and heat capacity are all temperature-independent and use the same values as employed in themore complex planetesimal models. As here we are onlyinterested in the effect of the disk temperature on plan-etesimal evolution, we do not consider radiogenic heating.The initial temperature of the planetesimal interior is thesame as in the previously described models.To simulate the influence of the disk temperature on theplanetesimal interior we prescribe the time-dependent disktemperature at 2 au as surface boundary condition, theorbital location with the highest temperatures where plan-etesimals form in the disk model (Fig. S9A). We model aprofile across the entire planetesimal diameter over time,covering the various evolution stages of the disk. To check Lichtenberg et al.
250 500 750 1000 1250 1500 1750
Maximum temperature, (K) D ep t h i n s i de p l ane t e s i m a l , ( k m ) (km)305070100300 B Varying planetesimal radius = 2 au 250 500 750 1000 1250 1500
Temperature, ( ) (K)
Time (Myr)0.30.60.81.05.0 C Change over time = 100 km, = 2 au0.1 0.2 0.3 0.5 0.7 1 2 3 5
Time in disk model, (Myr) ( K ) Orbital distance, (au)24 715 A Midplane temperature
Figure S9. Effects of external disk temperature on the thermal profile in planetesimals. ( A ) Midplane gas temperature for 2, 4, 7,and 15 au in the disk simulation (lines). ( B ) Resulting peak temperature profiles in planetesimals between 30–300 km in radius at 2 au.The dashed green line shows the temperature criterion for hydrous rock decomposition ( T decomp ). ( C ) Time evolution for a 100 km radiusplanetesimal at 2 au. Different line colors indicate different times during the disk lifetime. numerical convergence of this simplified code, we per-formed calculations with different resolutions ranging from1000 to 8000 grid points and tested timestep lengths rang-ing from 25 days to 1 yr. The results (Fig. S9) show thecode has converged and demonstrate that disk heatingof the planetesimals is limited to the near surface layers, thus only plays a volumetric role for small planetesimals,while the interiors of larger planetesimals remain mostlyunaffected. Maximum temperatures achieved at variousdepths inside differently sized planetesimals are plotted inFig. S9B and the time-evolution of temperature for the R P = 100 km case is displayed in Fig. S9C. These arethe maximum possible effects in the disk model. Internaltemperatures of bodies formed at more distant orbits (Fig.S9A) are even less affected. PLANETARY MIGRATION
To estimate the migration timescales for planetary em-bryos in Reservoir I and II (Fig. 1), we compute the torquefrom the disk on an embryo at each timestep and usethis to evolve the total angular momentum of an embryo of fixed mass with a first-order integrator. For type I mi-gration (see e.g. Baruteau et al. (2014) for a review) thetorque scales as Γ = q h Σ g r Ω , (S17)with the embryo-to-star mass ratio q P = M embryo /M (cid:63) , disk aspect ratio h , and Keplerian frequency Ω p . The ex-act magnitude and direction of the torque can be writtenas a sum of the so-called Lindblad or ‘wave’ torque – gen-erated by waves driven by the planet at resonant locationsin the disk – and the corotation torque – generated by gasthat is on average corotating with the planet. The wavetorque (Paardekooper et al. 2010) can be written as γ Γ wave / Γ = − . − . B + 0 . A, (S18)where γ is the ratio of specific heats, B is the negativepower law exponent of the local disk temperature profileand A is the negative power law exponent of the local diskdensity profile.The calculation of the corotation torque requires deter-mining how both diffusion within the gas and planetary ifurcation of planetary building blocks S11 mass affect the material in this region. At a given vis-cosity, the torque peaks in magnitude as the non-linear‘horseshoe drag’. With too little viscosity to replenish theangular momentum within this region, the horseshoe dragcan become saturated, leaving only the Lindblad torque.For sufficiently high viscosity, the torque is forced into alinear regime which is reduced in magnitude. We there-fore account for local disk conditions when computing theoverall magnitude and direction of the torque. The linearand non-linear versions of the barotropic component of thetorque (Paardekooper et al. 2010) are given by γ Γ lin , baro / Γ = 0 . / − A ) , (S19) γ Γ hs , baro / Γ = 1 . / − A ) (S20)respectively, while the entropy-related components of thelinear corotation torque and horseshoe drag are given by γ Γ lin , ent / Γ = (2 . − . /γ ) (cid:15) e , (S21) γ Γ hs , ent / Γ = 7 . (cid:15) e /γ, (S22)respectively, with (cid:15) e = B − ( γ − A being the negative ofthe local power law exponent of the specific entropy profile.To choose γ , we utilise the fact that the temperature ofthe disk is set only by radiative and viscous heating, andradiative cooling. We expect the radiative processes to bemuch faster than the dynamical timescale, so consider thedisk to be locally isothermal and use γ = 1 . This assump-tion implicitly adds infinite thermal diffusion to the gas,driving the entropy-related part of the corotation torqueinto the linear regime. For transitioning between the unsat-urated and saturated states of the torque, we simply switchvarious components of the corotation/horseshoe drag onand off as necessary – we consider the overall value of thetorque Γ tot to be defined by one of three regimes, largelygoverned by the behaviour of the barotropic component:• ‘Unsaturated’, in which the total torque is thesum of the wave torque, the linear entropy-related torque and the non-linear barotropic horseshoe drag( Γ tot = Γ wave + Γ lin , ent + Γ hs , baro ).• ‘Saturated’, in which the corotation region is as-sumed to be totally depleted of angular momentumsuch that the planet feels no torque from this re-gion and the only remaining component is the wavetorque ( Γ tot = Γ wave ).• ‘Linear’, in which the total torque is the sum of thewave torque and the linear forms of both the entropy-related and barotropic corotation torque ( Γ tot =Γ wave + Γ lin , ent + Γ lin , baro ).We must then choose one of these three regimes at eachtimestep. The maximal, non-linear, unsaturated horse-shoe drag is achieved when the disk viscosity parameter α v satisfies (Baruteau & Lin 2010) . q / h / < α v < . q / h / . (S23)Values of α v lower than the left-hand side of this criterionwill cause the torque to become saturated, while valueshigher than the right-hand side will push the torque intoits linear form. At each timestep we select a regime forthe torque based on which of the inequalities in Eq. S23is satisfied. We compute the total torque based on theregime and the value of Γ , using Eq. S17 and the diskproperties taken directly from the disk model. We linearlyinterpolate between different snapshots of this disk modelto obtain values for disk properties on timescales relevantfor migration. As a test of the validity of this model, wecompared the estimated migration timescales with those in two hydrodynamical simulations of a migrating Earth- mass planet using the two-dimensional FARGO-3D code(Ben´ıtez-Llambay & Masset 2016). In both tests the planetstarted at 1.8 au with α v = 1 × − , a surface densityprofile Σ = (1000g / cm ) r − and h ∝ r / . We createdthese initial conditions in both model setups. In the firsttest we used h = 0 . at 1 au, placing our model in the‘linear’ regime. The difference between the migration ratein our model and the FARGO-3D simulations was around7% over the 2000 yr simulation time. In the second test weused h = 0 . (the ‘unsaturated’ regime) and found thatthe migration rate in the simulation was 32% faster than inour model. The parameters of the latter test would – ac-cording to the model in Paardekooper et al. (2011) – par-tially saturate the corotation torque in our disc. Becausethe corotation torque promotes outward migration in thisdisc, the partial saturation would result in faster migrationin the hydrodynamical models as observed, so any differ-ence between our model and equivalent hydrodynamicalsimulations is due to this approximate treatment of satu-ration. We conclude that a more sophisticated approachwould not alter the conclusions of this paper. In addition to the type I migration timescales we inves-tigate whether a migrating protoplanet opens a gap in thedisk, which may slow down its migration. To determine ifand when this happens, we use the gap-opening criterion
P < (Crida et al. 2006), where P = 3 / hq / + 50 α v h q P . (S24)We consider any embryo which satisfies P < to be gapopening and calculate the threshold masses for each or-bital radius and timestep to estimate whether a protoplanetin Reservoir II may be halted in its migration within the con-text of our disk model. Lichtenberg et al. ( au )
100 10 1 100 10 1 50 10 150 A Reservoir II, giant planet cores
Core starting points0.75 Myr, 10 au1.3 Myr, 17 au2.5 Myr, 80 au2.0 Myr, 90 au O r b i t a l d i s t an c e , ( au ) S n o w l i n e B Reservoir I, rocky protoplanets
Terrestrial planet migration corridors0.4 Myr, 3.0 au, 0.01 0.1 0.6 Myr, 7.0 au, 0.1 0.2 2.0 Myr, 2.0 au, 0.01 0.2
Time in disk model, (Myr) ( au )
10 10 10 0.1 0.2 0.2 C Viscosity ratio sensitivity
Migration viscosity / = 1 / = 0.1 / = 0.01 Gap openingthreshold mass( )
Figure S10. Type I and type II orbital migration in the disk simulation. ( A ) Migration of massive giant planet cores from startinglocations and times (red, blue, yellow, and purple circles) located in Reservoir II. Lines show the migration path for a given mass (1,10, 50, 100, or 150 Earth masses). The colored background indicates the mass threshold for a giant planet core to open a gap in thedisk gas and transition from the type I to the type II migration regime (like the purple 150 M Earth core). ( B ) Type I migration for variousterrestrial planet embryo masses and starting locations and times (red, blue, and yellow circles) in Reservoir I. Lines and shades indicatethe migration paths for a given starting location and protoplanet mass between 0.01–0.1 M Earth (red), 0.1–0.2 M Earth (blue), or 0.01–0.2 M Earth (yellow). Planets that grow more massive during the disk lifetime would migrate inward of 0.5 au. ( C ) Disk viscosity sensitivityof the shown simulation tracks. For each starting location and time in Reservoir I (red circles) or II (blue circles), the disk viscosity thatdrives type I migration ( α mig ) is varied between α v (solid lines) and α t (dotted lines). ifurcation of planetary building blocks S13
Figs. 1 and S10 show migration tracks from this model.Fig. S10A shows the migration tracks of potentially early-formed giant planet cores in Reservoir II. In all but the 150Earth mass ( M Earth ) case, type I migration is too rapid toopen up a gap and transition to the (slower) type II regime.In more comprehensive simulation setups, migration andgrowth show a complex interplay (Bitsch et al. 2019a,b;Johansen et al. 2019), but our models illustrate that early-formed giant planet cores typically migrate rapidly towardthe star. Fig. S10B illustrates that also smaller-mass plan-ets, such as the embryos that accrete to form the innerterrestrial planets, are limited in their maximum mass ac-cretion during the disk phase. If they grow too fast andbecome too massive, they migrate rapidly inward. Laterformation beyond the snow line (to reach the present-dayorbits) would result in overly water-rich planets (Lichten-berg et al. 2019a; Bitsch et al. 2019b), inconsistent withthe present-day volatile inventory of the Solar System ter-restrial planets (Peslier et al. 2018). Fig. S10C displaysthe small sensitivity of the migration efficiency on the em-ployed midplane viscosity. Whether protoplanets in thetype I regime are migrating with α v or α t (see above) haslittle effect on the simulations.Our model of migration is simplified by necessity andthere are effects which cannot be accounted for in sucha model that may affect the dynamical evolution of mi-grating protoplanets (Kley 2019), such as the heatingtorque (Ben´ıtez-Llambay et al. 2015; Bitsch et al. 2019b),variations in angular momentum transport in wind- orturbulence-driven disks (Kimmig et al. 2020), the interplaybetween growth and migration (Johansen et al. 2019), anddynamical interactions between multiple migrating em-bryos (Hands et al. 2014; Hands & Alexander 2018). PLANETARY ACCRETION MODE
To determine the timescales for the onset of planetaryaccretion from an initially gravitationally collapsed popula- tion of planetesimals (Fig. 2), we employ theoretical pre-scriptions for planetesimal-planetesimal collisional growth(Ida & Lin 2004) and pebble accretion (Liu & Ormel 2018;Ormel & Liu 2018).In the scenario of planetary growth via ballistic accretionof planetesimals (Safronov 1972; Morbidelli et al. 2012),the growth rate of the planetary embryos is determinedby the velocity dispersion in the planetesimal population,which is equivalent to the eccentricity distribution. Ec-centricity excitation is damped by gas drag, but less sofor lower-mass planetesimals, which are excited by largerbodies. We use a previously-derived prescription for thegrowth timescale (i.e., the characteristic time for doublingthe initial mass of a planetary object) in this scenario(Kokubo & Ida 2002). In this context, the growth time scales as τ c, acc (cid:39) . × (cid:18) Σ d
10 g cm − (cid:19) − (cid:16) r (cid:17) / × (cid:18) M embryo M Earth (cid:19) / × (cid:18) M (cid:63) M (cid:12) (cid:19) − / × (cid:20) (cid:18) Σ g . × g cm − (cid:19) − / × (cid:16) r (cid:17) / × (cid:18) m P g (cid:19) / (cid:21) yr , (S25)with the solid (dust) surface density Σ d , semi-major axis r ,mass of the accreting planetary embryo M embryo , mass ofthe central star M (cid:63) , mass of the Sun M (cid:12) , gas surface den-sity Σ g , and mass of the accreting low-mass planetesimals m P . We determine the timescale for the onset of accretion from an initial planetesimal population produced from thestreaming instability, which may be the bottleneck for theaccretion of larger Mars-sized protoplanets within the disklifetime (Visser & Ormel 2016; Liu et al. 2019). Follow-ing estimates of the initial mass function of planetesimals(IMF) (Johansen et al. 2015; Abod et al. 2019), we choose M embryo = 3 . × kg and m P = 1 . × kg. Us-ing mean densities of ρ P = 2750 kg m − , this translatesto radii of R embryo = 300 km and R P = 50 km, respec-tively. All other parameters for the growth timescale arecalculated from the time-dependent evolution of the diskand dust coagulation model, as described above. Plan-etesimal growth becomes faster when more of the totalmass is distributed in smaller objects, as they are moreeasily dynamically excited. In the chosen IMF, the ma-jority of the planetesimal mass is in large objects, whichwe approximate by R P = 50 km. However, because ofthe / -exponent in Eq. S25, the simulation results areonly weakly affected by this choice. The timescale for thegrowth from pebble accretion (Ormel & Klahr 2010; Lam- brechts & Johansen 2012) is estimated as (Liu & Ormel2018; Ormel & Liu 2018) τ c , PA = M embryo / (cid:16) (cid:15) PA ˙ M pebble (cid:17) , (S26)with the pebble flux ˙ M pebbles and the pebble accretion ef-ficiency (cid:15) PA = f set (cid:15) set + (1 − f set ) (cid:15) bal , (S27)which approximates the effects of vertical turbulence onthe dust particles in the ballistic and settling regimes ofpebble accretion, depending on the settling fraction of in-dividual pebbles f set (Ormel & Liu 2018). The efficiencies Lichtenberg et al. G r o w t h t i m e sc a l e ( M y r) disk lifetime = 2 au = 15 au Collision-dominated growth A Larger / Smaller
Accreting planetesimals: = 50 km 0.002 (km)
Planetesimal accretion ( ) G r o w t h t i m e sc a l e ( M y r) = 2 au = 15 au Pebble-aided growth B Larger / ( )
Pebble accretion
Figure S11. Growth timescale sensitivity for pebble and planetesimal accretion in Reservoir I and Reservoir II.
The embryo radiiand masses (legends) indicated in A and B are equivalent for constant planet density. ( A ) Growth timescale in the planetesimal (collision)accretion scenario for different embryo radii (dashed, solid, dash-dotted, and dotted lines). Red and blue indicate growth at 2 or 15 au,respectively. The gray horizontal band (labeled ’ ≈ disk lifetime’) indicates the approximate lifetime of the circum-solar disk. The up- anddownward-pointing gray arrows indicate the trends for the growth timescale: increase for increasing embryo radius R embryo , decreasefor smaller planetesimals ( R P ) that accrete onto the embryo. ( B ) Same as in panel A , but for pebble accretion. The growth timescaledecreases for increasing protoplanet mass, but is highly sensitive to secular changes of the pebble flux in the disk simulation. in the settling regime are given by (cid:15) set = (cid:16) (cid:15) − , + (cid:15) − , (cid:17) − / , (S28) (cid:15) set , = A η disk (cid:114) q p τ s ∆ vv K f set , (S29) (cid:15) set , = A q p η disk h eff f , (S30)with the embryo-to-star mass ratio q P , the disk radial pres-sure gradient η disk , dimensionless stopping time τ s , rela-tive velocity between pebble and embryo ∆ v , Keplerianvelocity v K = (cid:113) G ( M (cid:63) + M embryo ) /a, (S31)effective pebble disk aspect ratio h eff , and numerical fitconstants A = 0 . and A = 0 . (Ormel & Liu 2018). The pebble accretion efficiencies in the ballistic regime are given by (cid:15) bal , = R embryo πτ s η disk a (1 − f set ) × (cid:115) q p aR embryo + (cid:18) ∆ vv K (cid:19) , (S32) (cid:15) bal , = 14 √ πη disk τ s h P (cid:0) − f (cid:1) × (cid:32) q p v K ∆ v R embryo a + R a ∆ vv K (cid:33) , (S33)with the pebble disk aspect ratio h P (Liu & Ormel 2018).The direct dependence on the disk aspect ratio in theabove formulation explains why pebble accretion in Fig. 2 ifurcation of planetary building blocks S15 is less efficient in the outer compared to the inner disk(Ormel & Liu 2018), which we discuss below. The pebbleflux ˙ M pebble is derived from the disk and dust coagula-tion model (Dra¸ ˙zkowska & Dullemond 2018). We assumethat the growing embryos have no inclination. The pebbleaccretion efficiency (cid:15) PA decreases with increasing inclina-tion, so the pebble accretion timescale would be longer foran accreting embryo with non-zero inclination. As a com-promise between the midplane and global viscosities ( α t and α v , respectively) in the disk model, we choose a ver-tical diffusion of α z = 1 × − (cf. Ormel & Liu 2018).Fig. S11 shows the sensitivity of the results in Fig.2 to these parameter choices. In the collisional growthmode (Fig. S11A), larger protoplanet mass increases thegrowth timescale, which means slower growth. In con-trast, in the pebble accretion mode, larger protoplanetmass decreases the growth timescale, which means fastergrowth. In the earliest disk stages, however, only for proto-planet masses substantially higher than the mass of Ceres( M Ceres , M Ceres ≈ . × − M Earth ) pebble accretionis faster than collisional accretion. Such high initial plan-etesimal masses are not produced in fluid dynamical sim-ulations of the streaming instability mechanism (Johansenet al. 2015; Simon et al. 2016, 2017; Abod et al. 2019;Li et al. 2019; Rucska & Wadsley 2021), nor supportedby asteroid belt size-frequency distribution reconstructions(Morbidelli et al. 2009; Delbo et al. 2017; Tsirvoulis et al.2018), or the Kuiper-belt impact cratering record (Singeret al. 2019). In Reservoir I at around 1 Myr after CAI for-mation, if sufficiently massive embryos are present, themass doubling timescale for pebble accretion onto Moon-like embryos is less than 1 Myr, which is the approximatetime interval of a high pebble flux in the inner Solar System(cf. Fig. 2).These results align with previous modeling resultsthat suggest a heterogeneous transition of the dominantgrowth mode from collisional to pebble accretion (Jo-hansen et al. 2015; Visser & Ormel 2016; Liu et al. 2019).
In Reservoir II the growth timescales for pebble accretiondisplay a larger deviation with embryo mass relative toReservoir I, however, similarly initial growth must proceedvia planetesimal collisions to reach masses susceptible forpebble accretion to commence. The increase in accretiontimescale for Reservoir I starting at ≈ OPEN SOURCE SOFTWARE PACKAGES
The numerical simulations and calculations in this workmade use of
SCIPY (Virtanen et al. 2020),
NUMPY (Harriset al. 2020),
PANDAS (McKinney 2010), and
ASTROPY (As-tropy Collaboration et al. 2018). Figures were producedusing
MATPLOTLIB (Hunter 2007) and
SEABORN (Waskom& the seaborn development team 2020).
SUPPLEMENTARY TEXT
EVIDENCE FOR COMPOSITIONAL DICHOTOMY &ACCRETION CHRONOLOGY
There are compositional and mass budget trends withheliocentric distance (Hayashi 1981; Gradie & Tedesco1982) for which various explanations have been proposed(Grimm & McSween 1993; Morbidelli et al. 2015, 2016).Mass budget calculations show ≈ M Earth in solids in theinner and ≈ M Earth = 1 . M Jupiter in solids and gasin the outer Solar System planets (Morbidelli et al. 2015;Nimmo et al. 2018). Water abundance increases with or-bital distance (DeMeo et al. 2015; Raymond et al. 2020).Hf-W dating of core segregation events (Kruijer et al. 2014)and Pb-Pb ages of chondrules (Bollard et al. 2017) indi-cate that the planet-forming embryos of Earth and Marsstarted their accretion early during the disk phase at (cid:46) (cid:38)
30 Myr for Earth and ≈ mechanism for the separation of the inner and outer So-lar System reservoirs is required if proto-Jupiter formedfar from the inner Solar System and late during the diskstage. This is consistent with Jupiter’s atmospheric volatilebudget ( ¨Oberg & Wordsworth 2019; Bosman et al. 2019)and the asymmetric distribution of its trojan asteroids (Pi-rani et al. 2019). The proposed filtering mechanism in theJupiter barrier scenario is sensitive to the dynamical con-figuration of the gas giants (Haugbølle et al. 2019) andthe dust destruction efficiency in the outer pressure bump(Dra¸ ˙zkowska et al. 2019). It is also unclear by what mech-anism gas accretion onto giant planets halts and if Jupiterwould stop to grow at sufficient mass during the disk life-time (Lambrechts et al. 2019a) if Jupiter’s core would havebeen formed (cid:46) Lichtenberg et al. formed Jupiter to explain the CC/NC dichotomy (Kruijeret al. 2017).The Hf-W time constraints on embryo growth indicateearly core-segregation events on the parent bodies of ironmeteorites in the inner Solar System reservoir from ≈ ≈ ≈ tiated planetary bodies, and their nitrogen and hydrogenisotope ratios suggest CCs are the principal source ofEarth’s volatile elements (Alexander et al. 2012; Marty2012), in addition to some potential nebular ingassing (Wuet al. 2018; Williams & Mukhopadhyay 2019). However,uncertainties persist in the timing of arrival of the volatilesource bodies and their nature and prior internal evolu-tion (Sarafian et al. 2014; Grewal et al. 2019; Sossi et al.2019; Piani et al. 2020). Deuterium abundances indicatean early contribution from interstellar volatile ices (Cleeveset al. 2014; Piani et al. 2018) or extensive iron oxidation(Sutton et al. 2017). Estimates of the initial water abun-dance in meteoritic parent bodies range from several percent (Marrocchi et al. 2018) to up to 50 per cent water by Table S1. Times of metal-silicate segregation and aque-ous alteration in meteorites.
Times relative to CAI formation(Connelly et al. 2012). CM, CR, and CI values are weightedmeans from several meteorites (Doyle et al. 2015; Jilly-Rehaket al. 2017). Used to generate Figs. 4A and 5A. References:(1) Kruijer et al. (2017), (2) Hunt et al. (2018), (3) Doyle et al.(2015), (4) Fujiya et al. (2013), (5) Fujiya et al. (2012), (6) Jilly-Rehak et al. (2017).Meteorite group Reservoir ∆ t CAI ± σ (Myr) Ref. Metal-silicate segregation
IC NC 0.3 ± ± ± ± ± ± ± ± ± ± ± Aqueous alteration
OC NC 2.4 +1 . − . (3)Tagish Lake CC 4.1 +1 . − . (4)CM CC 4.1 +0 . − . (5)CV CC 4.2 +0 . − . (3)CR CC 4.2 +1 . − . (6)CI CC 4.4 +0 . − . (4)CO CC 5.1 +0 . − . (3) mass (Lodders 2003; van Dishoeck et al. 2014). Traces offluid flow are widespread in the meteorite record, even inNC bodies (Lewis & Jones 2016; Lewis et al. 2018). Thecarbonaceous chondrites exhibit the widest range of water abundances (Alexander et al. 2012, 2018), but some NCchondrites also display substantial abundances with up to6000 ppm H O (McCubbin & Barnes 2019). The rangein measurements indicates that the delivery of water andother highly volatile elements to the inner terrestrial plan-ets was subject to the stochastic and evolving nature ofthe accretion process.
NUCLEOSYNTHETIC ISOTOPE VARIABILITY IN PLANETARYMATERIALS
Distinct nucleosynthetic compositions have been identi-fied for each meteorite parent body and planet, for whichsamples are available for analyses (i.e., Earth and Mars).These nucleosynthetic variations are caused by the het-erogeneous distribution of presolar dust, carrying a char- ifurcation of planetary building blocks
S17 Cr T i Non-carbonaceousreservoir (NC) Carbonaceousreservoir (CC)
COCK CV NWA 2976CM CR NWA 6704CBTafassasset CIAubrites ECsOCsHEDs AngritesMesosiderites PallasitesAcapulcoitesUreilites Earth MoonMars
CAIs
Variations in neutron-rich supernovae-derived isotopes
CC chondritesCC achondritesEarth, Moon, MarsNC chondritesNC achondrites A Mo Z r CAIsCV CRCB COCIOCECEarth deficit
Variations in -process isotopes
CAIsCC chondritesNC chondritesEarth B Figure S12. Nucleosynthetic isotope anomalies in Solar System planetary materials. ( A ) Nucleosynthetic Ti and Cr isotope datain ε notation (parts per 10,000 deviation of the ratio from the respective terrestrial standards). Average values are given for variouschondrite and achondrite groups (symbols). Tafassasset, NWA 2976 and NWA 6704 refer to individual meteorites. Earth, Moon and Marsreflect data from Earth’s mantle, lunar Apollo samples, and Martian meteorites for bulk compositions, respectively. HEDs (Howardite-Eucrite-Diogenites) are meteorites considered to originate from Vesta. CC = carbonaceous chondrites and related groups (blue), NC= non-carbonaceous groups (red), ECs = enstatite chondrites, OCs = ordinary chondrites. Data from (Trinquier et al. 2007; Leya et al.2008, 2009; Larsen et al. 2011; Sanborn et al. 2019; G¨opel et al. 2015). 2- σ uncertainties are indicated with red and blue lines on thesymbols. The dashed turquoise arrow points toward the location of CAIs in Ti-Cr space, off the displayed axes limits. Light gray horizontaland vertical lines center on the axes origin. ( B ) Nucleosynthetic Zr and Mo isotope data of various chondrite classes. Same notation andcolors as in panel A . Earth is most enriched in s -process isotopes compared to all other analyzed materials. Zr data from (Akram et al.2013, 2015), Mo data from (Burkhardt et al. 2011; Render et al. 2018). Both panels adapted from Mezger et al. (2020). acteristic isotopic composition that reflects its stellar for-mation site (Lugaro et al. 2018). This stardust contributed ≈
3% of the total dust in the protoplanetary disk (Hoppeet al. 2017). Stardust was incorporated to different extentsinto each planetary body and led to their distinct isotopecompositions (Ek et al. 2019). Two first order types ofanomalies can be distinguished (Fig. S12): variations in s -process isotopes (Ek et al. 2019; Render et al. 2018;Akram et al. 2015; Burkhardt et al. 2011) and variations in neutron-rich isotopes such as Ca, Ti and Cr (Trin- quier et al. 2007; Leya et al. 2008; Trinquier et al. 2009).The latter indicate nucleosynthesis in a supernova envi-ronment. The variability in neutron-rich isotopes shows adistinct difference between carbonaceous chondrites (CC)and the remaining Solar System material available for lab-oratory studies (i.e. non-carbonaceous chondrites, Marsand Earth, NC).The origin of this isotopic dichotomy has broadly beenattributed to two processes: (a) a change in compositionof the material infalling onto the disk because of a het-erogeneous distribution of stardust in the molecular cloud(Nanne et al. 2019; Jacquet et al. 2019) or (b) thermal pro-cessing and selective destruction of dust during the infallor in the disk itself. The latter process selectively removed specific carrier phases of nucleosynthetic anomalies (Trin-quier et al. 2009; Van Kooten et al. 2016) and assumesa homogeneous and well-mixed molecular cloud. In con-trast to the neutron-rich nucleosynthetic variability (e.g., in Ti or Cr), the s -process variations show a single ar-ray with no clear dichotomy gap (Fig. S12B) (Akram et al.2015; Ek et al. 2019), with the exception of Mo isotopes(Poole et al. 2017; Kruijer et al. 2017). Based on the s -process variations, dust grown in the interstellar medium (ISM) may have been destroyed by thermal processingleading to the s -process variations (Ek et al. 2019). ThisISM-dust contributes ≈
97% of the dust in the protoplan-etary disk (Hoppe et al. 2017) and has an average SolarSystem composition.In our proposed scenario of early Solar System ac-cretion (Fig. 6), a change in the composition of the in-falling material (model (a) ) can explain the initial establish-ment of the nucleosynthetic dichotomy (Nanne et al. 2019;Jacquet et al. 2019). Planetary objects in the ReservoirI region are dominated by material that is mainly fed tothe inner disk during later infall stages and thus influencesthe composition of these planetesimals, while Reservoir IIpreferentially samples older material from before the infallchanged composition. This older material, which includes Lichtenberg et al.
CAIs, is transported to the outer part of the disk duringthe Class I disk phase through viscous spreading and dif-fusional transport (Pignatale et al. 2018; Yang & Ciesla2012; Nanne et al. 2019; Jacquet et al. 2019). However, incontrast to previous suggestions (Kruijer et al. 2017; De-sch et al. 2018; Scott et al. 2018), our proposed scenariodoes not require the presence of proto-Jupiter to block theaddition of outer Solar System materials to Reservoir I inthe inner Solar System. Instead, both the chemical andisotopic dichotomy are an inherent consequence of thetwo subsequent episodes of planetesimal formation andthe pebble flux suppression due to the dust pile-up alongthe drifting snow line. Similarly, thermal processing of infalling dust (model (b) )is consistent with our proposed model and could have gen-erated or contributed to the observed nucleosynthetic het-erogeneities. However, selective destruction of dust in thedisk itself is less likely, because in our model most of thethermal processing of inner Solar System material occursinside of the planetesimals. This processing should notaffect the nucleosynthetic anomalies because these aremainly carried by refractory elements (Toth et al. 2020).REFERENCES
Abod, C. P., Simon, J. B., Li, R., et al. 2019,
ApJ , 883, 192,doi:
Akram, W., Sch¨onb¨achler, M., Bisterzo, S., & Gallino, R. 2015,
GeoCoA , 165, 484, doi:
Akram, W., Sch¨onb¨achler, M., Sprung, P., & Vogel, N. 2013,
ApJ , 777, 169, doi:
Alexander, C. M. O., Bowden, R., Fogel, M. L., et al. 2012,Science, 337, 721, doi:
Alexander, C. M. O., McKeegan, K. D., & Altwegg, K. 2018,
SSRv , 214, 36, doi:
Alibert, Y., Venturini, J., Helled, R., et al. 2018, Nat. Astron., 2,2397, doi:
Andrews, S. M., Wilner, D. J., Zhu, Z., et al. 2016,
ApJL , 820,L40, doi:
Andrews, S. M., Huang, J., P´erez, L. M., et al. 2018,
ApJL , 869,L41, doi:
Ansdell, M., Williams, J. P., van der Marel, N., et al. 2016,
ApJ ,828, 46, doi:
Astropy Collaboration, Price-Whelan, A. M., Sip˝ocz, B. M., et al.2018, AJ , 156, 123, doi: Aumatell, G., & Wurm, G. 2014,
MNRAS , 437, 690,doi:
Bagdassarov, N., Solferino, G., Golabek, G. J., & Schmidt,M. W. 2009, Earth Planet. Sci. Lett., 288, 84,doi:
Bai, X.-N. 2017,
ApJ , 845, 75, doi:
Bai, X.-N., Ye, J., Goodman, J., & Yuan, F. 2016,
ApJ , 818, 152,doi:
Balbus, S. A., & Hawley, J. F. 1991,
ApJ , 376, 214,doi:
Baruteau, C., & Lin, D. N. C. 2010,
ApJ , 709, 759,doi:
Baruteau, C., Crida, A., Paardekooper, S. J., et al. 2014, inProtostars and Planets VI, ed. H. Beuther, R. S. Klessen,C. P. Dullemond, & T. Henning, 667,doi:
Benedikt, M. R., Scherf, M., Lammer, H., et al. 2020,
Icarus ,347, 113772, doi:
Ben´ıtez-Llambay, P., Masset, F., Koenigsberger, G., & Szul´agyi,J. 2015,
Nature , 520, 63, doi:
Ben´ıtez-Llambay, P., & Masset, F. S. 2016,
ApJS , 223, 11,doi:
Birnstiel, T., Fang, M., & Johansen, A. 2016,
SSRv , 205, 41,doi:
Birnstiel, T., Klahr, H., & Ercolano, B. 2012,
A&A , 539, A148,doi:
Bitsch, B., Izidoro, A., Johansen, A., et al. 2019a,
A&A , 623,A88, doi:
Bitsch, B., Raymond, S. N., & Izidoro, A. 2019b,
A&A , 624,A109, doi:
Bollard, J., Connelly, J. N., Whitehouse, M. J., et al. 2017, Sci.Adv., 3, e1700407, doi:
Bosman, A. D., Cridland, A. J., & Miguel, Y. 2019,
A&A , 632,L11, doi:
Bouvier, A., & Wadhwa, M. 2010, Nat. Geosci., 3, 637,doi:
Brasser, R., & Mojzsis, S. J. 2020, Nat. Astron., 4, 492,doi:
Braukm¨uller, N., Wombacher, F., Funk, C., & M¨unker, C. 2019,Nat. Geosci., 12, 564, doi:
Burkhardt, C., Kleine, T., Oberli, F., et al. 2011, Earth Planet.Sci. Lett., 312, 390, doi:
Carrera, D., Gorti, U., Johansen, A., & Davies, M. B. 2017,
ApJ ,839, 16, doi:
Castillo-Rogez, J., Johnson, T. V., Lee, M. H., et al. 2009,
Icarus , 204, 658, doi: ifurcation of planetary building blocks
S19
Castillo-Rogez, J., & Young, E. D. 2017, in Planetesimals: EarlyDifferentiation and Consequences for Planets, ed. L. T.Elkins-Tanton & B. P. Weiss (Cambridge: CambridgeUniversity Press), 92–114Charnoz, S., Pignatale, F. C., Hyodo, R., et al. 2019,
A&A , 627,A50, doi:
Ciesla, F. J., & Cuzzi, J. N. 2006,
Icarus , 181, 178,doi:
Cleeves, L. I., Bergin, E. A., Alexand er, C. M. O. D., et al. 2014,Science, 345, 1590, doi:
Connelly, J. N., Bizzarro, M., Krot, A. N., et al. 2012, Science,338, 651, doi:
Costa, A., Caricchi, L., & Bagdassarov, N. 2009, Geochem.Geophys. Geosys., 10, Q03010,doi:
Crida, A., Morbidelli, A., & Masset, F. 2006,
Icarus , 181, 587,doi:
Cuzzi, J. N., & Zahnle, K. J. 2004,
ApJ , 614, 490,doi:
Dauphas, N., & Pourmand, A. 2011,
Nature , 473, 489,doi:
Davidson, J., Busemann, H., Nittler, L. R., et al. 2014,
GeoCoA ,139, 248, doi:
Delbo, M., Walsh, K., Bolin, B., Avdellidou, C., & Morbidelli, A.2017, Science, 357, 1026, doi:
DeMeo, F. E., Alexander, C. M. O., Walsh, K. J., Chapman,C. R., & Binzel, R. P. 2015, in Asteroids IV, University ofArizona Press, Tucson, ed. P. Michel, F. E. DeMeo, & W. F.Bottke, 13–41,doi:
Desch, S. J., Kalyaan, A., & O’D. Alexander, C. M. 2018,
ApJS ,238, 11, doi:
Doyle, P. M., Jogo, K., Nagashima, K., et al. 2015, Nat.Commun., 6, 7444, doi:
Dra¸ ˙zkowska, J., & Alibert, Y. 2017,
A&A , 608, A92,doi:
Dra¸ ˙zkowska, J., & Dullemond, C. P. 2018,
A&A , 614, A62,doi:
Dra¸ ˙zkowska, J., Li, S., Birnstiel, T., Stammler, S. M., & Li, H.2019,
ApJ , 885, 91, doi:
Dra¸ ˙zkowska, J., Alibert, Y., & Moore, B. 2016,
A&A , 594, A105,doi:
Dullemond, C. P., Apai, D., & Walch, S. 2006a,
ApJL , 640, L67,doi:
Dullemond, C. P., Natta, A., & Testi, L. 2006b,
ApJL , 645, L69,doi:
Ek, M., Hunt, A. C., Lugaro, M., & Sch¨onb¨achler, M. 2019, Nat.Astron., 4, 273, doi:
Elkins-Tanton, L. T. 2012, Annu. Rev. Earth Planet. Sci., 40,113, doi:
Ercolano, B., Jennings, J., Rosotti, G., & Birnstiel, T. 2017,
MNRAS , 472, 4117, doi:
Ercolano, B., & Pascucci, I. 2017, R. Soc. Open Sci., 4, 170114,doi:
Ferri`ere, K. M. 2001, Rev. Mod. Phys., 73, 1031,doi:
Flaherty, K., Hughes, A. M., Simon, J. B., et al. 2020,
ApJ , 895,109, doi:
Flaherty, K. M., Hughes, A. M., Rosenfeld, K. A., et al. 2015,
ApJ , 813, 99, doi:
Flaherty, K. M., Hughes, A. M., Rose, S. C., et al. 2017,
ApJ ,843, 150, doi:
Flock, M., Fromang, S., Gonz´alez, M., & Commerc¸on, B. 2013,
A&A , 560, A43, doi:
Fu, R. R., & Elkins-Tanton, L. T. 2014, Earth Planet. Sci. Lett.,390, 128, doi:
Fujiya, W., Sugiura, N., Hotta, H., Ichimura, K., & Sano, Y. 2012,Nat. Commun., 3, 627, doi:
Fujiya, W., Sugiura, N., Sano, Y., & Hiyagon, H. 2013, EarthPlanet. Sci. Lett., 362, 130, doi:
Gammie, C. F. 1996,
ApJ , 457, 355, doi:
Gerya, T. V. 2019, Introduction to numerical geodynamicmodelling (Cambridge University Press)Gerya, T. V., & Yuen, D. A. 2003, Phys. Earth Planet. Inter., 140,293, doi: —. 2007, Phys. Earth Planet. Inter., 163, 83,doi:
Golabek, G. J., Bourdon, B., & Gerya, T. V. 2014,
M&PS , 49,1083, doi:
G¨opel, C., Birck, J.-L., Galy, A., Barrat, J.-A., & Zanda, B. 2015,
GeoCoA , 156, 1, doi:
Gradie, J., & Tedesco, E. 1982, Science, 216, 1405,doi:
Greaves, J. S., & Rice, W. K. M. 2010,
MNRAS , 407, 1981,doi:
Grewal, D. S., Dasgupta, R., Sun, C., Tsuno, K., & Costin, G.2019, Sci. Adv., 5, eaau3669, doi:
Grimm, R. E., & McSween, H. Y. 1993, Science, 259, 653Grimm, R. E., & McSween, Jr., H. Y. 1989,
Icarus , 82, 244,doi:
Hands, T. O., & Alexander, R. D. 2018,
MNRAS , 474, 3998,doi:
Hands, T. O., Alexander, R. D., & Dehnen, W. 2014,
MNRAS ,445, 749, doi:
Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020,Nature, 585, 357–362, doi:
Harsono, D., Bjerkeli, P., van der Wiel, M. H. D., et al. 2018, Nat.Astron., 2, 646, doi:
Haugbølle, T., Weber, P., Wielandt, D. P., et al. 2019, AJ , 158,55, doi: Lichtenberg et al.
Hayashi, C. 1981, Progress of Theoretical Physics Supplement,70, 35, doi:
Hevey, P. J., & Sanders, I. S. 2006, Meteorit. Planet. Sci., 41,95, doi:
Hin, R. C., Coath, C. D., Carter, P. J., et al. 2017,
Nature , 549,511, doi:
Hoppe, P., Leitner, J., & Kodol´anyi, J. 2017, Nat. Astron., 1, 617,doi:
Hueso, R., & Guillot, T. 2005,
A&A , 442, 703,doi:
Hunt, A. C., Cook, D. L., Lichtenberg, T., et al. 2018, EarthPlanet. Sci. Lett., 482, 490, doi:
Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90Ida, S., & Guillot, T. 2016,
A&A , 596, L3,doi:
Ida, S., & Lin, D. N. C. 2004,
ApJ , 616, 567,doi:
Ida, S., Yamamura, T., & Okuzumi, S. 2019,
A&A , 624, A28,doi:
Jacquet, E., Pignatale, F. C., Chaussidon, M., & Charnoz, S.2019,
ApJ , 884, 32, doi:
Jilly-Rehak, C. E., Huss, G. R., & Nagashima, K. 2017,
GeoCoA , 201, 224, doi:
Johansen, A., Blum, J., Tanaka, H., et al. 2014, in Protostarsand Planets VI, ed. H. Beuther, R. S. Klessen, C. P.Dullemond, & T. Henning, 547,doi:
Johansen, A., Ida, S., & Brasser, R. 2019,
A&A , 622, A202,doi:
Johansen, A., Mac Low, M.-M., Lacerda, P., & Bizzarro, M.2015, Sci. Adv., 1, e1500109, doi:
Johansen, A., Oishi, J. S., Mac Low, M.-M., et al. 2007,
Nature ,448, 1022, doi:
Kimmig, C. N., Dullemond, C. P., & Kley, W. 2020,
A&A , 633,A4, doi:
Kita, N. T., Yin, Q.-Z., MacPherson, G. J., et al. 2013, Meteorit.Planet. Sci., 48, 1383, doi:
Kley, W. 2019, in Saas-Fee Advanced Course 45: FromProtoplanetary Disks to Planet Formation (Springer), 151Kokubo, E., & Ida, S. 2002,
ApJ , 581, 666, doi:
Krot, A. N., Nagashima, K., Alexander, C. M. O., et al. 2015, inAsteroids IV, ed. P. Michel, F. E. DeMeo, & W. F. Bottke,635–660, doi:
Kruijer, T. S., Burkhardt, C., Budde, G., & Kleine, T. 2017, Proc.Natl. Acad. Sci. U.S.A., 114, 6712,doi:
Kruijer, T. S., Touboul, M., Fischer-G¨odde, M., et al. 2014,Science, 344, 1150, doi:
Lambrechts, M., & Johansen, A. 2012,
A&A , 544, A32,doi:
Lambrechts, M., Lega, E., Nelson, R. P., Crida, A., & Morbidelli,A. 2019a,
A&A , 630, A82,doi:
Lambrechts, M., Morbidelli, A., Jacobson, S. A., et al. 2019b,
A&A , 627, A83, doi:
Larsen, K. K., Trinquier, A., Paton, C., et al. 2011,
ApJL , 735,L37, doi:
Lesur, G., Kunz, M. W., & Fromang, S. 2014,
A&A , 566, A56,doi:
Levison, H. F., Kretke, K. A., Walsh, K. J., & Bottke, W. F. 2015,Proc. Natl. Acad. Sci. U.S.A., 112, 14180,doi:
Lewis, J. A., & Jones, R. H. 2016, Meteorit. Planet. Sci., 51,1886, doi:
Lewis, J. A., Jones, R. H., & Garcea, S. C. 2018,
GeoCoA , 240,293, doi:
Leya, I., Sch¨onb¨achler, M., Kr¨ahenb¨uhl, U., & Halliday, A. N.2009,
ApJ , 702, 1118, doi:
Leya, I., Sch¨onb¨achler, M., Wiechert, U., Kr¨ahenb¨uhl, U., &Halliday, A. N. 2008, Earth Planet. Sci. Lett., 266, 233,doi:
Li, R., Youdin, A. N., & Simon, J. B. 2019,
ApJ , 885, 69,doi:
Lichtenberg, T., Golabek, G. J., Burn, R., et al. 2019a, Nat.Astron., 3, 307, doi:
Lichtenberg, T., Golabek, G. J., Dullemond, C. P., et al. 2018,
Icarus , 302, 27, doi:
Lichtenberg, T., Golabek, G. J., Gerya, T. V., & Meyer, M. R.2016,
Icarus , 274, 350, doi:
Lichtenberg, T., Keller, T., Katz, R. F., Golabek, G. J., & Gerya,T. V. 2019b, Earth Planet. Sci. Lett., 507, 154,doi:
Liu, B., & Ormel, C. W. 2018,
A&A , 615, A138,doi:
Liu, B., Ormel, C. W., & Johansen, A. 2019,
A&A , 624, A114,doi:
Lodders, K. 2003,
ApJ , 591, 1220, doi:
Lugaro, M., Ott, U., & Kereszturi, ´A. 2018, Prog. Part. Nucl.Phys., 102, 1, doi:
Manara, C. F., Morbidelli, A., & Guillot, T. 2018,
A&A , 618, L3,doi:
Marchi, S., Walker, R. J., & Canup, R. M. 2020, Sci. Adv., 6,eaay2338, doi:
Marrocchi, Y., Bekaert, D. V., & Piani, L. 2018, Earth Planet. Sci.Lett., 482, 23, doi:
Marty, B. 2012, Earth Planet. Sci. Lett., 313, 56,doi:
Maurel, C., Bryson, J. F. J., Lyons, R. J., et al. 2020, Sci. Adv., 6,eaba1303, doi: ifurcation of planetary building blocks
S21
McCubbin, F. M., & Barnes, J. J. 2019, Earth Planet. Sci. Lett.,526, 115771, doi:
McKinney, W. 2010, in Proceedings of the 9th Python in ScienceConference, ed. S. van der Walt & J. Millman, 51–56Mezger, K., Debaille, V., & Kleine, T. 2013,
SSRv , 174, 27,doi:
Mezger, K., Sch¨onb¨achler, M., & Bouvier, A. 2020, Space Sci.Rev., 216, 1Monteux, J., Golabek, G. J., Rubie, D. C., Tobie, G., & Young,E. D. 2018, Space Sci. Rev., 214, 39,doi:
Morbidelli, A., Bottke, W. F., Nesvorn´y, D., & Levison, H. F. 2009,
Icarus , 204, 558, doi:
Morbidelli, A., Lambrechts, M., Jacobson, S., & Bitsch, B. 2015,
Icarus , 258, 418, doi:
Morbidelli, A., Lunine, J. I., O’Brien, D. P., Raymond, S. N., &Walsh, K. J. 2012, Annu. Rev. Earth Planet. Sci., 40, 251,doi:
Morbidelli, A., Bitsch, B., Crida, A., et al. 2016,
Icarus , 267, 368,doi:
Najita, J. R., & Kenyon, S. J. 2014,
MNRAS , 445, 3315,doi:
Nakamura, T. 2006, Earth Planet. Sci. Lett., 242, 26,doi:
Nakato, A., Nakamura, T., Kitajima, F., & Noguchi, T. 2008,Earth Planets Space, 60, 855Nanne, J. A. M., Nimmo, F., Cuzzi, J. N., & Kleine, T. 2019, EarthPlanet. Sci. Lett., 511, 44, doi:
Nesvorn´y, D., Li, R., Youdin, A. N., Simon, J. B., & Grundy,W. M. 2019, Nat. Astron., 3, 808,doi:
Nimmo, F., & Kleine, T. 2007,
Icarus , 191, 497,doi:
Nimmo, F., Kretke, K., Ida, S., Matsumura, S., & Kleine, T. 2018,
SSRv , 214, 101, doi:
Norris, C. A., & Wood, B. J. 2017,
Nature , 549, 507,doi: ¨Oberg, K. I., & Wordsworth, R. 2019, AJ , 158, 194,doi: Ormel, C. W., & Klahr, H. H. 2010,
A&A , 520, A43,doi:
Ormel, C. W., & Liu, B. 2018,
A&A , 615, A178,doi:
Paardekooper, S. J., Baruteau, C., Crida, A., & Kley, W. 2010,
MNRAS , 401, 1950, doi:
Paardekooper, S. J., Baruteau, C., & Kley, W. 2011,
MNRAS ,410, 293, doi:
Pape, J., Mezger, K., Bouvier, A. S., & Baumgartner, L. P. 2019,
GeoCoA , 244, 416, doi:
Pascucci, I., & Tachibana, S. 2010, in Protoplanetary Dust:Astrophysical and Cosmochemical Perspectives, ed. D. A.Apai & D. S. Lauretta, 263–298Pascucci, I., Testi, L., Herczeg, G. J., et al. 2016,
ApJ , 831, 125,doi:
Patoˇcka, V., Calzavarini, E., & Tosi, N. 2020, Phys. Rev. Fluids,5, 114304, doi:
Peslier, A., Sch¨onb¨achler, M., Busemann, H., & Karato, S.-I.2018, Space Sci. Rev., 212, 743Piani, L., Marrocchi, Y., Rigaudier, T., et al. 2020, Science, 369,1110, doi:
Piani, L., Yurimoto, H., & Remusat, L. 2018, Nat. Astron., 2,317, doi:
Pignatale, F. C., Charnoz, S., Chaussidon, M., & Jacquet, E.2018,
ApJL , 867, L23, doi:
Pirani, S., Johansen, A., Bitsch, B., Mustill, A. J., & Turrini, D.2019,
A&A , 623, A169, doi:
Poole, G. M., Rehk¨amper, M., Coles, B. J., Goldberg, T., &Smith, C. L. 2017, Earth Planet. Sci. Lett., 473, 215,doi:
Raymond, S. N., Izidoro, A., & Morbidelli, A. 2020, in PlanetaryAstrobiology, ed. V. S. Meadows, G. N. Arnay, B. E. Schmidt,& D. J. Des Marais (University of Arizona Press), 287Render, J., Brennecka, G. A., Wang, S.-J., Wasylenki, L. E., &Kleine, T. 2018,
ApJ , 862, 26,doi:
Rucska, J. J., & Wadsley, J. W. 2021,
MNRAS , 500, 520,doi:
Rudge, J. F., Kleine, T., & Bourdon, B. 2010, Nat. Geosci., 3,439, doi:
Safronov, V. S. 1972, Evolution of the protoplanetary cloud andformation of the earth and planets.Sanborn, M. E., Wimpenny, J., Williams, C. D., et al. 2019,
GeoCoA , 245, 577, doi:
Sarafian, A. R., Nielsen, S. G., Marschall, H. R., McCubbin,F. M., & Monteleone, B. D. 2014, Science, 346, 623,doi:
Sarafian, A. R., Hauri, E. H., McCubbin, F. M., et al. 2017,Philos. Trans. R. Soc. A, 375, 20160209,doi:
Sato, T., Okuzumi, S., & Ida, S. 2016,
A&A , 589, A15,doi:
Schiller, M., Bizzarro, M., & Fernandes, V. A. 2018,
Nature ,555, 507, doi:
Schobert, B. N., Peeters, A. G., & Rath, F. 2019,
ApJ , 881, 56,doi:
Sch¨onb¨achler, M., Carlson, R. W., Horan, M. F., Mock, T. D., &Hauri, E. H. 2010, Science, 328, 884,doi: Lichtenberg et al.
Schoonenberg, D., & Ormel, C. W. 2017,
A&A , 602, A21,doi:
Scott, E. R. D., Krot, A. N., & Sanders, I. S. 2018,
ApJ , 854,164, doi:
Segura-Cox, D. M., Schmiedeke, A., Pineda, J. E., et al. 2020,
Nature , 586, 228, doi:
Shakura, N. I., & Sunyaev, R. A. 1973,
A&A , 24, 337Simon, J. B., Armitage, P. J., Li, R., & Youdin, A. N. 2016,
ApJ ,822, 55, doi:
Simon, J. B., Armitage, P. J., Youdin, A. N., & Li, R. 2017,
ApJL ,847, L12, doi:
Singer, K. N., McKinnon, W. B., Gladman, B., et al. 2019,Science, 363, 955, doi:
Sossi, P. A., Klemme, S., O’Neill, H. S. C., Berndt, J., & Moynier,F. 2019,
GeoCoA , 260, 204, doi:
Stevenson, D. J., & Lunine, J. I. 1988,
Icarus , 75, 146,doi:
Sutton, S., Alexander, C. M. O., Bryant, A., et al. 2017,
GeoCoA , 211, 115, doi:
Teague, R., Bae, J., & Bergin, E. A. 2019,
Nature , 574, 378,doi:
Teague, R., Guilloteau, S., Semenov, D., et al. 2016,
A&A , 592,A49, doi:
Toth, E. R., Fehr, M. A., Friebel, M., & Sch¨onb¨achler, M. 2020,
GeoCoA , 274, 286, doi:
Trinquier, A., Birck, J.-L., & All`egre, C. J. 2007,
ApJ , 655, 1179,doi:
Trinquier, A., Elliott, T., Ulfbeck, D., et al. 2009, Science, 324,374, doi:
Tsirvoulis, G., Morbidelli, A., Delbo, M., & Tsiganis, K. 2018,
Icarus , 304, 14, doi:
Turner, N. J., Fromang, S., Gammie, C., et al. 2014, inProtostars and Planets VI, ed. H. Beuther, R. S. Klessen,C. P. Dullemond, & T. Henning, 411,doi:
Tychoniec, Ł., Manara, C. F., Rosotti, G. P., et al. 2020,
A&A ,640, A19, doi: van Dishoeck, E. F., Bergin, E. A., Lis, D. C., & Lunine, J. I.2014, in Protostars and Planets VI, ed. H. Beuther, R. S.Klessen, C. P. Dullemond, & T. Henning, 835,doi:
Van Kooten, E. M. M. E., Wielandt, D., Schiller, M., et al. 2016,Proc. Natl. Acad. Sci. U.S.A., 113, 2011,doi:
Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat.Methods, 17, 261, doi:
Visser, R. G., & Ormel, C. W. 2016,
A&A , 586, A66,doi:
Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T.2009,
ApJ , 702, 1490, doi: —. 2011,
ApJ , 737, 36, doi:
Wang, H., Weiss, B. P., Bai, X.-N., et al. 2017, Science, 355,623, doi:
Warren, P. H. 2011, Earth Planet. Sci. Lett., 311, 93,doi:
Waskom, M., & the seaborn development team. 2020,mwaskom/seaborn, latest, Zenodo,doi:
Williams, C. D., & Mukhopadhyay, S. 2019,
Nature , 565, 78,doi:
Wu, J., Desch, S. J., Schaefer, L., et al. 2018, J. Geophys. Res.Planet., 123, 2691, doi:
Yang, L., & Ciesla, F. J. 2012, Meteorit. Planet. Sci., 47, 99,doi:
Youdin, A. N., & Goodman, J. 2005,
ApJ , 620, 459,doi: