© The Authors 2024

1 Introduction

Planet formation involves the growth from interstellar grains of sub-micron sizes to planets of thousands of kilometres in diameter, which is a process through at least 12 orders of magnitude in length scale. Details of the involved processes are still under ongoing research. Particularly, the formation of solid cores which subsequently accrete gas is a crucial yet still unclear step. This has been an active field of research for decades and requires further investigations.

Weidenschilling (1977) presented a classic problem in planet formation that, due to aerodynamic drag in protoplanetary discs, solids of 10 cm to 1 m in size typically have a radial drift timescale of ~100 yr, which is much shorter than the typical disc lifetime of 1–10 Myr. Furthermore, laboratory experiments of collisions (e.g. Wurm et al. 2005; Güttler et al. 2010) also show a general behaviour that millimetre-sized grains require extremely small relative velocities to grow, so that fragmentation and bouncing are avoided. These barriers of particle growth are often summarized as the ‘metre-size barrier’ in the literature. This implies that planetesimals of a kilometre in size have to form rapidly through the metre-sized scale from dust via an alternative process.

The Goldreich-Ward mechanism suggests the formation of planetesimals through gravitational collapse of a very dense dust disc as a result of dust settling (Goldreich & Ward 1973), where the dust disc needs to be ~10⁴ times thinner than the gas disc. However, Cuzzi et al. (1993) showed that this cannot occur in a protoplanetary disc. The dense dust disc at the midplane, along with the gas in it, rotates at the Keplerian velocity; however, the gas disc immediately above rotates at a sub-Keplerian velocity due to the radial pressure gradient. This results in a steep vertical velocity gradient at the dust-gas interface, which induces the Kelvin-Helmholtz instability, preventing the dust disc from settling and collapsing gravitationally.

However, settling a dust disc with a solid density comparable to the gas density is possible without triggering the Kelvin-Helmholtz instability. Analyses in multiple works (e.g. Youdin & Goodman 2005; Youdin & Lithwick 2007; Johansen et al. 2007, 2009; Bai & Stone 2010) suggest this can induce non-gravitational clumping of dusts due to disc turbulence or streaming instability. The over-densities of dust can subsequently collapse through gravity on an orbital timescale. Recent hydro-dynamic numerical simulations (e.g. Johansen et al. 2012, 2015; Simon et al. 2016, 2017) further show that dense filaments of solid particles undergo gravitational collapse and planetesimals up to about the size of Ceres are almost instantly formed. This process is a viable pathway for planetesimal formation.

The classical core accretion model of gas giant formation (Mizuno 1980; Pollack et al. 1996) requires a solid core of ~10 M_⊕. Beyond the critical mass, hydrostatic equilibrium in the gas envelope cannot be maintained, resulting in runaway gas accretion. The growth ends as the supply of gas is terminated due to gap opening in the disc or gas dispersal as the disc evolves.

Through N-body simulations, Kokubo & Ida (1998, 2000) showed that pairwise accretion of planetesimals results in runaway growth, where more massive bodies grow faster. As protoplanets grow massive enough to interact with each other gravitationally, their orbital separations remain larger than ~5 Hill radii and growth becomes oligarchic, where the growth rate is slower for more massive bodies. This results in a bimodal system of a few protoplanets and a population of small planetes-imals. Their extrapolation estimates that the growth timescale to reach 5–10 M_⊕ is of the order of 10–100 Myr beyond 5 au, which is much longer than the typical disc lifetime. Since a solid core of ~ 10 M_⊕ has to be formed before disc dispersal in order to accrete gas, a more efficient planetesimal growth mechanism is required.

Large populations of grains ranging from millimetres to tens of centimetres in radius, or pebbles, have been detected in pro-toplanetary discs by millimetre to centimetre observations (e.g. Testi et al. 2003; Wilner et al. 2005). These observations are consistent with the metre-size barrier mentioned above. The growth of these small particles is stalled and they remain throughout most of the lifetime of the discs (Cleeves et al. 2016). This lays the foundation for the notion of pebble accretion. In this scenario, a large population of pebbles, as leftover solids, co-exists with planetesimals, in contrast to the classical scenario where pebbles are neglected for the growth of planetesimals of the order of a kilometre and beyond. Planetesimals that are massive enough to gravitationally deflect pebbles from the gas streamline and have a long enough encounter time can accrete a significant fraction of the drifting pebbles. This emerges as a mechanism for efficient planetesimal growth commonly called ‘pebble accretion’ (Ormel & Klahr 2010; Lambrechts & Johansen 2012; Guillot et al. 2014; see Johansen & Lambrechts 2017; Ormel 2017, for review).

Kretke & Levison (2014) conducted a series of numerical simulations incorporating pebble accretion with an initial mass spectrum of ~10⁶ planetesimals. The Lagrangian Integrator for Planetary Accretion and Dynamics (LIPAD; Levison et al. 2012), an N-body code, was deployed, which utilizes statistical algorithms to follow a large number of particles represented by tracers. As a result of oligarchic growth, the simulations generally form hundreds of ~M_⊕ bodies at 4–10 au but further growth is halted due to gravitational scattering. The scattered oligarchs also pollute the inner Solar System with water and disrupt the outer Solar System.

To produce a Solar System analogue, the later work by Levison et al. (2015) modifies the pebble formation model that the pebble formation timescale is lengthened to ~1 Myr. This allows viscous stirring among planetesimals, which is on a shorter timescale compared to the growth timescale through pebble accretion. The less massive planetesimals are excited to orbit with higher inclinations. As the pebble density is lower farther away from the midplane of the disc, these inclined planetesimals are then starved of pebbles. This scenario yielded 1–4 planets at 5–15 au from the Sun without a stage of oligarchic growth. However, as noted in their work, gas accretion was cut off arbitrarily once the planet reaches the Jupiter mass M_J, instead of employing physical laws to stall the growth. Also, the embryos started to accrete gas in the simulations at around 8 Myr. The adopted gas accretion rate is likely unrealistically high as the disc has only ~4% of its initial surface density at this age in their model, which results in a generous gas accretion rate. Finally, planet migration, which puts a critical time constraint on planet formation, was not considered in the model either.

Matsumura et al. (2017), in turn, employed the Symplec-tic Massive Body Algorithm (SyMBA; Duncan et al. 1998), a direct N-body code, with modifications to include pebble accretion, planet migration and gas accretion. They explored the dependence on stellar metallicity, stellar accretion rate and the viscosity parameter of the disc. Without migration, 1–3 gas giants are formed at a few au in younger and less viscous discs. However, at the end of their 50 Myr simulations with migration, none of the results is consistent with the Solar System, as there are no giant planets left beyond 1 au. This shows that planet migration plays a crucial role in planet formation. Another major difference between the works by Levison et al. (2015) and Matsumura et al. (2017) is the number of particles simulated. Levison et al. (2015) use LIPAD, which simulates a large population of particles employing a statistical algorithm making viscous stirring among planetesimals possible. They also focused on growing gas giant analogous to the Solar System, and the domain of simulation is 4–15 au. In contrast, Matsumura et al. (2017) focus on the production of the observed exoplanetary systems, and the domain of simulation is 0.3–5 au instead.

More recently, Bitsch et al. (2019) adopt the slower migration prescription in the high-mass regime by Kanagawa et al. (2018). They employ the pebble and N-body code FLINT-STONE that also includes planet migration, eccentricity and inclination damping, as well as disc evolution. Their results show that with higher pebble mass flux and reduced planet migration rate, gas giants can indeed survive at wide orbits; with the final semimajor axes sensitive to the pebble mass flux and planet migration rate. Also, some of the resulting gas giants undergo scattering close to the Sun and end at a few au from the Sun. However, in these simulations, there are also other planets of a few to tens of M_⊕ that migrate into the inner disc with less than 1 au, in contrast to the Solar System. Similarly, Matsumura et al. (2021) is able to form cold giant planets but cannot simultaneously avoid massive planetary cores migrating into the inner Solar System.

These works incorporating pebble accretion into global N-body simulations show intriguing results that the formation of gas accreting cores is possible through pebble accretion. Yet, further investigations are required to produce results that are consistent with the Solar System. The present study aims at assembling the giant planets analogous to those in the Solar System. In contrast to previous N-body planet formation models (e.g. Matsumura et al. 2017, 2021; Bitsch et al. 2019) that focus on a small number of lunar-mass embryos, we assume an initial planetesimal disc with planetesimal sizes comparable to those formed via the gravitational collapse induced by streaming instability. This is made computationally possible by employing SyMBA parallelized (SyMBAp; Lau & Lee 2023), which is a parallelized version of SyMBA. In the following, Sect. 2 presents the methodology adopted in this work and the results are presented in Sect. 3. The discussion of the results, the implications and caveats are in Sect. 4.

2 Methods

We generally follow the model by Matsumura et al. (2017) where additional subroutines are coupled with the symplectic direct N-body algorithm SyMBA (Duncan et al. 1998) to study planet formation in a protoplanetary disc. To facilitate the integration of a self-gravitating planetesimal disc in this work, we instead employ SyMBAp (Lau & Lee 2023). Further improvements are also made on the models of pebble accretion, gas accretion and the transition to the high-mass regime of planet migration. The following includes a summary of various parts of the model and the modifications made in this work are described in detail.

Fig. 1 Time evolution of ${\dot{M}}_{*}$, with the initial age of the disc t0 = 0.5 Myr. The value of ${\dot{M}}_{*}$ is turned down linearly when it drops below 10−9 M⊙ yr−1 to mimic the effect of photoevaporation.

2.1 Disc model

We consider an axisymmetric protoplanetary disc around a Solar-type star of 1 M_⊙ in mass and 1 L_⊙ in luminosity undergoing steady gas accretion. The gas accretion rate can be expressed as $${\dot{M}}_{*}=3\pi {\text{Σ}}_{\text{g}}v$$(1)

with Σ_g the gas surface density. For the viscosity v, the Shakura & Sunyaev (1973) α-parametrization is adopted such that $$v={\alpha }_{\text{acc}}{c}_{\text{s}}{H}_{\text{g}}$$(2)

with the viscosity parameter α_acc = 10⁻³ set in this work. The isothermal sound speed is used and given by $c_{\text{s}}=\sqrt{k_{\text{B}}T/\mu }$ with the Boltzmann constant k_B, the disc midplane temperature T, the mean molecular weight of the gas µ = 2.34m_H, and the hydrogen mass m_H = 1.67 × 10⁻²⁷ kg. The gas disc scale height H_g is defined by H_g = c_s/Ω_ĸ, where the local Keplerian orbital frequency ${\text{Ω}}_{\text{K}}=\sqrt{G{M}_{*}/r^3}$ with the gravitational constant G, the mass of the central star M_*, and the distance from the star r. Following Hartmann et al. (1998), the evolution of the disc is propagated from the modulation of the stellar accretion rate by $$\mathrm{log\hspace{0.17em}}\left(\frac{{\dot{M}}_{*}}{{10}^{-8}{M}_{\odot }{\text{yr}}^{-1}}\right)=-1.4\mathrm{\hspace{0.17em}log\hspace{0.17em}}\left(\frac{t+t_0}{\text{Myr}}\right)$$(3)

with the time since the start of the simulation t and the initial age of the disc t₀ = 0.5 Myr. Figure 1 shows the time evolution of ${\dot{M}}_{*}$. When ${\dot{M}}_{*}$ drops below 10⁻⁹ M_⊙ yr⁻¹, ${\dot{M}}_{*}$ is linearly turned down to zero at t + t₀ = 5.5 Myr to mimic the effect of photoevaporation following Matsumura et al. (2017). With this setup, the initial stellar accretion rate is about 2.64 × 10⁻⁸ M_⊙ yr⁻¹ and reaches 10⁻⁹ M_⊙ yr⁻¹ when t ≈ 4.68 Myr.

In general, the inner part of the disc is dominated by viscous heating and the outer part is dominated by radiative heating. Since this work focuses on the formation of the giant planets in the Solar System, only radiative heating is considered for the disc, in contrast to the disc model in Matsumura et al. (2017, 2021) where viscous heating is also considered. The midplane temperature profile of the disc T is given by (Oka et al. 2011) $$T=150{\left(\frac{r}{\text{au}}\right)}^{-3/7} \text{K}.$$(4)

This setup yields the reduced disc scale height profile $${\widehat{h}}_{\text{g}}\equiv \frac{H_{\text{g}}}{r}\approx 0.024{\left(\frac{r}{\text{au}}\right)}^{2/7}.$$(5)

With Eq. (1) for the gas accretion rate, Eq. (2) for the α-parametrization, and Eq. (3) for the evolution of the stellar accretion rate, Eqs. (4) and (5) yield the gas surface density in the radiatively heated region $${\text{Σ}}_{\text{g}}\approx 2.7\times {10}^3{\left(\frac{{\alpha }_{\text{acc}}}{{10}^{-3}}\right)}^{-1}\frac{{\dot{M}}_{*}}{{10}^{-8}{M}_{\odot }{\text{yr}}^{-1}}{\left(\frac{r}{\text{au}}\right)}^{-15/14} \text{g} {\text{cm}}^{-2}$$(6)

This disc model yields a profile of the midplane pressure gradient parameter, where P is the midplane gas pressure, $$\begin{array}{cc}\eta & \equiv -\frac{{\widehat{h}}_{\text{g}}^2}{2}\frac{\partial \mathrm{\hspace{0.17em}ln\hspace{0.17em}}P}{\partial \mathrm{\hspace{0.17em}ln\hspace{0.17em}}r}\\ & \approx 8.02\times {10}^{-4}{\left(\frac{r}{\text{au}}\right)}^{4/7}.\end{array}$$(7)

2.2 Planetesimal disc

Instead of starting with lunar mass embryos as in Matsumura et al. (2017), a planetesimal disc is generated from 5 to 20 au initially with an initial mass function implemented in a manner similar to Lau et al. (2022) as summarized in the following. Plan-etesimals are drawn from the cumulative mass distribution in the work on planetesimal formation by Abod et al. (2019), which has the form of an exponentially truncated power law. The number fraction of planetesimals above mass m is given by $$\frac{N_{>m}}{N_{\text{ini}}}={\left(\frac{m}{m_{\min }}\right)}^{-0.3}\exp \left(\frac{m_{\min }-m}{0.3m_G}\right),$$(8)

for m ≥ m_min≥ with m_min being the minimum planetesimal mass considered, N_>m is the number of particles with a mass > m, N_ini is the initial number of particles, and m_G is a planetesimal gravitational mass. We have set m_min = 10⁻²m_G in this work, which is well below the peak of the distribution of the planetesimal mass in each logarithm mass bin as noted by Lau et al. (2022). The upper limit of m is also artificially set at 3m_G in the realization algorithm to avoid a mathematical singularity. This value is an order of magnitude larger than the characteristic mass of the initial mass function (0.3m_G), where Abod et al. (2019) also show that the maximum planetesimal mass is about an order of magnitude more massive than the characteristic mass. In this manner, only an insignificant number of massive planetesimals (~8 × 10⁻⁶N_ini) is lost. The form of the cumulative mass function is shown in Fig. 2.

For m_G, we adopt the critical mass for gravitational collapse of a dust clump in the presence of turbulent diffusion by Klahr & Schreiber (2020), which is given by $$\begin{array}{cc}{m}_{\text{G}}& =\frac{1}{9}{\left(\frac{\delta }{\text{St}}\right)}^{3/2}{\widehat{h}}_{\text{g}}^3{M}_{\odot }\\ & \approx 5.78\times {10}^{-4}{\left(\frac{\delta }{{10}^{-5}}\right)}^{3/2}{\left(\frac{\text{St}}{{10}^{-2}}\right)}^{-3/2}{\left(\frac{{\widehat{h}}_{\text{g}}}{0.038}\right)}^3{M}_{\oplus }\end{array}$$(9)

where δ is the small-scale diffusion parameter, which is independent of α_acc, and St is the Stokes number. In this work, we set δ = 10⁻⁵ and St = 10⁻² exclusively for planetesimal realization. While the strength of the small-scale diffusion is an active research topic in the field, the adopted value is motivated by the measurements of local diffusivity of dust particles in streaming instability presented in Schreiber & Klahr (2018).

In each simulation, the semimajor axis a of a new planetesimal is randomly drawn from 5 to 20 au, which implies a surface number density of planetesimals that scales with 1/r. The value of m_G is then evaluated with the local disc scale height. Afterwards, the mass m of this planetesimal is drawn from the mass function given by Eq. (8) with the chosen value of N_ini noted later in Sect. 2.6. Figure 3 shows the initial mass distributions of the realized planetesimal discs with one example shown for each of the chosen values of N_ini. The eccentricity e is randomly drawn from a Rayleigh distribution with the scale parameter 10^–6. The inclination i in radian is also drawn from a Rayleigh distribution but with the scale parameter 5 × 10^–7 instead. Other angles of the orbital elements are drawn randomly from 0 to 2π. The physical radius R_p is calculated by assuming an internal density ρ_s = 1.5 g cm^–3. The realization process repeats until the total number of planetesimals reaches the chosen value. The planetesimals are then evolved under full gravitational interactions between themselves and the central star, as well as additional effects of pebble accretion (Sect. 2.3), gas accretion (Sect. 2.4) and planet-disc interactions (Sect. 2.5).

Fig. 2 Adopted truncated power law initial planetesimal mass function as described by Eq. (8) based on Abod et al. (2019). It is presented in the unit of the planetesimal gravitational mass mG

Fig. 3 Initial mass distribution of the realized planetesimal discs. One example is shown for each of the chosen values of Nini. The width of each bin is 0.2 au.

Fig. 4 Time evolution of the pebble mass flux ${\dot{M}}_{\text{pf}}$ given by Eq. (11) with Z0 = 10–2 and αacc = 10–3 as set in this work.

2.3 Pebble accretion

We implement the ‘pebble formation front’ model (Lambrechts & Johansen 2014) to estimate the pebble mass flux ${\dot{m}}_{\text{peb}}$ As dust particles coagulate and grow into pebbles, their velocities are strongly influenced by the headwind. This causes a significantly inward drift of pebbles that provide a solid mass flux to the inner part of the disc. Since the dust growth timescale increases with radius in general, the source of the pebble mass flux, or the pebble formation front, evolves outwards in time. The location of the pebble formation front r_pſ is given by (Lambrechts & Johansen 2014) $$r_{\text{pf}}(t)={\left(\frac{3}{16}\right)}^{1/3}{\left(G{M}_{*}\right)}^{1/3}{\left({ϵ}_{\text{d}}{Z}_0\right)}^{2/3}{t}^{2/3}$$(10)

with the initial dust-to-gas ratio Z₀ and the particle growth parameter ∈_d = 0.05. The pebble mass flux ${\dot{M}}_{\text{pf}}$ is then calculated from the dust mass swept across by the pebble formation front per unit time, that is, $$\begin{array}{cc}{\dot{M}}_{\text{pf}}=& 2\pi r_{\text{pf}}{Z}_0{\text{Σ}}_{\text{g}}\left(r_{\text{pf}}\right){\dot{r}}_{\text{pf}}\\ \approx & \frac{{\dot{M}}_{\ast }}{{10}^{-8}{M}_{\odot }{\text{yr}}^{-1}}{\left(\frac{Z_0}{{10}^{-2}}\right)}^{5/3}{\left(\frac{{\alpha }_{\text{acc}}}{{10}^{-3}}\right)}^{-1}\\ & \times {\left(\frac{t+t_0}{\text{Myr}}\right)}^{-1/3}{10}^2{M}_{\oplus }{\text{Myr}}^{-1}.\end{array}$$(11)

A factor of$r_{\text{pf}}^{-1/14}$ is omitted for simplicity. We set Z₀ = 10^–2 in this work and Fig. 4 shows the time evolution of ${\dot{M}}_{\text{pf}}$ for the chosen parameters. We note that at 4.5 Myr, briefly before disc dispersal, r_pſ ≈ 350 au. This is comparable to the typical observed disc sizes, which is of the order of 100 au (e.g. Andrews et al. 2018; Long et al. 2018; Cieza et al. 2021). In Matsumura et al. (2017), ${\dot{M}}_{\text{pf}}$ is halved inside of the snow line. However, this treatment is not implemented in the present work as it focuses on the outer Solar System where particles are removed before they can reach the ice line in our model. The radial domain of this work is summarized later in Sect. 2.6. On the other hand, we follow Matsumura et al. (2021) and adopt the pebble disc scale height given by $$H_{\text{peb }}={\left(1+\frac{\text{St}}{{\alpha }_{\text{turb }}}\right)}^{-1/2}{H}_{\text{g}}$$(12)

with the Stokes number of pebble St. Following Ida et al. (2018), an α_turb parameter is introduced, which is about an order of magnitude smaller than α_acc as evaluated by Hasegawa et al. (2017). The latter is distinct from that in the classical α-parametrization, i.e. the α_acc parameter introduced in Sect. 2.1. In this work, we set α_turb/α_acc = 0.1. The α_turb parameter is also used for prescribing gas accretion (Sect. 2.4) and planet-disc interactions (Sect. 2.5) as described in the respective sub-sections.

Furthermore, the pebble flux available to each body is subtracted by the total pebble accretion rate of the superior bodies that are farther from the central star, if there are any. We define a pebble accretion efficiency ∈_PA such that the growth rate of a body i by pebble accretion is given by $${\dot{m}}_{\text{PA},i}={ϵ}_{\text{PA}}\max \left({\dot{M}}_{\text{pf}}-{\displaystyle \sum _{n=i+1}^N{\dot{m}}_{\text{PA},n}},0\right),$$(13)

where bodies (i + 1) to N are all the superior ones.

In this work, we also compare the pebble accretion efficiency of Ida et al. (2016) with modifications by Matsumura et al. (2021), ∈_IGM16, and that by Liu & Ormel (2018) and Ormel & Liu (2018), ∈_OL18. In the derivation of ∈_IGM16, the pebble-accreting body is assumed to be in a circular orbit as noted in Sect. 3.2 of Ida et al. (2016) and shown in Eq. (33) of their work regarding the pebble relative velocity. In contrast, Liu & Ormel (2018) and Ormel & Liu (2018) do not hold this assumption, and both the inclination and the eccentricity of the pebble-accreting body contribute to the pebble relative velocity. The modifications of ∈_IGM16 made by Matsumura et al. (2021) considered the inclination of the body. However, it only plays a role in the calculation of the pebble volume density as shown in Eq. (32) of their work but not in the calculation of the pebble relative velocity. The differences between the two pebble accretion prescriptions and the consequences are further discussed in Sect. 4.1.

When the planetesimals grow into massive cores, the process of pebble isolation occurs when they perturb the gas surface density profile and stop pebbles from reaching the planet itself as well as the inferior bodies that are closer to the central star, if there are any. We follow the assumption in Matsumura et al. (2017) that the required mass, which is often called the ‘pebble isolation mass’, is given by $$\begin{array}{cc}{m}_{\text{iso}}& =\frac{1}{2}{\widehat{h}}_{\text{g}}^3{M}_{\ast }\\ & \approx 9.14{\left(\frac{{\widehat{h}}_{\text{g}}}{0.038}\right)}^3{M}_{\oplus }\end{array}$$(14)

Once any planet reaches this mass, pebble accretion is stopped for this planet and all the inferior ones if there are any.

2.4 Gas accretion

When a massive core has formed and its solid accretion rate is low, gas can contract and form an envelope. We follow Ikoma et al. (2000) for the critical mass for runaway gas accretion, which is given by, for planet i, $$m_{\text{g},\text{crit}}=10{\left(\frac{{\dot{m}}_{\text{PA},i}}{{10}^{-6}{M}_{\oplus }{\text{yr}}^{-1}}\frac{\kappa }{1{\text{cm}}^2{\text{g}}^{-1}}\right)}^p{M}_{\oplus }.$$(15)

In this work, we set the parameter p = 0.25 (Ida & Lin 2004) and the envelope opacity k = 1 cm² g⁻¹. For cores that have reached this mass, we assume the gas envelope collapses on the Kelvin-Helmholtz timescale τ_KH given by (Ikoma et al. 2000; Ida & Lin 2004) $${\tau }_{\text{KH}}={10}^9{\left(\frac{m}{M_{\oplus }}\right)}^{-3}\left(\frac{\kappa }{1 {\text{cm}}^2 {\text{g}}^{-1}}\right)\text{yr}.$$(16)

There are two factors that limit the actual gas accretion rate considered in our model. First, the gas supply is limited by the stellar accretion rate as well as the gas accreted by the superior planets. Also, gap opening by the planet shall further limit the gas accretion rate. And, we assume gas accretion is exponentially cutoff when the planet’s Hill radius equals the local disc scale height, which is given by $m_{\text{Hill }}=3M_{*}{\widehat{h}}_{\text{g}}^3.$. These can be summarized as the expression for the gas accretion rate of planet i $$\begin{array}{l}{\dot{m}}_{\text{g},i}=\min [\frac{m}{{\tau }_{\text{KH}}},\\ \mathrm{max\hspace{0.17em}}\left({\dot{M}}_{*}-{\displaystyle \sum _{n=i+1}^N{\dot{m}}_{\text{g},n}},0\right)f_{\text{local }}\mathrm{exp\hspace{0.17em}}\left(-\frac{m}{m_{\text{Hill }}}\right)]\end{array}$$(17)

where planets (i + 1) to N are all the superior ones and the reduction factor f_local is given by (Ida et al. 2018) $$f_{\text{local }}=\frac{0.0308{\widehat{h}}_{\text{g}}^{-4}{\left(m/M_{*}\right)}^{4/3}{\alpha }_{\text{acc }}^{-1}}{1+0.04K}.$$(18)

The gap opening factor K is given by Eq. (24) in the next subsection (Sect. 2.5).

2.5 Planet-disc interactions

Other than the N-body gravitational interactions, the bodies also experience the torques due to the planet-disc interactions. We adopt the prescription based on dynamical friction by Ida et al. (2020) and the transition from the low-mass to the high-mass regime by Ida et al. (2018) based on the gap opening factor K by Kanagawa et al. (2015). The timescales for the non-isothermal case and finite inclination i, while $i<{\widehat{h}}_{\text{g}}$, (Appendix C and D of Ida et al. 2020 and Matsumura et al. 2021) are implemented. The evolution timescales of semimajor axis, eccentricity and inclination are defined respectively by $${\tau }_a\equiv -\frac{a}{ \text{d}a/\text{d}t},{\tau }_e\equiv -\frac{e}{ \text{d}e/\text{d}t},{\tau }_i\equiv -\frac{i}{ \text{d}i/\text{d}t}.$$(19)

These timescales are given by $$\begin{array}{l}{\tau }_a=\frac{t_{\text{wav}}^{\prime }}{2{\widehat{h}}_{\text{g}}^2}[\frac{{\text{Γ}}_L}{{\text{Γ}}_0}{\left(1-\frac{1}{C_{\text{M}}}\frac{{\text{Γ}}_L}{{\text{Γ}}_0}\sqrt{{\widehat{e}}^2+{\widehat{i}}^2}\right)}^{-1}\\ {+\frac{{\text{Γ}}_C}{{\text{Γ}}_0}\mathrm{exp\hspace{0.17em}}\left(-\frac{\sqrt{{\widehat{e}}^2+{\widehat{i}}^2}}{e_{\text{f}}}\right)]}^{-1},\end{array}$$(20) $${\tau }_e=1.282t_{\text{wav}}^{\prime }\left[1+\frac{{\left({\widehat{e}}^2+{\widehat{i}}^2\right)}^{3/2}}{15}\right],$$(21) $${\tau }_i=1.838t_{\text{wav}}^{\prime }\left[1+\frac{{\left({\widehat{e}}^2+{\widehat{i}}^2\right)}^{3/2}}{21.5}\right],$$(22)

where $\widehat{e}\equiv e/{\widehat{h}}_{\text{g}},\widehat{i}\equiv i/{\widehat{h}}_{\text{g}}$, and we follow Fendyke & Nelson (2014) for the factor $e_{\text{f}}=0.01+{\widehat{h}}_{\text{g}}/2$. The normalized Lindblad torque Γ_L/Γ₀ and corotation torque Γ_C/Γ₀ are described in detail by Paardekooper et al. (2011). The characteristic time including the transition to the high-mass regime $t_{\text{wav}}^{\prime }$ (Tanaka et al. 2002; Ida et al. 2018) is given by $$t_{\text{wav }}^{\prime }=\left(\frac{M_{*}}{m}\right)\left(\frac{M_{*}}{{\text{Σ}}_{\text{g}}{r}^2}\right)\left(\frac{{\widehat{h}}_{\text{g}}^4}{{\text{Ω}}_{\text{K}}}\right)(1+0.04K)$$(23)

with the gap opening factor K given by $$K={\left(\frac{m}{M_{*}}\right)}^2{\widehat{h}}_{\text{g}}^{-5}{\alpha }_{\text{turb }}^{-1}.$$(24)

As noted in Lau et al. (2022), it is more suitable to evaluate the value of Ω_K at the instantaneous distance from the star r of the body instead of its semimajor axis a in N-body simulations with large number of particles due to potential frequent encounters. We follow Ida et al. (2018) and introduce the α_turb parameter set to α_turb/α_acc = 0.1 as described in Sect. 2.3. The three timescales are applied to the equation of motion $$\mathit{a}=-\frac{v_{\text{K}}\cdot S_a}{2{\tau }_a}{\mathit{e}}_{\theta }-\frac{v_r}{{\tau }_e}{\mathit{e}}_r-\frac{v_{\theta }-v_{\text{K}}}{{\tau }_e}{\mathit{e}}_{\theta }-\frac{v_z}{{\tau }_i}{\mathit{e}}_z$$(25)

in the cylindrical coordinates (r, θ, z) with the velocity of the embryo υ = (υ_r, υ_θ, υ_z) and the local Keplerian velocity υ_K = rΩ_K. A switch for planet migration S_a is introduced to toggle the evolution of the semimajor axis, which is turned off and on respectively by setting S_a to 0 and 1 in this work.

2.6 Numerical setups

To explore the dependence on the total number of planetes-imals, three values of N_ini = {1000, 2000, 5000} are chosen. They translate respectively to a total planetesimal mass of about {0.02, 0.04, 0.1} M_⊕. We test two pebble accretion efficiency prescriptions ∈_PA = {∈_IGM16, ∈_OL18} described in Sect. 2.3 and the two states of S_a = {0,1} described in Sect. 2.5 that switches off or on the evolution of semimajor axis due to planet-disc interactions. Each simulation lasts for 6.5 Myr to allow for further dynamical evolution due to gravitational interactions after disc dispersal. Particles are removed if the heliocentric distance is less than 1 au or greater than 100 au. For each combination of the parameters, we conduct five simulations to sample the stochastic variations in the outcome. Thus a total of 60 simulations are conducted in this work and presented in the next section.

3 Results

The first part of this section (Sect. 3.1) presents the results with migration turned off, i.e. S_a = 0, followed by Sect. 3.2 where the results with migration turned on, i.e. S_a = 1, is presented.

3.1 Simulations without planet migration (S_a = 0)

3.1.1 Pebble accretion efficiency ∈_PA = ∈_IGM16

Figure 5 shows the results for N_ini = 1000, S_a = 0 and ∈_PA = ∈_IGM16. Each row presents a snapshot of the simulations at t = {0.10, 0.75, 2.50, 4.00, 6.50} Myr respectively. For the first three columns from the left, the total occurrences of particles across all five simulations are shown by heat maps. The left-most column shows the mass m in M_⊕ and the semimajor axis a. The next two columns to the right show the eccentricity e and inclination i against m respectively. The right-most column shows the differential mass distribution of the particles with each colour corresponds to one of the five simulations. Particles in one of the five simulations (blue) is also plotted with particles above 10⁻³ M_⊕ denoted by enlarged dots. For the last row (6.5 Myr), which shows the end results, particles above 10⁻³ M_⊕ in all simulations are shown individually (with a different colour for each simulation) without using heat maps.

The m–a plots show a rapid growth by pebble accretion in the inner part of the disc in the first 0.1 Myr of the simulations. Some planetesimals in the massive tail of the distribution have grown by more than 3 orders of magnitude dominantly by pebble accretion. The growth rate has a strong dependence on the distance from the star, and particles closer to the central star accrete pebble much faster, as predicted by Ida et al. (2016). This is also consistent with the analysis which includes both pebble and planetesimal accretion in Coleman (2021), though our simulations focus on the outer Solar System.

The e–m plots and the i–m plots show the early and fast growing bodies quickly heat up their neighbouring planetesimals from the beginning of the simulations to 0.75 Myr, increasing the eccentricities and inclinations of neighbouring planetesimals. The massive cores of ~M_⊕ stop further growth of the neighbouring smaller bodies by viscous stirring, with about 20 bodies having reached ~1–10 M_⊕ by 0.75 Myr. This effect of viscous stirring on pebble accretion is consistent with Levison et al. (2015) and further discussed in Sect. 4.1. The e and i of these cores are also damped and remain low in contrast to those of the smaller bodies, which allows these massive bodies to further increase in mass due to the proximity to the dense pebble disc. This effect is more noticeable from the differential mass distributions, i.e. the rightmost column, that only the particles in the massive tail of the initial planetesimal population can grow significantly while the rest remain about the same mass. The growth of these massive bodies is drastically different from the traditional oligarchic growth scenario, where the growth is slowed down by viscous heating that clears nearby planetesimals. Here, the more massive bodies can continue growth via pebble accretion until reaching the pebble isolation mass, which is a result of the perturbations to the gas disc.

As the simulations progress forward, the massive cores grow further by gas accretion and eject most of the small bodies from 0.75 to 4 Myr. At the end of the simulations, i.e. t = 6.50 Myr, some of the massive cores and gas giants (m > 10² M_⊕) formed have been ejected, and 1–4 gas giants remain but their locations vary greatly across the simulations. This indicates a strong stochastic behaviour due to dynamical instabilities that result from the formation of multiple gas giants in a short range of distance from the star. Also, ice giants (m ~ 10 M_⊕) do not survive in any of these simulations: they either became gas giants or were scattered out of the system by other giants. On the other hand, the results with N_ini = 2000 and 5000 do not show any qualitative difference from the presented results with N_ini = 1000 (Fig. 6).

Fig. 5 Results for the simulations with Nini = 1000, migration turned off (Sa = 0), and the pebble accretion efficiency ∈PA = ∈IGM16. Each row presents a snapshot of the simulations at the time indicated by the timestamp on the left. For the first three columns from the left, the total occurrences of particles across all five random simulations are shown by heat maps with 2 × 2 cells in each minor axis grid cell. The left-most column shows m and a. The next two columns to the right respectively shows the e and i, respectively, against a. The right-most column shows the differential mass distribution of all bodies with each colour corresponds to one of the five simulations. Particles in one of the five simulations (blue) is also plotted with enlarged dots denoting particles above 10−3 M⊕. For 6.5 Myr, i.e. the end of the simulations, particles above 10−3 M⊕ in all five simulations are shown individually without using heat maps. Further descriptions are in the text.

Fig. 6 End results for the simulations with Nini = 2000 and 5000, respectively, as indicated on the left, migration turned off (Sa = 0), and the pebble accretion efficiency ∈PA = ∈IGM16. The two columns here correspond to the left-most and the right most columns of Fig. 5, respectively. There is no qualitative difference in the end results among the simulations with the chosen set of Nini = {1000, 2000, 5000}.

3.1.2 Pebble accretion efficiency ∈_PA = ∈_OL18

Figure A.1 shows the results for N_ini = 1000, S_a = 0 and ∈_PA = ∈_OL18. Compared to the results for ∈_PA = ∈_IGM16, the growth by pebble accretion is generally slower, but still rapid. Some plan-etesimals grow by up to about 2 orders of magnitude in mass in the first 0.1 Myr and massive cores (m ~ M_⊕) are formed at 0.75 Myr. At 2.5 Myr, the massive cores in the inner part of the disc (~5–10 au) have reached the local pebble isolation mass and gas accretion begins with less than ~10 bodies having gained mass between the ~10 M_⊕ cores and the initial planetesimals. In the previous simulations (Fig. 5), this stage is reached at 0.75 Myr. This delay is caused by the change in the adopted pebble accretion efficiency ∈_PA, where∈_IGM16 is more efficient than ∈_OL18 as also shown in Matsumura et al. (2021). A comparison between the two efficiency prescriptions and the consequences are further discussed in Sect. 4.1. A more distinct dichotomy in mass is produced with ∈_PA = ∈_OL18 as shown by comparing the differential mass distribution in Fig. 5 for 0.75 Myr and that in Fig. A.1 for 2.50 Myr. A more significant number of planetesimals has reached ~ 10^-3 M_⊕ in the former case while a sharper cut near the upper end (~10^–4 M_⊕) of the initial distribution is shown in the latter case. At this stage, the intermediate-mass bodies between these two groups, which have mass of about 10^–5 – 10⁻¹ M_⊕, are generally dynamically colder, as shown by the e–m and i–m plots. As the simulations continue to 4.00 Myr, some bodies have become gas giants in the inner part of the disc, with some bodies of ~1–10 M_⊕ residing outside of 10 au, in contrast to the results shown in Fig. 5 at the same time.

At the end of the simulations, one to two gas giants and one to two ice giants are formed as well, which is the closest set of simulations in the work to reproduce the Solar System’s giant planets. A significant number of the initial planetesimals remain, especially in the outer part of the disc at around 20 au. This is distinct from the results with ∈_PA = ∈_IGM16, where no ice giants are formed and most of the initial planetesimals have been scattered at the end of the simulations, probably due to the higher number of gas giants. Nonetheless, the locations of the leftover bodies still vary greatly across the simulations, so that the stochastic nature of the system remains.

Figure A.2 shows the results for N_ini = 2000 instead with the same pebble efficiency prescription. Compared to Fig. A.1, the differential mass distribution shows that the massive tail extends for about twice as high in m. This leads to the formation of more massive cores in the subsequent evolution of the simulations. At the end of the simulations, more gas giants and fewer ice giants are formed in this case. Only two out of the five simulations has one to four ice giants, while this class of bodies is absent in the rest of the simulations. With N_ini = 5000, shown in Fig. A.3, only one simulation contains an ice giant at the end, which instead is located in the inner part of the disc at about 6 au. Here, we find a dependence on the value of N_ini, which is not present when ∈_PA = ∈_IGM16 (Sect. 3.1.1). This is likely caused by the difference in the rate of pebble accretion, which is further discussed in Sect. 4.1.

3.2 Simulations with planet migration (S_a = 1)

Figure A.4 shows the results for N_ini = 1000, with migration S_a = 1 and ∈_PA = ∈_IGM16. The snapshots of the m-a distribution show that once the cores reach ~M_⊕, they migrate inwards rapidly, even though α_turb/α_acc = 0.1. For the massive cores that grow from planetesimals in the inner part of the disc, they have moved out of the simulation domain before runaway gas accretion occurs. For the massive cores that remain by the end of the simulations, the depletion of the gas disc stops both the migration as well as gas accretion. As a result, only cores of a few M_⊕ are formed and survive in the simulations. A large fraction of the initial planetesimal population remains at the end as they are not scattered due to the absence of giant planets. Similarly, Fig. A.5 shows the results for ∈_PA = ∈_OL18 with S_a = 1 where only cores of a few M_⊕ are formed and survive. These cores are slightly less massive in this case compared to Fig. A.4. The results with N_ini = 2000 and 5000 do not show any qualitative difference from this results with migration in effect. Since the massive cores migrate rapidly and none reach the runaway gas accretion phase by the end of the simulation, the dependence on N_ini shown in the case without planet migration for ∈_PA = ∈_OL18 (Sect. 3.1.2) is no longer present in this case here.

4 Discussion

4.1 Pebble accretion efficiency

Pebble accretion has been shown by the results of our model (Sect. 3) to be a promising way to grow planetesimals efficiently such that massive cores of ~10 M_⊕ can form well before disc dispersal and accrete gas to become giant planets. Nonetheless, forming giant planets analogous to those in the Solar System still requires further modifications to the model. In the presented results without planet migration (S_a = 0), ice giants are formed only in the simulations with the pebble accretion efficiency prescription by Liu & Ormel (2018) and Ormel & Liu (2018), i.e. ∈_PA = ∈_OL18, as presented in Sect. 3.1.2. The ice giants in these simulations stop accreting gas because by the time they are massive enough to accrete a gaseous envelope the gas disc is dispersed. In contrast with ∈_PA = ∈_IGM16, as shown in Sect. 3.1.1, massive cores of ~10 M_⊕ are formed much earlier and the giant planets have enough time to reach the prescribed final mass (Sect. 2.4) before disc dispersal. This shows that the timing of the formation of the massive cores and the start of gas accretion plays an important role in the final architecture of the planetary system.

As noted by Matsumura et al. (2021), ∈_OL18 is generally a few times less efficient than ∈_IGM16 for the adopted value of α_turb. And, in the present work, the simulations begin with a mass spectrum of planetesimals which spans over two decades in mass, up to 10⁻⁴ M_⊕, instead of lunar-mass embryos. This demonstrates the effect of the pebble accretion onset mass and the effect of viscous stirring on pebble accretion efficiency more clearly as discussed in the following.

4.1.1 Pebble accretion onset mass

First, we focus on the limit that the eccentricity e of the pebble-accreting body is much lower than the midplane pressure gradient parameter η ~ 10⁻³. This is also an assumption held by Ida et al. (2016) in the derivation of the pebble accretion efficiency. Since we are considering the start of pebble accretion, the mass of the body is generally small and pebble accretion typically operates in the Bondi regime. In this case, the pebble relative velocity is determined by the headwind. For a high pebble relative velocity, the pebble encounter time is shortened so that pebbles may not be deflected enough from the gas streamline and not have enough time to settle onto the planetesimal. As such, the accretion is no longer in the settling regime. This reduction effect is captured in the pebble accretion efficiency prescription by Ida et al. (2016) as well as that by Liu & Ormel (2018) and Ormel & Liu (2018) but in slightly different manners.

Ida et al. (2016) adopt the reduction factor for the cross section in the settling regime of pebble accretion proposed by Ormel & Kobayashi (2012). This reduction factor is given by $${\kappa }_{\text{IGM}16}=\exp \left[-{\left(\frac{\text{St}}{\min \left(2,{\text{St}}_{\text{crit }}\right)}\right)}^{0.65}\right]$$(26)

with the critical Stokes number of pebble $${\text{St}}_{\text{crit }}=\frac{4m}{{\eta }^3{M}_{*}}.$$(27)

A similar reduction factor is also found in Liu & Ormel (2018), which is given by $${\kappa }_{\text{OL}18}=\exp \left[-\frac{1}{2}{\left(\frac{\text{Δ}v}{v_{\text{crit }}}\right)}^2\right]$$(28)

with the pebble relative velocity Δυ and the critical relative velocity $$v_{\text{crit }}={\left(\frac{m}{M_{*}\text{St}}\right)}^{1/3}{v}_K\text{. }$$(29)

In the head wind regime, Δυ = ηυ_K, and, with Eq. (28), the reduction factor can be expressed as $${\kappa }_{\text{OL18,hw }}\approx \exp \left[-{\left(\frac{\text{St}}{0.707\cdot {\text{St}}_{\text{crit }}}\right)}^{2/3}\right]$$(30)

for a more insightful comparison with κ_IGM16 in Eq. (26). By inspection, the dependence on the planetesimal mass m is virtually identical for both cases when m ≲ 2 × 10⁻⁴ M_⊕ for η = 10⁻³, while a factor of about 0.707 is multiplied to m for κ_OL18,hw. Figure A.6 shows the values of κ_IGM16 and κ_OL18,hw with an assumed St = 0.1 and r = 5 au in our disc model. For m < 10⁻⁵ M_⊕, κ_OL18,hw is generally a few times smaller than κ_IGM16. This is in agreement with the findings by Matsumura et al. (2021) and the early stage of the presented simulation results. When the bodies are still dynamically cold, the growth by pebble accretion with ∈_PA = ∈_OL18 is generally slower.

While restricting the discussion in the headwind regime with small e, a pebble accretion onset mass m_PA,hw can be defined (Visser & Ormel 2016; Ormel 2017) by setting Δυ = υ_crit, which yields $$m_{\text{PA},\text{hw}}=\mathrm{St}{\eta }^3{M}_{*}.$$(31)

For m = m_PA,hw, this means κ_OL18 ≈ 0.61 and κ_IGM16 ≈ 0.67. Figure A.7 shows a comparison of m_PA,hw and the planetesimal gravitational mass of the adopted initial planetesimal mass function, m_G, at different locations of the disc. The increase with r for m_PA,hw is steeper than that for m_G. This is in agreement with the results that the growth by pebble accretion is faster in the inner part of the disc. Also, m_G is about 5–10 times smaller than m_PA,hw from 5 to 20 au. This means the massive tail of the planetesimal population overlaps with the mass range for the sharp cut off in the values of the reduction factors for both prescriptions (κ_IGM16 & κ_OL18,hw) as shown in Fig. A.6. As a result, the randomness in the exact number of particles drawn near the top end of the distribution as well as that in their locations play a significant role to the final architecture of the modelled planetary systems.

This is more clearly shown by the difference in the results with N_ini = {1000, 2000, 5000} while all have S_a = 0 and ∈_RA = ∈_OL18. As the number of planetesimal increases, the largest drawn mass increases slightly as well due to the higher probability of getting at least one particle with such mass. This leads to an earlier formation of massive cores, which are more likely to become gas giants by the time of disc dispersal while fewer or no ice giants remain. Nonetheless, this effect is not observed with ∈_PA = ∈I_GM16 likely due to a generally more efficient pebble accretion such that gas accretion starts early for the massive cores with enough time to reach the mass of a gas giant even with N_ini = 1000.

Although our results show an apparent dependence on the initial number of particles N_ini, we emphasize that this can be a result of a statistical artefact. With the adopted initial mass function by Abod et al. (2019), as shown in Eq. (8), there is no upper limit on the planetesimal mass. Although an artificial upper limit of 10 times of the characteristic mass is imposed, this limit has a negligible effect on the actual realized planetesimal populations, where only a number fraction of planetesimals of ~8 × 10⁻⁶ is lost. Therefore, the massive tail of the initial planetesimal population drawn in this manner has a dependence on the number of particles, which sets the normalization constant of the initial mass function. This means a physical upper limit of planetesimal mass (e.g. Gerbig & Li 2023) is needed to remove this artefact for future investigations. Nonetheless, our results show the upper end of the initial planetesimal population plays the most important role in growth by pebble accretion while the rest of the small planetesimals do not affect their growth significantly.

We note that in Lambrechts & Johansen (2012), the transition mass of an embryo m_t is defined as the mass at which the Hill radius is comparable to the Bondi radius, i.e. $$m_t=\sqrt{\frac{1}{3}\frac{{\left(\eta {\nu }_K\right)}^3}{G{\text{Ω}}_K}\approx 0.578{\eta }^3{M}_{*}.}$$(32)

This mass is often adopted as the initial embryo mass in the works involving pebble accretion (e.g. Bitsch et al. 2015). The value of m_t is a few times larger than m_PA,hw from Eq. (31) for St = 0.1. This indicates these initial embryos can always grow by pebble accretion efficiently. The evolution from the point of planetesimal formation to the onset of pebble accretion is missing in this approach. We also note that the characteristic mass (0.3 M_G) of the adopted initial planetesimal mass function is about an order of magnitude less massive than m_PA,hw, which is comparable to the value adopted by Coleman (2021). However, in the expression for m_G in this work as shown by Eq. (9), the value of the small-scale diffusion parameter δ can be an order of magnitude larger or smaller than the adopted value (Schreiber & Klahr 2018). This translates to an even larger uncertainty in the initial planetesimal mass since m_G ∝ δ^3/2, which shall greatly change the results of our model and will require further investigations.

4.1.2 Pebble accretion and dynamical heating

However, as noted by Ida et al. (2016), the assumption of small e in the estimation of the pebble relative velocity only holds when e < η ~ 10⁻³. This condition breaks down quite early in the presented simulations, with a majority of the particles having e exceeding 10⁻³ by 0.75 Myr in all the presented simulations. Multiple planet formation models (e.g. Levison et al. 2015; Jang et al. 2022; Lau et al. 2022; Jiang & Ormel 2022) have shown the effect of increased pebble relative velocity due to dynamical heating on pebble accretion. Figure A.6 also includes the general form of κ_OL18 with Δv = max(0.76e, ηυ_K) (Liu & Ormel 2018) for e = 10⁻², where the curve is shifted towards higher m by more than two orders of magnitudes, i.e. a much larger m is required for efficient pebble accretion. Therefore, it is likely an important feature of a realistic model to consider the effect of pebble accretion being interrupted when the eccentricities grow, especially in the context of planet formation where massive cores and giant planets are formed among planetesimals. However, once planet migration is in effect and cores of ~1–10 M_⊕ are readily removed, they cannot continuously excite and eject the planetes-imals. Pebble accretion in this case is not severely interrupted by dynamical heating as shown in the results (Sect. 3.2), so that both pebble accretion prescriptions yield more similar results at the end of the simulations when migration and removal is present.

We note that there are other works on the initial planetesimal mass function (e.g. Simon et al. 2016, 2017; Schäfer et al. 2017; Gerbig & Li 2023), and this topic remains an active field of research. Meanwhile, the outcome of the subsequent growth of the initial planetesimals is sensitive to their initial mass as well as the distribution. Also, we assume a planetesimal disc as a part of the initial conditions, but its formation is not investigated in this work, which is also an active field of research (e.g. Drążkowska et al. 2016; Carrera et al. 2017; Schoonenberg et al. 2018; Lenz et al. 2019, 2020). These parts of the model concerning the initial planetesimals require further investigations for a more robust planet formation model.

4.2 Planet migration

When planet migration is turned off in our model, i.e. S_a = 0, the results with N_ini = 1000 and ∈_PA = ∈_OL18 (Sect. 3.1.2, Fig. A.1) show one to two gas giants and one to two ice giants beyond 6 au. This is in general agreement with Levison et al. (2012) in forming the giant planets in the Solar System without forming hundreds of massive cores in the process. In their work, planet migration is not considered either.

However, once planet migration is turned on in our model, i.e. S_a = 1, the results (Sect. 3.2) show that cores of ~1–10M_⊕ rapidly migrate towards the inner part of the disc and many leave the simulation domain as a migration trap is not implemented at the inner edge of the disc. This is in agreement with previous works on planet formation that include planet migration (e.g. Cossou et al. 2014; Coleman & Nelson 2016b; Matsumura et al. 2017; Jang et al. 2022). Although the migration timescale in the high-mass regime in this work is already lengthened by setting a turbulent-α parameter α_turb that is only one-tenth of the classical α parameter α_acc, it is still not enough to retain these massive cores at wide orbit in our model to form Solar-System-like giant planets. Further parameter search may be required to produce cold giant planets with planet migration in effect but the current results suggest that some massive cores are inevitably lost to the inner Solar System in the process as shown in other works (e.g. Bitsch et al. 2019; Matsumura et al. 2021).

Figure A.8 shows a heat map of the migration timescale τ_a in the m–r space at t = 0.5 Myr in our model. There is a region of rapid migration with τ_a ~ 10⁵ yr for m ~ 1–10 M_⊕ across the planetesimal disc. This is in agreement with the results that the massive cores have migrated significantly before runaway gas accretion can occur for them to enter the high-mass regime of migration where τ_a ~ 10⁶ yr. For the surviving cores, migration only stops as the gas surface density becomes very low that slows down migration but this also terminates gas accretion as shown in the results. Also, it seems to be a general result that multiple massive cores (~1–10 M_⊕) inevitably enter the inner Solar System with a smooth disc model where migration trap is not present except at the inner edge of the protoplanetary disc. In contrast, other works (e.g. Coleman & Nelson 2016a; Lau et al. 2022) have shown a possibility in retaining these cores at wide orbit due to the presence of a substructure in the gas disc. These findings and the recent observations of substructure in protoplanetary discs (e.g. Andrews et al. 2018; Long et al. 2018; Dullemond et al. 2018; Cieza et al. 2021) suggest that a substructure in the pro-toplanetary disc is a promising way to interrupt rapid migration and prevent the formation of super-Earths and hot Jupiters.

5 Conclusions

This work attempts to form the giant planets of the Solar System in a smooth protoplanetary disc. An initial planetesimal disc is simulated with the parallelized N-body code SyMBAp with additional subroutines to include the effects of pebble accretion, gas accretion, and planet-disc interactions with the protoplanetary disc.

Our model starts from planetestesimals (each with m < 10^–4 M_⊕) instead of planetary embryos (m ~ 10^–2 M_⊕). In this work, we demonstrate the difference between the pebble accretion prescription by Ida et al. (2016) and that by Liu & Ormel (2018) and Ormel & Liu (2018). In Ida et al. (2016), the pebble-accreting body is assumed to be in a circular orbit and the pebble relative velocity, which sets the pebble encounter time, is set by the headwind in the disc. In contrast, Liu & Ormel (2018) and Ormel & Liu (2018) do not hold this assumption and consider the relative velocity due to eccentricity and inclination. In the case that the number of embryos is small and they are well above the pebble accretion onset mass both prescriptions give similar results, as noted in Matsumura et al. (2021). However, in a planetesimal disc, viscous stirring becomes important and can effectively terminate growth by pebble accretion due to the increased pebble relative velocity and shortened pebble encounter time. This can occur when the inclinations of the bodies are small, and they are still well inside the pebble disc as also noted by Lau et al. (2022). Therefore, to realistically model planet formation via pebble accretion starting from planetesi-mals, it is crucial to consider the reduced pebble encounter time due to dynamical heating.

When planet migration is not considered, our model can reproduce one to two gas giants and one to two ice giants beyond 6 au, which is analogous to the giant planets in the Solar System. However, we also note that the results have a dependence on the initial number of planetesimals. Further studies on the processes involved in planetesimals formation is required to construct a more realistic model.

Once planet migration is in effect, massive cores of about 10 M_⊕ are readily removed as they migrate towards the inner boundary of the simulations. This shows that the formation of the giant planets in the Solar System requires an alternative and effective way to stop the migration of the first massive body formed before reaching the inner Solar System. Multiple works (e.g. Coleman & Nelson 2016a; Lau et al. 2022) have demonstrated that pressure bump in the disc can act as a migration trap while some other works (e.g. Jiang & Ormel 2022; Chrenko & Chametla 2023) do not support this scenario. Further investigations are required to characterize the disc conditions that can retain massive planetary cores and allow the formation of cold gas giants.

Acknowledgements

This work was supported in part by a postgraduate studentship (T.C.H.L.) and a seed fund for basic research (T.C.H.L. and M.H.L.) at the University of Hong Kong and Hong Kong RGC grant 17306720 (T.C.H.L. and M.H.L.). Parts of the computations were performed using research computing facilities offered by HKU Information Technology Services. T.C.H.L. also acknowledges funding from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under grant 325594231. This research was supported by the Munich Institute for Astro-, Particle and BioPhysics (MIAPbP) which is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’ s Excellence Strategy — EXC-2094 – 390783311.

Appendix A Additional figures

Fig. A.1 Same as Fig. 5 except ∈PA = ∈OL18.

Fig. A.2 Same as Fig. A.1 except Nini = 2000.

Fig. A.3 Same as Fig. A.1 and A.2 except Nini = 5000.

Fig. A.4 Same as Fig. 5 except with migration turned on (Sa = 1).

Fig. A.5 Same as Fig. A.1 except Sa = 1.

Fig. A.6 Reduction factor ĸ on the pebble accretion cross section in our disc model according to the prescription by Ida et al. (2016), ĸIGM16 (solid line), a similar reduction factor by Liu & Ormel (2018) in the head wind regime with small e and i ĸOL18,hw (dashed line), and that with e = 10–3 and small i ĸOL18,e = 10–2 (dotted line) as different mass m. The values of St = 0.1 and r = 5 au are set for this estimation.

Fig. A.7 Pebble accretion onset mass mPA,hw in the headwind regime (solid line) and the planetesimal gravitational mass mG (dashed line) at different locations of the disc r with the corresponding value of η shown in the upper x-axis. The values of mG are around the mass range for the steep cutoff in k shown in Fig. A.6.

Fig. A.8 Heat map of the migration timescale τa in the m–r space at t = 0.5 Myr in our model. A region of rapid migration (τa ~ 105 yr) is presence for m ~ 1 – 10 M⊕ across the planetesimal disc (r = 5 – 20 au).

References

All Figures

	Fig. 1 Time evolution of ${\dot{M}}_{*}$, with the initial age of the disc t0 = 0.5 Myr. The value of ${\dot{M}}_{*}$ is turned down linearly when it drops below 10−9 M⊙ yr−1 to mimic the effect of photoevaporation.
In the text

	Fig. 2 Adopted truncated power law initial planetesimal mass function as described by Eq. (8) based on Abod et al. (2019). It is presented in the unit of the planetesimal gravitational mass mG
In the text

	Fig. 3 Initial mass distribution of the realized planetesimal discs. One example is shown for each of the chosen values of Nini. The width of each bin is 0.2 au.
In the text

	Fig. 4 Time evolution of the pebble mass flux ${\dot{M}}_{\text{pf}}$ given by Eq. (11) with Z0 = 10–2 and αacc = 10–3 as set in this work.
In the text

Fig. 5 Results for the simulations with Nini = 1000, migration turned off (Sa = 0), and the pebble accretion efficiency ∈PA = ∈IGM16. Each row presents a snapshot of the simulations at the time indicated by the timestamp on the left. For the first three columns from the left, the total occurrences of particles across all five random simulations are shown by heat maps with 2 × 2 cells in each minor axis grid cell. The left-most column shows m and a. The next two columns to the right respectively shows the e and i, respectively, against a. The right-most column shows the differential mass distribution of all bodies with each colour corresponds to one of the five simulations. Particles in one of the five simulations (blue) is also plotted with enlarged dots denoting particles above 10−3 M⊕. For 6.5 Myr, i.e. the end of the simulations, particles above 10−3 M⊕ in all five simulations are shown individually without using heat maps. Further descriptions are in the text.

In the text

Fig. 6 End results for the simulations with Nini = 2000 and 5000, respectively, as indicated on the left, migration turned off (Sa = 0), and the pebble accretion efficiency ∈PA = ∈IGM16. The two columns here correspond to the left-most and the right most columns of Fig. 5, respectively. There is no qualitative difference in the end results among the simulations with the chosen set of Nini = {1000, 2000, 5000}.

In the text

	Fig. A.1 Same as Fig. 5 except ∈PA = ∈OL18.
In the text

	Fig. A.2 Same as Fig. A.1 except Nini = 2000.
In the text

	Fig. A.3 Same as Fig. A.1 and A.2 except Nini = 5000.
In the text

	Fig. A.4 Same as Fig. 5 except with migration turned on (Sa = 1).
In the text

	Fig. A.5 Same as Fig. A.1 except Sa = 1.
In the text

Fig. A.6 Reduction factor ĸ on the pebble accretion cross section in our disc model according to the prescription by Ida et al. (2016), ĸIGM16 (solid line), a similar reduction factor by Liu & Ormel (2018) in the head wind regime with small e and i ĸOL18,hw (dashed line), and that with e = 10–3 and small i ĸOL18,e = 10–2 (dotted line) as different mass m. The values of St = 0.1 and r = 5 au are set for this estimation.

In the text

	Fig. A.7 Pebble accretion onset mass mPA,hw in the headwind regime (solid line) and the planetesimal gravitational mass mG (dashed line) at different locations of the disc r with the corresponding value of η shown in the upper x-axis. The values of mG are around the mass range for the steep cutoff in k shown in Fig. A.6.
In the text

	Fig. A.8 Heat map of the migration timescale τa in the m–r space at t = 0.5 Myr in our model. A region of rapid migration (τa ~ 105 yr) is presence for m ~ 1 – 10 M⊕ across the planetesimal disc (r = 5 – 20 au).
In the text