The Effect of a Strong Pressure Bump in the Sun's Natal Disk: Terrestrial Planet Formation via Planetesimal Accretion Rather than Pebble Accretion

The gas disk surface density profile in our simulations is based on the minimum mass solar nebula model (MMSN; Weidenschilling 1977; Hayashi 1981). The central star is a solar-mass star. For simplicity, the gas disk surface density is kept constant over time in most of our simulations (Lenz et al. 2019). This is justified because we are mostly interested in the early evolution of the disk, namely, during its first 1 Myr. For completeness, in order to understand the impact of the gas surface density on our results we also present simulations considering an end-member scenario where the gas disk is initially very low mass. The gas disk extends from about 0.35 au to 350 au in all our simulations, which is consistent with the typical sizes of observed disks (Andrews & Williams 2007; Andrews et al. 2018).

We have performed simulations including and ignoring the effects of one pressure bump in the disk. To mimic the presence of the bump, we modify our original gas disk profile by invoking a simple rescaling function (Pinilla et al. 2012; Dullemond et al. 2018; Morbidelli 2020).

In disks with no pressure bumps, the vertically integrated gas density is represented by the so-called "Minimum MMSN" (Weidenschilling 1977; Hayashi 1981):

$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\mathrm{mmsn}}(r)=1700{\left(\displaystyle \frac{r}{1\mathrm{au}}\right)}^{-1.5}\displaystyle \frac{{\rm{g}}}{{\mathrm{cm}}^{2}},\end{eqnarray} \tag{ 1 }$

where r represents the heliocentric distance. In disks with a pressure bump, the power-law gas surface density profile is modified to mimic the local perturbation of the gas as

$\begin{eqnarray}&&\begin{array}{l}{{\rm{\Sigma }}}_{\mathrm{pb}}(r)={{\rm{\Sigma }}}_{\mathrm{mmsn}}(r)\\ \,\times \,\exp \left[-A\exp \left(-\displaystyle \frac{{\left(r-{r}_{\mathrm{pb}}\right)}^{2}}{w{\left({r}_{\mathrm{pb}}\right)}^{2}}\right)\right]\end{array}\end{eqnarray} \tag{ 2 }$

(Pinilla et al. 2012; Dullemond et al. 2018; Morbidelli 2020), where r_pb, A, and w represent the pressure bump radial location, amplitude, and width, respectively. We will perform simulations with and without pressure bumps, and in the latter cases assuming that it appears at different times during the disk's life.

In all our simulations, the pressure bump is considered to be a static structure, which is certainly very simplistic. However, as we are mostly interested in understanding how a strong pressure bump influences planet formation in the inner disk, we consider this approximation acceptable. We envision that our long-lived and azimuthally extended pressure bump is either caused by planets or disk physics and chemistry. Following Pinilla et al. (2012) and Dullemond et al. (2018) we set

$\begin{eqnarray}&&w({r}_{\mathrm{pb}})={{fH}}_{\mathrm{gas}}({r}_{\mathrm{pb}}),\end{eqnarray} \tag{ 3 }$

where f is a dimensionless parameter of order unity, and ${h}_{\mathrm{gas}}(r)={H}_{\mathrm{gas}}(r)/r=0.037{\left(\tfrac{r}{1\mathrm{au}}\right)}^{2/7}$ is the gas disk aspect ratio (e.g., Chiang & Goldreich 1997; Chambers 2009). H_gas is the disk scale height.

For simplicity, the radial gas disk temperature is also kept constant in our nominal simulations, ⁴ and it is given as

$\begin{eqnarray}&&T(r)={h}_{\mathrm{gas}}^{2}\displaystyle \frac{{{GM}}_{\odot }}{r}\displaystyle \frac{\mu }{{ \mathcal R }},\end{eqnarray} \tag{ 4 }$

where ${ \mathcal R }$ is the ideal gas constant, μ is the gas mean molecular weight set equal to 2.3 g mol⁻¹, G is the gravitational constant, and M_⊙ is the solar mass. In our nominal disk, the disk snow line is at about 5 au. Having the disk snow line at 5 au means that our disk is hotter than classical passive disks (Chiang & Goldreich 1997; Chambers 2009). Our disk temperature is more consistent with models of actively accreting young disks (e.g., Min et al. 2011; Bitsch et al. 2015).

The gas pressure at the disk's midplane is calculated in the locally isothermal limit as

$\begin{eqnarray}&&P(r)={c}_{{\rm{s}}}\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{x}}}}{{H}_{\mathrm{gas}}\sqrt{2\pi }},\end{eqnarray} \tag{ 5 }$

where the sound speed is c_s(r) = H _gas/Ω_k, and ${{\rm{\Omega }}}_{{\rm{k}}}(r)\,=\sqrt{{{GM}}_{\odot }/{r}^{3}}$ represents the Keplerian orbital frequency. Σ_x takes the form of one of our different gas disk profiles (e.g., Equation (1) or (2)).

Long-lived pressure bumps cannot be narrower than about the local disk pressure scale height (Yang & Menou 2010; Ono et al. 2016). A pressure bump narrower than the disk scale height would probably imply that the disk has not reached hydrostatic equilibrium yet. In this case, Rossby wave instability would be potentially triggered and the axial symmetry of the pressure bump would be lost (Dullemond et al. 2018). In light of these arguments, in our simulations, we restrict w to values where 1 ≤ f ≤ 3 (Pinilla et al. 2012). The amplitude of the pressure bump in our simulations is considered to be a free parameter, but we restrict the values of A to be roughly consistent with the density oscillations due to zonal flows driven by magnetohydrodynamical instabilities (Uribe et al. 2011) or gaps opened by giant planets (e.g., Crida et al. 2006).

The radial gas velocity is calculated via the following equation, which corresponds to the viscous evolution of the disk (Lynden-Bell & Pringle 1974):

$\begin{eqnarray}&&{v}_{{\rm{r}},\mathrm{gas}}=-\displaystyle \frac{3}{{{\rm{\Sigma }}}_{{\rm{x}}}\sqrt{r}}\displaystyle \frac{\partial }{\partial r}({{\rm{\Sigma }}}_{{\rm{x}}}\nu \sqrt{r}).\end{eqnarray} \tag{ 6 }$

In Equation (6), the gas disk viscosity is calculated within the α-viscosity paradigm (Shakura & Sunyaev 1973), where

$\begin{eqnarray}&&\nu =\alpha {c}_{s}{h}_{\mathrm{gas}},\end{eqnarray} \tag{ 7 }$

and α represents the Shakura–Sunyaev viscosity parameter. In order to calculate the radial velocity of the gas, we use finite-difference approximation between adjacent points of our radial grid to compute the derivative in Equation (6) (also in Equation (8)).

The azimuthal velocity of the gas is characterized by the pressure support parameter

$\begin{eqnarray}&&\eta =-\displaystyle \frac{1}{2}{\left(\displaystyle \frac{{c}_{{\rm{s}}}}{{v}_{{\rm{k}}}}\right)}^{2}\displaystyle \frac{\partial \mathrm{ln}P}{\partial \mathrm{ln}r},\end{eqnarray} \tag{ 8 }$

where v_k is the Keplerian velocity.

Our nominal gas disk model assumes from the beginning a solar dust-to-gas ratio everywhere in the disk, namely Z = 1.5% (e.g., Asplund et al. 2009). In our simulations, the initial dust surface density is represented by

$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\mathrm{dust}}=0.015{{\rm{\Sigma }}}_{\mathrm{mmsn}}.\end{eqnarray} \tag{ 9 }$

At the beginning of our nominal simulations, the total mass in dust in the inner (<5 au) and outer disk (>5 au) is about 18 and 180 M_Earth, respectively. Note that the solar system metallicity calculated by Asplund et al. (2009) includes several elements heavier than He and many species that do not condense inside the disk snow line (e.g., H₂O, CO, N₂, etc.). Thus, the very initial dust-to-gas ratio inside the snow line could be a factor of 2–3 smaller than that considered here. We explore the impact of the initial dust-to-gas ratios on our results in a limited number of simulations (see Section 3.1.2). These simulations indicate that invoking a lower initial dust-to-gas ratio inside the snow line does not affect the main, broad conclusions of our work.

To follow the radial evolution of the dust we solve the one-dimensional continuity equation for the column dust density (Birnstiel et al. 2010),

$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial {{\rm{\Sigma }}}_{\mathrm{dust}}}{\partial t}+\displaystyle \frac{1}{r}\displaystyle \frac{\partial }{\partial r}\left\{r\left[{\overline{v}}_{{\rm{r}},\mathrm{dust}}{{\rm{\Sigma }}}_{\mathrm{dust}}\right.\right.\\ \quad \left.\left.-{D}_{\mathrm{dust}}\displaystyle \frac{\partial }{\partial r}\left(\displaystyle \frac{{{\rm{\Sigma }}}_{\mathrm{dust}}}{{{\rm{\Sigma }}}_{{\rm{x}}}}\right){{\rm{\Sigma }}}_{{\rm{x}}}\right]\right\}\\ \qquad =\,\displaystyle \frac{\partial {{\rm{\Sigma }}}_{\mathrm{pla}}}{\partial t},\end{array}\end{eqnarray} \tag{ 10 }$

where the dust diffusivity

$\begin{eqnarray}&&{D}_{\mathrm{dust}}=\displaystyle \frac{\alpha {c}_{{\rm{s}}}{H}_{\mathrm{gas}}}{1+{{St}}^{2}},\end{eqnarray} \tag{ 11 }$

and ${\overline{v}}_{{\rm{r}},\mathrm{dust}}$, which is the mass-weighted radial velocity of dust, is calculated following Birnstiel et al. (2012, see their Equation (24)). Σ_pla represents the planetesimal surface density. We test the effects of different gas viscosities in our simulations, namely α = 10⁻⁴ and 10⁻³. The term on the right-hand side of Equation (10) accounts for the formation of planetesimals. We will describe the details of our model for planetesimal formation later in this section (see Section 2.4).

In Equation (10), a priori, we ignore the sublimation of ice dust grains that drift inwards eventually crossing the disk ice line (see in Morbidelli et al. (2016)). We will consider this effect later when modeling the growth of protoplanetary embryos. We numerically solve Equation (10) by invoking the flux-conserving implicit donor-cell algorithm with zero-density boundary conditions (Birnstiel et al. 2012). We use a radial logarithm grid with 600 cells (Desch et al. 2018).

The degree of coupling between the gas and dust motions is well defined by the so-called Stokes number (St). Dust grains with St ≪ 1 are very well coupled to the gas flow streaming lines, whereas dust grains with St ≫ 1 are well detached from the gas motion. In the Epstein regime, the dust grain size is smaller than the mean free path of the gas and the Stokes number takes the following form (Epstein 1924)L

$\begin{eqnarray}&&{St}=\displaystyle \frac{\pi a\rho }{2{{\rm{\Sigma }}}_{{\rm{x}}}},\end{eqnarray} \tag{ 12 }$

where a is the respective dust grain size, and ρ is the particle bulk density.

The radial velocity of dust particles is computed as (e.g., Okuzumi et al. 2012; Ueda et al. 2019)

$\begin{eqnarray}&&\begin{array}{l}{v}_{{\rm{r}},\mathrm{dust}}=-\displaystyle \frac{{St}}{{{St}}^{2}+{\left(1+Z^{\prime} \right)}^{2}}2\eta {v}_{{\rm{k}}}\\ \,+\,\displaystyle \frac{1+Z^{\prime} }{{{St}}^{2}+{\left(1+Z^{\prime} \right)}^{2}}{v}_{{\rm{r}},\mathrm{gas}},\end{array}\end{eqnarray} \tag{ 13 }$

where $Z^{\prime}$ represents the local gas-to-dust ratio.

In our simulations, the initial dust size is set equal to a₀ = 10⁻⁴ cm (Birnstiel et al. 2012; Dra̧żkowska et al. 2019). These dust particles grow on the collisional timescale to the following size:

$\begin{eqnarray}&&{a}_{\mathrm{grow}}(t)={a}_{0}\exp (t/{\tau }_{\mathrm{grow}}),\end{eqnarray} \tag{ 14 }$

where the growth timescale is defined as

$\begin{eqnarray}&&{\tau }_{\mathrm{grow}}=\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{x}}}}{{{\rm{\Sigma }}}_{\mathrm{dust}}{{\rm{\Omega }}}_{{\rm{k}}}}.\end{eqnarray} \tag{ 15 }$

The growth of dust particles is also controlled by drift and fragmentation. Overall, large dust grains tend to collide at higher relative velocities than small ones (Zsom et al. 2010). In the limit where the dust motion is dominated by the gas turbulent motion, one can define the largest size that dust grains can grow until they start to fragment due to high impact velocities. Following Birnstiel et al. (2012), this maximum size (a_frag) reads as

$\begin{eqnarray}&&{a}_{\mathrm{frag}}={f}_{{\rm{f}}}\displaystyle \frac{2{{\rm{\Sigma }}}_{{\rm{x}}}}{\pi \rho }\displaystyle \frac{{v}_{{\rm{f}}}^{2}}{3\alpha {c}_{{\rm{s}}}^{2}},\end{eqnarray} \tag{ 16 }$

where f_f = 0.37 is a fudge parameter, and v_f represents the fragmentation threshold velocity. Based on the results of laboratory experiments (e.g., Testi et al. 2014), we assume that the fragmentation threshold velocities of silicate and ice dust grains are different (Gundlach & Blum 2015). We assume that silicate dust grains inside the snow line have bulk densities equal to 3 g cm⁻³ and fragment at velocities of 1 m s⁻¹. Ice dust grains beyond the snow line—which are expected to be more sticky—have bulk densities equal to 1 g cm⁻³ and fragments at velocities of 10 m s⁻¹. In our nominal experiments, we model a smooth transition in the threshold fragmentation velocities and bulk densities at the snow line by invoking the following functions:

$\begin{eqnarray}&&{v}_{{\rm{f}}}(r)=\displaystyle \frac{10+\exp \left(\tfrac{{r}_{\mathrm{snow}}-r}{0.05{r}_{\mathrm{snow}}}\right)}{1+\exp \left(\tfrac{{r}_{\mathrm{snow}}-r}{0.05{r}_{\mathrm{snow}}}\right)}\,{\rm{m}}\,{{\rm{s}}}^{-1}\end{eqnarray} \tag{ 17 }$

$\begin{eqnarray}&&\rho (r)=\displaystyle \frac{1+3\exp \left(\tfrac{{r}_{\mathrm{snow}}-r}{0.05{r}_{\mathrm{snow}}}\right)}{1+\exp \left(\tfrac{{r}_{\mathrm{snow}}-r}{0.05{r}_{\mathrm{snow}}}\right)}\,{\rm{g}}\,{\mathrm{cm}}^{-3}.\end{eqnarray} \tag{ 18 }$

In all simulations we set r_pb = r_snow, where r_snow represents the location of the fixed disk water snow line in the disk. In our nominal simulations r_snow = 5 au, but we also performed simulations with a hotter disk where the snow line is at ∼10 au.

Dust fragmentation is also promoted by radial drift. The maximum dust grain size due to drift-induced fragmentation is given as (Birnstiel et al. 2012; Dra̧żkowska et al. 2019)

$\begin{eqnarray}&&{a}_{\mathrm{df}}=2{f}_{{\rm{f}}}\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{x}}}}{\rho \pi }\displaystyle \frac{{v}_{{\rm{f}}}}{\left|\eta \right|{v}_{{\rm{k}}}}.\end{eqnarray} \tag{ 19 }$

Although fragmentation is the major limiting effect on dust sizes, this may not be the case everywhere in the disk. Following Dra̧żkowska et al. (2019), we explicitly limit the size that dust particles grow in the disk before they drift. This effect is particularly important in the context of our two-dust-species approach to model correctly dust evolution in the outer part of the disk, where the dust coagulation timescale can be longer than the drift timescale (Dra̧żkowska et al. 2019). The drift-limited dust size is

$\begin{eqnarray}&&{a}_{\mathrm{drift}}=\displaystyle \frac{{f}_{{\rm{d}}}}{\pi \rho }\displaystyle \frac{{{\rm{\Sigma }}}_{\mathrm{dust}}}{2\left|\eta \right|},\end{eqnarray} \tag{ 20 }$

where f_d = 0.55 is a fudge factor used in Birnstiel et al. (2012) to better fit the results of the dual dust-species model to those following the evolution of multiple dust/pebble species.

The largest pebbles in the disk have sizes given by

$\begin{eqnarray}&&{a}_{\max }=\min ({a}_{\mathrm{grow}},{a}_{\mathrm{frag}},{a}_{\mathrm{df}},{a}_{\mathrm{St}=1},{a}_{\mathrm{drift}}),\end{eqnarray} \tag{ 21 }$

where a_St=1 represents the size of particles with St = 1.

Finally, we also impose that dust sizes cannot get smaller than our initial dust size:

$\begin{eqnarray}&&{a}_{\max }=\max ({a}_{\max },{a}_{0}).\end{eqnarray} \tag{ 22 }$

We will use the terms "dust" and "pebbles" to refer to our smallest and largest dust grains thereafter. As pebbles drift inwards due to gas drag they may be spontaneously concentrated into clumps via vortices (e.g., Lyra & Klahr 2011; Raettig et al. 2015) or zonal flows (e.g., Johansen et al. 2012; Dittrich et al. 2013; Bai & Stone 2014). If the local pebble density reaches a critical threshold the clump may collapse to form planetesimals. To mimic the formation of planetesimals in our simulations, we follow the approach of Lenz et al. (2019). As shown in Equation (10), our analytic treatment of planetesimal formation appears as a sink term in the advection equation for dust/pebble evolution, and it is given as

$\begin{eqnarray}&&\displaystyle \frac{\partial {{\rm{\Sigma }}}_{\mathrm{pla}}}{\partial t}=-\displaystyle \frac{\varepsilon }{d}| {v}_{{\rm{r}},\mathrm{peb}}| {{\rm{\Sigma }}}_{\mathrm{peb}}{ \mathcal H },\end{eqnarray} \tag{ 23 }$

where

$\begin{eqnarray}&&\displaystyle \frac{\partial {{\rm{\Sigma }}}_{\mathrm{pla}}}{\partial t}=\displaystyle \frac{\partial {{\rm{\Sigma }}}_{\mathrm{peb}}}{\partial t}.\end{eqnarray} \tag{ 24 }$

v_r,peb and Σ_peb represent the radial velocity and surface density of pebbles in the disk, respectively. We set the pebble trap distance due to zonal flows or vortices as d = 5H_gas (Lenz et al. 2019). ${ \mathcal H }$ is a Heaviside step function that controls planetesimal formation (Gerbig et al. 2019; Lenz et al. 2019), and ε is a dimensionless free parameter defining the planetesimal formation efficiency. In our simulations, we test several values for ε, from 10⁻³ to 0.8 in Gerbig et al. (2019, 2020), Lenz et al. (2019), and Voelkel et al. (2020). In this paper, we do not attempt to constrain ε. ε may be constrained by using our final disks as input of simulations of the late stage of accretion of terrestrial planets and identifying those disks that better produce terrestrial planet analogs. A very massive terrestrial disk (≲2 au) with total mass in solids larger than several Earth masses would probably fail to match the inner solar system terrestrial planets, in particular the low-mass Mars (Chambers 2001; Raymond et al. 2005). Meteorite ages are probably another constraint that can be used to bracket ε. Constraining plausible values of ε will be valuable in any future studies.

Local planetesimal formation is only triggered if the following condition is met:

$\begin{eqnarray}&&{\dot{M}}_{\mathrm{peb},\mathrm{crit}}\leqslant 2\pi r{{\rm{\Sigma }}}_{\mathrm{peb}}| {v}_{{\rm{r}},\mathrm{peb}}| ,\end{eqnarray} \tag{ 25 }$

where the critical pebble flux is

$\begin{eqnarray}&&{\dot{M}}_{\mathrm{peb},\mathrm{crit}}=\displaystyle \frac{{m}_{{\rm{c}}}}{\varepsilon {\tau }_{{\rm{c}}}}.\end{eqnarray} \tag{ 26 }$

The collapsing dust mass m_c is given by

$\begin{eqnarray}&&{m}_{{\rm{c}}}=\displaystyle \frac{4}{3}\pi {l}_{{\rm{c}}}^{3}{\rho }_{\mathrm{Hill}}.\end{eqnarray} \tag{ 27 }$

In Equation (26), τ_c represents the lifetime of traps that can concentrate pebbles to trigger planetesimal formation. We have set τ_c = 1000P_orb (Bai & Stone 2014), where P_orb is the orbital period at the clump-forming location (Lenz et al. 2020). This assumption is also consistent with particle concentration timescales observed, for instance, in numerical simulations of "pure" ⁵ streaming instability (e.g., Simon et al. 2016; Yang et al. 2017). In this work, we assume that planetesimal formation occurs only if St > 5 × 10⁻⁵, which is a condition stricter than that used in Voelkel et al. (2020; St > 0). Nevertheless, it may be considered too generous compared to the typical Stokes number required for planetesimal formation to take place in pure streaming instability (e.g., Yang et al. 2017). We will discuss the impact of this choice later in the paper (in Section 4). The Hill density (Gerbig et al. 2019) is defined as

$\begin{eqnarray}&&{\rho }_{\mathrm{Hill}}=\displaystyle \frac{9}{4\pi }\displaystyle \frac{{M}_{\odot }}{{r}^{3}}.\end{eqnarray} \tag{ 28 }$

By equating the diffusion and collapse timescale of the clump, one derives the critical length scale (Lenz et al. 2019)

$\begin{eqnarray}&&{l}_{{\rm{c}}}=\displaystyle \frac{2}{3}\sqrt{\displaystyle \frac{\alpha }{{St}}}{H}_{\mathrm{gas}}.\end{eqnarray} \tag{ 29 }$

Finally, if Equation (25) is satisfied, the surface density in planetesimals produced in a single time step dt is

$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\mathrm{pla}}=\displaystyle \frac{\varepsilon }{d}| {v}_{{\rm{r}},\mathrm{peb}}| {{\rm{\Sigma }}}_{\mathrm{peb}}{dt}.\end{eqnarray} \tag{ 30 }$

As planetesimals form in our simulations, we also model their growth via collisions with other planetesimals using a simple analytical prescription. We follow the growth of single planetesimals (thereafter we will refer to these bodies as protoplanetary embryos) at two specific locations of the disk: at 0.5 and 1.0 au. Our goal here is not to build an entire system of embryos in the terrestrial region but instead roughly infer the typical masses of embryos in the terrestrial region at the end of the gas disk phase and their growth mode (via planetesimal or pebble accretion). Modeling the growth and evolution of several embryos during the late stage of accretion of terrestrial planets is beyond the scope of this paper.

The orbits of planetesimals in a gaseous disk are affected by several processes such as gas-drag damping, gravitational interaction with other planetesimals, viscous stirring and dynamical friction due to interactions with protoplanetary embryos, and collisional damping (e.g., Ida & Makino 1993; Kokubo & Ida 2000; Thommes et al. 2003; Chambers 2006; Morishima et al. 2013). The interplay between these different processes typically leads to an equilibrium state where the planetesimals' random velocities are in the dispersion-dominated regime (e_pla, i_pla ≳ r_{m H}/(2^1/3 a_Emb)), where

$\begin{eqnarray}&&{r}_{\mathrm{mH}}\approx {a}_{\mathrm{Emb}}{\left(\displaystyle \frac{2{m}_{\mathrm{Emb}}}{3{M}_{\odot }}\right)}^{1/3}.\end{eqnarray} \tag{ 31 }$

r_mH represents the mutual Hill radius of adjacent embryos with identical mass m_Emb and mean semimajor axis a_Emb. In the limit where gravitational stirring and gas-drag damping are in equilibrium, the rms orbital eccentricities/inclinations of nearby planetesimals are defined as

$\begin{eqnarray}&&{e}_{\mathrm{pla}}=2{i}_{\mathrm{pla}}=2.7{\left(\displaystyle \frac{{R}_{\mathrm{pla}}{\rho }_{\mathrm{pla}}}{{\rm{\Delta }}{C}_{{\rm{d}}}{a}_{\mathrm{Emb}}{\rho }_{{\rm{x}}}}\right)}^{1/5}\displaystyle \frac{{r}_{\mathrm{mH}}}{{2}^{1/3}{a}_{\mathrm{Emb}}},\end{eqnarray} \tag{ 32 }$

where R_pla, ρ_pla, Δ, ρ_x, and C_d are the planetesimal radius, planetesimal bulk density, mutual separation of adjacent embryos in mutual hill radii (Kokubo & Ida 2000), the gas column density at the embryo's position, and gas-drag coefficient. We performed simulations with different sizes of planetesimals, R_pla varying from 10 to 500 km. We assumed a single planetesimal size in each simulation and ignored the effects of planetesimal fragmentation in our model (but see for instance Chambers 2006, 2008). The bulk density of planetesimals and embryos inside the snow line are set to 3 g cm⁻³. Following Kokubo & Ida (2000) we set Δ = 10, and C_d = 1 (Chambers 2006).

The growth of a protoplanetary embryo via accretion of planetesimals in the dispersion-dominated regime is well described by the particle-in-a-box approach. The growth rate of an embryo depends on embryo size, local planetesimal surface density, and the relative velocities of the embryo and planetesimals during close approaches. The accretion rate of an embryo is given by (e.g., Safronov 1972; Greenberg et al. 1978; Wetherill 1980; Spaute et al. 1991)

$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{{dM}}_{\mathrm{Emb}}}{{dt}}\approx F\displaystyle \frac{{{\rm{\Sigma }}}_{\mathrm{pla}}^{* }}{\sin ({i}_{\mathrm{pla}})a\sqrt{2}}\pi {R}_{\mathrm{Emb}}^{2}\\ \,\times \,\left(1+\displaystyle \frac{{v}_{\mathrm{esc}}^{2}}{{v}_{\mathrm{rel}}^{2}}\right){v}_{\mathrm{rel}},\end{array}\end{eqnarray} \tag{ 33 }$

where ${v}_{\mathrm{esc}}=\sqrt{\tfrac{2{{Gm}}_{\mathrm{Emb}}}{{R}_{\mathrm{Emb}}}}$ is the escape velocity at the embryo's surface, and v_rel the relative velocity between planetesimals and the planetary embryo. When calculating the planetary embryo's accretion rate, we update ${{\rm{\Sigma }}}_{\mathrm{pla}}^{* }$ at every time step to account for the reduction in the planetesimal surface density due to embryo–planetesimal accretion and also for the increase in surface density due to the formation of new planetesimals following our model. For simplicity, we assume that planetesimals and growing embryos have a negligible effect on the formation of new planetesimals. F is a fudge factor invoked to compensate for the underestimated accretion rate when calculating the relative velocities between planetesimals and embryos using the rms value of the planetesimals' random velocities. The deviation of the planetesimals' velocity relative to a Keplerian circular and planar orbit (embryo's orbit) is defined as (e.g., Kenyon & Luu 1998; Morbidelli et al. 2009)

$\begin{eqnarray}&&{v}_{\mathrm{rel}}={v}_{{\rm{k}}}\sqrt{\displaystyle \frac{5}{8}{e}_{\mathrm{pla}}^{2}+\displaystyle \frac{1}{2}{\left(\sin {i}_{\mathrm{pla}}\right)}^{2}}.\end{eqnarray} \tag{ 34 }$

We have found that F = 3 provides a good fit to the accretion rate of an embryo in N-body numerical simulations (see Figure 1; see also Greenzweig & Lissauer 1992).

Zoom In Zoom Out Reset image size

Pebble accretion onto protoplanetary embryos has been described in detail in several studies (e.g., Ormel & Klahr 2010; Lambrechts & Johansen 2012; Johansen et al. 2015; Levison et al. 2015; Liu & Ormel 2018; Lambrechts et al. 2019). Pebble accretion is well characterized by two regimes of accretion—namely the Bondi and Hill accretion regimes. In the Bondi regime, the relative velocity between a pebble and a protoplanetary embryo on circular and coplanar orbit is

$\begin{eqnarray}&&{\rm{\Delta }}v=\displaystyle \frac{\sqrt{4{{St}}^{2}+1}}{{{St}}^{2}+1}\eta {v}_{{\rm{k}}}.\end{eqnarray} \tag{ 35 }$

In the Hill regime, their relative velocity is higher and takes the form

$\begin{eqnarray}&&{v}_{{\rm{H}}}={{\rm{\Omega }}}_{{\rm{k}}}\displaystyle \frac{{r}_{\mathrm{mH}}}{{2}^{1/3}}.\end{eqnarray} \tag{ 36 }$

A smooth transition from the Bondi to the Hill accretion regimes may be obtained by solving the following combined velocity equation (Johansen et al. 2015) using an interactive method,

$\begin{eqnarray}&&{\hat{R}}_{\mathrm{acc}}=\sqrt{\displaystyle \frac{{t}_{{\rm{f}}}{{Gm}}_{\mathrm{Emb}}}{\delta v}},\end{eqnarray} \tag{ 37 }$

where δ v also depends on ${\hat{R}}_{\mathrm{acc}}$ as

$\begin{eqnarray}&&\delta v({\hat{R}}_{\mathrm{acc}})=0.25{\rm{\Delta }}v+0.25{{\rm{\Omega }}}_{{\rm{k}}}{\hat{R}}_{\mathrm{acc}},\end{eqnarray} \tag{ 38 }$

and

$\begin{eqnarray}&&{t}_{{\rm{f}}}=\displaystyle \frac{{St}}{{{\rm{\Omega }}}_{{\rm{k}}}}.\end{eqnarray} \tag{ 39 }$

We solve Equation (37) using the Newton–Raphson method for ${\hat{R}}_{\mathrm{acc}}$ by setting at first iteration δ v = Δv.

The pebble surface density is calculated in Equation (10) and the pebble surface density in the disk midplane is

$\begin{eqnarray}&&{\rho }_{\mathrm{peb},\mathrm{mid}}=\displaystyle \frac{{{\rm{\Sigma }}}_{\mathrm{peb}}}{{H}_{\mathrm{peb}}\sqrt{2\pi }}.\end{eqnarray} \tag{ 40 }$

The pebble disk scale height is (e.g., Johansen et al. 2015)

$\begin{eqnarray}&&{H}_{\mathrm{peb}}=H\sqrt{\displaystyle \frac{\alpha }{{St}}}.\end{eqnarray} \tag{ 41 }$

Finally, the pebble accretion rate of a protoplanetary embryo is

$\begin{eqnarray}&&\displaystyle \frac{{{dM}}_{\mathrm{Emb}}}{{dt}}=4\pi {R}_{\mathrm{acc}}^{2}{\rho }_{\mathrm{peb},\mathrm{mid}}\delta v({R}_{\mathrm{acc}}),\end{eqnarray} \tag{ 42 }$

where a good fit to the results of numerical integrations is obtained when

$\begin{eqnarray}&&{R}_{\mathrm{acc}}={\hat{R}}_{\mathrm{acc}}\exp (-0.4{\left({t}_{{\rm{f}}}/{t}_{\mathrm{pass}}\right)}^{0.65}).\end{eqnarray} \tag{ 43 }$

For a more detailed description of our pebble accretion prescription, we refer the reader to previous works (Johansen et al. 2015; Bitsch et al. 2019; Izidoro et al. 2021; Lambrechts et al. 2019). In Equation (43),

$\begin{eqnarray}&&{t}_{\mathrm{pass}}=\displaystyle \frac{{{GM}}_{\mathrm{Emb}}}{{\left({\rm{\Delta }}v+{v}_{{\rm{H}}}\right)}^{3}}.\end{eqnarray} \tag{ 44 }$

Finally, if R_acc < R_Emb and t_f > t _pass we set ${{dM}}_{{\rm{Emb}}}/{dt}=0$.

To model the growth of planetary embryos from planetesimal and pebble accretion, we use the output of our dust evolution simulations, viz. the evolution of planetesimal surface density and pebble flux in the inner disk. Planetary embryos in our simulations grow from planetesimal masses. The initial nominal radius of our protoplanetary embryos is 100 km, ⁷ which corresponds to a mass of about 3 × 10⁻⁶ M_⊕ for spherical planetesimal with a uniform density of 3 g cm⁻² . All individual planetesimals in the feeding zone of a growing embryo are assumed to have the initial radius and mass of the embryo and to not grow in time. As mentioned before, we will perform simulations assuming different initial radii for planetary embryos (see Section 3.1.1), ranging from 10 to 500 km to test how the initial planetesimal/embryo size impacts our results.

We do not model the growth of planetary embryos in the entire inner disk. We estimate the final masses and the growth mode of planetary embryos at 0.5 au and 1 au only. This approach is sufficient to provide information about the typical masses and growth mode of embryos in the terrestrial zone. During the embryos' growth, we ignore the effects of type I migration. We will discuss how it may impact our results later in the paper. In this section, we restrict ourselves to simulations with pressure bumps in the disk. Without pressure bumps, planetary embryos in the inner disk are expected to grow to masses beyond 1 M_⊕ (Lambrechts et al. 2019; Voelkel et al. 2020).

Figure 5 shows the growth of planetesimals in disks with α = 10⁻⁴ and a pressure bump (see bump parameters on the left-top panel) with different planetesimal formation efficiencies. The left panels show the embryo's mass evolution as a function of time. The right-side panels show the relative contribution of pebbles (via pebble accretion) to the final embryo's mass. The top panels of Figure 5 correspond to an embryo growing at 0.5 au whereas the bottom panels correspond to an embryo growing at 1 au. It can be seen that the final masses of the embryos may vary up to 3 orders of magnitude, depending on the planetesimal formation efficiency. However, embryos at both locations grow systematically from accreting other planetesimals rather than pebbles. In fact, the contribution of pebble accretion to embryos growing at 1 au is even lower, compared to an embryo at 0.5 au. This result does not depend much on planetesimal formation efficiency but is instead caused by the fast inward drift of pebbles in the inner disk. Embryos in the innermost parts of the disk grow at a faster rate by first accreting other planetesimals—because of shorter dynamical timescales—than an embryo at 1 au. The embryo at 1 au takes longer to reach masses where pebble accretion becomes relatively more efficient and pebbles in the inner disk are lost before it can accrete them. This result is the norm of almost all our simulations. Growing planetary embryos eventually stop growing via planetesimal accretion because as they grow, they increase the random velocities and reduce the density of planetesimals in their feeding zones. Figure 5 shows that the growth timescale of Earth-mass embryos is significantly shorter than Earth's accretion timescale, which is estimated to be between 30 and 150 Myr after solar system formation (Jacobsen 2005; Wood & Halliday 2005; Kleine & Walker 2017). In some scenarios, even Mars-mass embryos grow in timescales shorter than that of Mars', which is estimated to be about 1–10 Myr (Dauphas & Pourmand 2011). This is particularly true in our disks with high planetesimal formation efficiency. Our planetesimal–planetesimal accretion prescription, which is based on the oligarchic growth regime, remains particularly valid for planetary embryos with masses of about Earth mass or lower. Very massive planetesimal disks rapidly produce very massive planetary embryos—in particular in the inner regions (<1 au). Disks with massive planetary embryos transition into the subsequent chaotic growth regime earlier than disks with lower-mass embryos (Kokubo & Ida 2002; Kenyon & Bromley 2006; Walsh & Levison 2019). Therefore, in reality, in the late stages of our disks, planetary embryos in the inner parts of the disk could be even more massive than our final estimates.

Zoom In Zoom Out Reset image size

Figure 6 shows the evolution of embryos in our high-viscosity disk. In this scenario, because pebbles in the inner disk are even smaller than those in the low-viscosity disks, they are even less efficiently accreted by growing embryos. Therefore, in this case, the contribution of pebble accretion to the growth of embryos in the terrestrial region is less than 1%.

Zoom In Zoom Out Reset image size

3.1.1. Dependence on the Initial Planetesimal Size

Streaming instability simulations produce a distribution of planetesimal sizes, with the smallest ones having radii of about 50 km to the largest ones with radii of ∼500 km (e.g., Johansen et al. 2015; Simon et al. 2016). We have performed a set of simulations to test the dependence of our results on planetesimal sizes. Figure 7 shows the growth of an embryo at 0.5 au (top panels) and 1 au (bottom panels) in disks with different initial planetesimal sizes. Our planetesimals' radii vary from 10 to 500 km. We have found that disks with smaller planetesimals promote faster growth (e.g., Levison et al. 2010) because planetesimals accrete more efficiently from the beginning. However, even in this faster growth regime, the contribution of pebble accretion to the final embryo mass is always lower than 50%. This is also the case if we assume that 500 km planetesimals grow in a sea of 10 km planetesimals. Even if the embryos grow faster via accretion of small planetesimals, pebbles in the inner disk are lost before growing embryos can gain significant mass via pebble accretion. Again, the contribution of pebble accretion to embryos growing at 1 au is relatively smaller. We have also performed simulations for a disk with α = 10⁻³, and the results were qualitatively similar.

Zoom In Zoom Out Reset image size

3.1.2. Effects of the Gas Surface Density

In our simulations, for simplicity, we assume that the disk is not evolving with time because we are mostly interested in the early stages of disk evolution. However, in order to understand how the gas surface density may impact our results, we performed an additional set of simulations where we reduce the initial gas surface density by a factor of 10 (${{\rm{\Sigma }}}_{\mathrm{mmsn}}/10$). We keep the assumption that the gas surface density does not evolve in time. In our simulations invoking a low-mass disk scenario, we increase the initial disk dust-to-gas ratio by a factor of 10 in order to have initially the same amount of dust in the disk compared to our nominal disk. This increase in the dust-to-gas ratio may also facilitate planetesimal formation via streaming instability (Simon et al. 2016; Yang et al. 2017). Figure 8 shows that embryos growing both at 0.5 au and 1 au in a low-mass disk also grow mostly via planetesimal accretion. In Figure 8, the planetesimal radius (and initial embryo size) is 100 km.

Zoom In Zoom Out Reset image size

3.1.3. Simulations with Growing Pressure Bumps

So far in our model setup, the pressure bump in the disk was either absent or was assumed to exist since the beginning of the simulation. Here we test the effects of the pressure bump appearing at later times of the disk. In this case, pebbles from the outer disk drift into the inner solar system until the pressure bump forms. We assume that these pebbles also have NC-like isotopic compositions. If these pebbles have CC-like compositions, that would perhaps violate the solar system isotopic dichotomy (but see Schiller et al. 2018, 2020). In these simulations, we reduce the pebble flux by a factor of 2 at the snow line when computing the planetary embryo's growth.

In order to mimic the growth of the pressure bump with time, we use linear interpolation between our power-law gas disk profile (${{\rm{\Sigma }}}_{\mathrm{mmsn}}$) and nominal bumpy disk (Σ_pb). We test three scenarios. In the first one, the bump linearly grows from 0.1 to 0.2 Myr. In the second case, it grows from 0.2 to 0.3 Myr. In both scenarios the bump is completely formed before 0.5 Myr, otherwise most of the pebbles from the outer disk reservoir drift inwards, potentially making it difficult to grow the solar system giant-planet cores later due to the small number of pebbles left. Despite that, we did perform simulations where the bump forms even later, from 0.9 to 1 Myr. This is motivated by the fact that Kruijer et al. (2017) suggest that Jupiter has to form within 1 Myr to explain the meteorite isotopic dichotomy. The results found for this particular set of simulations are qualitatively similar to those presented in this section, in the scenario where the bump grows from 0.2 to 0.3 Myr.

Figure 9 shows the growth of embryos in simulations where the bump grows at different times. Here, we have performed simulations considering two planetesimal formation efficiencies, namely = 10⁻³ and 10⁻² and disks with different viscosities. Figure 9 shows that when α = 10⁻³, planetary embryos grow at most to Moon- to Mars-mass planetary embryos and that these embryos grow mainly via planetesimal accretion. The model result partially changes for low-viscosity disks. In this case, because pebbles are relatively larger, embryos at 0.5 au quickly grow and start to accrete pebbles before the bump forms and the pebble flux ceases. The contribution of pebble accretion for embryos growing in this setup is typically larger than 50%. However, we will show later that these embryos form so quickly that they should have migrated inwards, potentially to the disk inner edge (∼0.1 au). Therefore, they cannot correspond to the terrestrial planetary embryos that formed the solar system. Finally, even in this scenario, where the bump forms later and the inner solar system is invaded by a lot of pebbles from the outer disk, embryos growing at 1 au, again, grow relatively more slowly and mostly via planetesimal accretion. Because these embryos grow more slowly, they may escape large-scale inward migration (⪆1 au; we will address this issue later in Section 3.3).

Zoom In Zoom Out Reset image size

The threshold fragmentation velocity has a tremendous impact on the size of pebbles and their lifetime in the disk. Here we increase the pebbles sizes that can grow in the inner disk by increasing the threshold fragmentation velocity of pebbles inside the snow line from 1 to 3 m s⁻¹. We refer to this quantity as v_r,snow. Note that this increase in the fragmentation velocity inside the snow line is roughly equivalent to assuming that the gas disk midplane has a lower level of turbulence that controls pebble sizes (see Equation (17)). An increase in the threshold fragmentation velocity by a factor of 3 roughly correlates to a level of turbulence a factor of ∼10 lower in the disk midplane (although, in reality, things are a bit more complicated because the pebble disk scale height also depends on the turbulence level). This view is consistent with gas disk models showing a low level of turbulence in the gas disk midplane (Gammie 1996). Disk observations also show disks with very low (e.g., α < 10⁻⁵) levels of turbulence (e.g., Dullemond et al. 2018). Finally, although laboratory experiments suggest that silicate dust grains fragment at velocities of ∼1 m s⁻¹ (Blum 2018), dust coagulation models have considered velocities of up 10 m s⁻¹ (e.g., Birnstiel et al. 2010; Pinilla et al. 2012). So, we test this scenario for completeness. The fragmentation velocity of ice pebbles is the same as our nominal simulations, i.e., 10 m s⁻¹.

Figure 10 shows the results of two simulations where v_r,snow is higher than our nominal case. In the top-left panel we show the growth of embryos in a disk where the pressure bump is fully formed from the beginning of the simulation and the fragmentation velocity inside the snow line is v_r,snow = 3 m s⁻¹. In this case, planetary embryos grow to about Mars mass and most of their masses come from planetesimals, similar to our previous results. This is because pebbles in the inner disk are so large (1–10 cm) that they drift inwards very quickly. This effect is indicated on the right-top panel of Figure 10. First, pebble accretion dominates, but at ∼0.04 Myr, planetesimal accretion starts to dominate the growth of the embryos because pebbles in the inner region have already drifted inside 0.5 au. In the bottom panels we show the growth of embryos in a disk where the pressure bump starts 0.3 Myr after the beginning of our simulation. Here we also take into account the sublimation of inward-drifting ice pebbles crossing the snow line before the bumps form, assuming a 1:1 ice-to-silicate ratio. In this case, the embryo at 0.5 au grows to a mass larger than Earth mostly via pebble accretion. The embryo growing at 1 au grows to a mass larger than 0.5 M_⊕ in less than half a million years. We will show later in the paper that both embryos in this simulation would probably migrate inwards very quickly and become a hot planet/super-Earth (Izidoro et al. 2021; Lambrechts et al. 2019). Therefore, a system where embryos grow very quickly via pebble accretion in our disk model cannot correspond to the current solar system (e.g., Johansen et al. 2021).

Zoom In Zoom Out Reset image size

In order to quantify the effects of type I migration in our previous simulations, we performed several additional simulations to test how embryos with different masses would migrate in our disk. Here, we account for the disk evolution by assuming that gas surface density decays exponentially in an e-fold timescale of 1 Myr. We do not model the growth of embryos but instead we assume that they grow to a given mass and we release them to migrate at different times of the disk. Our goal here is not to precisely constrain the migration history of our embryos but only get a sense of their scale of radial migration and how it correlates with their formation time. The prescription of type I migration considered here is detailed in Izidoro et al. (2017, 2021) and follows Paardekooper et al. (2011). We also impose a planet trap at the disk inner edge, assumed to be at 0.1 au to mimic the effect of the disk inner edge (Romanova et al. 2003; Flock et al. 2019). Note that accounting simultaneously for the effects of type I migration, and pebble and planetesimal accretion in this work would require N-body numerical simulations, which would be computationally very expensive—due to the large number of planetesimals required—and impractical due to the large set of parameters tested in this study.

Figure 11 shows the final location of embryos released to migrate at different times. Each panel shows, from top to bottom, embryos with masses of 0.5 M_⊕, 0.1 M_⊕, and 0.01 M_⊕. Embryos are released to migrate from 0 to 5 Myr. Figure 11 shows that if an embryo of 0.5 M_⊕ forms at 0.5 au in less than 3.5 Myr, it should reach the disk inner edge. That would be the case of embryos growing mostly via pebble accretion in Figure 9 (second panel from top to bottom). Even at 1 au, such a massive embryo should not form earlier than 2 Myr. Mars-mass embryos migrate relatively slowly, so our simulations suggest that an embryo with this mass should not reach Mars mass before 1–2 Myr. Finally, Moon-mass planetary embryos migrate very little and could have formed at any stage of the disk. These results suggest that embryos in the terrestrial region in our disk could not have grown via a vigorous flux of (large) pebbles, otherwise they should have migrated too close to the Sun (Izidoro et al. 2021; Lambrechts et al. 2019). Figure 11 shows that avoiding large-scale inward migration of our embryos (e.g., those produced in Figures 5 and 6) requires a sufficiently low planetesimal formation efficiency (e.g., ≲ 10⁻²)). At the same time, a very small may form low-mass planetary embryos (and potentially a terrestrial disk with not enough mass to form all terrestrial planets). Thus, there seems to exist a region of the parameter space where the interplay between growth and migration is just right to make our current solar system. Identifying this parameter space and how it may change for different disk models remains a subject for future studies.

Zoom In Zoom Out Reset image size

One question is whether it is possible to suppress or slow down type I migration of these massive terrestrial embryos in our simulations? We discuss two different scenarios that we do not consider in this work. The first one is that—because our embryos, for simplicity, are assumed to be fully formed in simulations of this section—we do not include the effects of thermal torques due to accretion of solids (Benítez-Llambay et al. 2015; Masset 2017) when computing the embryos' final locations in the disk. Thermal torques may reverse the migration of a growing protoplanetary embryo (see Masset 2017). We have estimated the critical mass where the effects of thermal torques are not guaranteed to be efficient enough to promote significant outward migration (see Masset 2017). At 1 au, in our model disk, this corresponds to masses larger than about Moon mass (if the adiabatic index is assumed to be 1.4). Baumann & Bitsch (2020) and Guilera et al. (2019) studied the effects of thermal torques on growing protoplanetary embryos and found similar critical masses. The effects of disk winds, which trigger mass loss and angular momentum transports in the disk, have also been invoked to suppress inward type I migration (Ogihara et al. 2015a, 2018; Kimmig et al. 2020). Disk wind effects may induce gas disk surface density profiles with flat or positive slopes in the inner disk (Suzuki et al. 2016). Such profiles may help to halt or even reverse type I migration in the inner regions of the disk (≲1–10 au; Ogihara et al. 2015b). However, the very gas disk profile depends on the assumed disk wind model, with some scenarios not resulting in positive slopes (Bai 2016). Inward type I migration is only fully suppressed in disk wind simulations if the wind is sufficiently strong (positive gas surface density slope) and the corotational torque is assumed fully unsaturated. This latter assumption may not be correct for Earth-mass planets, and some level of inward migration for planets more massive than Mars may be difficult to avoid (Ogihara et al. 2018). Due to the uncertainties of the effects of disk winds on the migration of low-mass planetary embryos, we do not include these effects in our current work.

In this section, we analyze the planetesimal surface density profiles produced in our simulations with that envisioned in the minimum MMSN (Hayashi 1981). Here we perform an additional set of simulations where we assume that the disk surface density evolves with time (exponential decay with e-fold timescale of 1 Myr), which is certainly more realistic than our initial approach. Also, motivated by the results of 3D hydrodynamical simulations that show that (Bitsch et al. 2015) the inner parts of protoplanetary disks are not strongly flared, we impose that our disk aspect ratio is constant over radius as h(r) = 0.037 in another additional set of simulations. At 10 kyr, pebbles at 0.5 au in our nominal flared disk with α = 10⁻⁴ are a factor of ∼7 smaller than those at 4 au (see Figure 2). In our nonflared disk, this difference is smaller, only a factor of 2.5.

The distribution of solids in the MMSN yields a radial surface density profile as $7{\left(r/1\mathrm{au}\right)}^{-1.5}{\rm{g}}\,{\mathrm{cm}}^{-2}$. Figure 12 compares the surface density of planetesimals in our simulations with those in the MMSN disk. The disk parameters and pressure bump setup are shown at the top of the figure. The black solid line shows the final planetesimal surface density in our simulation of Figure 3 (left panels). The gray line shows the planetesimals' surface density in the simulation where the disk is not flared and evolves in time. As one can see, in both disks—flared and nonflared—our final radial distribution of planetesimals is much steeper than that of the MMSN model. Our simulations yield surface density of planetesimals as $\approx {\left(r/1\mathrm{au}\right)}^{-3.5}$ (see also Lenz et al. 2019). Steep surface density profiles may help to alleviate the so-called small-Mars problem of the classical model of the formation of terrestrial planets (Izidoro et al. 2015b). We do not model the late stage of accretion of terrestrial planets in this work, but this result provides an interesting avenue for future research.

Zoom In Zoom Out Reset image size