Making the Solar System

We consider a protoplanetary disk that evolves due to angular momentum lost to a disk wind. As a result, gas flows inwards with a mass flux F_gas. The disk is heated by gravitational energy released by the inflowing gas together with reprocessed radiation from the central star. Following Mori et al. (2021) and Kondo et al. (2022), we assume that heat generated by the disk wind is released in a layer at a particular optical depth, τ_wind, with respect to the disk surface.

The midplane temperature, T_mid, at a distance, a, from the star is

$$\begin{eqnarray}&&{T}_{\mathrm{mid}}^{4}={T}_{\mathrm{irr}}^{4}+\displaystyle \frac{{Q}_{\mathrm{wind}}}{2{\sigma }_{\mathrm{SB}}}\left[1+\displaystyle \frac{3{\tau }_{\mathrm{wind}}}{4}\right],\end{eqnarray} \tag{ 1 }$$

where σ_SB is the Stefan–Boltzmann constant, and Q_wind is the rate of energy release per unit disk area, given by

$$\begin{eqnarray}&&{Q}_{\mathrm{wind}}=\displaystyle \frac{{F}_{\mathrm{gas}}{{\rm{\Omega }}}^{2}}{4\pi },\end{eqnarray} \tag{ 2 }$$

where Ω is the angular velocity.

Following Chiang & Goldreich (1997), the temperature T_irr solely due to stellar irradiation, at a distance a is

$$\begin{eqnarray}&&{T}_{\mathrm{irr}}^{4}=\displaystyle \frac{2}{7}{\left(\displaystyle \frac{k}{{{GM}}_{* }\mu {m}_{H}}\right)}^{4/7}{\left(\displaystyle \frac{{L}_{* }}{4\pi {\sigma }_{\mathrm{SB}}}\right)}^{8/7}{a}^{-12/7},\end{eqnarray} \tag{ 3 }$$

where M_* and L_* are the stellar mass and luminosity, k is Boltzmann's constant, m_H is the mass of a hydrogen atom, and μ is the mean molecular weight of the gas. Following Chambers (2009), we assume that the luminosity of a solar-mass protostar at time t is given by

$$\begin{eqnarray}&&{L}_{* }=9.2{\left(1+\displaystyle \frac{t}{0.1\ \mathrm{Myr}}\right)}^{-2/3}{L}_{\odot }.\end{eqnarray} \tag{ 4 }$$

Following Suzuki et al. (2016), the inward radial velocity of the gas is given by

$$\begin{eqnarray}&&{u}_{\mathrm{gas}}={\alpha }_{\mathrm{wind}}{c}_{s},\end{eqnarray} \tag{ 5 }$$

where c_s is the sound speed of the gas and α_wind is a measure of the strength of the disk wind. The gas surface density is

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\mathrm{gas}}=\displaystyle \frac{{F}_{\mathrm{gas}}}{2\pi {{au}}_{\mathrm{gas}}}.\end{eqnarray} \tag{ 6 }$$

Following Hartmann & Bae (2018), we assume that the inward gas flux varies inversely with time, t, following an infall phase lasting for t_disk. The gas flux is given by

$$\begin{eqnarray}&&{F}_{\mathrm{gas}}={F}_{\mathrm{gas},0}{\left[1+\displaystyle \frac{t}{{t}_{\mathrm{disk}}}\right]}^{-1}.\end{eqnarray} \tag{ 7 }$$

We assume that pebbles are injected into the planet-forming region of the disk from larger distances with a mass flux

$$\begin{eqnarray}&&{F}_{\mathrm{peb},0}={{ZF}}_{\mathrm{gas}},\end{eqnarray} \tag{ 8 }$$

where Z is the dust-to-gas ratio, and we assume that the gas and pebble fluxes decline at the same rate, following Lambrechts et al. (2019).

The pebble flux, F_peb, at smaller distances is reduced with respect to F_peb,0 as a result of conversion of pebbles into planetesimals and pebble accretion onto embryos, which are calculated at each location in the disk (see below).

Following Dullemond et al. (2018), pebbles move radially inwards at a speed, u_peb, that depends on the inward advection of the gas, u_gas, and the drift rate, u_drift, with respect to the gas:

$$\begin{eqnarray}&&{u}_{\mathrm{peb}}={u}_{\mathrm{drift}}+\displaystyle \frac{{u}_{\mathrm{gas}}}{(1+{\mathrm{St}}^{2})},\end{eqnarray} \tag{ 9 }$$

where St is the pebble Stokes number, and

$$\begin{eqnarray}&&{u}_{\mathrm{drift}}=\displaystyle \frac{\mathrm{St}}{(1+{\mathrm{St}}^{2})}\displaystyle \frac{{c}_{s}^{2}}{{v}_{\mathrm{kep}}}\displaystyle \frac{d\mathrm{ln}P}{d\mathrm{ln}a},\end{eqnarray} \tag{ 10 }$$

where P is the gas pressure and v_kep = aΩ is the Keplerian velocity.

Thus, the pebble surface density at distance a is given by

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\mathrm{peb}}=\displaystyle \frac{{F}_{\mathrm{peb}}}{2\pi {{au}}_{\mathrm{peb}}}.\end{eqnarray} \tag{ 11 }$$

We assume that pebbles grow quickly by mutual collisions until they reach a maximum size at which collisions begin to cause erosion. Collisional velocities are assumed to be set mostly by turbulence in the gas, such that the pebble Stokes number is given by (Ormel & Cuzzi 2007)

$$\begin{eqnarray}&&\mathrm{St}=\displaystyle \frac{{v}_{\mathrm{frag}}^{2}}{3{\alpha }_{\mathrm{turb}}{c}_{s}^{2}},\end{eqnarray} \tag{ 12 }$$

where α_turb is a measure of the turbulence strength and v_frag is the collision velocity at which collisions become erosional. We allow for different material strengths inside and outside the ice line, assumed to lie at the distance where T_mid = 160 K, such that

$$\begin{eqnarray}\begin{array}{rcl}{v}_{\mathrm{frag}} & = & {v}_{\mathrm{rock}}\qquad T\gt 160\ {\rm{K}}\\ & = & {v}_{\mathrm{ice}}\qquad T\lt 160\ {\rm{K}}.\end{array}\end{eqnarray} \tag{ 13 }$$

We use a fixed value of v_rock = 100 cm s⁻¹, which is a typical experimental result for silicates (Guttler et al. 2010). The ice fragmentation speed is a model parameter.

In addition, we assume that the pebble flux decreases by a factor of 2 moving inwards across the ice line due to evaporation of the icy component.

Following Bitsch et al. (2018), an embryo with a mass exceeding the pebble-isolation mass, M_iso, will halt the inward flux of pebbles. As a result, pebble accretion and planetesimal formation will be halted at this distance and all distances closer to the star, where

$$\begin{eqnarray}&&{M}_{\mathrm{iso}}=\left(25+\displaystyle \frac{{\lambda }_{\mathrm{crit}}}{0.00476}\right){f}_{\mathrm{fit}}{M}_{\oplus },\end{eqnarray} \tag{ 14 }$$

where

$$\begin{eqnarray}&&{\lambda }_{\mathrm{crit}}=\displaystyle \frac{{\alpha }_{\mathrm{turb}}}{2\mathrm{St}}\end{eqnarray} \tag{ 15 }$$

and

$$\begin{eqnarray}\begin{array}{rcl}{f}_{\mathrm{fit}} & = & {\left(\displaystyle \frac{{H}_{\mathrm{gas}}}{0.05a}\right)}^{3}\left[0.34{\left(\displaystyle \frac{-3}{{\mathrm{log}}_{10}{\alpha }_{\mathrm{turb}}}\right)}^{4}+0.66\right]\\ & & \times \,\left[\displaystyle \frac{7}{12}-\displaystyle \frac{1}{6}\displaystyle \frac{d\mathrm{ln}P}{d\mathrm{ln}a}\right].\end{array}\end{eqnarray} \tag{ 16 }$$

We use the model of Lenz et al. (2019) for the formation of planetesimals and planetary embryos. Planetesimals are assumed to form continually from pebbles in most regions of the disk such that the planetesimal surface density, Σ_pml, increases at a rate

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{d{{\rm{\Sigma }}}_{\mathrm{pml}}}{{dt}} & = & \left(\displaystyle \frac{\epsilon {u}_{\mathrm{drift}}}{5{H}_{\mathrm{gas}}}\right){{\rm{\Sigma }}}_{\mathrm{peb}}\qquad \mathrm{St}\gt {\mathrm{St}}_{\min }\\ & = & 0\,\mathrm{St}\lt {\mathrm{St}}_{\min },\end{array}\end{eqnarray} \tag{ 17 }$$

where H_gas is the gas scale height and is an efficiency factor for the conversion of pebbles into planetesimals. Following Lenz et al. (2019), we assume that planetesimal formation only occurs for pebbles with a Stokes number above a critical value, ${\mathrm{St}}_{\min }$. Both and ${\mathrm{St}}_{\min }$ are parameters in our simulations.

Numerical simulations suggest that planetesimals form with a range of sizes, typically ∼100 km, with a decreasing probability at larger sizes (Johansen et al. 2015; Simon et al. 2016; Li et al. 2019). We model this by assuming that most planetesimals form with a radius R_pml = 50 km. However, when planetesimals are first able to form at a particular distance, we assume that a small number of larger planetary embryos form at this location. The embryos have initial mass M_emb,0.

We ignore the formation of new embryos in a particular region after the first ones have formed. Instead, we assume that existing embryos grow by accreting pebbles and planetesimals, and embryos occasionally merge such that their mean orbital spacing is maintained at b Hill radii, where b ≃ 10 (Kokubo & Ida 1998).

Following Ormel & Klahr (2010), a planetary embryo of mass M_emb accretes pebbles at a rate

$$\begin{eqnarray}&&\displaystyle \frac{{{dM}}_{\mathrm{emb}}}{{dt}}={{\rm{\Sigma }}}_{\mathrm{peb}}{v}_{\mathrm{rel}}\times \max \left[\displaystyle \frac{\pi {R}_{\mathrm{cap}}^{2}}{2{H}_{\mathrm{peb}}},2{R}_{\mathrm{cap}}\right],\end{eqnarray} \tag{ 18 }$$

where H_peb is the scale height of the pebbles, given by (Youdin & Lithwick 2007)

$$\begin{eqnarray}&&{H}_{\mathrm{peb}}={H}_{\mathrm{gas}}{\left(\displaystyle \frac{{\alpha }_{\mathrm{turb}}}{{\alpha }_{\mathrm{turb}}+\mathrm{St}}\right)}^{1/2},\end{eqnarray} \tag{ 19 }$$

and R_cap is the capture radius for pebble accretion, given by

$$\begin{eqnarray}\begin{array}{rcl}{R}_{\mathrm{cap}} & = & {r}_{H}\times \min \left[{\left(\displaystyle \frac{12{r}_{H}\mathrm{St}}{a\eta }\right)}^{1/2},2{\mathrm{St}}^{1/3}\right]\\ & & \times \,\exp \left[-{\left(\displaystyle \frac{\mathrm{St}}{{\mathrm{St}}_{\mathrm{crit}}}\right)}^{0.65}\right],\end{array}\end{eqnarray} \tag{ 20 }$$

where r_H is the embryo's Hill radius, and

$$\begin{eqnarray}&&\eta =\displaystyle \frac{1}{2}{\left(\displaystyle \frac{{c}_{s}}{{v}_{\mathrm{kep}}}\right)}^{2}\displaystyle \frac{d\mathrm{ln}P}{d\mathrm{ln}a}.\end{eqnarray} \tag{ 21 }$$

The exponential term in Equation (20) is an empirical factor that accounts for the reduced efficiency of pebble accretion for large pebbles, where

$$\begin{eqnarray}&&{\mathrm{St}}_{\mathrm{crit}}=\min \left[12{\left(\displaystyle \frac{{r}_{H}}{a\eta }\right)}^{3},\displaystyle \frac{2}{3}\right].\end{eqnarray} \tag{ 22 }$$

In addition, v_rel is the pebble-embryo relative velocity just before an encounter, given by

$$\begin{eqnarray}&&{v}_{\mathrm{rel}}={v}_{\mathrm{kep}}\times \max \left[\eta ,\displaystyle \frac{3{r}_{H}{\mathrm{St}}^{1/3}}{a}\right].\end{eqnarray} \tag{ 23 }$$

Following Chambers (2014), an embryo accretes planetesimals at a rate

$$\begin{eqnarray}&&\displaystyle \frac{{{dM}}_{\mathrm{emb}}}{{dt}}=\displaystyle \frac{\pi {R}_{\mathrm{emb}}^{2}{v}_{\mathrm{rel}}{{\rm{\Sigma }}}_{\mathrm{pml}}}{2{H}_{\mathrm{pml}}}\left[1+\displaystyle \frac{2{{GM}}_{\mathrm{emb}}}{{R}_{\mathrm{emb}}{v}_{\mathrm{rel}}^{2}}\right],\end{eqnarray} \tag{ 24 }$$

where R_emb is the radius of the embryo, Σ_pml is the surface density of planetesimals, v_rel is the relative velocity shortly before an encounter, given by

$$\begin{eqnarray}&&{v}_{\mathrm{rel}}={v}_{\mathrm{kep}}\times \max \left[e,\displaystyle \frac{{r}_{H}}{a}\right],\end{eqnarray} \tag{ 25 }$$

where e is the mean planetesimal eccentricity, and H_pml is the planetesimal scale height, given by

$$\begin{eqnarray}&&{H}_{\mathrm{pml}}={ai}\simeq \displaystyle \frac{{ae}}{2},\end{eqnarray} \tag{ 26 }$$

where i is the mean planetesimal inclination, and the last approximation is appropriate for a planetesimal in a dispersion-dominated regime.

We calculate the evolution of e taking into account viscous stirring from the embryos and damping due to gas drag. The viscous stirring rate is

$$\begin{eqnarray}&&\displaystyle \frac{{{de}}^{2}}{{dt}}=\displaystyle \frac{{M}_{\mathrm{emb}}{P}_{\mathrm{VS}}}{3{M}_{* }{bP}},\end{eqnarray} \tag{ 27 }$$

where P is the orbital period, b is the mean orbital spacing of embryos in Hill radii, and

$$\begin{eqnarray}&&{P}_{\mathrm{VS}}=73{\left[1+0.11\displaystyle \frac{{a}^{2}{e}^{2}}{{r}_{H}^{2}}\right]}^{-1}.\end{eqnarray} \tag{ 28 }$$

The gas drag damping rate is

$$\begin{eqnarray}&&\displaystyle \frac{{{de}}^{2}}{{dt}}=-\displaystyle \frac{2{e}^{2}}{{t}_{\mathrm{drag}}}\left[1+{\left(\displaystyle \frac{e}{\eta }\right)}^{2}\right],\end{eqnarray} \tag{ 29 }$$

where

$$\begin{eqnarray}&&{t}_{\mathrm{drag}}=\displaystyle \frac{6\rho {R}_{\mathrm{pml}}}{{\rho }_{\mathrm{gas}}\eta {v}_{\mathrm{kep}}},\end{eqnarray} \tag{ 30 }$$

where ρ is the bulk density and ρ_gas is the gas density.

We neglect radial scattering and radial drift of planetesimals.

We model the formation and evolution of systems like the solar system by considering the formation of planetesimals, and the accretion of pebbles and planetesimals onto planetary embryos in an evolving protoplanetary disk driven by a disk wind. We focus on the inner 30 au of the disk since this is where the Sun's planets are located. We are interested primarily on processes that occur while the disk is still present, so simulations are terminated after 3 Myr.

The disk is divided into 200 logarithmically spaced radial zones, spanning 0.1 to 30 au. Gas and pebbles are assumed to flow into this region from the outer edge, being supplied by more distant parts of the disk that we do not model here. The time-integrated gas and pebble fluxes are 0.086 M_⊙ and 290 M_⊕, respectively, over the entire simulation. In each zone, we track the inward flux of pebbles, the conversion of some pebbles into planetesimals, and the accretion of both populations by planetary embryos.

At each step in the evolution, we calculate the number of embryos (possibly <1) per radial zone. Embryos in the same zone are assumed to have identical masses and growth rates rather than following the growth of each embryo separately. We note that orbital migration and gas accretion onto embryos are not included for reasons described in Section 1. The embryos are assumed to maintain a constant spacing in Hill radii due to mergers and scattering, although embryo–embryo interactions are not calculated explicitly.

In this study, we concentrate on six parameters with values that are especially uncertain. These are (i) the strength of disk evolution driven by the disk wind, quantified by α_wind; (ii) the strength of turbulence in the gas, quantified by α_turb; (iii) the optical depth of the disk layer that is heated by the disk wind, measured from the surface, denoted by τ_wind; (iv) the fragmentation velocity of icy pebbles, denoted by v_ice; (v) the conversion efficiency from pebbles to planetesimals, quantified by ; and (vi) the minimum Stokes number of pebbles that can be converted into planetesimals, denoted by ${\mathrm{St}}_{\min }$.

These six parameters, together with the range of values that we explore, are listed in Table 1. The remaining model parameters, which we keep fixed, are listed in Table 2.

Table 1. Variable Model Parameters and Permitted Values

Parameter	Symbol	Value	Best Fit
(1)	(2)	(3)	(4)
Turbulence strength	αturb	10−5–10−3	6.75 ± 0.37 × 10−5
Disk wind strength	αwind	10−3–10−2	1.03 ± 0.04 × 10−3
Optical depth of heated layer	τwind	1–30	22.7 ± 1.8
Planetesimal formation efficiency		10−3–1	0.0407 ± 0.0003
Minimum St for planetesimals to form	${\mathrm{St}}_{\min }$	10−3–0.1	3.53 ± 0.29 × 10−3
Ice fragment speed	vice	20–500 cm s−1	186 ± 12 cm s−1

Download table as:  ASCIITypeset image

The combined mass of the inner planets of the solar system is approximately 2 Earth masses, and their orbits span roughly 0.4–1.5 au. It is unlikely that the inner planets had finished growing at the point when the solar nebula dispersed, but most of the material that would end up in these planets was probably present in this region at that time.

Thus, we require that the total mass, M_terr, of embryos and planetesimals interior to 1.5 au at the end of the simulation is about 2 Earth masses.

Most of the mass contained in the terrestrial planets exists in Venus and Earth, which occupy a relatively narrow region in terms of distance from the Sun. It is plausible that the region beyond Earth's orbit was depleted by processes associated with Jupiter (Walsh et al. 2011; Clement et al. 2018), which we do not attempt to model here. The relative lack of material interior to Venus requires a different explanation. Here, we assume that planetesimal formation and/or mass accretion onto planetary embryos was inefficient in this region, but not entirely prevented.

Specifically, we require that the total mass M_Merc of embryos and planetesimals interior to 0.5 au is 0.025–0.1 Earth masses, which is roughly within a factor of 2 of Mercury's mass.

The terrestrial planets contain very little water compared to the roughly solar composition of the solar nebula. The combined mass of water in Earth's oceans and interior is almost certainly <1% of the planet's mass, while the other inner planets are drier still (Marty 2012; Halliday 2013). While the origin of this depletion remains a matter of debate, here we assume it is the result of materials accreted by the terrestrial planets. We neglect the possibilty that water delivered by pebbles is not accreted by a planetary embryo due to evaporation in the envelope or refluxing back to the nebula (Chambers 2017; Ali-Dib & Thompson 2020; Johansen et al. 2021).

Thus, we require that the water mass fraction, f_ice, of the embryos and planetesimals interior to 1.5 au at the end of the simulation is <1%.

In the core-accretion model, Jupiter's growth begins with the formation of a solid core that grows too massive to sustain a hydrostatic atmosphere, thereby initiating runaway gas accretion. Simulations of core accretion typically find that runaway gas accretion requires a core of roughly 10 Earth masses (Pollack et al. 1996; Ikoma et al. 2000). We do not consider the possibility that ongoing pebble accretion may delay the onset of gas accretion.

Jupiter's core must form before the solar nebula dispersed, and early enough to allow time for gas accretion to occur. Here, we simply require that a 10-Earth-mass core forms at Jupiter's current location, but note that the pebble barrier constraint described below typically requires that this core forms early in the nebula's lifetime.

We require embryos at 5 au to have a mass M_Jup of at least 10 Earth masses at the end of the simulation.

The Nice model for the early evolution of the solar system indicates that the orbits of the giant planets migrated radially after the planets formed due to interactions with a disk of planetesimals at larger distances (Tsiganis et al. 2005). These planetesimals failed to form an additional planet, and the great majority were ultimately ejected from the solar system, with the remnants forming the modern Kuiper Belt. The Nice model suggests that roughly 15–20 Earth masses of planetesimals existed beyond Neptune's original orbit, occupying a region of roughly 20–30 au from the Sun (Nesvorny 2018).

Thus we require that the mass, M_KB, of planetesimals between 20 and 30 au is 15–20 Earth masses at the end of a simulation.

Meteorites can be divided into two distinct groups on the basis of the isotopic ratios of a number of elements (Kleine et al. 2020). This dichotomy spans meteorites from differentiated and undifferentiated parent bodies. The estimated ages for the two groups of parent bodies overlap, which suggests the meteorites sample planetesimals that formed or were modified in two isolated regions of the solar nebula. The origin of the dichotomy is uncertain, but meteorite data imply it was established by roughly 0.5 Myr after the first solids formed (Kleine et al. 2020). Here, we adopt a model in which the dichotomy is established and maintained when an embryo becomes massive enough to block the inward flow of pebbles (Kruijer et al. 2017).

Thus we require that the time, t_iso, needed for an embryo to reach the pebble-isolation mass somewhere in the disk is no more than 0.5 Myr.

Based on the constraints described above, we assign a score, S, to the output of a simulation as follows:

$$\begin{eqnarray}&&S={S}_{1}\times {S}_{2}\times {S}_{3}\times {S}_{4}\times {S}_{5}\times {S}_{6},\end{eqnarray} \tag{ 31 }$$

where S₁ is a score based on the combined mass of the terrestrial planets:

$$\begin{eqnarray}\begin{array}{rcl}{S}_{1} & = & \min \left[{x}_{1},\displaystyle \frac{1}{{x}_{1}}\right]\\ {x}_{1} & = & \displaystyle \frac{{M}_{\mathrm{tot}}(a\,\lt \,1.5\ \mathrm{au})}{2\ {M}_{\oplus }}.\end{array}\end{eqnarray} \tag{ 32 }$$

S₂ is a score based on the mass in the region occupied by Mercury:

$$\begin{eqnarray}\begin{array}{rcl}{S}_{2} & = & \min [1,4{x}_{2}]\qquad {x}_{2}\lt 1\\ & = & \displaystyle \frac{1}{{x}_{2}}\,{x}_{2}\geqslant 1\\ {x}_{2} & = & \displaystyle \frac{{M}_{\mathrm{tot}}(a\,\lt \,0.5\ \mathrm{au})}{0.1\ {M}_{\oplus }}.\end{array}\end{eqnarray} \tag{ 33 }$$

S₃ is a score based on the water mass fraction of the terrestrial planets:

$$\begin{eqnarray}\begin{array}{rcl}{S}_{3} & = & \min \left[1,\displaystyle \frac{1}{{x}_{3}}\right]\\ {x}_{3} & = & \displaystyle \frac{{M}_{\mathrm{tot},\mathrm{ice}}(a\,\lt \,1.5\ \mathrm{au})}{0.01{M}_{\mathrm{tot}}(a\,\lt \,1.5\ \mathrm{au})}.\end{array}\end{eqnarray} \tag{ 34 }$$

S₄ is a score based on the mass of Jupiter's core:

$$\begin{eqnarray}\begin{array}{rcl}{S}_{4} & = & \min [1,{x}_{4}]\\ {x}_{4} & = & \displaystyle \frac{{M}_{\mathrm{emb}}(a=5\ \mathrm{au})}{10\ {M}_{\oplus }}.\end{array}\end{eqnarray} \tag{ 35 }$$

S₅ is a score based on the mass of the primordial Kuiper Belt:

$$\begin{eqnarray}\begin{array}{rcl}{S}_{5} & = & \min \left[1,\displaystyle \frac{4{x}_{5}}{3}\right]\qquad {x}_{5}\lt 1\\ & = & \displaystyle \frac{1}{{x}_{5}}\,{x}_{5}\geqslant 1\\ {x}_{5} & = & \displaystyle \frac{{M}_{\mathrm{tot},\mathrm{pml}}(20\,\lt \,a\lt 30\ \mathrm{au})}{20\ {M}_{\oplus }}.\end{array}\end{eqnarray} \tag{ 36 }$$

And S₆ is a score based on the time t_iso at which the first embryo reaches the pebble-isolation mass:

$$\begin{eqnarray}\begin{array}{rcl}{S}_{6} & = & \min \left[1,\displaystyle \frac{1}{{x}_{6}}\right]\\ {x}_{6} & = & \displaystyle \frac{{t}_{\mathrm{iso}}}{0.5\ \mathrm{Myr}}.\end{array}\end{eqnarray} \tag{ 37 }$$

Using the planet formation model described in Section 2, we search for values of the six model parameters listed in Table 1 that yield systems that satisfy the solar system constraints listed above.

Following Chambers (2018), we use a particle-swarm optimization (PSO) scheme in which a swarm of particles explores the six-dimensional phase space of parameter values. Each particle has an instantaneous position x and velocity v in phase space. We also keep track of the best-fit set of parameters p found by each particle and the global best fit g found by any of the particles.

The PSO scheme undergoes a series of iterations. At each iteration, and for each particle i, we run a planet formation simulation using the set of parameters defined by x _i , and calculate the score for the output of the simulation using Equation (31). If the fit is better than the existing fit for the particle, or for the swarm as a whole, then the stored fit is updated. The particle location is then updated using

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{v}}}_{i} & = & {C}_{1}{{\boldsymbol{v}}}_{i}+{C}_{2}{z}_{1}({{\boldsymbol{p}}}_{i}-{{\boldsymbol{x}}}_{i})+{C}_{3}{z}_{1}({\boldsymbol{g}}-{{\boldsymbol{x}}}_{i})\\ {{\boldsymbol{x}}}_{i} & = & {{\boldsymbol{x}}}_{i}+{{\boldsymbol{v}}}_{i},\end{array}\end{eqnarray} \tag{ 38 }$$

where z₁ and z₂ are uniform random numbers between 0 and 1. If there is no improvement in the fit after 20 iterations, we nudge each particle a small distance in a random direction in phase space and continue.

We use 20 swarm particles. Following Chambers (2018), we choose C₁ = 0.83 and C₂ = C₃ = 1.65. Typically, a few thousand iterations are needed to find an excellent fit (score S > 0.98) to the solar system constraints. We performed 31 optimization runs, starting with different initial random particle positions x . In all but one case, a successful fit was found.

Figure 1 shows the distribution of parameter values for 30 fits found by the PSO scheme. These simulations produced excellent matches to the solar system constraints described in the previous section. In all cases, they scored at least 0.98, where 1 is the maximum possible score. The means and standard deviations for the best-fit parameters from the 30 simulations are given in Table 1.

Zoom In Zoom Out Reset image size

For each of the six parameters, the best-fit cases cluster around a particular value with a range much narrower than the range of values that were explored. However, some spread is observed. Thus, simulations that match the solar system constraints occupy a small but finite region of phase space.

The distributions for five of the parameters lie away from the edges of the permitted region of phase space, which implies that these boundaries did not bias the results. (The values of the optical depth parameter, τ_wind, lie close to the highest permitted value, which is 30, but do not exceed 25.) The disk-wind parameter, α_wind, is an exception, with the best-fit cases clustered against 10⁻³, which is the lowest value considered here.

This suggests there is an additional region of phase space, at smaller values of α_wind, that satisfies the solar system constraints. However, such values may correspond to unrealistically massive disks. When α_wind = 10⁻³, the initial disk mass interior to 30 au is 0.026 M_⊙. Here, we assume that this region is continually supplied by material flowing inwards from more distant parts of the disk. Thus, the majority of the mass should lie outside 30 au initially. The total disk mass is likely to be several times higher than 0.026 M_⊙, which puts the disk at the high end for observed protoplanetary disk masses (Tychoniec et al. 2020). The disk surface density is proportional to 1/α_wind, so values of α_wind < 10⁻³ would correspond to even more massive disks. For this reason, we do not explore smaller values of α_wind.

Some of the parameter values seen in Figure 1 can be readily understood. For example, the temperature structure of the inner disk is determined by τ_wind. If τ_wind is too small, we would expect the water ice fraction of the terrestrial planets to increase above the 1% constraint we imposed in Section 3. More specifically, we require that the ice line remains beyond the outer boundary of the terrestrial-planet region (1.5 au) until a more distant embryo becomes massive enough to stop the influx of pebbles. In the best-fit case, the first embryo reaches this pebble-isolation mass at about 0.5 Myr. Figure 2 shows the evolution of the ice line in this case, and we see that it does indeed cross 1.5 au at roughly 0.5 Myr.

Zoom In Zoom Out Reset image size

The preferred value of the planetesimal formation probability also has a simple explanation, at least to first order. Provided that the pebble Stokes number exceeds ${\mathrm{St}}_{\min }$, the mass of planetesimals that form in the Kuiper Belt depends mainly on . The requirement that 15–20 Earth masses of planetesimals form in the Kuiper Belt constrains to a fairly narrow range of values.

We can get additional insights into the best-fit parameter values by looking at the evolution in detail, as we do below, and by seeing the effect of varying the parameters, which is explored in Section 5.

Figure 3 shows the state of a typical successful simulation at four epochs in time. The green and blue curves in the figure show the pebble and planetesimal surface densities, respectively, as a function of distance. The red curves show the embryo masses at each distance.

Zoom In Zoom Out Reset image size

Early in the simulation, planetesimals and embryos do not form within about 0.5 au of the Sun. In this region, the pebble Stokes number St is below the threshhold ${\mathrm{St}}_{\min }$ for planetesimals to form in the model used here (Lenz et al. 2019). The Stokes number is set by turbulent fragmentation, with St ∝ 1/T_mid. As a result, St increases with radial distance (see Figure 4), so planetesimals and embryos can form beyond a certain distance in the disk where $\mathrm{St}={\mathrm{St}}_{\min }$. The Stokes number also increases over time as the disk cools, so the inner edge of the planetesimal- and embryo-forming region moves inwards. The location and evolution of this inner edge are the main factors that control the amount of solid material in the region now occupied by Mercury.

Zoom In Zoom Out Reset image size

For almost the first 0.5 Myr, pebbles are present throughout the disk, with a surface density decreasing with distance such that Σ_peb ∝ F_peb/(au_peb). The pebble flux, F_peb, declines over time and so does Σ_peb. There is also a small jump at the ice line caused by ice evaporation and the change in the fragmentation velocity from v_rock to v_ice.

The pebble radial velocity, u_peb, has two components: (i) the drift of pebbles with respect to the gas u_drift ∝ a^1/2, and (ii) the advection of the gas ${u}_{\mathrm{gas}}\propto {T}_{\mathrm{mid}}^{1/2}$, where the scalings are appropriate for the model used here. These components are shown in Figure 4. The drift component dominates outside the ice line due to the higher fragmentation strength of ice than rock in the best-fit case, which corresponds to higher St and faster drift. Inside the ice line, however, the radial motion of the gas is typically more important. The difference explains the change in slope of Σ_peb at the ice line seen in the first panel in Figure 3.

As the simulation progresses, the pebble surface density profile flattens somewhat. This is because pebble accretion becomes more efficient as embryos grow larger, so the flux of pebbles reaching the inner disk is reduced due to the loss of pebbles accreted at larger distances.

At about 0.5 Myr, embryos located near 2.5 au reach the pebble-isolation mass. The flux of pebbles to the inner disk is halted at this point, which can be seen in the last two panels of Figure 3. Over time, embryos at progressively larger distances also reach the pebble-isolation mass (see Figure 5), so that an increasingly large fraction of the disk ceases to receive pebbles.

Zoom In Zoom Out Reset image size

The blue curves in Figure 3 show the planetesimal surface density Σ_pml. Planetesimals form rapidly, especially in the inner disk, except for the region where $\mathrm{St}\lt {\mathrm{St}}_{\min }$. By 0.1 Myr, the surface density Σ_pml exceeds the pebble surface density across most of the planet-forming region.

Planetesimals form at a rate $\propto {{\rm{\Sigma }}}_{\mathrm{peb}}/({{aT}}_{\mathrm{mid}}^{1/2})$. The denominator in this expression increases with radial distance, which means the surface density profile for the planetesimals is steeper than that for the pebbles.

At early times, planetesimals form at a faster rate than they are accreted by embryos throughout the disk. Planetesimal accretion rates increase over time since they depend on the cumulative mass of planetesimals that have formed in a region. Conversely, planetesimal formation rates decline as the pebble flux decreases. The two rates become equal at about 0.4 Myr at 1 au (see Figure 6), after which the planetesimal surface density declines. Shortly after this, the flux of pebbles into the inner disk ceases, planetesimal formation stops, and Σ_pml declines more rapidly.

Zoom In Zoom Out Reset image size

In the outer regions of the disk, planetesimal formation is less rapid than in the inner disk, but planetesimal accretion is even more inefficient. Few planetesimals are removed, and the surface density increases steadily over time. By 3 Myr, Σ_pml at 25 au is only a few times smaller than at 1 au, as shown in Figure 6.

The prolonged nature of planetesimal formation in the outer disk explains why a massive proto-Kuiper Belt is able to form despite the fact that planetesimal formation is inefficient at large distances in the model of Lenz et al. (2019). Conversely, planetesimals form rapidly in the inner disk, but formation stops at an early stage due to pebble isolation, and planetesimals are also swept up quite rapidly by embryos in this region.

The red curves in Figure 3 show the embryo mass as a function of distance. Embryo growth rates are fastest just beyond the ice line. Embryos grow more slowly inside the ice line, and more slowly still in the outer parts of the disk. By 3 Myr, embryos between 2.5 and 12 au have reached the pebble-isolation mass and stopped growing apart from a gradual increase due to planetesimal accretion.

The time required to reach the pebble-isolation mass increases with distance from the Sun. The solar nebula also cools over time. As a result, the final embryo mass increases more slowly with distance than the instantaneous pebble-isolation mass (see Figure 5), and is very roughly 10 Earth masses over a range of distances.

Figure 7 shows the growth histories of embryos at four particular locations: 0.4, 1, 5, and 15 au from the Sun. Growth is initially fastest at 1 au, with embryos reaching lunar mass by 0.5 Myr. Growth slows noticeably at this point as pebble accretion ceases.

Zoom In Zoom Out Reset image size

At 5 au, embryos initially grow quite slowly. The growth rate increases rapidly, becoming especially fast between 0.2 and 0.3 Myr. These objects reach an Earth mass by 0.44 Myr and 10 Earth masses by 0.8 Myr. At this point, pebble accretion ceases and growth essentially stalls. At 15 au, embryos follow a similar growth trajectory to those at 5 au, but on a timescale that is an order of magnitude longer.

At 0.4 au, unlike the other locations, embryo formation is delayed until 0.4 Myr. Once formed, embryos grow rapidly due to both pebble and planetesimal accretion. Pebble accretion ceases at 0.5 Myr, but the embryos continue to grow quickly due to planetesimal accretion alone. The reservoir of planetesimals is soon depleted, however, and with no source of new planetesimals the growth rate slows dramatically after 1 Myr.

A key feature of the solar system is the small masses of the terrestrial planets compared to the giant planets. This is reproduced in the growth histories of the embryos at 1 and 5 au in Figure 7, an effect that would be magnified further if gas accretion was included.

We can understand one of the main reasons for this difference by examining the rate of pebble accretion using the expressions in Section 2. For most of the growth history considered here, embryos undergo pebble accretion in the three-dimensional (3D) growth regime. In this case, we can express the growth rate as

$$\begin{eqnarray}&&\displaystyle \frac{{{dM}}_{\mathrm{emb}}}{{dt}}=\displaystyle \frac{{M}_{\mathrm{emb}}}{{t}_{\mathrm{grow}}},\end{eqnarray} \tag{ 39 }$$

where

$$\begin{eqnarray}\begin{array}{rcl}{t}_{\mathrm{grow}} & = & 4.3{\alpha }_{\mathrm{turb}}\left(\displaystyle \frac{{M}_{* }}{{F}_{\mathrm{peb}}}\right)\left(\displaystyle \frac{{c}_{s}^{4}}{{v}_{\mathrm{kep}}^{3}{v}_{\mathrm{frag}}}\right)\left(\displaystyle \frac{{u}_{\mathrm{peb}}}{{u}_{\mathrm{drift}}}\right)\\ & & \times \,\exp \left[2{\left(\displaystyle \frac{\mathrm{St}}{{\mathrm{St}}_{\mathrm{crit}}}\right)}^{0.65}\right]\end{array}\end{eqnarray} \tag{ 40 }$$

and

$$\begin{eqnarray}&&{\mathrm{St}}_{\mathrm{crit}}=\min \left[\displaystyle \frac{4{M}_{\mathrm{emb}}}{{M}_{* }{\eta }^{3}},\displaystyle \frac{2}{3}\right].\end{eqnarray} \tag{ 41 }$$

Figure 8 shows t_grow at 1 au and 5 au as a function of the embryo mass. Here we have assumed that T_mid = 230 and 90 K, respectively, at these locations, which are the values at 0.3 Myr in the best-fit simulations. The pebble fluxes are taken to be 5 × 10⁻⁵ M_⊕ yr⁻¹ and 10⁻⁴ M_⊕ yr⁻¹ at 1 au and 5 au, respectively, where the difference is caused by ice evaporation at the ice line.

Zoom In Zoom Out Reset image size

At 1 au, the growth timescale is roughly 2 × 10⁵ yr, almost independent of embryo mass. At 5 au, t_grow is much longer than this for the initial embryo mass of 10⁻⁵ M_⊕. However, t_grow declines rapidly as M_emb increases. For M_emb > 10⁻³ M_⊕, growth is 4 times faster at 5 au than at 1 au.

The strong dependence of t_grow on M_emb at 5 au is caused by the exponential term in Equation (40), and the fact that the critical Stokes number St_crit for effective pebble accretion depends on the embryo mass. When the simulation begins, St ≫ St_crit at 5 au, and growth is very slow. As the embryos grow, St_crit increases until the exponential term is no longer significant.

At 1 au, the embryos are always massive enough that St < St_crit, and the exponential term in Equation (40) is less important. The notable difference in behavior from that at 5 au is attributable to the strong radial dependence of the 1/η³ factor in Equation (41).

Several factors in Equation (40) favor a shorter growth timescale at 5 au than 1 au when the exponential term can be neglected: (i) the pebble flux is twice as high at 5 au, (ii) the fragmentation speed is also about 2 times higher, which helps explain why the optimization scheme favors v_ice > v_rock, (iii) the quantity u_peb/u_drift is about 3 times smaller (this can be seen in Figure 4), and (iv) the sound speed is also smaller. The combination of all these factors outweighs the dependence on v_kep, which favors more rapid growth at 1 au.

We note that an important part of the reason why embryos grow faster at 5 au than at 1 au is the contribution of inward gas advection to the radial velocity of pebbles at 1 au. Pebbles at 1 au are swept inwards more rapidly by the gas than they would be solely due to drift caused by the pressure gradient. As a result, the pebble surface density at 1 au is substantially reduced, and growth is correspondingly slower.

This effect can be seen in Figure 9, which shows simulations that include (red curves) or ignore (green curves) the contribution of gas advection to the pebble radial velocity. All other model parameters are the same as the case shown in Figure 3. When gas advection is accounted for, embryos inside the ice line are roughly an order of magnitude less massive than they would be otherwise.

Zoom In Zoom Out Reset image size

An additional effect that limits the growth of embryos in the inner disk is the removal of inflowing pebbles before they reach this region. This "pebble filtering factor" is partly the result of pebble accretion by embryos (Guillot et al. 2014). A second, less appreciated aspect is the filtering of pebbles due to the ongoing conversion of pebbles into planetesimals.

We can see the effects of pebble filtering in Figure 10, where the solid curves show the pebble flux as a function of distance at 0.1, 0.3, and 0.6 Myr (blue, green, and red curves, respectively).

Zoom In Zoom Out Reset image size

At 0.1 Myr, embryos remain small throughout the disk, and pebble accretion has little effect on the pebble flux. The flux declines slowly with decreasing distance, which is almost entirely due to planetesimal formation. There is also an abrupt decrease by a factor of 2 at the ice line. The flux becomes constant inside 0.5 au due to the absence of embryos and planetesimal formation.

At 0.3 Myr, the behavior outside 5 au is similar to the earlier epoch. The only difference is that the pebble flux entering the planet-forming region has declined by about 40%. Between 5 and 2 au, however, the situation changes dramatically. The pebble flux declines by a factor of 3, due largely to efficient pebble accretion by embryos with masses comparable to Earth. Note that this pebble filtering factor acts to starve the inner disk of pebbles, and reduces the pebble accretion rates in the terrestrial-planet region beyond the factors considered in the previous section.

The red curve in Figure 10 shows the pebble flux at 0.6 Myr. At this stage, embryos at about 3 au have reached the pebble-isolation mass, truncating the flow of pebbles to regions closer to the Sun. At larger distances, the pebble flux resembles the situation at 0.3 Myr, albeit with a further reduction due to the decreased overall pebble flux.

For comparison, the dotted curves in Figure 10 show a simulation in which there is no planetesimal formation. In this case, the pebble flux is essentially independent of distance apart from the discontinuity at the ice line, and the region just outside the ice line where efficient pebble accretion takes place at 0.3 and 0.6 Myr.

Pebble accretion and planetesimal accretion onto embryos both play a role in planet formation but their relative importance varies with time and location. For example, Figure 11 shows the cumulative mass fraction accreted in the form of pebbles by embryos. The green and red curves show the pebble contributions after 0.5 and 3 Myr, respectively.

Zoom In Zoom Out Reset image size

At both epochs, the disk can be divided into three regions based on the pebble fraction of the embryos. The embryos interior to about 1.5 au have accreted less than half their mass from pebbles, although the pebble contribution is significant, typically at least 10%. Beyond 1.5 au, the pebble contribution rises steeply until pebbles dominate the embryo mass budget. Further out, the pebble contribution falls rapidly, at a distance that increases over time. Thus, planetesimals are the dominant component in the outermost disk, although we note that embryos in this region have undergone little growth overall (see Figure 3).

The rapid increase in pebble fraction at 1.5 au roughly coincides with the location of the ice line at 0.5 Myr. This is also the time when the inward flux of pebbles is first halted by the presence of large embryos further out. Pebble accretion is much more effective beyond the ice line, at this stage, for reasons we described above. Conversely, planetesimal formation tends to become less effective with increasing distance. The transition is smoothed somewhat due to the inward motion of the ice line at earlier epochs. Nonetheless, the transition occurs over a short range of distances.

At 3 Myr, the transition from mostly planetesimal accretion to pebble-dominated accretion occurs at nearly the same location as at 0.5 Myr. However, the planetesimal contribution has increased for all embryos out to about 7 au. This is the result of ongoing planetesimal accretion during the interval, whereas pebble accretion has stopped due to the absence of inward-drifting pebbles. Embryos in the region 12–20 au are still growing at this point, having not yet reached the pebble-isolation mass. Pebbles are the dominant component in this region.

We note that the pebble contribution is likely to decrease if residual planetesimals are swept up after the solar nebula has dispersed. This is likely to be particularly important for the terrestrial-planet region. For example, if we assume that all remaining planetesimals in this region are accreted, then the total pebble contribution falls below about 5% for the region between 1 and 1.5 au.

We can see the role of planetesimal accretion by comparing the red and green curves in Figure 12, which shows two simulations using identical parameters with and without planetesimal accretion, respectively. In the latter case, the embryos can grow only by pebble accretion.

Zoom In Zoom Out Reset image size

As expected, embryo growth is slower when planetesimal accretion is turned off. The difference varies with distance but is apparent throughout the disk. It makes sense that planetesimal accretion should lead to faster embryo growth in the inner disk where planetesimal accretion is quite efficient. However, embryos are also much larger, or grow much faster in the outer disk.

The reason for this behavior becomes clearer by examining the green curve in the final panel of Figure 12. Typically embryos beyond the ice line have either reached the pebble-isolation mass (roughly 2–6 au), or they remain close to their initial mass (>6 au).

As we saw earlier, pebble accretion is largely ineffective when St ≫ St_crit. This is the case for all initial embryos outside the ice line. Initially, therefore, growth by pebble accretion is very slow, and it remains so until the embryos have grown large enough that St_crit ∼ St. (Recall from Equation (41) that St_crit ∝ M_emb.) Once this happens, embryos grow very rapidly by pebble accretion until they reach the pebble-isolation mass.

When planetesimal accretion is included, significant early growth occurs despite the ineffectiveness of pebble accretion. This demonstrates an important role of planetesimal accretion outside the ice line, namely to give an early boost to embryo masses, allowing them to reach the point where rapid pebble accretion can take over. This allows embryos to reach the pebble-isolation mass over a much larger region of the disk by the end of the simulation. Thus, early growth due to the accretion of planetesimals leads to much more effective formation of giant-planet cores by triggering rapid pebble accretion.

Figure 13 shows the effect of varying the turbulence strength, α_turb, on the values of the solar system constraints. All other parameters have the same values as the case shown in Figure 3. Changing the turbulence strength affects the evolution by altering pebble accretion in two ways. First, larger α_turb leads to more energetic pebble collisions, and thus smaller pebbles and smaller St. Second, stronger turbulence increases the pebble scale height, H_peb, reducing the space density of pebbles.

Zoom In Zoom Out Reset image size

Interestingly, the embryo mass at 5 au, M_Jup, is largely unaffected by changes in α_turb for values <4 × 10⁻⁴. These embryos all reach the pebble-isolation mass, M_iso, within 3 Myr and essentially stop growing. M_iso itself depends very weakly on α_turb, so the final embryo mass at 5 au is roughly constant.

However, raising α_turb typically results in slower growth at 5 au, leading to larger values of the pebble-isolation time, t_iso, at least for α_turb > 5 × 10⁻⁵. When α_turb > 5 × 10⁻⁴, embryos take more than 3 Myr to reach the pebble-isolation mass. For these embryos, the growth timescale ${t}_{\mathrm{grow}}\propto {\alpha }_{\mathrm{turb}}^{2}$ in Equation (40) because gas advection dominates the radial motion of pebbles, so that u_peb/u_drift ∝ 1/α_turb. As a result, M_Jup declines steeply with further increases in α_turb.

We note that for α_turb < 5 × 10⁻⁵, the growth rate at 5 au is nearly constant. In this case, the weak turbulence and large values of St lead to very small values of H_peb. Much of the growth of embryos at 5 au takes place in the two-dimensional (2D) growth regime as a result, which is slower than the equivalent 3D growth rate would be for R_cap < H_peb. This effect counteracts the nominally faster growth caused by larger St.

The timing of pebble isolation also determines the ice fraction, f_ice, of the terrestrial planets. When α_turb is small, pebble isolation occurs before the ice line enters the terrestrial-planet region, leading to f_ice = 0. Stronger turbulence delays the onset of pebble isolation so that icy pebbles can enter the terrestrial-planet region, leading to larger f_ice. When α_turb > 2 × 10⁻⁴, solid material in the terrestrial-planet region becomes almost 50% ice by mass.

The total mass of material, M_terr, inside 1.5 au depends on α_turb in a complicated manner. For the smallest values of α_turb, we have $\mathrm{St}\gt {\mathrm{St}}_{\min }$ throughout the disk so that embryos can form everywhere. As turbulence increases, the region where embryos can form shrinks, reducing the extent of the terrestrial-planet region. This is reflected in the decline in the mass, M_Merc, interior to 0.5 au.

In addition, larger α_turb reduces the efficiency of pebble accretion in the inner disk, reducing embryo growth. Growth in the terrestrial region slows down with the onset of pebble isolation because only planetesimal accretion can continue at this point. As we saw above, the timing of pebble isolation, t_iso, is nearly constant for α_turb < 5 × 10⁻⁵. Thus, the net effect of raising α_turb is to reduce the mass of the terrestrial planets with little change for embryos at 5 au.

When α_turb > 10⁻⁴, we see that M_terr decreases rapidly. At this point, the inner edge of the region where embryos can form is approaching the outer edge of the terrestrial-planet region at 1.5 au. The situation changes abruptly when α_turb exceeds 2 × 10⁻⁴. Now, the increase in t_iso means that icy pebbles can enter the region inside 1.5 au. These pebbles have larger St, allowing embryos to form here, albeit quite late in the evolution. Further increases in α_turb delay pebble isolation even more, allowing M_terr to increase despite the reduced efficiency of pebble accretion.

Finally, we note that the mass, M_KB, of planetesimals forming in the proto-Kuiper Belt is essentially independent of α_turb. The planetesimal formation rate depends only on the pebble flux and H_gas provided that u_peb ≃ u_drift. This is true in the outermost disk for all but the largest values of α_turb considered here.

The ice fragmentation speed, v_ice, is one of the two main factors that determine the pebble Stokes number beyond the ice line, the other being the turbulence strength. Thus v_ice affects whether embryos and planetesimals can form (which requires $\mathrm{St}\gt {\mathrm{St}}_{\min }$), and also affects the pebble accretion rate onto embryos. Figure 14 shows the effect of varying v_ice, with the other parameters held fixed.

Zoom In Zoom Out Reset image size

Small values of v_ice lead to very small St, so that embryos cannot form in many parts of the disk. In extreme cases, when v_ice < 25 cm s⁻¹, we find that embryos are unable to form in the Kuiper Belt region (here taken to be >20 au). For v_ice < 40 cm s⁻¹ embryos cannot form at 5 au. This is the reason why M_KB and M_Jup, respectively, fall to zero when v_ice is small.

As the ice fragmentation speed increases beyond about 70 cm s⁻¹, larger pebbles are able to exist, and this leads to more effective pebble accretion. As a result, the pebble-isolation time decreases. Once v_ice exceeds about 230 cm s⁻¹, however, t_iso starts to increase again. At this point, the pebble scale height has become small enough that much of the growth of embryos beyond the ice line takes place in the 2D regime, which is slower than the equivalent 3D growth rate for R_cap < H_peb. From this point on, the reduced pebble surface density due to higher drift rates becomes the more important factor, and growth rates decline with increasing v_ice. There is a broad "sweet spot" for low values of t_iso with v_ice ∼ 200 cm s⁻¹, which explains why the best-fit simulations tend to lie in this part of phase space.

Inside the ice line, the formation of embryos and planetesimals is controlled by v_rock rather than v_ice, so bodies can form in this region even if they are prevented from forming beyond the ice line. However, growth in the terrestrial-planet region is affected by v_ice indirectly because of its influence on the timing of pebble isolation t_iso.

For v_ice > 65 cm s⁻¹, the solid mass M_terr in the terrestrial region is clearly controlled by t_iso. Larger t_iso allows more time for embryos in the terrestrial zone to grow by pebble accretion, leading to a larger total mass in the this region. The mass M_Merc interior to 0.5 au closely tracks the behavior of M_terr since it is controlled by the same factors.

For v_ice < 65 cm s⁻¹, M_terr is essentially constant. In this case, some embryos in the terrestrial-planet region are able to reach the pebble-isolation mass (roughly 2.5 M_⊕ at 1 au) before the end of the simulation. This shuts off the supply of pebbles, halting pebble accretion and preventing any more planetesimals from forming. The remaining planetesimals can still be accreted, but this does not change the total mass in the region.

The fact that embryos inside the ice line are able to reach pebble isolation explains why t_iso becomes roughly constant for small v_ice despite large variations in the behavior beyond the ice line. In this case, the first (and often only) embryos to reach pebble isolation lie inside the ice line.

The planetesimal formation efficiency parameter, , controls the rate at which pebbles are converted into planetesimals provided that $\mathrm{St}\geqslant {\mathrm{St}}_{\min }$. It also determines the fraction of the pebble flux that does not form planetesimals and is thus available for pebble accretion by embryos closer to the Sun. Figure 15 shows the effect of varying while the other parameters remain fixed.

Zoom In Zoom Out Reset image size

When is small, it has little effect on embryo growth outside the ice line. Only a small fraction of pebbles are converted into planetesimals, so planetesimal accretion is negligible in this region and pebble accretion rates are nearly unaffected. Both the embryo mass M_Jup at 5 au and the time t_iso of the onset of pebble isolation are roughly constant as a result.

Larger values of lead to longer times until pebble isolation begins. Although more planetesimals are available for planetesimal accretion this is more than offset by the reduction in the pebble flux so that t_iso increases. When t_iso approaches the simulation length of 3 Myr, the embryos at 5 au also fail to reach the pebble-isolation mass, and M_Jup starts to decline rapidly. This behavior is similar to the steep drop off in M_Jup for large α_turb seen in Figure 13, and for the same reason.

Once t_iso exceeds about 0.6 Myr, icy pebbles are able to enter the terrestrial-planet region before embryos at larger distances reach pebble isolation. This leads to a significant ice fraction of the solid material in the terrestrial-planet region.

For < 0.1, the mass of solid material in the terrestrial region and the region occupied by Mercury both increase with increasing . This occurs for two reasons. First, increased planetesimal formation directly adds to the reservoir of solid material in these regions, even if the planetesimals take time to be accreted by embryos. Second, increasing delays pebble isolation at larger distances, as we noted above. This allows longer for planetesimal formation and pebble accretion to operate in the inner disk.

However, when becomes larger than about 0.15, both these factors become outweighed by the reduced flux of pebbles reaching the inner disk due to planetesimal formation further out. At this point, both M_Merc and M_terr begin to decline rapidly with increasing .

Finally, as expected, the mass M_KB of planetesimals that form in the proto-Kuiper Belt depends linearly on .

In our model, planetesimals and planetary embryos are able to form only when the pebble stokes number St exceeds ${\mathrm{St}}_{\min }$. Figure 16 shows the effect of varying ${\mathrm{St}}_{\min }$ while the other parameters retain their values from the simulation shown in Figure 3.

Zoom In Zoom Out Reset image size

Beyond the ice line, the main effect of ${\mathrm{St}}_{\min }$ is to determine where and when embryos can form. Both M_Jup and t_iso are nearly constant for ${\mathrm{St}}_{\min }\lt 0.03$. In this case, embryos are able to form at the start of the simulation and growth proceeds normally.

For larger values of ${\mathrm{St}}_{\min }$, embryos cannot form initially at 2.5–5 au because the disk is too hot. The resulting turbulent velocities are high enough, and pebble collisions are energetic enough, that $\mathrm{St}\lt {\mathrm{St}}_{\min }$. Later, when the disk cools, St increases and embryos are able to form here. However, the delay means that embryos take longer to reach the pebble-isolation mass, or do not reach this mass at all before the simulation ends. This leads to larger t_iso and smaller M_Jup. When ${\mathrm{St}}_{\min }$ exceeds about 0.07, embryos are unable to form at 5 au at any time.

Inside the ice line, the pattern is similar except that embryo and planetesimal formation is prevented at much smaller values of ${\mathrm{St}}_{\min }$. For example, we require ${\mathrm{St}}_{\min }\lt 0.0035$ in order for planetesimals and embryos to form in the region occupied by Mercury, otherwise planet formation is prevented here entirely.

In principle, planetesimal formation in the proto-Kuiper Belt will also be affected by ${\mathrm{St}}_{\min }$. However, in all cases considered here, we find that $\mathrm{St}\gt {\mathrm{St}}_{\min }$ at the start of the simulation or soon afterwards. As a result, the mass of planetesimals in this region is essentially independent of ${\mathrm{St}}_{\min }$.

Figure 17 shows the effect of varying the disk-wind parameter, α_wind, while the other parameters have the same values as the simulation shown in Figure 3.

Zoom In Zoom Out Reset image size

The changes caused by varying α_wind are more subtle than for some of the model parameters considered above. The main role of α_wind is to determine the inward radial velocity, u_gas, of the gas, which is a component of the inward velocity, u_peb, of the pebbles. This in turn affects the surface density of pebbles, Σ_peb. Larger values of α_wind lead to smaller values of Σ_peb, and thus lower rates of planetesimal formation and pebble accretion.

As we saw in Figure 4, u_peb in the outer disk is dominated by radial drift rather than gas advection. For this reason, α_wind has only a modest effect on M_Jup and M_KB. The effect on the timing, t_iso, of pebble isolation is significant, however, with t_iso increasing by about a factor of 2 for the range of α_wind values considered here. As we noted earlier, when t_iso exceeds about 0.5 Myr, icy pebbles are able to enter the terrestrial-planet region before pebble isolation begins, and this is reflected in the nonzero values of f_ice at larger values of α_wind.

One might expect that the delayed onset of pebble isolation would also lead to an increase in the solid mass, M_terr, in the terrestrial-planet region. However, this effect is more than offset by the decreased efficiency of pebble accretion and planetesimal formation due to the reduction in Σ_peb. As a result, M_terr declines with increasing disk-wind strength.

The disk wind is assumed to heat the disk in a layer with an optical depth of τ_wind from the disk surface. This parameter helps determine the midplane temperature, particularly in the inner part of the disk. It is less important in the outer disk where the temperature is determined mainly by stellar irradiation. Figure 18 shows the effect of varying τ_wind while other parameter values are held fixed.

Zoom In Zoom Out Reset image size

An increase in τ_wind leads to a hotter inner disk, reducing the amount of time that icy pebbles can exist here. As a result, the ice fraction of the inner planets declines for larger values of τ_wind.

When τ_wind is small, a significant amount of the terrestrial-planet region receives icy pebbles before the supply is terminated by pebble isolation. Icy pebbles are stronger than rocky ones for the best-fit value of v_ice. The pebbles are larger as a result, and pebble accretion proceeds more rapidly. More solid material is available due to the presence of ice. Together, these factors lead to faster growth in the inner disk when τ_wind is small, and M_terr is also larger as a result.

Because small values of τ_wind lead to larger values of St, we find that planetesimals and embryos are able to form closer to the Sun in this case. As a result, the mass, M_Merc, of solid material in the region occupied by Mercury increases substantially, exceeding an Earth mass when τ_wind = 1.

A striking feature of the solar system is the large difference in mass between the inner and outer planets. Even when the gaseous envelopes of the outer planets are removed, they are still an order of magnitude more massive than the terrestrial planets.

Morbidelli et al. (2015) used simple growth models to show that planetesimal accretion alone cannot explain the mass difference between the Mars-mass bodies that probably formed in the inner solar nebula and the 10-Earth-mass cores of the giant planets. Instead, this difference could have arisen if these objects formed by pebble accretion provided that the pebble Stokes number dropped by a factor of 10 inside the ice line. They argued that the evaporation of water ice from pebbles drifting inwards across the ice line would have released much smaller rocky grains, leading to a sharp drop in pebble size.

The results in Sections 4 and 5 support this conclusion to some extent. However, we find that additional factors are likely to play a role in determining planetary sizes.

The rocky grains released by ice evaporation will presumably collide and stick together, quickly forming larger and larger aggregates until collisions become too energetic. Thus, the pebble size inside the ice line is controlled by the velocity at which collisions change from growth to fragmentation, just as it is outside the ice line. Unfortunately, the fragmentation velocity of water ice under nebular conditions is not known. Some experiments suggest that icy particles can stick at much higher speeds than rocky ones, but other work suggests this is not the case (Gundlach & Blum 2015; Musiolik & Wurm 2019). Planet formation studies can help constrain v_ice. For example, the best-fit simulations in Section 4 predict an ice fragmentation speed about twice that of rock.

This means that pebbles should be roughly 4 times smaller inside the ice line than outside, in our model, compared to a factor of 10 found by Morbidelli et al. (2015). Despite this difference in pebble size, the best-fit simulations in Section 4 yield essentially the same result as Morbidelli et al. (2015): bodies inside and outside the ice line have masses of roughly 0.1 and 10 Earth masses, respectively. Since pebble accretion typically becomes more effective for larger pebbles, this implies that pebble accretion is inherently less efficient in our model, at least inside the ice line. This can be attributed to the inclusion of the gas advection term in Equation (9) for the pebble radial velocity, which was not considered by Morbidelli et al. (2015). Including gas advection substantially reduces the pebble surface density in the inner disk , lowering the pebble accretion rate in the same way that the presence of very small rocky pebbles would (see Figure 9).

Morbidelli et al. (2015) considered embryos with initial masses M_emb,0, similar to the Moon. However, models and observations of small-body populations suggest that planetesimals were smaller than this when they first formed—a few hundred kilometers at most (Morbidelli et al. 2009; Simon et al. 2016). With this in mind, we used an initial embryo mass of 10⁻⁵ M_⊕, equivalent to a diameter of 350 km, while the other planetesimals are 100 km in diameter.

For embryos this small pebble accretion can be very slow because encounters with an embryo do not last long enough for gas drag to be effective. This is the reason for the exponential term in Equation (40). Since the growth rate increases exponentially with mass, embryos that are somewhat larger than M_emb,0 can grow very fast. This implies embryos will either grow to the pebble-isolation mass or remain near M_emb,0. Intermediate-mass bodies like Mars should be very rare. A similar argument was made by Brasser & Mojzsis (2020).

However, the picture changes substantially when planetesimal accretion is included. Planetesimal accretion leads to steady growth interior to the ice line. Beyond the ice line, planetesimal accretion is slow, but it is faster than pebble accretion for small embryos. Thus, planetesimal accretion can provide the boost needed for rapid pebble accretion to occur later in the evolution. As shown in Figure 12, planetesimal accretion allows multiple Mars-sized objects to form in the terrestrial-planet region, and it expands the region where giant-planet cores can form in the outer disk.

Izidoro et al. (2021a) have also examined the relative contributions of planetesimal and pebble accretion to the formation of the Sun's terrestrial planets. They looked at embryo growth at 0.5 and 1 au in a disk with an imposed pressure bump at 5 au that formed in the first 0.3 Myr of the disk's history. Such a bump could have arisen due to a large embryo or been a long-lived feature of the solar nebula. These authors found that planetesimal accretion dominated growth at 1 au. This was also true at 0.5 au for disks with strong turbulence. However, for weak turbulence, similar to the best-fit case considered here, Izidoro et al. (2021a) found rapid pebble accretion at 0.5 au, forming Earth-mass planets in a few hundred thousand years. This difference with the current work can be attributed to the larger minimum pebble Stokes number needed for planetesimal formation used here, which delays the formation of embryos at 0.5 au.

The planetesimal formation model used here was also used by Lenz et al. (2020) to calculate the initial planetesimal surface density, Σ_pml, in the solar nebula. These authors compared Σ_pml with the mass needed to form the observed planets and small-body populations at seven locations in the solar system. This allowed the authors to constrain uncertain parameters such as the planetesimal formation efficiency , the initial disk mass M_disk,0, and the fragmentation velocity v_frag, which was assumed to be the same for ice and rock.

Lenz et al. (2020) found that a wide range of parameter values produced a good match to the solar system. Like the simulations described here, they also found that planet formation close to the Sun was slowed or prevented due to small values of the pebble Stokes number. The best-fit case found by Lenz et al. (2020) has = 0.05, M_disk,0 = 0.1 solar masses, and v_frag = 200 cm s⁻¹. These are similar to the values found in this study if we take v_ice for the fragmentation velocity. Their preferred value of the turbulence strength α_turb = 3 × 10⁻⁴ is larger than in this study, but the authors noted that a reasonable fit can be obtained for values as small as α_turb = 10⁻⁵. The minimum Stokes number ${\mathrm{St}}_{\min }$ required for planetesimal formation was also similar to the one found here, allowing for the fact that Lenz et al. (2020) used a smooth decrease in formation efficiency rather than an abrupt transition.

These similarities are encouraging. However, we note that the current study provides much stronger constraints on the formation of the solar system. Although Lenz et al. (2020) considered the mass of planetesimals needed to form the planets, they did not allow for the time required to grow the planets by accreting the planetesimals. They also did not consider pebble accretion, which was probably the dominant growth mode in the outer disk. These factors were included in the current study, and this explains why we find that good fits occupy a much narrower region of phase space (see Figure 1).

Recently, it has become apparent that meteorites, and by extension their parent bodies, can be divided into two distinct groups on the basis of their isotopic compositions (Kleine et al. 2020). These are typically referred to as the carbonaceous chondrite (CC) and noncarbonaceous (NC) groups. The two groups of parent bodies have a range of ages that overlap each other. This indicates the parent bodies or their precursors formed in two regions of the solar nebula that coexisted but did not exchange much material. Kruijer et al. (2017) have proposed that the formation of Jupiter's core generated a barrier between the inner and outer solar system when it reached the pebble-isolation mass. This formed two isolated reservoirs, one inside the barrier and one outside, that gave rise to NC and CC meteorite groups, respectively.

Conversely, Brasser & Mojzsis (2020) calculated that Jupiter's growth by pebble accretion was too slow to form an effective barrier at the time required. They estimate that Jupiter's core took roughly 1 Myr to reach the pebble-isolation mass. During this time, enough pebbles drifted past proto-Jupiter to contribute roughly 50% of the total bulk of the terrestrial planets. This appears to be at odds with models that suggest that CC material makes up no more than a few percent of Earth by mass (Marty 2012).

In the simulations in Section 4, pebbles only contribute 10%–20% of the mass of embryos inside the ice line at 3 Myr (see Figure 11). The rest of the embryo mass comes from planetesimal accretion, which was not considered by Brasser & Mojzsis (2020). This fraction falls further if we allow for the sweep up of residual planetesimals after the solar nebula has dispersed. In this case, the pebble contribution becomes <5% for embryos in the region 1–1.5 au.

However, this debate misses an important point. If planetesimals form continually from inward-drifting pebbles, it makes no sense to distinguish between the contributions from pebble and planetesimal accretion. For example, pebbles that were accreted at 0.4 Myr would have had an identical composition to planetesimals that formed at 0.4 Myr. If these planetesimals were also accreted by an embryo, their effect on the embryo's composition would be the same as the pebbles themselves.

A better explanation for the meteorite dichotomy is that the composition of inward-drifting pebbles changed over time. The parent bodies of the NC meteorites, together with the terrestrial planets, sampled pebbles and planetesimals that formed early, before the barrier was established. The parent bodies of the CC meteorites sampled pebbles and planetesimals that formed over a wider range of times, possibly dominated by material that formed after the barrier formed.

Brasser & Mojzsis (2020) have pointed out that the evolution of the ice line plays an important role in determining the masses of the terrestrial planets because it controls whether pebbles in this region are rocky or icy. This in turn determines the fragmentation velocity and Stokes number of the pebbles and the rate of pebble accretion.

We extend this conclusion by noting that the location of the ice line at the time of pebble isolation is the key factor in determining the characteristics of the inner planets. A similar prediction was made by Morbidelli et al. (2016). In the absence of strong viscous heating, the ice line is likely to enter the terrestrial-planet region at an early stage. Thus, the low masses of the terrestrial planets and their nearly ice-free compositions both require that pebble isolation occurred early in the evolution of the solar nebula.

In the previous sections we have examined the early stages of planetary growth in the solar system using a simple model of planet formation. This avoided physical processes that remain highly uncertain, and allowed us the study a wide range of parameter space. However, the simplifications mean the results should be treated with caution. Here, we briefly note some modeling caveats.

The model uses a specified, time-varying flux of fully formed pebbles injected into the outer edge of the model grid at 30 au. A detailed model of pebble growth is not included, and the pebble size at each location is assumed to be the maximum size that pebbles can attain before collisions cause fragmentation rather than sticking.

We chose not to use a detailed dust growth model because existing calculations tend to predict the rapid loss of solids from the outer parts of disks (Brauer et al. 2007), in disagreement with observations of dust-rich mature disks (Ansdell et al. 2016). Instead, we assume that a significant reservoir of solid material survives until late times, but the pebble flux declines over time in common with the gas flux. We note that these two fluxes do not necessarily have a fixed ratio as assumed here.

The assumption that pebbles grow rapidly with most of the mass concentrated near a maximum size is supported by detailed growth models (Brauer et al. 2008). The maximum size can be set by either fragmentation or the radial drift rate (Birnstiel et al. 2012). The latter typically applies in the outermost regions of a disk or when pebble Stokes numbers are very large. Fragmentation-limited growth is more relevant for the cases considered here.

In the model used here, the inward flux of pebbles is halted when a planetary embryo reaches the pebble-isolation mass. In reality, some small particles may be able to cross the partial gap in the disk opened by such an embryo (Weber et al. 2018). Thus, pebble isolation is unlikely to completely shut off the inward flow of solid material as we have assumed.

The orbits of planets can migrate radially due to tidal interactions with disk material. However, planetary migration is sensitive to the thermodynamics and flow of material close to the planet. Various modeling efforts have arrived at different conclusions, so that the rate and direction of migration are uncertain in many cases (Paardekooper et al. 2011; Lega et al. 2014; Benitez-Llambay et al. 2015; Guilera et al. 2021). For this reason we chose to follow some other studies such as Morbidelli et al. (2015) and neglect migration.

Orbital migration is likely to be important in at least some systems. For example, it is a plausible explanation for the short orbital periods of many super-Earths (Izidoro et al. 2021b). It is less clear how important migration was in the solar system. In the terrestrial-planet region planetary embryos may have been too small to undergo much migration while the solar nebula was present. Migration may have been a significant factor for the giant planets, although studies suggest the inward migration of Jupiter was limited by the presence of Saturn (Masset & Snellgrove 2001; Morbidelli & Crida 2007). We note that Jupiter's core could have formed somewhat further out than its current location and still had time to accrete gas within the model used here. A different embryo, located closer to the Sun, at say 2.5 au, could have been responsible for shutting off the pebble supply to the inner disk and causing the compositional dichotomy seen in meteorites, without necessarily growing into a giant planet. However, the fate of such an embryo is unclear and not addressed by this study.

Finally, we note that gas accretion onto planetary embryos is not included in the model. Currently, there are large uncertainties in the gas accretion rate for massive planets due to 3D hydrodynamical effects that are still being studied (Szulagyi et al. 2014; Lambrechts et al. 2019). For this reason, we chose to consider only whether embryos grow large enough to initiate gas accretion rather than estimating their final masses. We focus on Jupiter alone since the presence of Jupiter may have strongly influenced core formation and gas accretion rates for planets that formed later and further from the Sun.