Natural separation of two primordial planetary reservoirs in an expanding solar protoplanetary disk.

Meteorites display an isotopic composition dichotomy between noncarbonaceous (NC) and carbonaceous (CC) groups, indicating that planetesimal formation in the solar protoplanetary disk occurred in two distinct reservoirs. The prevailing view is that a rapidly formed Jupiter acted as a barrier between these reservoirs. We show a fundamental inconsistency in this model: If Jupiter is an efficient blocker of drifting pebbles, then the interior NC reservoir is depleted by radial drift within a few hundred thousand years. If Jupiter lets material pass it, then the two reservoirs will be mixed. Instead, we demonstrate that the arrival of the CC pebbles in the inner disk is delayed for several million years by the viscous expansion of the protoplanetary disk. Our results support the hypothesis that Jupiter formed in the outer disk (>10 astronomical units) and allowed a considerable amount of CC material to pass it and become accreted by the terrestrial planets.

INTRODUCTION

The chronological ages of meteorites are inferred either by dating the formation of their main components (e.g., chondrules) or by dating specific physicochemical processes (e.g., core-mantle segregation or aqueous alteration) to derive the parent body accretion time through thermal modeling. Such studies indicate that the accretion of meteorite parent bodies in the solar system started already within a few hundred thousand years from the formation of the calcium-aluminium-rich inclusions (CAIs) (generally taken to be the starting time of solar system formation) () and lasted for several million years (–). The difference in composition between the noncarbonaceous (NC) group, originating in the inner solar system, and the carbonaceous (CC) group, originating in the outer solar system, needs to be present in the very early disk phase, within a few 105 years to explain the formation of primitive achondrites and iron meteorites of NC composition. These two reservoirs must remain radially separated at least for 2 to 3 Ma for the primitive NC and CC parent bodies to form around that time with distinct isotopic compositions ().

Disk thermal processing has been proposed to explain the origin of these two disk reservoirs (, ). In this scenario, solar system dust is made up of refractory presolar grains and more volatile mantles grown in the interstellar medium (ISM) (). The presolar grains are s-process dominated and likely originate from asymptotic giant branch stars, whereas the ISM grains have a near-solar isotopic composition (CI-like) condensed from the gas that is already well mixed from different isotopic sources (e.g., supernovae). The NC composition is then a result of the selective thermal destruction of either ISM-grown mantles or fragile presolar carriers (, ). During the early infall phase, the episodic accretion associated with an FU Orionis–type outburst temporally raises the disk temperature and moves the water ice line to the 20 – to 30–astronomical unit (AU) disk region (). This would have naturally resulted in the separation of the NC component from the bulk CI component out to those distances, provided that the supernova-derived material resides chiefly in the fragile mantle of ice and organics. Alternatively, the initial two distinct disk reservoirs could also be inherited from isotopically heterogeneous accretion of molecular cloud gas onto the disk (, ).

Previous studies have suggested that an early and rapid growth of Jupiter is necessary to explain the subsequent 2- to 3-Ma separation of these two distinct reservoirs (, , ). Here, we find a fundamental inconsistency with the Jupiter barrier scenario, which requires an alternative model. We demonstrate that, after the formation of the protoplanetary disk, the combination of viscous spreading and radial drift of dust particles can consistently explain the nucleosynthetic isotopic compositions in both the asteroid and planetary record—without invoking Jupiter to block the inward flow of solid material.

RESULTS

We provide a simple estimation of the solid particles’ radial drift time scale. Their radial drift velocity is given by (see Materials and Methods)where vK ≡ ΩKr is the Keplerian velocity, r is the orbital distance, ΩK is the angular velocity, is the headwind prefactor (a measure of how large the gas rotation velocity deviates from the Keplerian velocity), P is the gas pressure, Hg is the disk gas scale height, and ug is the gas velocity at the disk midplane, which is proportional to the midplane turbulent viscous parameter αt. The Stokes number of the particle is defined as τs ≡ ΩKtstop, where tstop is the stopping time (a measure of how fast the particle adjusts its momentum to the surrounding gas).

The particle’s Stokes number in the fragmentation-limit growth regime iswhere vF is the fragmentation threshold velocity. In the inner disk regions with a moderate midplane turbulence (e.g., at 5 AU), αt ≲ τs,F ≪ 1, and the drift time scale can approximately expressed as

The above time scale strongly depends on the fragmentation threshold vF. For a high fragmentation threshold (vF ≫ 1 m/s), the larger particle sizes lead to more rapid orbital decay of drifting particles, which results in a faster mixing between NC and CC reservoirs. On the other hand, for disks with a high turbulent level αt ≫ τs,F, solid particles are well coupled to the gas, and their radial drift is determined by the gas accretion flow.

We construct a one-dimensional (1D) model for the gas and solid components in a protoplanetary disk. The disk is evolved using the conventional α-viscosity approach. Most of the gas in the expanding disk is accreted onto the central star, while the rest spreads outward to conserve the angular momentum (). The solid particles are allowed to grow to their fragmentation-limited size (), yielding particle sizes that are smaller further from the star and later in time. The effects of the gas spreading and dust radial drift are included in the calculations.

We show in Fig. 1A that the particles initially at r ≫ rexp move outward together with the expanding gas, where rexp is the separation radius between the gas inflow and outflow. Their motion then becomes inward when the particles’ radial drift overwhelms the outward spreading of the gas at later times. We find that dust particles starting at 1 to 10 AU drift into the inner disk within 1 Ma, whereas those initially located in the 20- to 30-AU region sweep the asteroid belt 2 to 4 Ma later, ending up with a millimeter characteristic size. This is also the case when we consider a disk with less vigorous turbulence in the midplane than the upper surface layer. In this so-called dead zone model, the local midplane turbulence further depends on the column density of the gas in the protoplanetary disk (see Materials and Methods). We additionally account for the early disk building phase in Fig. 1C, where the angular momentum is dominated by the infall process with gas only being accreted inward. After a typical infall time of 0.3 Ma, the disk viscosity dominates angular momentum transport and leads the outward spreading of the gas. Our result implies that the inner disk can preserve its NC nature until the CC dust particles drift through, as long as the initial boundary of NC/CC reservoirs (rCC) is larger than the original viscous expansion radius of gas (rexp).

Fig. 1.

Radial drift of dust particles in a viscously spreading disk.

Three panels give (A) uniform turbulent model, (B) dead zone model, and (C) dead zone + infall model. The black curves trace the radial drift of dust particles whose size is limited in the fragmentation regime, while the gray dashed lines indicate the gas trajectories. The background color refers to the size of the particle at the corresponding time and disk location, with the magenta dash-dot lines marking pebbles of 0.1, 1, and 10 mm in size. The disk parameters are vF = 1 m/s, rd = 10 AU, Cd = 0.14 M⊙, Tirr0 = 130 K, and α = 5 × 10−3 in the uniform turbulent model and αA = 10−2 and αD = 5 × 10−4 in the dead zone model (fiducial run). In the dead zone + infall model, the early infall period where gas only accretes inward lasts for 0.3 Ma, and after that, the disk gas transports angular momentum through viscous spreading. The outward flow of gas in the expanding disk brings along with it the solids from the outer regions of the disk. As the solids detach from the gas flow, they drift into the inner disk regions with the arrival time determined by their initial distance from the star.

In the prevailing view, the early and rapid formation of Jupiter has been invoked to explain the NC-CC separation (, , ). When the core of Jupiter grows to approximately 10 to 20 M⊕ (), it starts carving a gap in the gas disk, preventing the outer CC solids from drifting inward and mixing with the inner NC reservoir. This model requires proto-Jupiter to form early within 0.5 to 1 Ma to explain the formation of differentiated planetesimals in both reservoirs. We identify an inconsistency in this model in Fig. 2A: If Jupiter appears early in the 5-AU region and totally blocks the flow of solids, then the solids of NC composition interior of Jupiter would drift into the Sun within a few hundred thousand years (Fig. 2A). This would then preclude the formation of the primitive NC parent bodies after 2 to 3 Ma.

Fig. 2.

Radial drift of pebbles and the pebble flux entering the inner terrestrial planet formation region.

The trajectories of the NC and CC dust particles are shown in magenta and blue, while the growth of Jupiter is traced by the black curve with 1, 3, 10, and 100 M⊕ and full mass in circles. The gas moves only inward in (A) and (B), whereas the gas outward viscous spreading is included in (C). In (A), we assume that Jupiter formed early and that it entirely blocked the inward drifting solid particles after 0.5 Ma, while in (B), Jupiter is assumed to let all dust particles pass through. In the first case, the inner (NC) reservoir is quickly depleted, while in the second case, the CC reservoir quickly takes over the inner disk. In (C), we show that viscous expansion can keep the two disk reservoirs separated within 2 to 3 Ma. Jupiter is here assumed to form at a wide orbit after 1 Ma. In top panels, the black dashed lines represent that Jupiter is not an efficient dust filter, while the black solid line assumes that solid particles are all blocked by Jupiter. The vertical dashed lines in the bottom panels mark when the NC and CC reservoirs get truncated. For example, in the left bottom panel, the NC reservoir is depleted at t = 0.45 Ma, while Jupiter blocks the CC reservoir at t = 0.5 Ma. The disk parameters are vF = 1.5 m/s, rd = 8 AU, Cd = 0.14 M⊙, Tirr0 = 130 K, κ0 = 10−2, αA = 8 × 10−3, and αD = 10−3.

Thus, a second scenario is considered in Fig. 2B where Jupiter allows material to permeate past the gap and replenish the inner reservoir. However, in this circumstance, the CC reservoir quickly overtakes the inner regions of the protoplanetary disk, and there would be no compositional dichotomy preserved in the primitive NC/CC parent bodies formed after 2 to 3 Ma. Reducing or stopping the radial drift of solids interior of Jupiter could solve this inconsistency (, ), but this suggestion comes at the expense of ignoring or reducing the radial drift of solids—a fundamental physical process operating in protoplanetary disks—to separate the reservoirs.

In this context, it has then been suggested that the NC/CC reservoirs were isolated because of the trapping of pebbles in speculative long-lived pressure bumps in the inner disk () or through continuous and efficient planetesimal formation around the water ice line (, ). The latter model requires that the inner disk experiences either intensive irradiation or strong viscous heating to keep the water ice line away from the terrestrial planet zone and additionally extremely high planetesimal formation efficiencies that convert hundreds of Earth masses of CC material to planetesimals. This raises the concern that a large fraction of those planetesimals would subsequently become implanted in the asteroid belt due to the migration of Jupiter ().

Here, we instead propose that the viscous expansion of the gas naturally separates the two disk reservoirs (Fig. 2C). As such, there is no need to invoke very early Jupiter formation to explain the NC-CC separation. In this case, the core of proto-Jupiter forms in a wide orbit and then undergoes inward migration. Adopting a state-of-the-art pebble accretion and planet migration model (), we provide an illustration of Jupiter’s growth and migration track in Fig. 2C. In our model, Jupiter finishes the growth of its 10 M⊕ core at a time of t ≃ 3 Ma in the 10-AU disk region and reaches its full mass (by gas accretion) and the current position at t ≃ 4.2 Ma.

These two contrasting models for the NC/CC dichotomy in the context of the known constraints on parent body formation timings are illustrated in Fig. 3 (see also Materials and Methods). The expansion of the protoplanetary disk is a key physical mechanism needed to understand the asteroid and planetary record of the solar system. As shown in Fig. 2C, with viscous disk spreading, the outer CC material arrives in the inner disk region only after 2 to 3 Ma. This expansion allows the parent bodies of NC differentiated meteorites, achondrites, enstatite, and ordinary chondrites to form sequentially during the first 2 to 3 Ma in the inner disk region (Fig. 3C). The later formed ordinary chondrites reveal a modest compositional deviation from earliest formed ureilites toward CC chondrites, implying a continuous admixing between the original inner NC disk reservoir and outer CI-like material (e.g., ureilites and CI chondrites as end members for the mixing) (). This is a natural consequence of the planetesimals forming at the time where the NC and CC reservoirs partly mix, as illustrated by the shaded region in Fig. 3C. The small (<100-km) NC asteroids are protected against further accretion of substantial amounts of drifting CC pebbles. The pebble accretion efficiency depends strongly on the inclinations and eccentricities of the planetesimals. The internal gravitational stirring of the primitive NC asteroids will quickly raise their eccentricities to a degree where subsequent accretion of CC material is insignificant (see Materials and Methods). A similar early excitation of the orbit of Mars has been proposed to explain how Mars terminated its accretion earlier than Earth and thus incorporated a lower fraction of CC material ().

Fig. 3.

Meteorite dating records and a schematic illustration of the two solar system formation scenarios.

(A) Symbols give the peak accretion ages of meteorites with uncertainties shown in horizontal lines. (B and C) The magenta and blue lines show example trajectories of NC and CC dust particles, while the dashed line in the viscous spreading scenario is the separation radius for gas inflow/outflow. For the early Jupiter scenario, the core grows in situ within 1 Ma and blocks all particles drifting from the outer part of the disk. The inner disk drains out solids quickly, which conflicts with the later formation of ordinary/enstatite chondrites and the incorporation of CC material into the terrestrial planets. For the viscous spreading scenario, Jupiter’s core forms beyond (and slightly later than) the CC differentiated parent bodies. The parent bodies of primitive NC meteorites (ordinary chondrites and enstatite chondrites) form at 2 to 3 Ma, while the inner disk region is changing composition. The primitive CC parent bodies form at t ∼ 3 to 4 Ma and are scattered into the asteroid belt together with the CC differentiated parent bodies by Jupiter’s migration.

The differentiated meteorite parent bodies of CC composition form at relatively large orbital distances (r > rCC) within 1 Ma, whereas the parent bodies of CC chondrites form at t ≃ 3 to 4 Ma exterior to the orbit of Jupiter and are later scattered into the asteroid belt by the combined gravity of Saturn and Jupiter (). The terrestrial planets—proto-Earth and Mars—accrete considerable CC material and complete their main accretion within 4 to 5 Ma (, , ). Our model allows Jupiter’ core to finalize its growth at a relatively late time (e.g., t ∼ 2 to 3 Ma) and at large orbital distances, such as r > 10 AU, in agreement with planet formation models that include inward migration due to gravitational torques exerted by the protoplanetary gas disk on the planet (, ). Large-scale inward migration of proto-Jupiter plays a decisive role in sweeping through the formation zone of CC differentiated planetesimals. The orbital crossings and close encounters between Jupiter and those planetesimals result in the CC iron meteorite parent bodies to be frequently scattered into the inner asteroid belt region ().

DISCUSSION

We first discuss the role of Jupiter on dust filtration. Hydrodynamic studies find that, while large millimeter-to-centimeter size particles are stopped at the pressure maximum exterior to the planet, small grains that are well coupled to the disk gas can pass the gap (–). How the above size-dependent filtering affects the isotope anomalies is not consistently accounted for in the early Jupiter scenario (). We show in Fig. 4 the efficiency of Jupiter as a dust filter, taking into account a size distribution of solids. Two different viscosity prescriptions are considered: a constant α everywhere in the disk (blue) and a higher value exterior to the gap (red). The latter assumes an elevated viscosity due to the generation of vortices by the Rossby wave instability associated with gap opening in the inner regions of the disk where the background turbulence is inherently low (). The critical Stokes number of particles to cross the gap decreases as Jupiter’s mass increases, whereas the maximum Stokes number of particles in the fragmentation regime remains constant at the planet location. When the mass of Jupiter is low, all particles are below the critical Stokes number and can drift through the gap. Jupiter begins to block particles only when the critical Stokes number is lower than the fragmentation-limited particle size. Although the maximum size particles dominate the total mass of the whole dust population, a nonnegligible fraction of them would unavoidably fragment into smaller particles. We assume that the dust size distribution follows a power law, n(a) ∝ a, where p = − 3.5, amin = 0.1 μm, and amax is the maximum size in the fragmentation limit. We find that a dust fraction of (acrit/amax)4− still diffuse through the gap region when amax > acrit. As shown in Fig. 4E, when Jupiter grows rapidly, 24 and 65% of CC dust would drift through Jupiter’s gap after 5 Ma at αgap = 0.001 and 0.003. This dust filtering works much less effectively for a 3-Ma gas accretion time scale in Fig. 4F.

Fig. 4.

Evolution of Jupiter’s mass, particle Stokes numbers, and CC dust mass interior to Jupiter after its core formation.

(A and B) The mass of Jupiter starting from 20 M⊕ at t = 1 Ma, with a characteristic growth time scale of 0.5 (left) and 3 Ma (right). (C and D) Critical (maximum) Stokes number τs, crit that particles can diffuse through the gap versus disk particles’ Stokes number in the fragmentation regime for αgap = 0.001 (blue) and 0.003 (red). The blue dots and red triangles refer to τs,crit(Mp) at two different αgap, where the data from the low masses and Jupiter’s full mass are adopted from () and setup I from (), respectively. From (), we calculate the critical sizes when the dust-to-gas ratio interior of the planet drops below 50% of the exterior value and convert those into Stokes numbers. (E and F) The mass of CC dust particles that drift across the Jupiter’s gap for αgap = 0.001 (blue) and 0.003 (red). The gray line shows the CC dust drift exterior to the orbit of Jupiter for comparison. The other disk parameters are α = 10−3, Tirr0 = 130 K, Cd = 0.14 M⊙, vF = 1 m/s, and rd = 10 AU. The growing Jupiter lets a significant fraction of the dust mass past its orbit, unless it grows very fast and has a low turbulent diffusion outside the gap.

We note that larger centimeter-sized CAIs, condensed in the inner regions of the protoplanetary disk and likely transported to the outer disk by winds and outflows, will be efficiently trapped by the growing Jupiter and their igneous nature will prevent the formation of CAI fragments that could pass the gap, and hence, this can explain how the ordinary chondrites avoided accreting CAIs (, ).

Our scenario supports that the disk angular momentum is transported outward by the turbulent viscosity, for instance, via the magnetorotational instability (MRI) (). Although MRI turbulence is likely to be suppressed due to nonideal magnetohydrodynamic (MHD) effects in the inner 10-AU disk region where disk winds could drive a laminar accretion flow, it can operate at large disk radii where the CC dust particles reside initially. In the outer disk region where ambipolar diffusion dominates, the turbulent level depends on the strength of the net magnetic field, dust abundance, and extreme ultraviolet/cosmic ray penetration depth (). The disk turbulence can be inferred from the spectral energy distribution or molecular line velocity measurements. Such observations show that some disks display a relatively high α ≳ 10−3 − 0.1 (, ), compatible with measurements from nonideal MHD disk simulations (). Moreover, a recent study found that, even in the magnetized wind-driven disk where the MRI is fully quenched, the magnetic torque can also drive the expansion of the gas in the outer disk region, analogous to the viscous spreading effect (). This implies that our expanding disk scenario can be applicable for a broad range of circumstances. Despite uncertainties in measuring disk size, the observed sizes of class II disks are plausibly larger than those of class I disks, in alignment with the above theoretical predictions ().

In summary, we propose that a generic outward-then-inward radial motion of pebbles, driven first by the viscous expansion of the disk gas and thereafter by the radial drift of pebbles, is the key ingredient that explains the isotopic compositions of the entire known meteoritic and planetary record. In this view, the observed NC/CC dichotomy is set by a secular change to the composition of the disk and not by an early spatial separation by an in situ–formed Jupiter (, ). We show that the gap opening by Jupiter will not stop all inward-drifting solid particles and that Jupiter lets considerable amounts of CC dust and pebble fragments enter the inner disk to contribute to the mass budgets of the terrestrial planets and to the parent bodies of the ordinary chondrites.

Our work implicitly supports a distant origin of Jupiter, which is consistent with orbital migration theory (, , ). Many heavy elements including ultravolatile noble gases are found to be enriched in Jupiter’s atmosphere compared to the solar photospheric values. This feature cannot be explained by in situ accretion and therefore strongly favors that Jupiter’s core formed in a distant, cold orbit (Fig. 3C) (). In this view, heavy elements can be accreted as solids that would later become diluted into the envelope, resulting in the observed enrichment pattern.

MATERIALS AND METHODS

Chronology of meteorites

We summarize the literature chronological meteorite studies in Fig. 3A. The ages of chondritic parent bodies can be deduced from the time when chondrules were last melted and crystallized. This age is determined from either the long-lived U-corrected Pb-Pb dating or short-lived Al-Mg and Hf-W datings. In the Pb-Pb dating approach, recent studies proposed that chondrules form simultaneously with CAIs and last for 3 to 4 Ma, whereas the age peaks around 0 to 1 Ma (, ). In the Al-Mg dating approach, chondrules typically record a lower initial Al isotopic ratio compared to that of CAIs, the latter of which is treated as the canonical value for the onset of the solar system (26Al/27Al =5.5 × 10−5). This is interpreted as a 1- to 4-Ma–protracted formation of chondrules after CAIs, with a major age peak at t = 2 to 3 Ma (–). Note that Pb-Pb dating of single chondrules shows different ages from that of the pooled chondrule separates, the latter of which seems more aligned with the results obtained from the Al-Mg dating (). On the other hand, the short-lived Al-Mg dating relies on the assumption of an initial homogeneous distribution of 26Al throughout the inner solar system, the validity of which is subject to debate in the community (, ).

On the other hand, the age of meteorites can instead be inferred from the isotope dating of specific physical/chemical events occurred in their history. For iron meteorites (differentiated planetesimals), Hf-W dating probes the ages when their parent bodies’ cores became segregated from mantles. Then, the formation time of these bodies can be traced back via thermal modeling of their internal evolution based on the assumption that 26Al is the predominant heat source. Similar approaches have been used for achondrites and ordinary and CC chondrites, with different dating events such as magma crystallization, formation of minerals like carbonates, and secondary fatalities. In such a way, Vesta and the parent bodies of ureilites and angrites are deduced to form within 0.5 Ma after the CAI formation (). The parent bodies of iron meteorites are found to form earlier than 0.5 Ma for the NC population and earlier than 1 Ma for the CC population (). The parent body accretion of ordinary chondrites takes place at t ∼ 1.6 to 1.8 Ma (), while those of CC chondrites occur at ∼2.5 to 4 Ma after CAIs, spanning from CO, CV, and CM types (, ).

Furthermore, the accretion history of terrestrial planets can be constrained from their isotopic compositions. Many isotopes such as 48Ca, 50Ti, and 54Cr of Earth, Mars, and Moon are found to be different from either NC or CC asteroid populations (, ). The terrestrial planets originally grew in the NC environment. One expects that these NC/CC disk reservoirs are temporarily separated early on and mixed together later through drifting of CC dust particles. As such, to account for their isotopic compositions, Earth and Mars need to accrete a nonnegligible fraction of CC material (an upper limit of 30 to 40% depending on different estimated isotopes) during the protoplanetary disk phase. Inferred from Hf-W chronometry, Mars attains its half mass within 3 Ma (), while the formation of Earth takes much longer, approximately at 30 to 100 Ma (). The above Earth age derivation, however, could be substantially affected by the Moon-forming impact. It is plausible that proto-Earth had completed the main accretion within the gas disk lifetime and the last Moon-forming giant impact happened a few tens of million years later (). The early and rapid formation of proto-Earth is more consistent with iron isotope evidence (). Johansen et al. () proposed a pebble accretion terrestrial planet formation model that explains their architecture, bulk compositions, and isotopic variability. They suggested that Mars and proto-Earth accrete CC materials 3 to 4 Ma after the CAI formation.

Spreading disk model

We want to track the 1D evolution of the gas and solid components in the protoplanetary disk. A conventional disk model is constructed, including both the gas spreading and dust radial drift. We focus on the stage when the star completed its main accretion and disk has already formed. Thus, the disk-cloud interaction and infall of envelope material onto the disk are neglected for our fiducial setup. The influence of the early infall is further addressed in the “Early disk evolution” section. Although the disk formation and evolution are found to be much more complicated in zoom-in 3D hydrodynamic simulations (), the advantage of our simplified 1D treatment is that it captures the essential physical processes, at the same time with the feasibility of testing the influence of model parameters at a considerable detail. Such explorations can inspire and guide future sophisticated numerical modelings.

We recapitulate the classical theory of viscous spreading disk (, ). On the basis of the assumption that viscosity increases linearly with radial distance ν = νd(r/rd), the characteristic time scale of viscous accretion/spreading is expressed aswhere rd is the characteristic disk radius at t = 0 years and νd is the corresponding viscosity at rd. Here, we define the dimensionless time as τν = 1 + t/tν. The gas surface density and disk accretion rate from the self-similarity solution are written as ()andwhere Md0 is the initial gas disk mass and r is the radial distance from the central star. The gas velocity is given by

Hereafter, we consider a realistic two-component disk structure based on different heating sources. The inner disk region is heated by viscous dissipation, while the outer disk region is dominated by stellar irradiation (). The standard α prescription is adopted for viscosity ν = αRgT/(μΩK) (), where Rg is the gas constant, T is the disk temperature, μ is the gas mean molecular weight, is the angular velocity, and M⋆ is the stellar mass. The self-similarity solution is applicable for the outer stellar irradiated disk region. The gas temperature and surface density are given bywhere ϕ is the grazing angle where the stellar light hits the disk surface, L⋆ is the stellar luminosity, Cd is the equivalent disk mass when the disk is purely stellar irradiated from Eq. 5, σSB is the Stefan-Boltzmann constant, and AU is the mean distance between Earth and Sun. Note that ϕ is much lower than the unit and should principally depend on r in a sophisticated disk model (leading to Tirr ∝ r−3/7) (). Here, we simply adopt it to be a constant. For the fiducial disk model, when ϕ = 0.01 and L⋆ = 1 L⊙, the irradiation temperature Tirr0 at 1 AU is 130 K. Because of the variation of ϕ (related to disk opacity) and more luminous young star with a higher L⋆ during its premain sequence, Tirr0 could be varied from 130 to 300 K.

For the inner viscously heated disk region, the viscosity depends on both the disk surface density and radial distance, and hence, the self-similarity solution is not exactly applicable. We assume that, in the inner disk region, the gas flow is maintained by Eq. 6 such that (Σgν)vis = (Σgν)irr. In the inner viscously heated region, the energy balance between heating and cooling provides thatwhere Te is the effective gas temperature at the disk surface. The heat that is generated from viscous dissipation in the midplane radiates vertically toward to the surface layer. For an optically thick radiative disk, the midplane temperature gives and the optical depth is τ = κΣg. We follow a simple disk opacity law κ = κ0T (), and then, the energy equation readsWe obtainandwhereand ΩK0 is the Keplerian angular velocity at 1 AU.

The transition radius rtran between these two disk regions can be solved by equating Tvis and Tirr

The initial disk mass isand the outer edge of the disk (or disk size, rout) defines that includes 95% of its mass. The above two quantities can be solved numerically. For the fiducial setup in Table 1, the initial disk mass Md0 is 0.065 M⊙ for the uniform turbulence model and 0.06 M⊙ for the dead zone model. Figure 5 illustrates the evolution of gas temperature and surface density for our fiducial setup.

Table 1.

Model setups for the parameter study.

The radius of particles is chosen to be either 0.1 (fiducial), 0.3, or 1 mm in the bouncing-limited growth regime or calculated from their Stokes number of in the fragmentation-limited growth regime where vF = 1 (fiducial), 0.3, or 3 m/s.

Run	Dead zone		Uniform turbulence α	T irr0 (K)	C d (M⊙)	v F (m/s)	R peb (mm)	r d (AU)	ρ• (g cm−3)
	α( = αA)	αD
Fiducial	10−2	5 × 10−4	5 × 10−3	130	0.14	1	0.3	10	1.5
Low α	10−3	5 × 10−4	10−3	130	0.14	1	0.3	10	1.5
Hot disk	10−2	5 × 10−4	5 × 10−3	280	0.14	1	0.3	10	1.5
Low − mass disk	10−2	5 × 10−4	5 × 10−3	130	0.014	1	0.3	10	1.5
Large size disk	10−2	5 × 10−4	5 × 10−3	130	0.14	1	0.3	20	1.5
Porous particle	10−2	5 × 10−4	5 × 10−3	130	0.14	1	0.3	10	0.15
High vF	10−2	5 × 10−4	5 × 10−3	130	0.14	3	∖	10	1.5
Low vF	10−2	5 × 10−4	5 × 10−3	130	0.14	0.3	∖	10	1.5
Large particle	10−2	5 × 10−4	5 × 10−3	130	0.14	∖	1	10	1.5
Small particle	10−2	5 × 10−4	5 × 10−3	130	0.14	∖	0.1	10	1.5

Fig. 5.

Two-component structure for a spreading disk.

(A) Disk temperature as a function of radial distance. (B) Gas surface density as a function of radial distance. The thickness of line represents different time of t = 0 years and 0.5, 1, and 2 Ma, respectively. (C) Time evolution of the transition radius. The parameters are adopted from the fiducial run with a uniform turbulence model in Table 1.

We note that the disk angular momentum is transported from the inner region of the disk to the outer region through the viscous stress. As such, most of the disk gas moves inward, while the remaining gas expands to larger radii to conserve angular momentum. On the basis of Eq. 7, the radius that separates mass inflow and outflow is expressed as

At the disk radius larger than rexp, the mass flux moves outward and ug > 0; at the disk radius smaller than rexp, the mass flux moves inward and ug < 0. Inserting Eq. 17 into Eq. 7, we rewrite the gas velocity as

The radial drift velocity of solid particle is given by ()where vK ≡ ΩKr is the Keplerian velocity. The Stokes number of the particle τs ≡ ΩKtstop, where tstop is the stopping time that measures how fast the particle adjusts its momentum to the surrounding gas. The headwind prefactor that measures the disk gas pressure gradient is given bywhere hg = Hg/r is the disk aspect ratio, Hg is the gas scale height, and Cη is an order unit factor related to the disk structure. Combining Eqs. 18 and 19, the outward movement of dust particles requires

The final size that particles can reach depends on the growth pattern. We consider two circumstances: Particles are either in fragmentation-limited or bouncing-limited regimes (). In the fragmentation-limited growth regime, the highest Stokes number of particles is given by , where cs = HgΩK is the gas sound speed and vF = 1 m/s is the adopted fragmentation threshold speed, related to the composition and size of the monomer grains (, ). The turbulent speed of the dust particles is given by . Inserting it into Eq. 21, we obtain the criterion of particles’ outward drift in the fragmentation regime

As can be seen from the above equation, a high α and hot disk promotes the outward movement of dust particles.

On the other hand, some laboratory experiments find that the growth of the particles is stalled at around 0.1 to 1 mm when the bouncing effect becomes dominant (, ). In this work, we assume that particles in the bouncing-limited growth regime reach a radius of 0.3 mm. The particles’ Stokes number readswhere ρ• is the internal density of the dust particles and R0 is the radius of the bouncing particles. Inserting Eq. 23 into Eq. 21, we obtain the criterion of particles’ outward drift in the bouncing regime

We assume that dust particles have already reached the maximum sizes at the beginning of numerical integration for the above two cases. This is valid as long as the growth time scale is shorter than the viscous spreading time scale. The growth time scale from micrometer-sized grains to millimeter-sized particles that initiate substantial radial drift is given by tgrow ∼ ξ(ZΩK)−1, where ξ ∼ 10 and Z = 10−2 is the canonical dust-to-gas ratio, while the viscous time scale is . It is easy to calculate that tgrow ≪ tν for standard disk parameters.

The turbulent velocity is isotropic when the MRI is the predominant disk turbulence source (). The MRI is not able to operate when the ionization level is low. Hence, a layered accretion model has also been proposed, where the turbulence in the midplane is strongly damped and accretion is mainly through an ionized, active surface layer (). We account for both these two models: One assumes a uniform disk turbulence (constant α everywhere), and the other incorporates the effect of the dead zone. In the latter model, α—which sets the global viscous evolution of the disk—is much higher than the turbulent level αt in the quiescent midplane layer. The ionization relates to the column density of the disk gas, which increases with orbital distance. We adopt a simple prescription for αt in the dead zone modelwhere αD and αA = α are the low and high boundary of turbulent viscosity parameter and rin is the inner disk radius where the disk temperature is 1000 K.

Parameter study

We investigate how the radial drift of dust particles influenced by different parameters in both the uniform turbulence and dead zone models. The particles are either in the fragmentation-limited or in the bouncing-limited growth regimes. The parameter setups are provided in Table 1, while the simulation results are depicted in Figs. 6 and 7. For the size of the particles limited by the fragmentation threshold, we find in Fig. 6 that particles maintain their outward movement for a less time when the disk is less turbulent, the temperature is lower, the initial disk size is larger, and the fragmentation threshold velocity is higher. Meanwhile, the outward-then-inward drifting particles that eventually enter the asteroid belt have a larger size when the disk is less turbulent, the temperature is lower, the initial disk size is larger, the disk is more massive, the fragmentation threshold velocity is higher, and/or the particles are more porous with a lower internal density.

Fig. 6.

Parameter study on the radial drift of dust particles with fragmentation-limited sizes.

The uniform turbulence (left) and dead zone (right) disk models. The black curves trace the radial drift of dust particles, while the gray dashed lines indicate the gas spreading at the same initial disk locations. The color refers to the size of the particle at the corresponding time and disk location. The model setups are presented in Table 1.

Fig. 7.

Parameter study on the radial drift of dust particles with bouncing-limited sizes.

The uniform turbulence (left) and dead zone (right) disk models. The black curves trace the radial drift of dust particles, while the gray dashed lines indicate the gas spreading at the same initial disk locations. The color refers to the Stokes number of the particle at the corresponding time and disk location. The model setups are presented in Table 1.

We note that the fragmentation velocity is a crucial factor in determining the duration of the NC-CC separation, because the Stokes number of the particles directly correlates to . We find that very large pebbles form when the fragmentation threshold is high. Such particles quickly decouple from the diffusing gas and drift into the inner disk within less than 1 Ma (see the high vF case in Fig. 6). It is therefore difficult to generate the NC-CC separation if the fragmentation velocity is very high.

For the size of the particles limited by the bouncing barrier, we find in Fig. 7 that particles beyond certain disk radii (≳20 AU) drift in a convergent manner. Furthermore, out particles have high Stoker numbers that already decouple from the disk gas to move inward, while particles with slightly inner orbits still couple with the outwardly spreading gas. The bouncing-limited size is also varied at 0.1, 0.3, and 1 mm. When the dust particles reach a smaller size, they can couple to the expanding gas for a longer time, resulting in a longer NC-CC separation. Similarly, porous dust particles with a lower internal density have a lower Stokes number compared to those with a compact structure. Therefore, the porous dust particles can maintain the outward expansion for a longer time. To summarize, to separate the two disk reservoirs for more than 2 to 3 Ma, a lower bouncing-limited size or particles with a lower internal density are preferred.

To analyze the results statistically, we also perform a high number of Monte Carlo sampling simulations by varying two disk parameters at each time. The CC material injection radius rCC is uniformly chosen from 5 to 35 AU, the initial disk size rd is uniformly chosen from 10 to 40 AU, the disk temperature at 1 AU ranges from 120 to 400 K, and α is log uniformly selected from 10−3 to 10−2 in the uniform turbulent model or αD from 10−4 to 10−3 in the dead zone model. The particles in the bouncing-limited regime are assumed to be 0.6 mm in diameter. The results are given in Fig. 8.

Fig. 8.

Monte Carlo simulation outcomes as a function of disk parameters.

The dead zone and uniform turbulence models with particles in the fragmentation- and bouncing-limited regimes are illustrated from top to bottom. The color corresponds to the time of CC dust particles entering the inner disk region of r < 1 AU. The circles with black edge color refer to the systems that CC particles entering the terrestrial planet formation region during 3 to 4 Ma after the CAI formation. We fix the other model parameters and vary the initial CC material injection radius versus midplane α (left), initial disk size (middle), and Tirr0 (right), respectively. The default parameters are αD = 5 × 10−4, αA = 10−2 (for the dead zone model), or α = 5 × 10−3 (for the uniform turbulence model); rd = 10 AU; vF = 1 m/s; Tirr0 = 130 K; and Cd = 0.14 M⊙, respectively.

The duration of the NC-CC separation in the expanding disk model is strongly influenced by the properties of the protoplanetary disk. Figure 8 shows parameter dependencies of the arrival time of CC particles in the inner disk. For the explored range of parameters, the time for CC particles reaching the inner 1-AU disk region spans from 105 to 107 years, with a later entering time for a larger rCC. Overall, the inner boundary of CC material larger than 20 AU is needed for the corresponding drift time of greater than a few million years. The drift time scale of particles beyond rCC sensitively depends on the initial disk size, because it directly correlates to places where disk gas can diffuse inward/outward (rexp).

Early disk evolution

Our study separates the disk-forming stage from the disk evolution stage. We define the former phase as the period in which gas from the parent molecular cloud accretes onto the central protostellar core and the surrounding circumstellar disk, while in the latter phase, this process has come to an end. The disk is then approximated as a closed and isolated system where gas accretion onto the host star is paired with an expanding outer disk, in line with conventional disk theory and observations (). In the “Results” section, we show that this disk evolution phase is sufficient to explain the observed isotopic dichotomy. There is, however, a significant degree of uncertainty on the earlier disk building phase (, , ).

Here, we show that our assumptions on the disk formation stage are in line with theoretical and observational constraints. First, disks likely form rapidly: A conventional Bonnor-Ebert sphere collapses within a free-fall time scale of Ma, where we take the core mass to be a solar mass and a critical radial extend of rBE to be ≈0.1 pc. Although some early condensates like CAIs may form within this early period, most solids form in the disk evolution phase (Fig. 3). This estimated free-fall time scale is compatible with the range of disk formation time scales measured in the state-of-the-art magnetohydrodynamical simulations (, , ).

It is less clear how exactly mass is deposited onto the disk during this disk-forming phase. In the classical picture, not considering magnetic fields, material will collapse inside out when the infalling material conserves its angular momentum. The cloud-collapsing front expands outward, with high angular momentum material deposited onto the disk at a later time, e.g., ().

Current MHD simulations similarly argue for inside-out formation without a strong imprint of magnetic braking: Ideal MHD simulations typically show a strong and changing misalignment between the total angular momentum and magnetic field vectors of young disks, when taking into account the turbulent nature of the giant molecular cloud (). Therefore, disks undergo relatively isotropic accretion—with strong local and temporal variation—of gas at ∼50 AU that results in a radial mass flow in the disk midplane (, ). We therefore assume that, within the infall phase, the gas only accretes inward. After tinfall, the disk evolution switches to the nominal situation described in the “Spreading disk model” section. Recently, it has been argued that magnetic braking and vertical accretion play a larger role (, ). However, a statistical analysis of these results, starting from the molecular cloud level, shows that disk properties for ideal and nonideal MHD (including ambipolar diffusion) simulations appear to be largely similar, but, in the latter case, disks tend to be less massive (<0.1 M⋆) (). Regardless, the current mass and radial extent of the solar system points toward it not forming from such a potential population of small, low-mass disks, but instead from an intermediate-sized disk of size >50 AU.Moreover, outer solar system solids with a low degree of thermal processing, like comets and CI chondrites, were more likely delivered radially from the outer regions of the collapsing cloud. This is in line with the model presented here with an outer, more pristine, CC component. Delivering the bulk material closer to the host star would give a potentially excessive degree of thermal processing.

Jupiter’s growth and migration

We adopt the state-of-the-art pebble-driven planet formation model (). The key equations for Jupiter’s growth and migration are recapitulated as follows.

The pebble accretion efficiency in the 2D and 3D limits readand

We take into account of the gas flow velocity in the above equations, where Mp and M⋆ are the masses of the planet and the star, respectively, and Δv is the relative velocity between the planet and the pebble.

The pebble mass flux () is assumed to be attached to the gas flow with a constant mass flux ratio such that . Particles of very low Stokes number are well coupled to the disk gas. The gas and dust particles drift together with the same velocities, and therefore, ξ remains equal to the initial disk metallicity. When the Stokes number is high, pebbles drift faster than the disk gas. In this case, to maintain a constant flux ratio, Σpeb/Σg becomes lower than the nominal disk metallicity, and the pebble surface density gets reduced. The mass accretion onto the planetary core can be written as

The core mass growth is terminated when the planet reaches the pebble isolation mass

After reaching the pebble isolation mass, the planet starts rapid gas accretion

The first term in Eq. 30 corresponds to the Kelvin-Helmholtz contractionwhere κenv is the opacity in of the planetary envelope.

The second term in Eq. 30 indicates how much gas the in the planet Hill is able to get accreted by the planetwhere RH = (Mp/3M⋆)1/3r is the Hill radius of the planet and Σg and Σgap are the unperturbed gas surface density and the surface density at the bottom of the gap, respectively. The relative depth of the gap can be given by

We also define the gap opening mass Mgap when Σgap is reduced by 50% compared to the unperturbed value.

Briefly, the Kelvin-Helmholtz contraction is the dominant gas accretion channel when the planet mass is low. When (dMpg/dt)KH > (dMpg/dt)Hill, the gas accretion is limited by the amount of gas within the planet’s Hill sphere. Once the planet becomes sufficiently massive, the available gas is eventually determined by the global disk accretion flow. In this work, κenv and facc are the two variables that determine the time scale of gas accretion.

The planet interacts with sounding disk gas, inducing orbital migration. Such a planetary migration rate is given by

This formula summarizes for both the low-mass planet in the type I and massive planet in the type II torque regimes. The type I migration rate readswhere fI = 2 is adopted from equation 8 in () approximately for the migration prefactor in a local isothermal disk.

We start the growth track of a protoplanet with 1 M⊕, and t0, r0, ξ, facc, and κenv are the parameters that can be varied to generate the current Jupiter’s mass and location. One example is illustrated in Fig. 2C, where t0 = 105 years, r0 = 45 AU, ξ = 0.04, κenv = 0.05 g/cm2, and facc = 0.9. The other parameters are vF = 1.5 m/s, rd = 8 AU, Cd = 0.14 M⊙, Tirr0 = 130 K, αA = 8 × 10−3, and αD = 10−3 (dead zone model).

Pebble accretion for asteroids on eccentric and inclined orbits

We explain how the inclination and eccentricity affect an asteroid’s growth by pebble accretion. The orbits of the planetesimals are excited through their mutual gravitational interactions. The random velocities of massive embryos remain low, whereas the eccentricities/inclinations of small asteroids become high. There are two consequences. First, the orbits of asteroids are lifted away from the dust midplane layer due to increased inclinations, resulting in fewer pebbles able to be accreted. Second, the pebble-asteroid encounter time becomes short because of their high relative velocity. When the encounter time is shorter than the stopping time of pebbles, the pebble accretion efficiency drops greatly. Both effects are more profound for lower-mass planetesimals.

Figure 9 shows the mass growth of asteroids varying with their eccentricities and inclinations. The parent bodies of ordinary and enstatite chondrites form in the NC reservoir at t = 2 Ma. However, although 100 M⊕ CC dust particles drift across their orbits over the next 1 Ma, these asteroids accrete an insignificant amount to reverse their isotopic signatures, as long as their eccentricities/inclinations are beyond a few times 10−3.

Fig. 9.

Planetesimal mass growth as a function of inclination.

The initial sizes of planetesimals are 100 (thick lines) and 20 km (thin lines), whereas their orbits are either circular (solid lines) or with eccentricities two times of inclinations (dashed lines). The disk is adopted to be the fiducial setup of the uniform turbulence model, and the size of particles is bouncing-limited to be 1 ml. The planetesimals are assumed to form at t = 2 Ma, and the following mass growth by pebble accretion is accounted for when 100 M⊕ pebbles drift bypass their orbits. The pebble accretion prescription is adopted from (, ).