© The Authors 2025

1 Introduction

Most protoplanetary discs fully dissipate within 3–5 Myr, with hardly any discs surviving past 10 Myr (Mamajek 2009; Ribas et al. 2015; Li & Xiao 2016). Classical theories of planet formation struggle to explain the formation of large planets, especially gas giants such as Jupiter and Saturn, within this short time-frame. For example, classical core formation through runaway accretion of planetesimals (km-sized or larger) takes longer than 10 Myr beyond 5 au from the star because of the low planetesimal number densities at these distances (Goldreich et al. 2004; Levison et al. 2010).

Pebble accretion (PA) is a proposed solution to this problem, whereby planetary cores quickly grow by efficiently accreting mm- to cm-sized solids (Lambrechts & Johansen 2012). The theory proposes that planetary seeds¹ can efficiently accrete aerodynamically small particles called pebbles, due to an interplay between gravitational and dissipative forces (Ormel & Klahr 2010). This leads to high growth rates, even for planets further out in the disc. The theory gained significant traction over the past decade and is supported by the detection of large reservoirs of centimetre-sized particles in protoplanetary discs (Testi et al. 2003; Wilner et al. 2005; Ricci et al. 2010).

In a scenario without gas, pebbles only accrete onto the planet if their trajectory directly collides with the planetary surface. This accretion scenario is called the ballistic regime, characterised by short interaction times and low accretion rates because of the high relative velocities and small collision cross-sections involved (Ormel 2017). The accretion of large particles (e.g. planetesimals) is always ballistic, since planetesimals are too large to be significantly influenced by gas.

In the presence of a protoplanetary gas disc, however, aerodynamically small pebbles lose energy due to drag from the headwind. As a result, the pebbles settle into the gravitational field of nearby planetesimals or planetary embryos and slowly spiral inwards until they accrete onto the planet. In this settling regime, the pebble accretion rate no longer depends on the physical radius of the planetesimal, but on the size of its gravitational field and thus on its mass (Ormel 2017). The accretion cross-section, the region from within which material is accreted onto the planet, is significantly enhanced compared to the gas-free scenario (Ormel & Klahr 2010), potentially becoming as large as the planet’s Hill sphere (Lambrechts & Johansen 2012), resulting in rapid growth. Moreover, unlike planetesimals, pebbles are highly mobile, drifting from the outer disc to the inner disc due to drag (Weidenschilling 1977). As a result, the pebble reservoir is constantly replenished, allowing for further growth.

Within the Solar System, PA has been used in varying degrees of success, to explain the formation of gas giants (Levison et al. 2015; Matsumura et al. 2017, 2021; Raorane et al. 2024; Lau et al. 2024), as well as the masses and orbits of Venus, Earth, and Mars (Johansen et al. 2021), the size distribution of asteroids (Johansen et al. 2015), and the prograde spin preference of the large bodies in the Solar System (Visser et al. 2019; Takaoka et al. 2023; Yzer et al. 2023).

The formation of exoplanetary systems through PA has been studied significantly less. Schoonenberg et al. (2019) studied the formation of the TRAPPIST-1 planetary system through pebble and planetesimal accretion, proposing a specific solution in which the planets formed at the snowline and migrated inwards sequentially. This model explains the mass and the water contents of the TRAPPIST-1 planets. However, Schoonenberg et al. (2019) started with large planetary embryos with radii of 1200 km (~1.8×10⁻³ M_E), to limit the number of particles in the numerical simulation. In fact, most PA studies start with large embryos because of computational constraints (see e.g. Morbidelli et al. 2015; Matsumura et al. 2021; Johansen et al. 2021). Nevertheless, the streaming instability model of planetesimal formation predicts planetesimals start with masses that are one to two orders of magnitude smaller (Simon et al. 2016).

The streaming instability is the consequence of the conservation of momentum in pebble accretion models. As the pebbles are slowed down by the aforementioned headwind and drift inwards from the outer disc to the inner disc, they impart their momentum to the gas and speed it up. This reaction, called the back-reaction, locally speeds up the gas and reduces the head-wind. This, in turn, reduces the rate at which the pebbles drift inwards, allowing the pebbles to locally pile up. As more pebbles pile up, they impart more momentum onto the gas, further reducing the headwind. Because of this positive feedback loop, a small pebble over-density can grow into a massive cloud of pebbles. Once the mass of the pebble cloud exceeds a critical value, the cloud collapses into planetesimals due to self-gravity (Youdin & Goodman 2005; Johansen & Youdin 2007; Youdin & Johansen 2007; Bai & Stone 2010a,b). Planetesimals formed through the streaming instability have radii between 50 and 450 km (Simon et al. 2016). This means that the planetesimals must have already significantly grown before they reach the embryo mass used by Schoonenberg et al. (2019).

In this study, we analyse the formation of terrestrial planets through pebble accretion, starting with around 400 planetesimals with radii between 175 and 450 km. Using numerical N-body simulations that include pebble and planetesimal accretion and type I and II migration (Goldreich & Tremaine 1979; Tanaka et al. 2002; Paardekooper et al. 2010; Paardekooper et al. 2011), we aim to answer the question if PA dominated growth can explain the formation of planetary systems close to the star. Aside from the Sun and TRAPPIST-1 (0.09 M_⊙), we tested the theory of PA for an M-dwarf star (0.20 M_⊙), a star at the edge between M-dwarfs and K-dwarfs (0.49 M_⊙), and a K-dwarf star (0.70 M_⊙) (Habets & Heintze 1981). Altogether, these stars represent the most abundant stellar types in the Milky Way.

The structure of this paper is as follows. Section 2 introduces models describing the disc (Sect. 2.1), pebble accretion (Sect. 2.2), and planetary migration (Sect. 2.3). Section 3 discusses the simulation set-up and parameters. In Section 4 the results of the analytical model with only PA and planet migration are presented. The results of the full N-body simulations are presented in Sect. 5 and further results are given in Sect. 6. Finally, our main conclusions are summarised in Sect. 7.

2 Models

In this study, N-body simulations were performed with around 400 planetesimals with radii between 175 and 450 km and densities of 3 g cm⁻¹. These planetesimals were represented by physical, gravitating particles in the simulation to account for planetesimal accretion, and to track the dynamics and stability of the planetary systems that form. Pebbles were not included as physical particles. Instead, the pebble flux was calculated based on the disc conditions and the mass accretion rate of pebbles on planetesimals and planets was calculated from analytical equations for the accretion efficiency from Ida et al. (2016) and Ormel & Liu (2018).

2.1 Disc model

The disc used in this study was assumed to be a steady accretion disc, meaning that the rate at which its material moves inwards (either gas or pebbles) is independent of the distance to the star. This is generally a good approximation in the inner regions of the disc that we are interested in.

Following Chambers (2009), and Ida et al. (2016), we assumed that the disc consists of two regimes governed by different types of heating. The inner disc is dominated by viscous heating, while stellar irradiation heats the outer disc (Hueso & Guillot 2005; Oka et al. 2011). The mid-plane temperature, T, in these two regimes is given by (Ida et al. 2016): $$T_{\mathrm{v}\mathrm{i}\mathrm{s}\mathrm{c}}\simeq T_{\mathrm{v}\mathrm{i}\mathrm{s}\mathrm{c},0}{M}_{\ast 0}^{3/10}{\alpha }_3^{-1/5}{\dot{M}}_{\ast 8}^{2/5}{\left(\frac{r}{\mathrm{a}\mathrm{u}}\right)}^{-9/10}\mathrm{K},$$(1a)

for the viscous regime, and (see also Chiang & Goldreich 1997): $$T_{\mathrm{i}\mathrm{r}\mathrm{r}}\simeq T_{\mathrm{i}\mathrm{r}\mathrm{r},0}{L}_{\ast 0}^{2/7}{M}_{\ast 0}^{-1/7}{\left(\frac{r}{\mathrm{a}\mathrm{u}}\right)}^{-3/7}\mathrm{K},$$(1b)

for the irradiative regime. In these equations, M_*0 and L_*0 are the mass and luminosity of the central star, normalised by the solar values M_⊙ and L_⊙. Moreover, ${\dot{M}}_{\ast 8}$ and α₃ are renormalisations of the gas accretion rate onto the star, ${\dot{M}}_{\ast }$, and the α-viscosity parameter, α_acc, from Shakura & Sunyaev (1973), respectively, using typical values at the beginning of the accretion phase (see e.g. Ida et al. 2016), such that: $${\dot{M}}_{\ast 8}\equiv \frac{{\dot{M}}_{\ast }}{{10}^{-8}{M}_{\odot }/\mathrm{y}\mathrm{r}}\text{\hspace{0.17em}and\hspace{0.17em}}{\alpha }_3\equiv \frac{{\alpha }_{\mathrm{a}\mathrm{c}\mathrm{c}}}{{10}^{-3}}.$$(2)

Finally, we assumed that T_visc,0 = 200 K and T_irr,0 = 150 K (Chiang & Goldreich 1997; Ida et al. 2016).

The gas scale height, H ≡ c_s/Ω_k, is related to these midplane temperatures through the orbital frequency, ${\mathrm{\Omega }}_{\mathrm{k}}=\sqrt{G{M}_{\ast }/r^3}$, and the sound speed, $c_{\mathrm{s}}=\sqrt{\gamma k_{\mathrm{b}}T/\left(\mu m_{\mathrm{p}}\right)}$, with k_b the Boltzmann constant, and m_p the proton mass. Assuming a heat capacity ratio, γ = 7/5, and a mean molecular weight, μ = 2.34, the scale height is given by: $$H_{\mathrm{g},\mathrm{v}\mathrm{i}\mathrm{s}\mathrm{c}}\simeq C_{\mathrm{H}}{T}_{\mathrm{v}\mathrm{i}\mathrm{s}\mathrm{c},0}^{1/2}{M}_{\ast 0}^{-7/20}{\alpha }_3^{-1/10}{\dot{M}}_{\ast 8}^{1/5}{\left(\frac{r}{\mathrm{a}\mathrm{u}}\right)}^{21/20}\mathrm{a}\mathrm{u},$$(3a) $$H_{\mathrm{g},\mathrm{i}\mathrm{r}\mathrm{r}}\simeq C_{\mathrm{H}}{T}_{\mathrm{i}\mathrm{r}\mathrm{r},0}^{1/2}{L}_{\ast 0}^{1/7}{M}_{\ast 0}^{-4/7}{\left(\frac{r}{\mathrm{a}\mathrm{u}}\right)}^{9/7}\mathrm{a}\mathrm{u},$$(3b)

in which the subscripts visc and irr correspond to the viscous and irradiative regime, respectively, and C_H ≈ 1.9949 × 10⁻³ au K^−1/2. Finally, following the steady accretion assumption, and the relation between ${\dot{M}}_{\ast },H_{\mathrm{g}}$, and Σ_g, $${\dot{M}}_{\ast }=3\pi {\alpha }_{\mathrm{a}\mathrm{c}\mathrm{c}}{H}_{\mathrm{g}}^2{\mathrm{\Sigma }}_{\mathrm{g}}{\mathrm{\Omega }}_{\mathrm{k}},$$(4)

the gas surface density is given by $${\mathrm{\Sigma }}_{\mathrm{g},\mathrm{v}\mathrm{i}\mathrm{s}\mathrm{c}}\simeq C_{\mathrm{\Sigma }}{T}_{\mathrm{v}\mathrm{i}\mathrm{s}\mathrm{c},0}^{-1}{M}_{\ast 0}^{1/5}{\alpha }_3^{-4/5}{\dot{M}}_{\ast 8}^{3/5}{\left(\frac{r}{\mathrm{a}\mathrm{u}}\right)}^{-3/5}\mathrm{g}{\mathrm{c}\mathrm{m}}^{-2},$$(5a) $${\mathrm{\Sigma }}_{\mathrm{g},\mathrm{i}\mathrm{r}\mathrm{r}}\simeq C_{\mathrm{\Sigma }}{T}_{\mathrm{i}\mathrm{r}\mathrm{r},0}^{-1}{L}_{\ast 0}^{-2/7}{M}_{\ast 0}^{9/14}{\alpha }_3^{-1}{\dot{M}}_{\ast 8}{\left(\frac{r}{\mathrm{a}\mathrm{u}}\right)}^{-15/14}\mathrm{g}{\mathrm{c}\mathrm{m}}^{-2},$$(5b)

with C_Σ ≈ 3.7708 × 10⁵ g cm⁻² K.

The boundary between the viscous and radiative regimes lies at the orbital radius, r, where T_visc(r) = T_irr(r), and is given by (see also Chambers 2009; Ida et al. 2016): $$r_{\mathrm{v}\mathrm{i}\mathrm{s}\mathrm{c}-\mathrm{i}\mathrm{r}\mathrm{r}}\simeq {\left(\frac{T_{\mathrm{v}\mathrm{i}\mathrm{s}\mathrm{c},0}}{T_{\mathrm{i}\mathrm{r}\mathrm{r},0}}\right)}^{70/33}{L}_{\ast 0}^{-20/33}{M}_{\ast 0}^{31/33}{\alpha }_3^{-14/33}{\dot{M}}_{\ast 8}^{28/33}\mathrm{a}\mathrm{u}.$$(6)

This boundary shifts radially inwards as the disc evolves, due to the decreasing value of ${\dot{M}}_{\ast }$ as the disc is being depleted.

In fact, all of the time evolution of the disc parameters was modelled by the time dependence of ${\dot{M}}_{\ast (8)}$. In the steady accretion state, this stellar gas accretion rate is given by (Hartmann et al. 1998): $$\begin{array}{rl}{\dot{M}}_{\ast }& =\frac{M_{\mathrm{d},0}}{2t_{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}(2-\xi )}{(\frac{t}{t_{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}}+1)}^{-(5/2-\xi )/(2-\xi )}\\ & \simeq \frac{7}{13}\frac{M_{\mathrm{d},0}}{t_{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}}{(\frac{t}{t_{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}}+1)}^{-20/13},\end{array}$$(7)

where m_d,0 is the initial disc mass, t_diff is the diffusion time, and $\xi =-\frac{\mathrm{d}\ln {\mathrm{\Sigma }}_{\text{g,irr}}}{\mathrm{d}\ln r}=15/14$ is the negative slope of the gas surface density power law in the irradiative regime (see Eq. (5b)). The diffusion time, t_diff, is, in turn, related to the accretion α-viscosity parameter, α_acc, which is assumed to be constant in a steady accretion disc, given by (Hartmann et al. 1998): $${\alpha }_{\mathrm{a}\mathrm{c}\mathrm{c}}=\frac{h_{\mathrm{g},\mathrm{D}}^{-2}}{6\pi (2-\xi {)}^2}\frac{t_{\mathrm{o}\mathrm{r}\mathrm{b},\mathrm{D}}}{t_{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}},$$(8)

where h_g,D ≡ H_g,D/r_D is the aspect ratio of the disc and t_orb,D = 2π/Ω_k,D is the orbital period of a circular orbit, both measured at the outer edge of the disc.

Following Matsumura et al. (2021), we assumed the initial values for the disc mass, M_d,0, the diffusion timescale, t_diff, and the outer radius of the disc, r_D, to fix the disc evolution to a specific model. Especially the latter parameter is important for pebble accretion since the pebble flux rapidly decreases once the pebble formation front reaches the edge of the disc (Sato et al. 2016). The initial conditions are shown in Table 1, and discussed in Sect. 3.

2.2 Pebble accretion model

The growth rate of a planet undergoing pebble accretion can be parameterised as $${\dot{M}}_{\mathrm{p}}=ϵ{\dot{M}}_{\mathrm{F}},$$(9)

where ϵ is the pebble accretion efficiency and ${\dot{M}}_{\mathrm{F}}$ is the pebble mass flux available to the planet. The accretion efficiency strongly depends on the size of the pebbles and the mass of the planet. Meanwhile, the pebble mass flux depends on the location of the pebble formation front, which is the region of the disc in which dust coagulates into pebbles.

2.2.1 Pebble radius and mass flux

This study used semi-analytical expressions for the monodisperse (i.e. single-sized) pebble radius, R_p, and Stokes number, τ_s, based on the work of, among others, Ida et al. (2016). The Stokes number, τ_s, describes how the orbit of a pebble is influenced by the gas of the disc; in particular, by drag from the headwind discussed in the introduction. It depends on both the pebble radius and the local gas conditions.

We assumed that the pebble radius is drift-limited, which is to say that pebbles grow in situ until their drift timescale, t_drift becomes shorter than their growth timescale, t_grow. This is a fair approximation when the location at which most pebbles are being formed, the pebble formation front, is in the outer disc (Ida et al. 2016). The effects of the bouncing and fragmentation barriers were assumed to be negligible, which is valid for the icy grains outside the snowline for discs with α_turb ≲ 10⁻² (Ida et al. 2016). For a discussion about the bouncing or fragmentation barriers for pebbles inside the snowline, we refer to Appendix A.

In this study, the planetesimals were located at orbital radii smaller than about 3 au. Using the assumption of a drift-limited pebble radius in combination with our disc conditions, we find that the Stokes number of pebbles in these inner regions of the disc typically ranges between 0.1 and 2 (see, e.g., Fig. A.1). For the full analytical expressions for the pebble radius, we refer to Ida et al. (2016).

The pebble mass flux, ${\dot{M}}_{\mathrm{F}}$, is determined by the amount of dust that is kicked up by the pebble formation front per unit of time. Following the models of Lambrechts & Johansen (2014); Ida et al. (2016); Sato et al. (2016); Ida et al. (2019), we assumed the pebble mass flux is given by: $$\begin{array}{rl}{\dot{M}}_{\mathrm{F}}=& 3.099\times \frac{\ln {10}^4}{\ln \frac{R_{\mathrm{p}\mathrm{e}\mathrm{b}}}{R_0}}{T}_{\mathrm{i}\mathrm{r}\mathrm{r},0}^{-1}{L}_{\ast 0}^{-2/7}{M}_{\ast 0}^{8/7}{\alpha }_3^{-1}{\dot{M}}_{\ast 8}{\left(\frac{{\mathrm{\Sigma }}_{\mathrm{p}\mathrm{g},0}}{0.01}\right)}^2\\ & \times {\left(\frac{r_{\mathrm{D}}}{\mathrm{a}\mathrm{u}}\right)}^{-4/7}{\left(\frac{t}{t_{\mathrm{p}\mathrm{f}}}\right)}^{-8/21}{(1+\frac{t}{t_{\mathrm{p}\mathrm{f}}})}^{-{\gamma }_{\mathrm{p}\mathrm{f}}}{\mathrm{M}}_{\oplus }/\mathrm{y}\mathrm{r}.\end{array}$$(10)

In this expression, R_peb/R₀ is the typical ratio between pebble and dust particle radius (~10⁴), Σ_pg,0 ≡ Σ_p,0/Σ_g,0 is the initial ratio between the pebble and gas surface density, r_D is the outer radius of the disc, and t_pf is the time it takes the pebble formation front to reach this outer radius. When t > t_pf, the pebble mass flux rapidly decreases due to the factor (1 + t/t_pf)^−γ_pf, introduced by Sato et al. (2016). Here, γ_pf = 1 + γ_pf,2 × (300 au/r_D with γ_pf,2 ~ 0.15 a fit parameter.

Figure 1 shows the pebble mass flux, and cumulative pebble mass as a function of time for the five stellar masses used in this study. There is a clear change in the slope of the pebble mass flux after the pebble formation front reaches the outer edge of the disc, after approximately 10⁵ yr.

Since pebbles contain significant amounts of water ice, the snowline, located where T_disc = 170 K, is an important boundary in analytical pebble prescriptions. Once the pebbles drift inwards past the snowline, their ice sublimates, changing the radius and density of the pebble, and reducing the total flux of pebble mass by 50% (Lodders 2003). The location of the snowline is given by r_snow ~ max(r_snow,visc, r_snow,irr), whereby $$r_{\text{snow},\text{visc\hspace{0.17em}}}\simeq {\left(\frac{T_{\text{snow\hspace{0.17em}}}}{T_{\text{visc},0}}\right)}^{-10/9}{M}_{\ast 0}^{1/3}{\alpha }_3^{-2/9}{\dot{M}}_{\ast 8}^{4/9}\mathrm{a}\mathrm{u},$$(11a)

is the snowline calculated in the viscous regime, and $$r_{\text{snow},\mathrm{i}\mathrm{r}\mathrm{r}}\simeq {\left(\frac{T_{\mathrm{s}\mathrm{n}\mathrm{o}\mathrm{w}}}{T_{\mathrm{i}\mathrm{r}\mathrm{r},0}}\right)}^{-7/3}{L}_{\ast 0}^{2/3}{M}_{\ast 0}^{-1/3}\mathrm{a}\mathrm{u},$$(11b)

the snowline calculated in the irradiative regime. For a detailed discussion of the influence of sublimation on the pebble radius, we refer to Appendix A.

We did not consider shading effects, which can complicate the behaviour of the snowline in the irradiative regime; nor did we take into account the slow outward movement of the snowline due to a decrease in vapour pressure resulting from a diminishing influx of icy pebbles over time (Schoonenberg et al. 2019). We only considered the inwards movement of the snowline when it is located in the viscously heated part of the disc, due to the decrease in ${\dot{M}}_{\ast 8}$ as the disc depletes.

Fig. 1 Radial pebble mass flux (solid lines, left-hand axis) and cumulative pebble mass (dashed lines, right-hand axis) in the inner disc as a function of time for the five stellar masses used in this study. After around 0.1 Myr, the pebble formation front reaches the outer edge of the disc, after which the slope of the pebble mass flux becomes significantly more negative. The disc around the 0.09 M⊙ star is weighted as 10% of the stellar mass, while all others are weighted as 5% (see the discussions in Sects. 3.1 and 4).

2.2.2 Pebble accretion efficiency

In this study, we considered two different prescriptions for the accretion efficiency: ϵ_IGM16 by Ida et al. (2016), and ϵ_OL18 by Ormel & Liu (2018). The model of Ida et al. (2016) is valid in the settling regime for planets on circular orbits, and contains terms for both 2D and 3D accretion. In the settling regime, the relative velocity between the planet and the pebbles is small enough, and therefore the encounter time long enough, for the pebbles to settle in the gravitational field of the planet, and spiral inwards, leading to an accretion cross-section that is orders of magnitude larger than the geometric cross-section of the planet. According to Ida et al. (2016), the accretion efficiency in this regime is given by $${ϵ}_{\mathrm{I}\mathrm{G}\mathrm{M}16}=min(1,\frac{C_{ϵ,\mathrm{i}}{C}_{ϵ}{b}^2}{4\sqrt{2\pi }{h}_{\mathrm{p}}}\frac{\sqrt{1+4{\tau }_{\mathrm{s}}^2}}{{\tau }_{\mathrm{s}}}(1+\frac{3b}{2\chi \eta })).$$(12)

The parameters in this equation are $$\begin{array}{rl}& \chi =\frac{\sqrt{1+4{\tau }_{\mathrm{s}}^2}}{1+{\tau }_{\mathrm{s}}^2},\\ & b\equiv \frac{B}{r}=min(1,\sqrt{3{\tau }_{\mathrm{s}}^{1/3}{r}_{\mathrm{H}}/(\chi \eta )})\times 2\kappa {\tau }_{\mathrm{s}}^{1/3}{r}_{\mathrm{H}},\\ & \kappa =\exp (-{\left(\frac{{\tau }_{\mathrm{s}}}{min(2,{\tau }_{\mathrm{s}}^{\ast })}\right)}^{0.65}),\\ & {\tau }_{\mathrm{s}}^{\ast }=4\left(M_{\mathrm{p}}/M_{\ast }\right)/{\eta }^3\\ & C_{ϵ}=min(\sqrt{\frac{8}{\pi }}\frac{h_{\mathrm{p}}}{b},1),\\ & h_{\mathrm{p}}\simeq {(1+\frac{{\tau }_{\mathrm{s}}}{{\alpha }_{\text{turb\hspace{0.17em}}}})}^{-1/2}{h}_{\mathrm{g}}\simeq {\left(\frac{{\tau }_{\mathrm{s}}}{{\alpha }_{\text{turb\hspace{0.17em}}}}\right)}^{-1/2}{h}_{\mathrm{g}},\\ & C_{ϵ,\mathrm{i}}=\frac{1}{2}\frac{(\mathrm{erf}\left(\frac{z+B}{\sqrt{2}{H}_{\mathrm{p}}}\right)-\mathrm{erf}\left(\frac{z-B}{\sqrt{2}{H}_{\mathrm{p}}}\right))}{\mathrm{erf}\left(\frac{B}{\sqrt{2}{H}_{\mathrm{p}}}\right)}.\end{array}$$

The reduction factor C_ϵ,i was proposed by Matsumura et al. (2021) to account for the effect of the orbital inclination of the planet on the amount of pebbles it encounters. Finally, r_H ≡ R_H/r, where R_H = r(M_p/3M_*)^1/3 is the Hill radius of the planet with mass, M_p.

The second accretion efficiency model, proposed by Ormel & Liu (2018), is based on 2D (Liu & Ormel 2018) and 3D simulations of pebble accretion, and includes the orbit-averaged influence of eccentricity and inclination of the planet’s orbit, as well as disc turbulence. In the 2D limit, the accretion efficiency is given by (Liu & Ormel 2018): $${ϵ}_{\mathrm{O}\mathrm{L}18,2\mathrm{D}}=\frac{A_2}{\eta }\sqrt{\frac{M_{\mathrm{p}}}{M_{\ast }}\frac{\mathrm{\Delta }{v}_{\mathrm{y}}}{{\tau }_{\mathrm{s}}{v}_{\mathrm{K}}}}{f}_{\mathrm{s}\mathrm{e}\mathrm{t}},$$(13)

and in the 3D limit by (Ormel 2017; Liu & Ormel 2018; Ormel & Liu 2018): $${ϵ}_{\mathrm{O}\mathrm{L}18,3\mathrm{D}}=\frac{A_3}{\eta h_{\mathrm{p},\mathrm{e}\mathrm{f}\mathrm{f}}}\left(\frac{M_{\mathrm{p}}}{M_{\ast }}\right)f_{\mathrm{s}\mathrm{e}\mathrm{t}}^2.$$(14)

In these equations, A₂ = 0.32 and A₃ = 0.39 are fitting constants. Furthermore, Δv_y is the azimuthal approach velocity, h_p,eff is the effective pebble aspect ratio, and f_set is the settling fraction (see Ormel & Liu 2018). The effective pebble aspect ratio includes a correction for the inclination, i_p, of the planet’s orbit. The settling fraction, f_set, dependents on both the inclination, and the eccentricity of the planetary orbit, through the vertical and azimuthal approach velocities, Δv_y and Δv_z, respectively.

If these approach velocities become larger than the critical settling velocity, v_*, given by (Ormel & Klahr 2010; Liu & Ormel 2018): $$v_{\ast }={\left(\frac{M_{\mathrm{p}}}{M_{\ast }{\tau }_{\mathrm{s}}}\right)}^{1/3}{v}_{\mathrm{K}},$$(15)

the pebbles have too little time to settle and spiral inwards. Accretion then enters the inefficient ballistic regime, in which only pebbles that are on a direct collision course with the planet are accreted. The expressions for the ballistic regime are provided by Liu & Ormel (2018).

Herein lies the core difference between ϵ_IGM16 and ϵ_OL18. Whereas the IGM16 model assumes that all planets are on circular orbits, the OL18 model considers the fact that for planets on significantly excited orbits, the relative velocities between the planet and the pebbles become too large for settling, leading to a rapid reduction in PA efficiency.

This effect is demonstrated in Fig. 2. This figure shows the orbit-averaged pebble accretion rate for planets with different semi-major axes, and masses, as a function of eccentricity (e), and inclination (i). These values have been calculated for our standard disc model (see Table 1 for initial conditions) around a solar-mass star, 0.01 Myr after the formation of the disc. The results for varying e were calculated with i = 0, and vice versa.

For the smallest planetesimals, an eccentricity between 10⁻³ and 10⁻² is already sufficient to reduce the growth rate by two orders of magnitude in the OL18 model. Meanwhile, in the IGM16 model, only for eccentricities of order 0.1 and above does the orbit-averaged accretion rate change, and by less than a factor of 2, primarily due to the fact that the orbit starts crossing the snowline (1.37 au for these conditions), leading to an increase (for 0.5r_snow < a < r_snow) or decrease (for a > r_snow) in the orbit-averaged encountered pebble flux.

Similarly, an induced inclination leads to a much steeper decline in accretion rate in OL18, than in IGM16, since the former considers both the reduced pebble density at high altitudes, and the increased relative velocity between the planet and the pebbles, while IGM16 only includes the reduction factor C_ϵ,i for the encountered pebble density from Matsumura et al. (2021).

For the N-body simulations, this means that if the planetesimals are quickly excited, PA might come to a halt in OL18, while the planets in IGM16 continue growing, leading to more and more massive planets in the latter simulations.

Fig. 2 Orbit-averaged pebble accretion rate as a function of eccentricity (left panels) and inclination (right panels) for planets around a 1.0 M⊙ star, 0.01 Myr after the formation of the disc. The top plots show the accretion rate for a 10−3 ME planetesimal at different orbital radii, while the bottom plots shows results for planets of different masses at 1 au. The solid lines were calculated with OL18, the dashed lines with IGM16. The OL18 prescription contains explicit expressions for the influence of the eccentricity and inclination on the accretion rate, while the IGM16 model only includes corrections for part of the influence of the inclination. For large e, the change in rp along a single orbit becomes significant enough for the variations in the encountered disc conditions to influence the accretion rate, both in the OL18 and the IGM16 model. The y-axis of the top-left plot changes from a logarithmic to a linear scale at 1 × 10−7 ME/yr to highlight the effect.

2.2.3 Pebble isolation mass

Pebble accretion ceases once the planet’s mass exceeds the pebble isolation mass, M_iso (Lambrechts & Johansen 2014). At this point, the planet significantly perturbs the gas disc, creating a local pressure bump just outside its orbit. Since pebbles drift against the pressure gradient, the local maximum traps the pebbles, preventing them from drifting further inwards. This not only halts pebble accretion for the planet in question, but also for all planets interior to it.

There are multiple prescriptions for M_iso. In this study, we used the prescription by Ataiee et al. (2018), who used 2D hydrodynamical simulations with gas and dust with turbulence parameters in the range of α_turb = [5 × 10⁻⁴, 1 × 10⁻²] to determine the dependence of M_iso on the disc aspect ratio, pressure gradient, and (turbulent) viscosity. They proposed $$\begin{array}{rl}{M}_{\mathrm{i}\mathrm{s}\mathrm{o}}\simeq & h_{\mathrm{g}}^3\sqrt{37.3{\alpha }_{\mathrm{t}\mathrm{u}\mathrm{r}\mathrm{b}}+0.01}\\ & \times (1+0.2{\left(\frac{\sqrt{{\alpha }_{\mathrm{t}\mathrm{u}\mathrm{r}\mathrm{b}}}}{h_{\mathrm{g}}}\sqrt{\frac{1}{{\tau }_{\mathrm{s}}^2}+4}\right)}^{0.7})M_{\ast }.\end{array}$$(16)

In the 2D low-viscosity limit (α_turb ~ 10⁻⁴) used in our study, the prescription of Ataiee et al. (2018) agrees well with the results of the main competing model by Bitsch et al. (2018) (see Fig. 8 of Ataiee et al. 2018). In the high-viscosity limit (α_turb ~ 10⁻²), the two models do significantly differ, by a factor of 3, but this range is not relevant to this study.

2.3 Planetary migration model

In this study, the equations of motion of planetesimals were solved numerically in order to analyse the dynamical evolution of the system. Aside from the gravitational force terms from the central star and the other planetesimals in the system, migration through planet-disc interactions enters into a planet’s equation of motion.

We used the equations of motion for migrating planets proposed by Ida et al. (2020). This prescription predicts planet-disc interactions well, both in the subsonic (Tanaka & Ward 2004) and the supersonic (Muto et al. 2011) regimes. The migration component of the equation of motion is given by $${\left(\frac{\mathrm{d}\mathbf{v}}{\mathrm{d}\mathrm{t}}\right)}_{\text{migr\hspace{0.17em}}}=-\frac{v_{\mathrm{K}}}{2{\tau }_{\mathrm{a}}}{\mathbf{e}}_{\theta }-\frac{v_{\mathrm{r}}}{{\tau }_{\mathrm{e}}}{\mathbf{e}}_{\mathrm{r}}-\frac{v_{\theta }-v_{\mathrm{K}}}{{\tau }_{\mathrm{e}}}-\frac{v_{\mathrm{z}}}{{\tau }_{\mathrm{i}}}{\mathbf{e}}_{\mathrm{z}},$$(17)

where v_K is the Keplerian orbital speed at radius r; v_r, v_θ and v_z are the velocity components in the radial, azimuthal and vertical direction, respectively; e_r, e_θ and e_z are the corresponding unit vectors; and τ_a, τ_e and τ_i are the characteristic timescales for the evolution of the semi-major axis, eccentricity and inclination, respectively. The timescales for type I migration are discussed in detail by Ida et al. (2020).

As planets grow more massive and transition from type I to type II migration, their inwards motion slows down. According to Kanagawa et al. (2018), type II migration is the same as type I migration, but with a reduced surface density due to the planet gap that has formed. They showed that the timescale of semi-major axis decay can be written as $${\tau }_{\mathrm{a}}^{\mathrm{\prime }}\simeq (1+0.04K){\tau }_{\mathrm{a}},$$(18)

with $$K={\left(\frac{M_{\mathrm{p}}}{M_{\ast }}\right)}^2{\left(\frac{H_{\mathrm{g}}}{a}\right)}^{-5}{\alpha }_{\mathrm{t}\mathrm{u}\mathrm{r}\mathrm{b}}^{-1},$$(19)

which is valid for both type I and type II migration. However, there is no consensus on the evolution of the eccentricity and inclination during type II migration. Following Matsumura et al. (2021), we assumed that τ_e and τ_i evolve similarly to τ_a, such that $${\tau }_{\mathrm{e}}^{\mathrm{\prime }}\simeq (1+0.04K){\tau }_{\mathrm{e}},\text{\hspace{0.17em}and\hspace{0.17em}}{\tau }_{\mathrm{i}}^{\mathrm{\prime }}\simeq (1+0.04K){\tau }_{\mathrm{i}}.$$(20)

These adjusted timescales were used in the equation of motion of Eq. (17). The back reactions of the planets on the gas were ignored.

3 Simulation parameters and initial conditions

In this study, we used a modified version of the N-body simulator SyMBA (Duncan et al. 1998) to analyse the planetary systems that form around different low-mass stars. This version of SyMBA was parallelised by Lau & Lee (2023), and has been augmented to include pebble accretion, type I and II migration, and eccentricity and inclination damping (Matsumura et al. 2017, 2021). The pebbles were not included as actual particles in the N-body simulations. Instead, we analysed two different pebble accretion models which describe the mass accretion rate of the planets, given the disc conditions and the Stokes number of the pebbles: the model of Ida et al. (2016) (IGM16), and the one of Ormel & Liu (2018) (OL18). We performed simulations for five different stellar masses: 0.09 M_⊙ (M-dwarf, TRAPPIST-1), 0.20 M_⊙ (M-dwarf), 0.49 M_⊙ (M-dwarf/K-dwarf), 0.70 M_⊙ (K-dwarf), and 1.00 M_⊙ (G-dwarf) (Habets & Heintze 1981). These stellar types make up the bulk of the main sequence stars in the galaxy. Moreover, because these stars are relatively long-lived, and have a habitable zone at small orbital radii where terrestrial planets are expected to form, planets around these stars are prime candidates in the search for life.

3.1 Disc parameters

The initial mass of a protoplanetary disc typically ranges between 1% and 10% of the mass of the central star (Pascucci et al. 2016). For the 0.09 M_⊙ stars, we focussed on discs at the high end of this range. These discs allow for the formation of larger planets due to the increased solid mass, which compensates for the fact that in our simulations of these systems, we include fewer planetesimals than the simulations with higher mass stars, because of computational constraints following from the required step size in the 0.09 M_⊙ simulations, which is discussed below. For all other stars, the initial disc mass was assumed to be 5% of the stellar mass. Furthermore, we assumed r_D = 100 au, based on the observed typical size of protoplanetary discs (Andrews et al. 2018; Andrews & Williams 2007; Vicente & Alves 2005), and we took t_diff = 0.5 Myr, which fits our assumption that the disc has a lifetime of about 5 Myr. Finally, we assumed α_turb = 10⁻⁴.

The other disc parameters follow from the equations presented in Sect. 2.1 and by Ida et al. (2016). In these equations, we assumed the mass-luminosity relation for low-mass protostars to be approximately L_*/L_⊙ ≃ M_*/M_⊙, following the fits with power laws of Baraffe et al. (2015), which are close to linear.

Two important control parameters that are not directly derived from the equations are the inner radius of the disc, r_in, and the inner truncation radius, r_trunc. The gas of the disc does not extend all the way to the surface of the star, most likely due to interactions with the star’s magnetosphere. The star’s magnetic field disrupts the disc and creates a cavity out to the radius where the magnetic energy density is equal to the kinetic energy density of the gas (Long et al. 2005). For a typical T Tauri star, the magnetospheric boundary lies between 0.05 and 0.1 au (Romanova & Lovelace 2006; Romanova et al. 2019). Interior to this inner radius, the surface density is assumed to rapidly decrease, such that (Brasser et al. 2018) $${\mathrm{\Sigma }}_{\mathrm{g}}(r<r_{\text{in\hspace{0.17em}}})={\mathrm{\Sigma }}_{\mathrm{g}}\left(r_{\text{in\hspace{0.17em}}}\right)\times \tanh \left(\frac{r-0.95r_{\text{in\hspace{0.17em}}}}{H_{\mathrm{g}}(r)}\right).$$(21)

The resulting gas cavity stops type I migration, creating a trap for migrating planets between 0.95–1.00 r_in. This is because the decreased surface density prevents planets from forming the density waves responsible for the torques that cause planetary migration. Without an inner edge of the disc, or any other (artificial) pressure bump stopping migration, all planets with a mass comparable to Earth or higher, drift into the central star within a few hundred thousand years of their formation, significantly limiting the odds of finding stars with (exo)planets, which is not in line with observations.

Particles that drifted even further inwards and crossed the inner truncation radius, r_trunc, were assumed to have accreted onto the star and were removed from the simulation. This r_trunc is not necessarily a physical parameter, but a computational constraint. Since SyMBA uses a single, fixed time step size, Δt, for all particles, the total number of steps in the simulation, and therefore the total runtime, is determined by the smallest necessary time step. For symplectic integrators like SyMBA, this time step should typically be around 1/20 th of the orbital period at the truncation radius (Wisdom & Holman 1991).

However, since the formation and evolution region of terrestrial planets lies so close to the star, r_trunc must be small, compared to simulations focussing on gas giant formation. As a result, these simulations take many months to complete. Given our limited computing time at the Snellius supercomputer (10⁶ CPU hours), we decided to use slightly longer time steps. For all stellar masses except 0.09 M_⊙, we set r_trunc to match an orbit with a period of 0.01 year, with the time step equalling 1/15th of this value. To allow for planets to be pushed to stable orbits interior to r_in by other massive planets, without immediately being removed from the simulation, we set r_in to an orbit with a period of 0.03 years for most stars, which corresponds to about 0.1 au for a solar-mass star.

These values are not appropriate for the 0.09 M_⊙ star, however. The most famous 0.09 M_⊙-system, TRAPPIST-1, has five planets with periods shorter than 0.03 years, three of which lie within the habitable zone. Two planets have periods shorter than 0.01 year. For this star, we therefore took r_in at the radius with an orbital period of 0.01 yr, so that all planets in the habitable zone were exterior to it, and r_trunc at an orbital period of 0.005 yr. Given our computational constraints, we had to employ a step size of 1/10 th of the orbital period at r_trunc, and reduce the total number of particles in the simulation.

3.2 Initial planetesimal distribution

As discussed in Sect. 1, a probable source for planetesimals in protoplanetary discs with significant pebble reservoirs is the streaming instability. In short, a positive feedback loop in the back-reaction of the pebbles on the gas can locally eliminate the headwind, as a result of which pebbles stop drifting inwards and start piling up. Once these dense pebble filaments exceed the threshold of ρ_peb ~ ρ_g, they collapse under their own gravity and form planetesimals with typical radii between about 100 and 400 km (Simon et al. 2016).

In this study, we assumed that planetesimals have an initial radius between 175 and 450 km. Since the number distribution of planetesimals is given by a power law, including planetesimals smaller than 175 km rapidly increases the number of particles in the simulation, and therefore the computational load, without significantly influencing the total mass of the planetesimal ring, which, for a shallow slope, is dominated by the more massive particles. To maintain a reasonable planetesimal ring mass, while keeping the number of particles low, we therefore ignored the smallest planetesimals from Simon et al. (2016) and focussed on slightly more massive ones. Nevertheless, these planetesimals are still significantly smaller than those used by, for instance, Schoonenberg et al. (2019).

For the size distribution of the planetesimals, we sampled a truncated Pareto distribution, such that $$R_{\mathrm{p}\mathrm{l}}=R_{\mathrm{p}\mathrm{l},min}{\zeta }^{-1/\beta },$$(22)

where R_min = 175 km, ζ is a randomly generated number from a uniform distribution between (0, 1) and β is the slope of the Pareto distribution. We assumed β = 2.5, which follows from the collision equilibrium (Dohnanyi 1969). Only planets with R_pl ≤ 450 km were accepted into the simulation. The initial mass of the planetesimals was estimated by assuming they are homogeneous spheres with a bulk density of 3000 kg m⁻³, and ranges between 10⁻⁵ and 2 × 10⁻⁴ M_E.

For the initial semi-major axis, we tested two different models. In the first model, the planetesimals were uniformly spread over a wide planetesimal ring. These simulations studied a scenario in which pebble accretion dominates over planetesimal accretion because of the relatively low probability of encountering other planetesimals.

The inner radius of the planetesimal ring in these simulations was set at 0.2 au, so that in every simulation, the most massive planets were able to migrate inwards by some distance before reaching the inner edge of the disc. The outer edge was set at 0.8 au around 0.09 and 0.20 M_⊙ stars, 1.0 au around 0.49 M_⊙ stars, 2.0 au around 0.70 M_⊙ stars, and 2.2 au around 1.00 M_⊙ stars. These outer edges were motivated by initial estimates of pebble accretion efficiencies of planets with different initial masses around the respective stars (see Sect. 4). Based on these models, planets with radii smaller than 450 km around 0.09 or 0.20 M_⊙ stars do not accrete any pebbles at orbital radii larger than about 0.5 to 1.0 au. Around 0.70 and 1.00 M_⊙ stars, 450 km planetesimals can efficiently accrete pebbles at orbital distances larger than 2.0 au. However, if the semi-major axis range for planetesimals around these stars is expanded even further, the distance between them becomes so large that they are virtually isolated and have a very low probability of encountering other planetesimals. Given our limited number of planetesimals, we decided to limit the outer edge of the planetesimal ring to between 2 and 2.2 au for these stars.

In the second model, the planetesimals were released in a narrow annulus around the snowline. This model was motivated by the fact that the streaming instability requires a locally enhanced pebble density (Carrera et al. 2015; Yang & Johansen 2014; Yang et al. 2017, 2018), which the snowline could provide. Volatiles that evaporate from the pebbles inside the snowline could diffuse back across the snowline where the vapour pressure is low, and re-deposited onto the solids, enhancing the solid-to-gas ratio just outside the snowline (Schoonenberg & Ormel 2017; Drążkowska & Alibert 2017; Liu et al. 2019).

The resulting dense filament of pebbles that can collapse into planetesimals has a typical width of Δr ≃ ηr_snow (Yang & Johansen 2016; Li & Xiao 2016; Liu et al. 2019), which is approximately 0.0066 au for a solar-mass star using the disc conditions of Liu et al. (2019). However, Liu et al. (2019) assumed a fully irradiated disc with much higher temperatures than we did, which means that both their η and r_snow were larger than ours. Moreover, their simulations allowed for the planetesimals to be injected one by one, so that the simulations have time to create a semi-stable system. SyMBA injects all planets at once, because of which, initialising simulations with such narrow planetesimal rings is highly unstable. Therefore, we started these simulations with a ring that is slightly further expanded. Liu et al. (2019) found that within a few thousand years, the planetesimal ring has expanded to a width of about 0.1 au, due to gravitational interactions between the planets. We initialised our planetesimals in these simulations with a semi-major axis uniformly taken from between (0.85r_snow 1.15r_snow), with values shown in Table 2.

For both planetesimal distributions, the initial eccentricity and inclination are uniformly generated from the ranges (0, 0.01) and (0°, 0.3°), respectively. Though these initial ranges are small, the smaller planetesimals are rapidly excited to much larger values due to gravitational interactions with more massive planetesimals. The longitude of the ascending note, argument of periapsis and mean anomaly are all uniformly chosen from the range (0°, 360°).

Finally, we took an initial planetesimal ring mass of 0.010 M_⊕ for the 0.09 M_⊙ stars, and 0.015 M_⊕ for all other stars. This translates into approximately 275 massive particles in the 0.09 M_⊙ simulations and 400 in the other simulations. Liu et al. (2019) found that for their disc model, the planetesimal ring around a solar-mass star weighs about 0.039 M_⊕, using the fact that the total solid mass available to build planetesimals is $$M_{\text{avail\hspace{0.17em}}}=2\pi r_{\text{snow\hspace{0.17em}}}\mathrm{\Delta }r{\mathrm{\Sigma }}_{\mathrm{p}}\left(r_{\text{snow\hspace{0.17em}}}\right)=2\pi r_{\text{snow\hspace{0.17em}}}\mathrm{\Delta }r{\mathrm{\Sigma }}_{\mathrm{g}}\left(r_{\text{snow\hspace{0.17em}}}\right)\frac{H_{\mathrm{p}}\left(r_{\text{snow\hspace{0.17em}}}\right)}{H_{\mathrm{g}}\left(r_{\text{snow\hspace{0.17em}}}\right)},$$

of which approximately 50% is converted into planetesimals (Simon et al. 2016). However, using this equation together with our disc conditions results in unreasonably small planetesimal ring masses (10⁻³ M_⊕ for 1.00 M_⊙, 2.8 × 10⁻⁴ M_⊕ for 0.09 M_⊙). As previously mentioned, the core difference is that Liu et al. (2019) assumed a disc that is significantly hotter and more turbulent (α_turb = 10⁻³), and contains pebbles whose Stokes number is not drift limited, but fixed to 0.1. These three model differences lead to a difference of over an order of magnitude in the estimated initial planetesimal ring mass.

Assuming our disc’s turbulence is enhanced right after the disc’s formation, and the pebbles’ Stokes number has not yet reached the drift equilibrium and is instead 0.1 at the moment the planetesimals are being formed, our initial ratio H_p/H_g is comparable with that of Liu et al. (2019). In this case, we find an initial planetesimal ring mass of 0.013 M_⊕ for a solar-mass star, which matches reasonably well with the value of 0.015 M_⊕ we assumed.

4 Semi-analytical results: PA as a function of orbital radius and initial planetesimal mass

To gain a preliminary insight into the role of PA in planet formation, and the general trends expected in the full N-body simulation results in Sect. 5, the equations of Sect. 2 were evaluated using a semi-analytical algorithm. This algorithm determined the rate at which solitary planets on idealised orbits accrete pebbles. It did not include planetary migration, nor the influence of planets on each other’s orbit or accretion rate. In this section, the influence of a planet’s orbital radius and initial mass on its final mass is discussed. In Appendix A, the semianalytical algorithm is used to compare different pebble-sized models inside the snowline.

Figure 3 shows the relative growth (M_final/M_initial, M/M₀ for short) of these idealised, solitary planets as a function of orbital radius on the x-axis and initial mass on the y-axis. These results were calculated using the PA prescription of Ormel & Liu (2018) (OL18). The results for the prescription of Ida et al. (2016) (IGM16) are shown in Fig. 4.

The orbits of the planets in these calculations are circular and uninclined. The black contours indicate the regions of the parameter space within which planets grow more massive than a Mars mass (M_M), an Earth mass (M_E) and two Earth masses. The contour lines for planets more massive than 5 and 10 Earth masses are never reached in these inner regions of the disc due to the pebble isolation mass. The red dashed contour line indicates the boundary above which the planets have reached the pebble isolation mass and cannot grow any further.

Based on the results in Figs. 3 and 4, neither the OL18 nor the IGM16 PA model seems able to produce Earth-mass planets around 0.09 and 0.20 M_⊙ stars, at least not for the initial mass range of planetesimals in the SyMBA simulations, even though the isolation mass around these stars does allow for the formation of Earth-like planets. The IGM16 prescription allows for slightly more efficient growth, but even in this scenario, the planets do not grow much heavier than a Mars mass. In fact, in calculations with a 0.09 M_⊙ star with a 0.05 M_* disc, planets only reach up to about a lunar mass, being limited primarily by the low pebble mass flux, totalling only a few Earth masses (see Fig. 1 for the pebble mass flux in a 0.10 M_* disc for the 0.09 M_⊙ star). We predict that it is highly improbable for Earth-mass planets to form in an N-body simulation with such conditions. Therefore, we focussed on a more massive disc of 0.10 M_* around 0.09 M_⊙ stars, which contains twice the pebble flux a 0.05 M_* disc does. Nevertheless, we predict planetesimal accretion and mergers between several high mass embryos must play an important role for Earth-like planets and TRAPPIST-1-like systems to emerge from the N-body simulations with these low-mass stars.

A final, noteworthy feature in the 0.09 and 0.20 M_⊙ results of Figs. 3 and 4 is the discontinuity in M/M₀ at an orbital radius shorter than 1 au, the exact location of which depends on the mass of the central star. This feature is caused by the snowline, interior of which the pebble mass flux is halved. It is most visible in the 0.09 and 0.20 M_⊙ stars calculations. For the more massive stars, the feature is less apparent, since most of the planets in these calculations are not limited by the available pebble flux and accretion efficiency, but by the pebble isolation mass.

Around the 0.49, 0.70, and 1.00 M_⊙ stars, planetesimals systematically grow into Earth-mass planets for a wide range of orbital separations, for both OL18 and IGM16. However, unlike what Figs. 3 and 4 might suggest, only a few planetesimals are expected to actually grow to Earth-like masses in the full SyMBA simulations of Sect. 5, because the planets are not solitary, and need to share the total available pebble flux. Moreover, if a planet relatively far out in the disc reaches the isolation mass early on, it blocks the pebble flux to all planets interior to it, halting their growth.

The time at which the planets in Fig. 4 reach their isolation mass (M_iso) is shown in Fig. 5 for IGM16. The situation for OL18 is similar, except for the fact that the isolation mass is never reached around 0.20 M_⊙ stars. The OL18 results are therefore not separately shown. Around 0.49, 0.70, and 1.00 M_⊙ stars, the first planets reach M_iso within the first 50 000 to 100 000 years of the simulation. The range of orbital radii for which M_iso is rapidly reached, is wider, and the time in which the planets do so is shorter, the more massive the central star is. This could limit the growth of planets closer to the star prematurely, especially for planets inside the snowline, where growth is slower due to the reduced pebble mass flux. This also suggests that the initial disc conditions are more important than how the conditions evolve over time, at least for the largest planets in the system.

Fig. 3 Final mass divided by the initial mass of planets, calculated using only the semi-analytical PA model OL18, shown as a function of orbital radius on the x-axis and initial mass (M0) on the y-axis. The eccentricity and inclination are zero. Black contours indicate final masses higher than a Mars mass (MM), Earth mass (ME) and two Earth masses. Planets above the red dashed line have reached the local pebble isolation mass. Around 0.49, 0.70, and 1.00 M⊙ stars, PA produces Earth-mass planets in a wide range of orbital separations. Around 0.09 and 0.20 M⊙, PA does not form planets larger than a few Mars masses, even though the pebble isolation mass allows for the formation of Earth-sized planets. For 0.09 M⊙ stars, planets do not grow larger than about a lunar mass (ML) in the default disc of 5% of the stellar mass, which is why a more massive disc of 10% of the stellar mass is used for this star.

Fig. 4 Same plots as Fig. 3, but calculated using the semi-analytical PA model IGM16. Pebble accretion using IGM16 is more efficient than using OL18, even for the circular, uninclined orbits assumed here. This is most clearly visible from the contours in the 0.09 and 0.20 M⊙ subplots.

Fig. 5 Time after which planets reach the pebble isolation mass tiso ≡ t(M=Miso), as a function of orbital radius and initial mass, using IGM16. White indicates the isolation mass has not been reached. The contours showing the final mass are the same as in Fig. 4. Around 0.49, 0.70, and 1.00 M⊙ stars, the first planets reach Miso within the first 50 000 to 100 000 years of the simulation, which could halt pebble accretion for planets close to the star prematurely.

5 Full N-body simulation results

In this section, the results from the SyMBA N-body simulations are presented. These simulations include pebble accretion, planetesimal accretion, and type I and type II migration. Simulations were performed for 0.09, 0.20, 0.49, 0.70, and 1.00 M_⊙ stars, and for four different models: two models based on the PA prescription of Ormel & Liu (2018) (OL18 and OL18-Ring), and two for the prescription of Ida et al. (2016) (IGM16 and IGM16-Ring). In the ‘Ring’ simulations, the planetesimals were initialised on a narrow ring around the snowline. In the other simulations, the planetesimals were released over a wider range of orbital radii. All other disc conditions were identical for the different models (see Sect. 3 for more details).

For each combination of stellar mass and model, eight simulations were performed with randomly generated planetesimals. The simulations were run for 5 Myr of evolution. During the first ~2.5–3 Myr of growth, all bodies in the simulation were self-gravitating. During the final stages of the simulations, only planetesimals more massive than 10⁻³ M_E were considered self-gravitating, in order to reduce the computational load caused by small planetesimals. As will be shown in the sections below, nearly all planets form within the first million years. The influence of the smallest planetesimals during the final stages of the simulations is negligible.

Each simulation² consisted of ~400 planetesimals, with radii between 175 and 450 km, following a truncated Pareto distribution (see Eq. (22)). Planets that moved interior to the truncation radius, defined as the location around the star with an orbital period of 0.01 yr, were removed from the simulation. The size of the time step was 6.7 × 10⁻⁴ yr. In total, the simulations took approximately 6 months to complete, using 4 CPU cores per simulation.

In Sect. 5.1, the simulation results for a solar-mass star are presented, and two specific simulations are discussed in detail. The results for the other stars are provided in Sect. 5.4. Finally, the long-term stability of the systems, and the influence of gas accretion, which was not included in the standard results, are analysed in Sects. 5.6 and 5.7.

5.1 Overview of the simulation results for a solar-mass star

The resulting planetary systems after 5 Myr of simulation are shown in Fig. 6. This figure shows all simulations for all four models around a solar-mass star. The first thing that stands out is that, as predicted in Sect. 4, the IGM16 PA-prescription generates planets in greater number, and of significantly higher mass than the OL18 prescription. Nevertheless, regardless of the PA model, every system produces one or more Earth-like planets, which we define as planets with masses between 0.67 and 1.5 M_E (indicated by the black edge around the markers). Earth-like planets appear in comparable quantities in the general, wide planetesimal disc, spanning from 0.2 to 2.2 au, as in the narrow planetesimal ring around the snowline, spanning from about 1.2 to 1.6 au, even though in the latter, the probability of planetesimals coming together is larger. It is therefore unlikely that in our model, significant planetesimal accretion is a requirement for the formation of Earth-like planets, and even for these relatively small planetesimals of up to 450 km in radius, PA is sufficient for the formation of large planets.

Nearly all of these large planets have migrated to the innermost regions of the disc, especially in the OL18 simulation. The inner edge of the gas disc r_in is indicated by the vertical dashed line at about 0.1 au. This inner edge acts as a trap, preventing the planets from migrating further inwards because of the rapidly decaying gas surface density, which prevents the planets from producing the density waves that generate the torques required for planet migration. Without this boundary, all planets would have continued drifting into the central star.

The current understanding of the mechanisms preventing planets from migrating too close to the star is still incomplete, mainly because of uncertainties in the influence of torques near the disc’s inner edge (see e.g. Brasser et al. 2018). The use of a gas cavity within r_in is a simple, yet general solution, motivated by the star’s magnetosphere disrupting the innermost regions of the disc, causing the gas surface density to rapidly drop at around 0.05 to 0.1 au (Long et al. 2005; Romanova & Lovelace 2006; Romanova et al. 2019).

Nevertheless, many planets have migrated significantly further inwards than r_in. This is due to mean motion resonances (MMRs; Terquem & Papaloizou 2007). Migrating protoplanets often get captured in MMRs, forming chains of low-mass planets with orbital periods that are in resonance, migrating through the disc together. As the first planet reaches the inner edge of the disc, the entire chain slows down and stalls close to r_in. However, the resonance chain also transfers part of the torque that is experienced by the planets that are still in the gas disc, to the planets that are in the gas-free region, pushing them further inwards.

If the chain becomes too long and massive, it become dynamically unstable, leading to orbital crossings, particles being ejected, and giant impacts (Izidoro et al. 2017; Izidoro et al. 2021), which break the resonances. This typically happens when the gas disc disperses, or shortly thereafter. Since the disc in this study is exponentially drained with a diffusion time of 0.5 Myr, there is not a specific moment at which the gas of the disc has fully dissipated, but generally, most collisions happen within the first ~2 Myr (Ogihara et al. 2015; Zawadzki et al. 2021; Hatalova et al. 2023), though late dynamic instabilities can occur for up to a 100 Myr after the formation of the disc (not modelled in this study).

Evidence of giant impacts is seen in the IGM16 results in Fig. 6. Many of the planets in the innermost regions of these discs have masses exceeding 5 or even 10 M_E, far above the pebble isolation mass, which is closer to 1 or 2 M_E (see Figs. 3 and 4). The IGM16 models produce so many large planets, that the systems become unstable, causing the many protoplanets to merge into giant planets when the MMR chain reaches r_in.

Fig. 6 Simulated planetary systems around a 1.00 M⊙ star for the four PA models. The size of the markers represents the mass of the planet. The simulations are grouped by their model, and each horizontal line represents a simulation. The vertical dashed line at ~10−1 au is the inner radius of the gas disc, within which the gas density rapidly drops and type I and type II migration halts. Any planet that enters the gray region, interior to the inner truncation radius, is removed from the simulation. The cyan shaded regions represent the conservative (darker shaded) and optimistic (lighter shaded) habitable zone. The colour of the planets indicates their AMD stability, discussed in Sect. 5.6. Only planets more massive than Mars are shown.

5.2 Planets in the habitable zone around solar-mass stars

A consequence of the rapid planetary migration is that very few planets remain in the habitable zone (HZ), at least in simulations with the OL18 model. The location of the conservative and optimistic HZ are shown using the darker and lighter cyan shading in Fig. 6. The HZ was calculated using the algorithms from Kopparapu et al. (2013, 2014), and assumed solar values for the effective temperature and luminosity, and a planet mass of 1 M_E. The planet mass influences the expected thickness of the atmospheres, and therefore the maximum strength of the greenhouse effect.

The normal OL18 simulations produced only one Earth-like planet in the HZ, and one Mars-like planet, in the eight different realisations of the system, with two other simulations with Mars-like planets close to the HZ. The OL18-Ring simulations produced even fewer planets in or close to the HZ, even though all planetesimals started in the HZ. This shows that forming an Earth-like planet in or around the HZ using PA is not as challenging as keeping it there, given the rapid migration.

The IGM16 simulations produce far more planets in the habitable zone, even more in the ‘Ring’ configuration than in the normal planetesimal distribution. This is because the IGM16 model is more efficient than the OL18 model at accreting pebbles onto small planetesimals (see Sect. 4), especially onto those that have been excited (see Fig. 2), since the reduction in accretion efficiency for planets on eccentric orbits is ignored. The IGM16 prescription is therefore far more likely to produce planets at late stages of the disc evolution than the OL18 prescription. The gas of the disc dissipates before these late-forming planets have time to migrate significantly inwards, allowing them to remain in the HZ.

5.3 Dynamical evolution of solar-like systems

To get a closer look into the general growth track of planets, and the formation of planets in the HZ in particular, we present an analysis of the dynamic evolution of OL18 simulation 2, which produced the Earth-like planet in the HZ, in Fig. 7. This simulation shows both the general trends observed in OL18 simulations, and the specific sequence of events that led to the formation of an Earth-like planet in the HZ.

5.3.1 The dynamical evolution of OL18 systems

As can be seen in Fig. 7, the first massive planets form within the first 0.1 Myr of the simulation, as predicted by Fig. 5. Most of these planets formed from planetesimals at the high end of the size distribution. Higher mass seeds have an advantage over lower-mass seeds, since they can efficiently accrete pebbles for higher values of e and i, and have a higher chance of starting their growth early on, when there is still little competition for the pebble flux. However, the initial mass is not the primary limiting factor to the growth of planets.

This is demonstrated by the blue planet, which grew to be the biggest in the system, despite its seed being only 3.125×10⁻⁵ M_E (245 km in radius). One of the reasons this planet was able to grow so large is that it was initialised close to the outer edge of the planetesimal disc, at 2.02 au from the star. With no large seeds exterior to it, and with the interior seeds being far enough not to perturb its orbit, the blue planet was allowed to develop undisturbed.

The blue planet’s exponential growth started equally abruptly as that of the brown and pink protoplanets (shown dashed and transparent, for they later merged with other planets), which were initialised around the same region as the blue planet, but with significantly higher mass. The blue planet did not need to first gradually grow to a specific mass threshold before entering the rapid accretion regime. Instead, it had to rid itself of its initial eccentricity and inclination through disc interactions.

With a starting eccentricity and inclination of e = 2.75 × 10⁻³ and i = 3.25 × 10⁻², the blue planetesimal was initialised in the inefficient ballistic regime (see Fig. 2 for reference). Only when its eccentricity dropped below 10⁻³, at around 20 000 yr, did the planet enter the settling regime, and start runaway pebble accretion. By the time its eccentricity increased again, as a result of a close encounter with the red planet at around 50 000 yr, the blue planet was sufficiently massive for it to remain in the settling regime, despite its high eccentricity.

After about 100 000 yr, the blue planet reached the pebble isolation mass. The orange planet reached the isolation mass earlier, which promptly stopped the growth of the purple and gray (dashed) planets, but since the blue planet was exterior to the orange one, the blue planet could continue to grow. After the blue planet reached its isolation mass, it started migrating inwards, joining in an MMR chain with the brown (dashed), pink (dashed) and orange planet. Together, the chain migrated to the inner edge of the gas disc at 0.09 au in about 0.5 Myr.

Several dynamic instabilities, one of which being caused by the purple planet migrating inwards and joining the chain, led the brown and pink planets to collide with the blue and orange planets, respectively, and forced the latter two into the gas cavity. For a more detailed discussion about different mean motion resonance structures and their dynamic evolution, we refer to, for example, Brasser et al. (2022), and Hatalova et al. (2023).

The evolution of the planets discussed above is very typical for all OL18(-Ring) systems. The first planets start growing within 1000–10 000 yr of the planetesimal formation. After about 50 000–100 000 yr, the first planets in the outer disc reach the isolation mass, halting the growth of the protoplanets interior to them. Within about 0.5 Myr, the most massive planets in the outer disc migrate to the inner edge of the disc, forming a resonance chain with the other large planets they drag with them. When these large planets cross the orbits of the smaller protoplanets closer to the star, they excite the orbits of the smaller protoplanets. If the mass of these smaller protoplanets is great enough for them to quickly lose the induced eccentricity and inclination through disc torques (e.g. Matsumura et al. 2021), they resume their growth. Within the next ~1 Myr, this second group of planets grows to the isolation mass, and migrates to the inner edge of the disc, joining the MMR chain. This often leads to multiple dynamic instabilities, which push the innermost planets into the gas cavity, and cause other planets to merge or be ejected.

Fig. 7 Dynamical evolution tracks (Mp, a, e and i) of all large planets in OL18 simulation 7 around a solar-mass star. The different coloured solid lines represent different planets, the red one being the planet that is currently in the habitable zone (cyan shaded region). The transparent dashed lines represent large planetesimals that merged with the planets. Sudden stepwise increases in planet masses signify these mergers. Only planets with masses >0.1 ME are included; however, for this simulation, there are no other objects with masses of ≳10−3 ME.

5.3.2 The dynamical evolution of OL18 planets in the habitable zone

The formation of Earth-like planets in the habitable zone (HZ) in OL18 simulations is less trivial, and therefore far less common, than the general formation of Earth-like planets in these systems. Since most Earth-mass planets are formed in just a few ten thousand years, but rapidly migrate to the inner edge of the disc within the next million years, a very specific sequence of events is required for an Earth-like planet to remain in the HZ. The only realisation for this case is the red planet in Fig. 7.

This planetesimal was initialised relatively far out in the disc at 1.85 au, which is a prerequisite for any planet in our simulations to end up in the HZ, since planets in our model generally do not migrate outwards. Like the other planets, the red planet entered runaway pebble accretion once its eccentricity and inclination dipped far enough for pebbles to settle into its gravitational field, at about 35 000 yr into the simulation. However, its growth was stopped at a very specific moment in its evolution, which allowed it to end up in the habitable zone. This happened due to two consecutive close encounters, first with the brown (dashed) planet and then with the blue planet. These encounters kicked the red planet even further out in the disc and increased its eccentricity and inclination by about two orders of magnitude, which halted its growth. The key element of these events is that the planet had exactly the right mass to eventually end up in the HZ. Had it been more massive, then it would have dampened its eccentricity and inclination faster, reached an Earth-like mass sooner, and it would have had time to migrate to the inner disc, as the purple planet did. Had it been less massive, then it might not have lost its eccentricity in time, and would have remained a sub-lunar object, like all the other planetesimals in its vicinity.

In fact, the mass of the red planet at the time of its excitation might have already been slightly too high for it to remain in the HZ. After all, at the end of the simulation, the planet is still migrating inwards. This migration might quickly cease due to the gas disc becoming too thin, especially if we assume photoevaporation kicks in, and blows away the remainder of the gas. Nevertheless, it might also be that this planet drifts out of the HZ because of its mass, just as all the others, if the simulations were continued for another million years.

On the other hand, the planet required 4.5 Myr to reach its final mass, which means its mass should not have been much lower either at the time of its excitation, for it would not have had enough time to grow then. This only demonstrates how specific the conditions need to be for an Earth-like planet to remain in the HZ when using a simple migration model.

Nevertheless, since the planet formed outside the snowline and accreted all of its mass from pebbles there, it consists for up to 50% of water. Even if a large fraction of this water is lost to evaporation due to the heat from the planet’s formation, there should still be more than enough left for the planet to possibly be suitable for life.

Fig. 8 Dynamical evolution tracks (Mp and a) of all large planets in IGM16 simulation 1 around a solar-mass star. The different coloured solid lines represent different planets. The transparent dashed lines represent large planetesimals that merged with the planets. Sudden stepwise increases in planet masses signify these mergers. Only planets with masses >0.1 ME are included. The cyan shaded region represents the habitable zone. The thick black ellipses in the lefthand panel highlight the first, second, and third generation of planets. IGM16 produces far more planets than OL18. Those formed as part of the third generation, after 1–2 Myr, remain in the HZ because they have insufficient time for migration.

5.3.3 The dynamical evolution of IGM16 systems

Figure 8 shows the dynamic evolution of a single system with the IGM16 PA-prescription. The evolution of the largest planets follows the same patterns as of those in the OL18 systems. The main difference is that IGM16(-Ring) produces far more planets than OL18(-Ring). In fact, the IGM16 PA-prescription is so efficient that it is able to create a third generation of planets. This is most likely due to the fact that the influence of the eccentricity and the ballistic regime are not included in IGM16, which were the main limiting factors for growth in the OL18 simulations. As the number of massive planets in the OL18 simulation grows, the other planetesimals become increasingly excited, making it less and less likely that additional planets form. By ignoring this negative feedback loop, the IGM16 model likely significantly overestimates the probability of planets arising from the planetesimal disc.

Either way, similarly to the OL18 simulations, the first generation of planets forms far out in the disc, reaches the isolation mass within the first 50 000–100 000 yr, and then migrates inwards. The second generation of planets forms when the first generation moves interior to it, and reaches the isolation mass at about 1 Myr into the simulation. However, unlike in the OL18 model, as the second generation migrates inwards, a third generation of planets has time to grow, starting its final growth phase between 1 and 2 Myr into the simulation. These planets do have time to grow to Earth-like masses, but not to migrate significantly inwards, which is why IGM16 simulations produce significantly more planets further out in the disc, in particular in the habitable zone.

Furthermore, since the system is oversaturated with Earth-mass planets, it is highly unstable, leading to a lot of mergers between the innermost massive planet, and the other planets joining the MMR chain. As a result, the innermost planet absorbs many of its neighbours, growing to over 10 M_E. Gas accretion was not included in these simulations, and because of the low pebble isolation mass in these regions of the disc, planets generally do not grow massive enough to start gas accretion, for which cores of masses between 5 and 10 M_E are required (Mizuno et al. 1978; Stevenson 1982; Bodenheimer & Pollack 1986; Hubickyj et al. 2005). However, due to these many mergers between massive planets, the IGM16 model is capable of creating gas giants close to the star. These are discussed in Sect. 5.7.

On a final note, both Figs. 7 and 8 suggest a clear preference for planets to form in the outer regions of the planetesimal disc (r >~ 0.8 au), even though the planetesimals in these simulations are evenly distributed between 0.2 and 2.2 au. An important factor is likely that the pebble mass flux inside the snowline, which starts at about 1.4 au, is reduced by 50%, giving planets outside the snowline a significant advantage. Another reason could be that the Stokes number in the inner disc is too high for the small planetesimals to efficiently accrete pebbles due to the high relative velocity (see Fig. A.1).

5.4 Planet formation for different stellar masses

5.4.1 0.49 and 0.70 M_⊙ stars

The results of the simulations around 0.49 and 0.70 M_⊙ stars are shown in Fig. 9. The results are very similar to those around a solar-mass star. Earth-like planets are systematically formed, generally more than one per system. This is due to the fact that the pebble isolation mass to which nearly all planets grow in our simulations lies at around an Earth-mass for our assumed disc conditions, limiting the formation of larger cores.

Similar to in the solar-mass simulations, nearly all planets in the OL18(-Ring) simulations migrate to the inner edge of the disc, where they form mean motion resonance (MMR) chains, which push some planets even further in. Ultimately, no planet remains in the habitable zone. To avoid making assumptions about the effective temperature of the stars, the inner and outer edge of the habitable zone in these simulations were calculated using the simpler relationships: $$r_{\mathrm{i}}=\sqrt{\frac{L/L_{\odot }}{1.1}}\text{\hspace{0.17em}and\hspace{0.17em}}{r}_{\mathrm{o}}=\sqrt{\frac{L/L_{\odot }}{0.53}},$$

in which 1.1 and 0.53 are constants representing the stellar flux at the inner and outer edge of the HZ, respectively (Kasting et al. 1993; Whitmire & Reynolds 1996). Furthermore, L is the absolute bolometric luminosity of the star on the main sequence, which for M and K dwarfs is approximately given by (Cuntz & Wang 2018): $$\frac{L}{L_{\odot }}={\left(\frac{M_{\ast }}{M_{\odot }}\right)}^{n(M)}(0.20{\mathrm{M}}_{\odot }\le M\le 0.85{\mathrm{M}}_{\odot }),$$

with the fitted exponent $$n(M)=-141.7M^4+232.4M^3-129.1M^2+33.29M+0.215.$$

This estimate of the HZ is generally much more narrow than both the conservative and optimistic HZ of the more sophisticated model of Kopparapu et al. (2013, 2014), which was used to calculate the HZ around the 0.09 and 1.00 M_⊙ stars, which have been modelled to TRAPPIST-1 and the Sun, respectively. The HZ around the 0.20, 0.49, and 0.70 M_⊙ might therefore be underestimated.

For the IGM16(-Ring) simulations, relatively fewer Earth-like planets remain in the HZ than in the 1.00 M_⊙ simulations. This is because the HZ lies significantly closer inwards than the snowline, which means that the third generation no longer forms inside the HZ, but needs to migrate to it. Moreover, since growth around these stars is slightly slower, the third generation has less time to reach Earth-like masses, and especially for the 0.49 M_⊙ star, most late-forming planets grow only to a few Mars masses an example of which can be seen in Fig. B.1 in Appendix B.

Fig. 9 Simulated planetary systems around the 0.49 (left) and 0.70 (right) M⊙ stars for the four PA models. The size of the markers represents the mass of the planet, while their colour represents the planet’s AMD stability (see Sect. 5.6). The layout of the figure is explained in the caption of Fig. 6.

5.4.2 0.09 and 0.20 M_⊙ stars

The results of the 0.09 and 0.20 M_⊙ simulations are shown in Fig. 10. The most important observation is that no Earth-like planets form around 0.09 stars, irrespective of the model that is used, according to the predictions in Figs. 3 and 4. In half of the IGM16 simulations with a broad planetesimal disc around a 0.20 M_⊙ star, a single Earth-like planet managed to form through the merger of the many sub-Earth planets that arise in this model. However, the highest mass planets produced by the OL18(-Ring) model around 0.20 M_⊙ stars are a few Mars masses, suggesting that the OL18 PA model is unable to form Earth-like planets around these low-mass stars. This is in stark contrast with the results of Ormel et al. (2017), and Schoonenberg et al. (2019). The latter authors find that for each combination of their parameters, multiple Earth-mass planets form around a 0.09 M_⊙ star, even though they use the same OL18 PA-prescription that in our 0.09 M_⊙ simulations cannot produce planets more massive than a few lunar masses, regardless of the initial width of the planetesimal disc.

There are two core differences between our models and those of Schoonenberg et al. (2019). They assumed that a very massive disc of planetesimals forms from the streaming instability, weighing around 1.4 M_E for most of their models, with each planetesimal having a mass of 3.6×10⁻³ M_E (1200 km) in most simulations. For their fiducial model, the total mass in planetesimals formed by the streaming instability was about 7.5% of the total solid mass in their disc. We, on the other hand, used a planetesimal disc of only 0.010 M_E, which is comparable to the mass predicted by Liu et al. (2019).

Secondly, Schoonenberg et al. (2019) did not consider the evolution of the disc itself. In our model, the gas accretion rate onto the star decreases over time, as a result of which many of the other parameters, such as the location of the snowline, and the gas surface density, change as well. Moreover, in our model, once the pebble formation front reaches the outer edge of the disc, after about 0.2 Myr, the pebble flux rapidly declines (Sato et al. 2016). Growth during the first few hundred thousand years is therefore essential, and around 0.09 and 0.20 M_⊙ stars, this growth is too slow, at least for the small planetesimals that are expected from the streaming instability (Simon et al. 2016). Figure 3 suggests that for embryos larger than 10⁻³ M_E, planets of at least a few Mars masses should be able to form.

However, since these large planetesimals are not expected to form during the streaming instability, our results suggest that it is not possible to explain the formation of TRAPPIST-1-like systems using a simple evolving disc model with pebble accretion as the primary mechanism for growth. Further studies are needed to determine if an increased pebble mass flux could produce Earth-mass planets around these stars, and if such an increase can be justified.

Fig. 10 Simulated planetary systems around the 0.09 (left) and 0.20 (right) M⊙ stars for the four PA models. The size of the markers represents the mass of the planet, while their colour represents the planet’s AMD stability (see Sect. 5.6). The layout of the figure is explained in the caption of Fig. 6.

5.5 General overview and statistics

Figure 11 provides an overview of all planets that have formed over the different OL18 and IGM16 simulations, showing their mass, eccentricity, and inclination as a function of orbital radius. The results for the OL18-Ring and IGM16-Ring simulations are shown in Fig. B.2 in Appendix B, and are generally very similar.

In the OL18 simulations, nearly all planets have masses below 3 M_E, and no planet has a mass higher than 5 M_E, meaning that the planets are too small for gas accretion to play any significant role in their formation (Mizuno et al. 1978; Stevenson 1982; Bodenheimer & Pollack 1986; Hubickyj et al. 2005).

The IGM16 model generates planets in greater number, and of significantly higher mass, with several exceeding 10 M_E. Barring the planets around 0.09 and 0.20 M_⊙ stars (shown in blue and orange), whose mass was limited by their own growth rate, rather than the pebble isolation mass or the interference from more massive planets, most of the planets in the OL18(-Ring) simulations are clustered in a narrow region of the parameter space, with masses between 0.5–2 M_E, and orbital radii between 0.05–0.2 au. Meanwhile, the IGM16(-Ring) model produces planetary objects over the entire range from 10⁻⁵ to 10 M_E, and from 0.05 to 3 au.

The core difference between the models is seen in the eccentricity and inclination results, which points to the problem in the IGM16 model, which was already pointed at in Sects. 2.2.2, and 5.3.3. By ignoring the rapid decrease in accretion efficiency due to planets on eccentric or inclined orbits entering the ballistic regime, excited planetesimals in the IGM16 model continue growing.

The result is a clear correlation in e and i as a function of M and a in the IGM16(-Ring) model. After the disc is stirred up by the first generation of planets reaching their isolation mass and migrating inwards within the first ten to a hundred thousand years, a disc of small, heavily excited planetesimals remains. What follows is a continuous stream of small planets with high e and i, starting their growth at orbital radii between ~0.8 and 2.2 au, growing massive enough to dampen their e and i, and to migrate inwards, and finally joining the MMR chain in the inner disc, where they are excited again by dynamic instabilities. If the excitation-induced ballistic regime were included, most of these planets would not have formed, and the mass of the most massive planets would be significantly less, since they would have fewer other planets to merge with.

Figure 12 provides a few core statistics describing the systems that have formed, grouped by stellar mass and model. All quantities have been averaged over the eight simulations per model. The mean total mass in planetesimals and maximum planet mass follow very similar trends, with a strong increase in value from the 0.09 to the 0.20 to the 0.49 M_⊙ star, after which the values for the OL18 model flatten out, while they continue growing for the IGM16 model. Even if the entire planetesimal disc around the 0.09 M_⊙ star merged together, the resulting planet would not be larger than a few Mars masses in the OL18 model, and would be a little short of 1 M_E for the IGM16 model. This makes it highly unlikely that additional simulations, with slightly different initial conditions and minor tweaks to the model, could suddenly produce TRAPPIST-1-like systems with multiple Earth-mass planets. Significant model changes, such as a drastic increase in the pebble mass flux, a further disc edge, or other modifications that lead to far more rapid growth in the initial disc phases, are required to explain the formation of TRAPPIST-1. However, these changes are currently not supported by the existing literature.

The total mass in planetesimals remains fairly constant between the 0.49, 0.70, and 1.00 M_⊙ OL18 results, as does the maximum planet mass and number of planets, even though more massive stars have more massive discs and therefore more solids available for growth. This constant trend arises because, in all OL18(-Ring) simulations, the dynamical growth profile has the same shape, and when the first planets have reached the isolation mass, these planetesimal discs have all accreted more or less the same amount of mass. Most of the planets can no longer grow after this point, due to their excited orbits. This is irrespective of the time at which the point was reached and the amount of solids that are left in the disc. Since the IGM16 models do not have this limitation, the trend of increased planetesimal disc mass, maximum planet mass, and number of planets continues here.

Fig. 11 Overview of the mass, eccentricity, and inclination of all planets formed in the OL18 and IGM16 simulations. Earth-like planets have been highlighted using a black edge around the marker. In the e and i plots, the size of the marker is proportional to its mass.

Fig. 12 Core mean quantities of all simulations, grouped by stellar mass and model. All quantities are calculated per simulation, and then averaged over the 8 simulations per group. (A) Mean total mass of the planetesimals and planets in the disc. (B) Mean maximum planet mass. (C) Mean number of Earth-like (0.67–1.5 ME) planets. The number of Mars-like planets (0.1–0.67 ME), and massive planets (>1.5 ME) are stacked on top with increasing transparency.

5.6 System structure and long-term stability

As mentioned at the start of Sect. 5, the simulations in this study were terminated after 5 Myr. This end time does not necessarily mark the end of the dynamic evolution of the systems. Since there is no sharp cut-off time after which we assume all gas of the disc has dissipated, there is still a little gas present at the end of these simulations. As the rest of the gas slowly disappears, or is rapidly blown away by photoevaporation from the igniting star, further dynamic instabilities are expected to occur.

To make a first-order assessment of the long-term orbital stability of the planets in our systems, without resorting to computationally expensive simulations, we make use of two metrics: the Hill stability criterion (Chambers et al. 1996), and the angular momentum deficit (AMD; Laskar 1997, 2000; Laskar & Petit 2017; Petit et al. 2017).

The Hill criterion predicts that orbit crossings will occur if the separation between two planets is less than the critical value of $2\sqrt{3}$ times their mutual Hill radius. These orbit crossings result in system instabilities and collisions. The mutual Hill radius is given by $$R_{\mathrm{H}}=\frac{a_{\mathrm{p}1}+a_{\mathrm{p}2}}{2}{\left(\frac{M_{\mathrm{p}1}+M_{\mathrm{p}2}}{M_{\ast }}\right)}^{1/3},$$

where the subscripts p1 and p2 indicate the two planets. The separation of the planets is often represented by the dynamic spacing, Δ, which is the ratio between the separation and the mutual Hill radius.

For the evaluation of the system stability through the Hill criterion, we only consider the likelihood of orbit crossings between the larger planets in the system, since collisions of large planets with small planetesimals or planetary embryos do not significantly alter the system’s architecture. For the 0.49, 0.70, and 1.00 M_⊙ simulations, everything more massive than 0.1 M_E is considered a planet. Since the 0.09 and 0.20 M_⊙ simulations only contain a few planets (or none) in this mass range, the definition of a planet here is any object more massive than 0.01 M_E, which is the typical mass of a planetary embryo (e.g. Woo et al. 2021, 2022; Voelkel et al. 2021). The results are shown in Fig. 13.

None of the simulated planet pairs are Hill-unstable. Still, about 30% of the planet pairs in the OL18 simulation have a dynamic spacing Δ < 10, which suggests they are less likely to remain stable in the long-term (gigayear-scale) (Pu & Wu 2015). For the IGM16 simulations, this is 19% of the planet pairs.

However, the Hill stability should be interpreted with caution, especially when trying to estimate long-term stability. Planets with dynamic spacings Δ < 10 can still be very stable, especially if they are in orbital resonance with each other, as is the case for most of our simulated planets. Moreover, the Hill stability criterion is intended for two-planet scenarios, since it does not consider mutual inclinations between orbits, which are known to influence the evolution and long-term stability of planetary systems significantly.

For the multi-planet systems we observe, the AMD is a more appropriate measure of the stability, since it accounts for both eccentricities and inclinations. The AMD is a measure of the excitation of a system, indicating its deviation from a perfectly coplanar and circular system. If the total AMD is below a critical value, collisions between planets cannot occur, and the system is stable. The total AMD is given by (Laskar 1997, 2000) $$C=\sum _{k=1}^n{\mathrm{\Lambda }}_k(1-\sqrt{1-e_k^2}\cos i_k),$$(23)

with ${\mathrm{\Lambda }}_k=m_k\sqrt{\mu a_k}$ the angular momentum of planet k if its orbit were circular. The AMD stability coefficient is then given by $$\beta =\frac{C}{{\mathrm{\Lambda }}^{\mathrm{\prime }}{C}_{\mathrm{c}}},$$(24)

in which Λ′ is the circular momentum of the outer planet of the pair, and C_c = C_c(α, γ) is the critical AMD of the pair, which depends on the ratio of the semi-major axes, α = a_in/a_out, and the ratio of the masses, γ = m_in/m_out (see Laskar & Petit 2017). A planet pair is stable if β < 1, or equivalently, if log β < 0. A system is stable if all planet pairs are stable.

The AMD stability of the planets in the simulations is presented through the colour scales in Figs. 6, 9, and 10. For the calculation of the total AMD of the system, all planetesimals, and planets were included in the sum of Eq. (23). However, similarly to in the Hill stability calculations, we only consider the stability between orbits of planets more massive than 0.1 M_Earth (0.01 M_Earth in the 0.09 and 0.20 M_⊙ simulations). These are the planets shown in Figs. 6, 9, and 10. For the innermost planet, the stability is calculated with respect to the central star, in which case C_c = 1.

For the OL18(-Ring) simulations, the systems around the 0.49 and 0.70 M_sun stars have all reached stability, and the 2 out of 53 planets in the 1.00 M_⊙ simulations that are still unstable, are close to stability. This means that for the OL18(-Ring) simulations, planets settle into a long-term stable system very quickly, in less than 5 Myr. An exception to this is seen in the OL18 systems around 0.20 M_⊙ stars, most of which have not yet reached stability. Approximately 25% of the planet pairs in these OL18 simulations are unstable, and 10% in the OL18-Ring simulations, meaning that these systems would likely experience mergers if the simulations were extended beyond the 5 Myr, possibly leading to larger planets. However, since most of these unstable planet pairs involve a Moon and a Mars mass object, or in the most extreme case two Mars mass planets, it is highly unlikely that Earth-mass objects would suddenly arise if the simulations were continued for another few million years.

The IGM16(-Ring) simulations are a different story. At the end of our simulations, these systems are still highly unstable, especially for more massive stars. These systems are dynamically active for far longer than the OL18 systems, still forming massive planets after 2 Myr of evolution. Moreover, since these systems are packed with massive planets, there have been more close encounters, leading to greater average eccentricity and inclination (see Figs. 11 and B.2), and the MMR chains are rendered more unstable. To find the long-term architecture of these systems, the simulations must be extended for at least a few more Myr, but likely until 40 Myr or even beyond 100 Myr (Hatalova et al. 2023).

Fig. 13 Dynamic spacing Δ of all planet pairs in all simulations, plotted against their period ratio Pout/Pin. The planet pairs are grouped per PA model, with their colour indicating the stellar mass. The horizontal lines indicate two important values of the dynamic spacing: the critical value $\mathrm{\Delta }=2\sqrt{3}$ (solid), and the value Δ = 10 (dashed).

5.7 The influence of gas accretion on the simulated systems

As previously mentioned, gas accretion was not included in the simulations. For most simulations, gas accretion should not play any significant role, since cores with masses between 5 and 10 M_E are required to start runaway gas accretion (Mizuno et al. 1978; Stevenson 1982; Bodenheimer & Pollack 1986; Ikoma et al. 2000; Hubickyj et al. 2005). Only in the IGM16 simulations around 0.70 and 1.00 M_⊙ stars are there occasionally planets within this mass range.

To validate that gas accretion did not play a significant role in the OL18 simulations, and to get a first-order idea of the influence of gas accretion on the IGM16 results, the simulations were repeated, starting from the snapshots in which the first generation of planets have just formed. The gas accretion rate in these simulations is given by the equations provided in Matsumura et al. (2021). To conserve computing resources, the minimum mass for self-gravitating particles was set at 0.005 and 0.01 M_E for the OL18 and IGM16 models, respectively. This threshold is too high for all the dynamics of the system to be included, since the gravitational pull of objects slightly smaller than the Moon can influence the orbits of other objects, which is not included in this model. Moreover, the threshold was imposed early in the disc, at a time when sub-lunar objects could still grow through pebble accretion, although in the case of the OL18 model, they hardly ever do. These results therefore do not supersede the results presented in the previous sections. They merely serve to provide an indication of whether gas accretion could significantly alter the conclusions of this study.

The results for a 0.70 and 1.00 M_⊙ star are presented in Fig. 14. The other stars are not included in the figure, since gas accretion did not play a significant part in these cases. As expected, gas accretion has an insignificant influence on planet formation in the inner disc for the OL18 model. The pebble isolation mass is too low for planets to grow to the required mass for runaway gas accretion through PA, and too few Earth-like planets form for the MMR chain to collapse into massive planetary cores.

In the IGM16 simulations, however, multiple gas giants form through the mergers of the many planetary embryos that are created. The IGM16 simulations produce more gas giants than the IGM16-Ring simulations, with some systems around the solar-mass star producing two ~100 M_E gas giants in the same simulation. However, since the number of planets formed in the IGM16(-Ring) simulations is most likely a gross overestimation, so is the number of gas giants. Future studies should confirm or deny this.

6 Discussion

This study further cements the idea that pebble accretion could explain the formation of Earth-mass planets around low-mass stars of ≳0.5 M_⊙. Around ≲0.20 M_⊙ stars, no Earth-like planets were observed, except for in the 0.20 M_⊙ IGM16 model, which managed to produce a single Earth-like planet in half of the simulations.

Fig. 14 Simulated planetary systems around the 0.70 and 1.00 M⊙ stars for the four PA models including gas accretion. The mass threshold for particles to become self-gravitating was increased compared to the simulations without gas accretion to conserve computing resources, meaning that these simulations are less accurate than the others. The purpose of this figure is therefore only to indicate whether gas accretion could significantly influence the results from the normal simulations. For the OL18 simulations, gas accretion does not play a significant part. In the IGM16 simulations, on the other hand, multiple gas giants form. The layout of the figure is explained in the caption of Fig. 6.

6.1 OL18 vs IGM16

In this study, two different prescriptions for the pebble accretion efficiency ϵ were tested, one by Ida et al. (2016, IGM16), and one by Ormel & Liu (2018, OL18). A few key measurements that show the strong difference between these models are presented in Table 3. In general, the IGM16 model produces many more planets, with a significantly higher maximum mass, than the OL18 model. The IGM16 models are also much more excited than the OL18 models, with total angular momentum deficits (AMD_total; a measure of the deviation from a completely circular system) that are orders of magnitude higher, especially for the most massive stars in the simulation, as a result of which most IGM16 simulations are still highly unstable after the 5 Myr of simulations.

Both models have their strengths and weaknesses. For example, the OL18 model uses orbit-averaged corrections for the eccentricity and inclination dependence of the accretion rate, but these do not take into account the influence of the argument of periapsis, nor the variations in encountered disc conditions along eccentric orbits. Nevertheless, these small inaccuracies will likely not have a significant influence on the overall accretion of planets, since these parameters only become relevant for highly eccentric and inclined orbits, in which virtually no accretion takes place.

The shortcomings of IGM16, on the other hand, are very significant. As has extensively been discussed in Sects. 2.2.2, 5.3.3, and 5.5, the core issue of IGM16 is that it assumes that all planets are on circular orbits and that PA is always in the settling regime. The OL18 simulations show, however, that the main limitations to the growth of planetesimals are their eccentricity and inclination, which can cause them to enter the highly inefficient ballistic pebble accretion regime. This effect is especially of importance given that, as the first planets reach around an Earth-mass and migrate to the inner edge of the disc within a few hundred thousand years of initialisation, they stir up the rest of the disc, inciting high eccentricities and inclinations in all other planetesimals. This increases the relative velocity, Δv, between the planetesimals and the pebbles so much that the pebbles no longer have time to settle in the gravitational field of the planetesimals, which significantly decreases the accretion efficiency.

In a realistic scenario, only the planetesimals that had already significantly grown before their excitation, are able to rid themselves of their eccentricity and inclination through interactions with the disc, and resume their growth in the later stages of the disc. This is what is observed in the OL18 simulations. However, by ignoring the effects of the excitation and assuming planets can continue growing uninteruptedly, IGM16 grossly overestimates the accretion rate of planetesimals after the formation of the first planets, and produces a large second, and even a third generation of planets. The many Earth-mass planets in the second generation migrate to the inner disc edge, where they cause dynamic instabilities in the mean motion resonance chains and merge with the first generation of planets, producing multiple massive super-Earths and even gas giants. The third generation, on the other hand, does not have time to migrate inwards and fills the entire disc with Earth-like planets.

The IGM16 model most likely overestimates the number of planets that are formed, the final mass of the planets at the disc’s inner edge, and the probability of encountering Earth-like planets in the habitable zone. A comparison with observations of exoplanetary systems should be performed to fully validate this. The IGM16 model could be very powerful and useful, since its analytical approach and general intuitiveness allow for a very detailed analysis of the influence of different disc conditions on planet formation. However, until its expressions for Δv are augmented to include the influence of the eccentricity and inclination, and equations for the accretion rate in the ballistic regime are added, the PA prescription from Ida et al. (2016) is not appropriate for N-body simulations. In the comparison with other works in the following section, therefore, we only considered the OL18 results.

6.2 Comparison with other studies

In Sect. 5.4.2, we compare our results for a 0.09 and 0.20 M_⊙ star to those of Schoonenberg et al. (2019). We find that without their assumption of a massive planetesimal disc and with the inclusion of disc evolution, Earth-like planets cannot form around these low-mass stars.

In this section, we compare our results to two other recent studies: the analysis from Lambrechts et al. (2019) on the influence of the pebble mass flux on the type of planetary system that forms and that of Johansen et al. (2021), who used more specific PA simulation set-ups in an attempt to explain the mass, orbits, and compositions of the terrestrial planets in the Solar System.

The claims of Johansen et al. (2021) are still very much disputed (see e.g. the recent work of Morbidelli et al. 2025 and the response by Johansen et al. 2024). The analysis of currently available meteorite samples shows a dichotomy in the isotopic composition between the inner and outer Solar System (Warren 2011), which is inconsistent with our current understanding of PA. After all, most of the accreted pebble mass originates in the outer disc, so if PA played a significant part in the formation of the terrestrial planets around the Sun, their isotopic composition should be equal to that in the outer Solar System (Mah et al. 2022). In this section, we, therefore, do not compare our results to those of Johansen et al. (2021) to try to explain the formation of the Earth specifically, but to analyse the formation of Earth-mass planets at significant orbital radius (in particular within the habitable zone) around Sun-like stars through PA in general. We refer to these systems as ‘Solar System-like’, neglecting the isotopic dichotomy of the actual Solar System.

We show that the artificial introduction of a few 10⁻³ M_E planetary embryos, as done by Johansen et al. (2021), is not necessary for PA to form Earth-mass planets and that the planetesimal size distribution following the streaming instability contains sufficiently massive objects for rapid growth. However, with our general model, we are unable to explain the formation of the aforementioned Solar System-like planetary systems. Because of the rapid growth and migration in our model, there is a low probability of planets remaining in the habitable zone. More importantly, in every simulation, there is a first generation of Earth-mass planets that has migrated to the inner edge of the disc. These planets are not present in the Solar System.

For Solar System-like systems to form, some force must be present to prevent planetary embryos from growing during the first 0.5–1 Myr or there must be pressure bumps in the disc preventing migration to the inner edge of the disc (Chambers 2023). Johansen et al. (2021) focussed on preventing early growth by assuming all initial planetesimals have a radius of 100 km (which is too small for them to growth through pebble accretion) and by introducing the actual protoplanets of 10⁻³ M_E representing Venus, Earth, Theia, and Mars after t = 0.66, 0.93, 1.50, and 2.67 Myr, respectively. Though these starting times were based on the expected growth time to 10⁻³ M_E from an initial mass function of planetesimals that peaks at M ~ 10⁻⁷ − 10⁻⁶ M_E, they seem inconsistent with our systems that also include a few planetesimals at the higher end of the expected mass range (Simon et al. 2016). The ~400 km planets in our simulations grow into massive planets much earlier. Even if relatively more very low-mass planetesimals were included in our model (which could coagulate over a longer period of time to form embryos later on), the first generation of planets that grew directly from the most massive planetesimals from the streaming instability would still be present in the innermost regions of the systems (see e.g. also Mah et al. 2022).

An alternative explanation for our observed systems and why they differ from the Solar System is presented by Lambrechts et al. (2019). They propose that there are two modes of growth, one that produces super-Earths – planets with masses between those of Earth and Neptune – at the inner edge of the disc, and one which produces Earth-like planets spread throughout the disc. The growth mode the system experiences is determined by the radial pebble mass flux. For a low pebble flux, such that the total pebble mass that reaches the inner disc is less than ≈110 M_E around a solar-mass star, the embryos within the snowline grow slowly without significant migration. The resulting widely spaced population of approximately Mars-mass embryos becomes unstable when the gas disc fully dissipates, which they assume is after 3 Myr. Collisions between these Mars-mass embryos create a small number of terrestrial planets, with masses of at most five Earth-masses. For high pebble fluxes, with a total pebble mass of more than ≈190 M_E around a solar-mass star, the embryos within the snowline rapidly grow sufficiently massive to migrate to the inner edge of the gas disc, where they continue accreting pebbles, and merge as a result of dynamic instabilities, forming a system of closely spaced super-Earths with masses between five and twenty Earth-masses.

This latter growth mode is in some aspects very akin to what we see in our 0.49, 0.70, and 1.00 M_⊙ systems. This would be consistent given that the total pebble mass in the 1.00 M_⊙ simulation is ≈215 M_E (see Fig. 1), which corresponds to 1.3% of the total disc mass, placing the simulation in the super-Earth growth mode. However, in many aspects, our simulations are also significantly different. For instance, we find that most planets form outside the snowline. Moreover, we find that the pebble isolation mass is closer to 1 M_E (Ataiee et al. 2018), instead of 10 M_E, as assumed by Lambrechts et al. (2019). This, combined with the fact that only few planets form, means that the objects at the inner edge of the disc are not in the 5–20 M_E range, but in the 1–3 M_E range.

Moreover, Lambrechts et al. (2019) filled their discs with 25 embryos of 10⁻² M_E each, which could grow into Mars-mass protoplanets for low pebble fluxes. It is not yet clear if the small planetesimals used in this study would be able to grow fast enough and in sufficiently large numbers under these conditions.

That being said, the 0.20 M_⊙ simulations in Fig. 10, for which the total pebble mass (calculated by integrating Eqs. (7) and (10)) is approximately 0.4% of the total disc mass, do somewhat resemble the situation sketched in the terrestrial growth mode of Lambrechts et al. (2019), suggesting that there might be a sweet spot in our simulations as well, for the formation of many Mars-like protoplanets which can form Earth-like planets with large orbital radii after the complete dissipation of the disc. However, if there is such a pebble flux sweet spot, it will most likely be extremely narrow, given the fact that the planetesimals in the 0.20 M_⊙ OL18 simulations have not yet accreted enough material to form a single Earth-like planet if they were all brought together, while the protoplanets in the 0.20 M_⊙ OL18-Ring simulations have already migrated to the inner edge of the disc, despite not having accreted sufficient material either. Further research, using different pebble mass fluxes and disc dissipation conditions, is needed to determine whether Solar System-like planetary systems can be created using our general PA model.

6.3 Limitations of this study

Most of the caveats and assumptions of this study have been discussed in Sects. 2 through 4. These include the approximations that go into the OL18 and IGM16 PA prescriptions in Sect. 2.2, the planetesimal number and distribution in Sect. 3.2, and the consequence of the Stokes number being ≫1 in the innermost regions of the disc, discussed in Appendix A.

Other aspects of the disc and planets that might influence PA and migration, but that have been ignored in this study, include gas accretion, as discussed in Sect. 5.7, disc winds, which could limit type I migration in the inner disc (Ogihara et al. 2018; Chambers 2019), the formation of a primordial atmosphere around the planetary embryos (Brouwers & Ormel 2020; Takaoka et al. 2023; Yzer et al. 2023), which could alter the PA efficiency and the planetary composition, and the early formation of gas giants further out in the disc, which could limit the pebble flux to the inner disc. However, these are all higher-order effects that are beyond the scope of this research.

There are three model caveats based on recent developments in the field of PA that we think warrant a slightly more detailed discussion than the list above. These concepts have not yet been widely adopted in numerical PA simulations, but they are worth considering in future work.

Firstly, the migration model used in this study did not consider thermal torques. Recent work shows that in non-isothermal discs with low viscosities, such as the discs used in this study, thermal torques can slow down, or even reverse planet migration (Masset 2017; Guilera et al. 2019; Guilera et al. 2021), possibly providing a means to limiting the rapid inwards migration seen in most of our simulations. The thermal torque consists of a cold component, originating from two cold and dense lobes that appear on either side of the planet’s orbit due to thermal diffusion of the disc, and a heating component, which has the opposite sign of the cold component, and is caused by the heat released by the planet from accreting solids, which creates hot and low-density lobes in its surroundings (Masset 2017; Guilera et al. 2019).

However, as noted by Guilera et al. (2019), the thermal torques of Masset (2017) are incompatible with the torques from Paardekooper et al. (2011) used in this study. While Masset (2017) expanded upon the work of Jiménez & Masset (2017), deriving the thermal torque by considering the effects of thermal diffusion in the disc which had previously been ignored, Paardekooper et al. (2011) derived their migration torques based on 2D radiative hydrodynamics simulations which already included thermal diffusion. The torques from Paardekooper et al. (2011) might, therefore, at least partially incorporate the effects of the cold torque. The heating component of the thermal torque has not been included in the models of Paardekooper et al. (2011), since they do not consider the planet’s luminosity. However, adding only this heating component without the full cold component is not appropriate either, since it is the balance between the cold and heating components that determine the influence of the thermal torque. The full migration model used in this study should therefore be replaced by that of Masset (2017) for the thermal torque to be considered. We leave this for our future studies.

Secondly, like most PA studies, our models considered monodisperse (single-sized) pebble accretion. The pebble size does depend on the local disc conditions, and is physically motivated by the drift barrier, but it does not include a size distribution. Lyra et al. (2023) show that in polydisperse (multi-species) PA, the mass of onset of PA lies 1–2 orders of magnitude lower, and the PA efficiency of small planets in the Bondi regime lies 1–2 orders of magnitude higher than in monodisperse PA. In these scenarios, the most efficiently accreted pebble size can be different from the size that dominates the pebble mass flux. Polydisperse PA could promote the formation of larger planets in the 0.09 and 0.20 M_⊙ simulations and might therefore help explain the formation of systems like TRAPPIST-1. We therefore considered employing polydisperse PA models in future studies.

The final model caveat is related to the pebble isolation mass (PIM) of Ataiee et al. (2018) used in this study. This semi-analytical expression for the PIM depends only on the disc conditions at the location of the planet. However, using numerical simulations, Chametla et al. (2022) show that the PIM also depends on the eccentricity of the planetary orbit. For eccentricities lower than the disc’s aspect ratio h_g, an increase in eccentricity leads to a decrease in PIM. For even higher eccentricities, the PIM increases monotonically. Chametla et al. (2022) find that for low-viscosity discs (α_turb ≤ 10⁻³), like the ones used in our study, the PIM is well-defined, and its eccentricity dependence is best reproduced by the expression $$\begin{array}{rl}{M}_{\text{iso\hspace{0.17em}},2}(\hat{e})=& \frac{M_{\text{iso\hspace{0.17em}}}(\hat{e}=0)}{M_{\mathrm{E}}}\\ & \times \{1+(400{\alpha }_{\text{turb\hspace{0.17em}}}+0.2)\hat{e}-\frac{5\hat{e}}{[8+{\hat{e}}^3]}\}{M}_{\mathrm{E}},\end{array}$$(25)

in which $M_{\text{iso}}(\hat{e}=0)$ is the PIM of Ataiee et al. (2018) in Eq. (16), and $\hat{e}=e/h_{\mathrm{g}}$.

To get a first idea of how much this new pebble isolation mass would have influenced the planets in our simulation, we used the planets evolutionary data, shown in e.g. Fig. 7, to trace the ratio between the mass (M_p) and PIM (M_iso,2) of each planet throughout their pebble-accreting growth phase. The maximum values are reported in Fig. 15.

The influence of the eccentricity dependence of the PIM is relatively small for most planets. None of the planets in the OL18 − 0.09 M_⊙ simulations have reached M_iso or M_iso,2. For the OL18–1.00 M_⊙ simulations, 80% (94%) of the planets do not exceed M_iso,2 by more than 5% (7.5%), with the highest excess being 15%. The planets around 0.49 M_⊙ stars seem to be influenced most, with 42% of the planets exceeding M_iso,2 by more than 5%, and a maximum excess of 36%. Similar values are found for the IGM16 simulations, though these results are slightly skewed by the large sample of third generation planets that never come close to M_iso.

It is important to note, however, that these results give an upper limit to the effects using the PIM of Chametla et al. (2022) in our simulations could have had. A maximum M/M_iso,2 of 1.1 does not necessarily mean an overestimation of the planet’s mass by 10%. There is a balance between the growth time scale and the eccentricity demping timescale. If a planet reaches M_iso,2 and stops growing, it could continue to circularise its orbit, assuming that it is not close to another large planet. This in turn brings M_iso,2 closer to M_iso, an effect which is not considered in the calculation of Fig. 15. Moreover, though only parts of the planetary evolution with a non-zero pebble accretion rate are included in the calculation of the maximum M/M_iso,2, the magnitude of this accretion rate is not considered in the estimations. For most slow-growing planets, such as those around 0.20 M_⊙ stars, the planets do not grow bigger than M_iso,2, as much as M_iso,2 shrinks below their mass due to orbital excitations at later stages of the disc evolution.

The effects of an increase in isolation mass due to high eccentricity cannot be traced from our data, because there are too many scenarios other than M_iso,2 > M_iso that lead to a maximum M/M_iso,2 < 1, such as low PA rates, late formation, or exterior planets reaching the isolation mass. However, especially for the OL18 simulations, it is unlikely that planets would grow bigger when using M_iso,2 instead of M_iso, since planets with high eccentricities cannot efficiently accrete pebbles.

All in all, it is unlikely that the implementation of the pebble isolation mass of Chametla et al. (2022) in our simulation would influence the final planet masses by more than a few percent. Nevertheless, it is worth including in all future studies, since it is a more complete model of the pebble isolation mass.

Aside from the three modelling limitations discussed above, there are a few caveats following from computational constraints. The available computing time governed the number of planetesimals, step size, and inner truncation radius. As explained in Sect. 3.1, the step size used in this study is slightly larger than advised by Wisdom & Holman (1991). However, Hatalova et al. (2023), who used an even larger step size, showed that decreasing their step size by a factor of four had no significant influence on the outcome of the simulations. Moreover, whereas they studied planetary growth through planetesimal accretion, for which the precise dynamics between all planetesimals is essential, our main growth method is pebble accretion, which is much less sensitive to the step size. Naturally, a sufficiently small step size is important for our simulations as well, in order to properly detect close encounters³, and to model the dynamics of the systems that form at the inner edge of the disc. However, if the planetesimal accretion simulations of Hatalova et al. (2023) are not significantly influenced by the chosen step size, then it is highly unlikely that our pebble accretion simulations would be.

Another important caveat comes in the form of the aforementioned inner truncation radius, interior to which planets are assumed to have merged with the star and are removed from the simulation. This is a necessary computational constraint that follows directly from the step size, as explained in Sect. 3.1. However, this radius is not a physical boundary, and there are stars with planets on stable orbits interior to this boundary. As can be seen in Figs. 6 and 9, there are several systems with planets that have been pushed close to r_trunc. In fact, several massive particles have been removed from the simulations because they crossed the boundary, primarily in the IGM16 runs. A smaller value of r_trunc is required to determine the actual stability of these ejected planets. However, this requires a much smaller step size, which was not possible given the available computing resources.

Fig. 15 Cumulative distribution function of the maximum ratio between planet mass and eccentricity-dependent PIM for the planets in the different simulations. Only growth phases with non-zero pebble accretion rates for planets with a final mass above 0.01 ME are considered in the calculation of the maximum.

6.4 Recommendations for future research

A useful next step for this research would be to perform an in-depth comparison between the simulated systems, and the observed exoplanetary systems, to determine which types of systems are most similar to the ones we predict. Other important steps include improving the IGM16 model to include expressions for the ballistic regime, and the influence of the eccentricity, inclination, and argument of periapsis on the relative velocity between the planetesimals and the pebbles. Furthermore, chemical models for the evolution of the pebble size inside the snowline should be performed, to test our proposed sublimation model.

However, the most relevant next step in our opinion would be to develop a general model that can explain the formation of Earth-like planets in the habitable zone systematically. This means either limiting the initial growth rate of the planets, or decreasing their migration speed. As discussed in Sect. 6.2, the pebble mass flux and the disc dissipation conditions are the most promising parameters for limiting the planetary growth rate, based on the work of Lambrechts et al. (2019). Performing a parameter study of these two variables is, therefore, a very important continuation of this research. Limiting migration could be done by introducing a pressure maximum close to the habitable zone, for example created by disc winds in the inner regions of the disc, or sublimation of volatiles from pebbles at the snowline. Moreover, as discussed in sect. 6.3, the inclusion of thermal torques from Masset (2017) might reduce the rate of inwards migration. Further research is required to determine how realistic these scenarios are, and how they are best combined.

7 Conclusion

We studied the formation of terrestrial planets around low-mass stars using a version of the N-body integrator SyMBA (Duncan et al. 1998), modified to include pebble accretion (PA), type I and II migration, and eccentricity and inclination damping (Matsumura et al. 2017, 2021). The main findings from these simulations are as follows:

• Earth-like planets are consistently formed around 0.49, 0.70, and 1.00 M_⊙ stars, irrespective of the model that is used;

• Around 0.09 and 0.20 M_⊙ stars, no Earth-mass planets are formed. Even if the final mass of all planetesimals in the disc would somehow concentrate onto one planet, the planet would not be more massive than a few Mars masses;

• In the 0.49, 0.70, and 1.00 M_⊙ simulations, the first planets reach the pebble isolation mass within a few hundred thousand years. As these planets rapidly migrate inwards, they excite the orbits of the other planetesimals so much that PA becomes highly inefficent, limiting planetary growth at later times;

• The IGM16 model does not consider the fact that planetesimals with high eccentricities and inclinations enter the ballistic regime. Therefore, this model unrealistically generates far more Earth-like and higher mass planets than the OL18 model;

Overall, PA has a high tendency to create Earth-like planets around low-mass stars of about 0.5 M_⊙ and higher. However, for a planet to remain in the HZ, a series of very specific events must occur. The planet must start its formation relatively slowly and far out in the disc. Its growth must then be stopped at precisely the right moment by interference from a massive planet around it. If its mass is too high, it will rid itself of the induced eccentricity and inclination too soon, grow to the isolation mass too fast, and migrate to the inner edge of the disc like all other planets. If its mass is too low, it holds on to its orbital excitations for too long and will not have time to grow to an Earth-like size, due to the decaying pebble flux.

If other processes in the disc, such as the sublimation of volatiles at the snowline (or disc winds) could create a pressure maximum within the habitable zone, which would act as a sufficiently strong migration trap, the probability of finding Earth-like planets in the HZ would significantly increase. A lower initial growth rate due to a reduced pebble flux could also increase this likelihood. Further research is needed to determine whether these scenarios are realistic.

Acknowledgements

We would like to thank Prof. Dr. C. Dominik for the insights he provided during discussions about this work. We also thank Prof. Dr. W. Lyra for his comment after reading a pre-print version of our paper. Finally, we thank the anonymous referee for their valuable commentary, making this paper more concise and clear. This work used the Dutch national e-infrastructure with the support of the SURF Cooperative and the Dutch Research Council (NWO) using grant no. EINF-7075. M.J. Yzer further acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program under grant agreement NO. 805445.

Appendix A The pebble radius within the snowline

Figure A.1 shows the pebble radius R_peb (top panel) and Stokes number τ_s (bottom panel) as a function of orbital radius for three different times. As discussed in sec 2.2.1, the pebble radius in this study is drift-limited, which is to say, the pebble grows until it reaches the critical Stokes number τ_crit,1 for which t_grow = t_drift. The focus of Fig. A.1 is on what happens once the pebbles cross the snowline. In the original model, here referred to as ‘no sublimation’, the equations of Ida et al. (2016) are smoothly extrapolated to within the snowline, so that the pebble radius and composition are unaltered. However, this assumption is highly unrealistic and inconsistent with the assumption that sublimation of volatiles causes the pebble mass flux to be halved (e.g. Morbidelli et al. 2015).

An alternative model is that of fragmentation. This model assumes that due to the sublimation of volatiles, pebbles break apart at the snowline into mm-sized silicate grains, and remain at this size due to the bouncing and fragmentation barrier (Morbidelli et al. 2015; Ida et al. 2016; Matsumura et al. 2017). However, the mechanics behind fragmentation, and whether it occurs in the first place, are still uncertain. Moreover, because in the fragmentation model τ_s ⋘ 1, other model assumptions, such as the assumption that the radial viscous diffusion velocity term u_ν is negligible, might become invalid.

We therefore propose a third model, which we refer to as the ‘sublimation’ model. Similarly to the fragmentation model, we assume that pebbles disintegrate into dust at the snowline due to sublimation of their icy contents. However, this dust then recoagulates to a new radius, with a new density dominated by rock. We assume ρ_s ~ 1.0 g cm⁻³ outside the snowline, and ρ_s ~ 2.5 g cm⁻³ inside the snowline. This, combined with the 50% reduction in the pebble mass flux ${\dot{M}}_{\mathrm{F}}$, results in a change in pebble radius in both the Epstein and the Stokes regime, as well as a shift in the location of the of boundary r_ES between the two regimes.

In fact, sublimation could introduce a second and third Epstein-Stokes boundary, as can be seen in the 1.0 Myr line (shown in orange) of R_peb in Fig. A.1. Initially in the outer disc, the pebble is in the Epstein regime and grows as it drifts inwards. Just outside the snowline, the pebble transitions from the Epstein regime to the Stokes regime, signified by the circle marker, after which it no longer grows. However, once the pebble crosses the snowline and falls apart, it grows back to the new value of τ_s,crit1 (see Ida et al. 2016, for the full expressions), which is again in the Epstein regime. As the pebble drifts further inwards, it continues growing, until it transitions to the Stokes regime one final time.

The pebble radius in the Stokes regime still has a few caveats, though. For instance, there is a discontinuity in the pebble radius of the green solid and dashed lines (5 Myr) of Fig. A.1, which is caused by a transition from the irradiative to the viscous regime, marked by the green triangle at 0.1 au. This discontinuity follows directly from the analytical expressions, which contain a change in parameter dependencies and slopes at the viscous-irradiative boundary.

When implementing this change in parameters in the analytical expressions for the pebble size of Ida et al. (2016), the pebble radius increases once the pebbles cross r_visc–irr. However, τ_s,crit1 inside r_visc–irr remains smaller than τ_s just outside the boundary, meaning that sudden in-situ growth is not expected. Much is still unclear about what happens with R_peb in the Stokes regime within the snowline, and further studies with dust coagulation models within the snowline are required to create a more exact model. For simplicity, we assume the analytical expressions from Ida et al. (2016) remain valid within the snowline and within the viscous-irradiative boundary, albeit with a different value for ρ_s and ${\dot{M}}_{\mathrm{F}}$. We acknowledge, however, that a model in which R_peb can only grow if τ_s(R_peb) < τ_s,crit1, might be more self-consistent.

Fig. A.1 Influence of sublimation and fragmentation on the pebble radius (top) and Stokes number (bottom) within the snowline of a 1.00 M⊙ star. The solid lines represent the sublimation model, in which pebbles fragment and immediately recoagulate with a solid-dominated internal density. The dotted lines represent the fragmentation model without recoagulation. The dashed lines show the default model, in which the pebbles’ size and composition remain unaltered as they cross the snowline. The colours of the lines signify different times of evaluation. The circles, triangles, and squares represent crossings of the Epstein-Stokes boundary, the viscous-irradiative boundary and the snowline respectively. At the snowline, Rpeb instantaneously decreases due to sublimation, and even more so due to fragmentation, which leads to a local reduction in τs.

Fig. A.2 Influence of sublimation and fragmentation on the accretion efficiency from OL18 (top) and IGM16 (bottom) within the snowline of a 1.00 M⊙ star. The accretion efficiencies are calculated for a 10−3 ME planetesimal. The line styles match those in Fig. A.1. Fragmentation causes a sudden decrease in accretion efficiency at the snowline, especially after 0.01 and 1.0 Myr, due to the τs becoming so small that pebbles couple to the gas. Without fragmentation, the efficiency rapidly decreases in the inner regions of the disc (≲ 0.2 au), due to τs becoming too large for pebbles to settle. The sudden and sharp increase in ϵIGM16 in the innermost part of the disc is an unintended modelling effect from the accretion impact parameter B becoming smaller than the planet radius Rpl.

The influence of the different pebble radius models within the snowline on the accretion efficiency⁴ is shown in Fig. A.2. These results have been calculated for a 10⁻³ M_E planetesimal, which is a typical mass for planets in the early stages of rapid PA. Fragmentation at the snowline causes a rapid drop in accretion efficiency, especially at early times in the disc. This reduction might be even stronger if aerodynamic deflection of tiny pebbles coupled to the gas is considered (Visser & Ormel 2016). This deflection depends on the considered gas flow model around the planet, which is a level of detail that is beyond the scope of this study.

Sublimation, on the other hand, causes an increase in accretion efficiency compared to the non-sublimation model, due to the high Stokes number being slightly reduced, allowing for the pebble to be slowed down slightly faster. Nevertheless, in the innermost regions of the disc, the Stokes number becomes so large that pebbles are no longer efficiently slowed down and do not settle in the planet’s gravitational field. This causes the rapid drop in accretion efficiency in both OL18 and IGM16.

In IGM16, the efficiency reduction for τ_s ≫ 1 is explicitly modelled in the accretion cross-section B through the exponential reduction factor κ from Ormel & Kobayashi (2012) (see the parameters of Eq. (12)). However, this approach has a problem of its own. Unlike OL18, the version of IGM16 used in this study does not include expressions for the ballistic regime, which occurs when the accretion impact parameter B, which is to say the largest impact parameter for which pebbles can accrete, becomes less than the geometric limit R_pl. In this regime, pebbles no longer accrete because they settle in the gravitational field at a distance B away from the planet, but because their trajectory directly intersects the planet’s surface. The planet’s surface thus becomes the impact parameter B for accretion. As τ_s rapidly increases in the innermost regions of the disc, κ drastically decreases, so that B ⋘ R_pl. In our version of IGM16, the cross-section cannot decrease below R_pl. However, using the geometric limit in the settling efficiency equations leads to a serious overestimation of the accretion efficiency, as can be seen from the sharp increase in ϵ_IGM16 in Fig. A.2. This is because ϵ_IGM16 scales approximately with $b^3\equiv B^3{r}^{-3}=R_{\text{pl\hspace{0.17em}}}^3{r}^{-3}$, the slope of which matches exactly with what is observed in the figure.

Allowing B to asymptotically decrease to 0 due to κ and ignoring the ballistic regime altogether, has less significant consequences on the growth rate of planets in these inner regions of the disc, since, as can be seen from the OL18 results, growth in the ballistic regime is small.

Either way, this shortcoming of IGM16 has little to no effect on the full SyMBA simulations, since it applies to a region far closer to the star than the planetesimal disc. Planets can only reach these regions through migration, and planets generally reach the pebble isolation mass in a much shorter timescale than they migrate.

The final caveat of the sublimation model is that in the inner regions of the disc, τ_s ≫ 1. As a result, the drift slows down significantly, since $$t_{\mathrm{d}\mathrm{r}\mathrm{i}\mathrm{f}\mathrm{t}}\propto \frac{1+{\tau }_{\mathrm{s}}^2}{{\tau }_{\mathrm{s}}}\approx {\tau }_{\mathrm{s}},\text{\hspace{0.17em}for\hspace{0.17em}}{\tau }_{\mathrm{s}}>1.$$

As τ_s continues to increase, the drift timescale again becomes longer than the growth timescale. This occurs when τ_s exceeds (Ida et al. 2016) $$\begin{array}{rl}{\tau }_{\mathrm{s},\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}2}& \simeq {\left(\frac{3\sqrt{3\pi }}{16|3-0.5\gamma |}\frac{{\alpha }_{\mathrm{a}\mathrm{c}\mathrm{c}}}{\eta }\frac{{\dot{M}}_{\mathrm{F}}}{{\dot{M}}_{\ast }}\right)}^{-1/2}\\ & \simeq 13{\left({\alpha }_3^{-1}{\dot{M}}_{\ast 8}{\dot{M}}_{\mathrm{F}4}^{-1}\right)}^{1/2}{\left(\frac{\eta }{{10}^{-3}}\right)}^{1/2}\propto r^q.\end{array}$$(A.1)

In theory, this could result in runaway coagulation (Okuzumi et al. 2012), since τ_s ∝ R². This means that once τ_s > τ_s,crit2, the condition will always be satisfied and pebble drift will come to a complete halt.

Figure A.3 shows that for the sublimation model around a solar-mass star, τ_s indeed exceeds τ_crit,2, suggesting the pebbles could enter this phase of runaway coagulation. However, when comparing the growth timescale to the drift timescale, as is done in Fig. A.4, it turns out that in the region where t_grow < t_drift, t_drift is of the order of a few years, and t_grow is only a few orders shorter. Moreover, when the compact silicate pebbles within the snowline come together, they do not perfectly stick to one another (Morbidelli et al. 2015; Ida et al. 2016), which significantly limits the growth rate, meaning that the growth timescale is underestimated. Therefore, it is probable that pebbles drift into the central star before they have time to significantly grow. For smaller stars, τ_s always remains within the confines of τ_crit.

Some might argue that the fragmentation model should be preferred, since these small pebbles never exceed τ_s,crit2. However, these small pebbles have Stokes numbers significantly smaller than τ_s,crit1, meaning they are also prone to rapid growth. In fact, pebbles in the fragmentation model stay at the snowline for tens of thousands of years due to their long drift timescale. Meanwhile, their growth timescale is of the order of a year or less, as can be seen in Fig. A.4. We believe that it is, therefore, far more likely that the mm-sized silicate grains of the fragmentation model would start growing, than that the pebbles in the sublimation model would. Moreover, in the fragmentation model, the growth of the pebbles would be around the snowline, and would therefore significantly impact the growth rate of the planetesimals. The potential runaway growth of pebbles in the sublimation model would occur far closer to the star, where, as argued before, they could only encounter planets that have already reached the pebble isolation mass.

We therefore favour the sublimation model over the fragmentation model, and use this model for the full SyMBA simulations. However, detailed dust coagulation computations using the disc models of Ida et al. (2016) are required to further validate the assumptions mentioned above. This is beyond the scope of this study.

Fig. A.3 Comparison between τcrit and τs in both the sublimation and fragmentation model for a 1.00 M⊙ star. In the sublimation model, τs > τcrit,2 for r ≲ 0.2 au, which could lead to runaway coagulation. As the disc evolves, τs decreases, and remains within the bounds up to smaller orbital radii. In the fragmentation model, τs ≪ τcrit,1, which suggests these pebbles no longer drift and should therefore grow in situ as well. The values for τcrit,1 and τcrit,2 were calculated by numerically solving tgrow = tdrift.

Fig. A.4 Comparison between tdrift and tgrow in both the sublimation and fragmentation model for a 1.00 M⊙ star. For both models, tgrow < tdrift in the inner disc, suggesting (runaway) coagulation could occur. However, while the fragmented pebbles remain stuck around the snowline for tens of thousands of years, four orders of magnitude longer than their growth timescale, pebbles in the sublimation model drift into the star in only a few years, giving them significantly less time to grow. It is therefore unlikely that runaway coagulation would be a serious problem in the sublimation model.

Appendix B Supplementary figures

This appendix presents additional figures showing more data or alternative visualisations complementary to Sect. 5.

Fig. B.1 The dynamical evolution tracks (Mp and a) of all large planets in IGM16 simulation 1 around a 0.49 M⊙ star. The coloured solid lines represent different planets. The transparent dashed lines represent large planetesimals that merged with the planets. Only planets with masses > 0.01 ME are included. The cyan shaded region indicates the habitable zone. In the 0.49 M⊙ simulations, the planets that form at a late stage in the disc (t ≳ 1 Myr) form outside the HZ and remain too small to migrate into it.

Fig. B.2 Overview of the mass, eccentricity, and inclination of all planets formed in the OL18-Ring and IGM16-Ring simulations. Earth-like planets have been highlighted using a black edge around the marker. In the e and i plots, the size of the marker is proportional to its mass. Overall, the results are very similar to those of the normal runs, presented in Fig. 11.

References

All Tables

All Figures

Fig. 1 Radial pebble mass flux (solid lines, left-hand axis) and cumulative pebble mass (dashed lines, right-hand axis) in the inner disc as a function of time for the five stellar masses used in this study. After around 0.1 Myr, the pebble formation front reaches the outer edge of the disc, after which the slope of the pebble mass flux becomes significantly more negative. The disc around the 0.09 M⊙ star is weighted as 10% of the stellar mass, while all others are weighted as 5% (see the discussions in Sects. 3.1 and 4).

In the text

Fig. 2 Orbit-averaged pebble accretion rate as a function of eccentricity (left panels) and inclination (right panels) for planets around a 1.0 M⊙ star, 0.01 Myr after the formation of the disc. The top plots show the accretion rate for a 10−3 ME planetesimal at different orbital radii, while the bottom plots shows results for planets of different masses at 1 au. The solid lines were calculated with OL18, the dashed lines with IGM16. The OL18 prescription contains explicit expressions for the influence of the eccentricity and inclination on the accretion rate, while the IGM16 model only includes corrections for part of the influence of the inclination. For large e, the change in rp along a single orbit becomes significant enough for the variations in the encountered disc conditions to influence the accretion rate, both in the OL18 and the IGM16 model. The y-axis of the top-left plot changes from a logarithmic to a linear scale at 1 × 10−7 ME/yr to highlight the effect.

In the text

Fig. 3 Final mass divided by the initial mass of planets, calculated using only the semi-analytical PA model OL18, shown as a function of orbital radius on the x-axis and initial mass (M0) on the y-axis. The eccentricity and inclination are zero. Black contours indicate final masses higher than a Mars mass (MM), Earth mass (ME) and two Earth masses. Planets above the red dashed line have reached the local pebble isolation mass. Around 0.49, 0.70, and 1.00 M⊙ stars, PA produces Earth-mass planets in a wide range of orbital separations. Around 0.09 and 0.20 M⊙, PA does not form planets larger than a few Mars masses, even though the pebble isolation mass allows for the formation of Earth-sized planets. For 0.09 M⊙ stars, planets do not grow larger than about a lunar mass (ML) in the default disc of 5% of the stellar mass, which is why a more massive disc of 10% of the stellar mass is used for this star.

In the text

	Fig. 4 Same plots as Fig. 3, but calculated using the semi-analytical PA model IGM16. Pebble accretion using IGM16 is more efficient than using OL18, even for the circular, uninclined orbits assumed here. This is most clearly visible from the contours in the 0.09 and 0.20 M⊙ subplots.
In the text

Fig. 5 Time after which planets reach the pebble isolation mass tiso ≡ t(M=Miso), as a function of orbital radius and initial mass, using IGM16. White indicates the isolation mass has not been reached. The contours showing the final mass are the same as in Fig. 4. Around 0.49, 0.70, and 1.00 M⊙ stars, the first planets reach Miso within the first 50 000 to 100 000 years of the simulation, which could halt pebble accretion for planets close to the star prematurely.

In the text

Fig. 6 Simulated planetary systems around a 1.00 M⊙ star for the four PA models. The size of the markers represents the mass of the planet. The simulations are grouped by their model, and each horizontal line represents a simulation. The vertical dashed line at ~10−1 au is the inner radius of the gas disc, within which the gas density rapidly drops and type I and type II migration halts. Any planet that enters the gray region, interior to the inner truncation radius, is removed from the simulation. The cyan shaded regions represent the conservative (darker shaded) and optimistic (lighter shaded) habitable zone. The colour of the planets indicates their AMD stability, discussed in Sect. 5.6. Only planets more massive than Mars are shown.

In the text

Fig. 7 Dynamical evolution tracks (Mp, a, e and i) of all large planets in OL18 simulation 7 around a solar-mass star. The different coloured solid lines represent different planets, the red one being the planet that is currently in the habitable zone (cyan shaded region). The transparent dashed lines represent large planetesimals that merged with the planets. Sudden stepwise increases in planet masses signify these mergers. Only planets with masses >0.1 ME are included; however, for this simulation, there are no other objects with masses of ≳10−3 ME.

In the text

Fig. 8 Dynamical evolution tracks (Mp and a) of all large planets in IGM16 simulation 1 around a solar-mass star. The different coloured solid lines represent different planets. The transparent dashed lines represent large planetesimals that merged with the planets. Sudden stepwise increases in planet masses signify these mergers. Only planets with masses >0.1 ME are included. The cyan shaded region represents the habitable zone. The thick black ellipses in the lefthand panel highlight the first, second, and third generation of planets. IGM16 produces far more planets than OL18. Those formed as part of the third generation, after 1–2 Myr, remain in the HZ because they have insufficient time for migration.

In the text

	Fig. 9 Simulated planetary systems around the 0.49 (left) and 0.70 (right) M⊙ stars for the four PA models. The size of the markers represents the mass of the planet, while their colour represents the planet’s AMD stability (see Sect. 5.6). The layout of the figure is explained in the caption of Fig. 6.
In the text

	Fig. 10 Simulated planetary systems around the 0.09 (left) and 0.20 (right) M⊙ stars for the four PA models. The size of the markers represents the mass of the planet, while their colour represents the planet’s AMD stability (see Sect. 5.6). The layout of the figure is explained in the caption of Fig. 6.
In the text

	Fig. 11 Overview of the mass, eccentricity, and inclination of all planets formed in the OL18 and IGM16 simulations. Earth-like planets have been highlighted using a black edge around the marker. In the e and i plots, the size of the marker is proportional to its mass.
In the text

Fig. 12 Core mean quantities of all simulations, grouped by stellar mass and model. All quantities are calculated per simulation, and then averaged over the 8 simulations per group. (A) Mean total mass of the planetesimals and planets in the disc. (B) Mean maximum planet mass. (C) Mean number of Earth-like (0.67–1.5 ME) planets. The number of Mars-like planets (0.1–0.67 ME), and massive planets (>1.5 ME) are stacked on top with increasing transparency.

In the text

	Fig. 13 Dynamic spacing Δ of all planet pairs in all simulations, plotted against their period ratio Pout/Pin. The planet pairs are grouped per PA model, with their colour indicating the stellar mass. The horizontal lines indicate two important values of the dynamic spacing: the critical value $\mathrm{\Delta }=2\sqrt{3}$ (solid), and the value Δ = 10 (dashed).
In the text

Fig. 14 Simulated planetary systems around the 0.70 and 1.00 M⊙ stars for the four PA models including gas accretion. The mass threshold for particles to become self-gravitating was increased compared to the simulations without gas accretion to conserve computing resources, meaning that these simulations are less accurate than the others. The purpose of this figure is therefore only to indicate whether gas accretion could significantly influence the results from the normal simulations. For the OL18 simulations, gas accretion does not play a significant part. In the IGM16 simulations, on the other hand, multiple gas giants form. The layout of the figure is explained in the caption of Fig. 6.

In the text

	Fig. 15 Cumulative distribution function of the maximum ratio between planet mass and eccentricity-dependent PIM for the planets in the different simulations. Only growth phases with non-zero pebble accretion rates for planets with a final mass above 0.01 ME are considered in the calculation of the maximum.
In the text

Fig. A.1 Influence of sublimation and fragmentation on the pebble radius (top) and Stokes number (bottom) within the snowline of a 1.00 M⊙ star. The solid lines represent the sublimation model, in which pebbles fragment and immediately recoagulate with a solid-dominated internal density. The dotted lines represent the fragmentation model without recoagulation. The dashed lines show the default model, in which the pebbles’ size and composition remain unaltered as they cross the snowline. The colours of the lines signify different times of evaluation. The circles, triangles, and squares represent crossings of the Epstein-Stokes boundary, the viscous-irradiative boundary and the snowline respectively. At the snowline, Rpeb instantaneously decreases due to sublimation, and even more so due to fragmentation, which leads to a local reduction in τs.

In the text

Fig. A.2 Influence of sublimation and fragmentation on the accretion efficiency from OL18 (top) and IGM16 (bottom) within the snowline of a 1.00 M⊙ star. The accretion efficiencies are calculated for a 10−3 ME planetesimal. The line styles match those in Fig. A.1. Fragmentation causes a sudden decrease in accretion efficiency at the snowline, especially after 0.01 and 1.0 Myr, due to the τs becoming so small that pebbles couple to the gas. Without fragmentation, the efficiency rapidly decreases in the inner regions of the disc (≲ 0.2 au), due to τs becoming too large for pebbles to settle. The sudden and sharp increase in ϵIGM16 in the innermost part of the disc is an unintended modelling effect from the accretion impact parameter B becoming smaller than the planet radius Rpl.

In the text

Fig. A.3 Comparison between τcrit and τs in both the sublimation and fragmentation model for a 1.00 M⊙ star. In the sublimation model, τs > τcrit,2 for r ≲ 0.2 au, which could lead to runaway coagulation. As the disc evolves, τs decreases, and remains within the bounds up to smaller orbital radii. In the fragmentation model, τs ≪ τcrit,1, which suggests these pebbles no longer drift and should therefore grow in situ as well. The values for τcrit,1 and τcrit,2 were calculated by numerically solving tgrow = tdrift.

In the text

Fig. A.4 Comparison between tdrift and tgrow in both the sublimation and fragmentation model for a 1.00 M⊙ star. For both models, tgrow < tdrift in the inner disc, suggesting (runaway) coagulation could occur. However, while the fragmented pebbles remain stuck around the snowline for tens of thousands of years, four orders of magnitude longer than their growth timescale, pebbles in the sublimation model drift into the star in only a few years, giving them significantly less time to grow. It is therefore unlikely that runaway coagulation would be a serious problem in the sublimation model.

In the text

Fig. B.1 The dynamical evolution tracks (Mp and a) of all large planets in IGM16 simulation 1 around a 0.49 M⊙ star. The coloured solid lines represent different planets. The transparent dashed lines represent large planetesimals that merged with the planets. Only planets with masses > 0.01 ME are included. The cyan shaded region indicates the habitable zone. In the 0.49 M⊙ simulations, the planets that form at a late stage in the disc (t ≳ 1 Myr) form outside the HZ and remain too small to migrate into it.

In the text

	Fig. B.2 Overview of the mass, eccentricity, and inclination of all planets formed in the OL18-Ring and IGM16-Ring simulations. Earth-like planets have been highlighted using a black edge around the marker. In the e and i plots, the size of the marker is proportional to its mass. Overall, the results are very similar to those of the normal runs, presented in Fig. 11.
In the text