EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 688 (August 2024) A&A, 688 (2024) A22 Full HTML
Open Access
Issue
A&A
Volume 688, August 2024
Article Number A22
Number of page(s) 12
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/202450464
Published online 31 July 2024

© The Authors 2024

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1 Introduction

Planet formation is a multi-step process spanning over 40 orders of magnitude in mass. In recent years, there has been significant progress in the understanding of this process driven by the numerous discoveries of exoplanets and observations of protoplanetary discs (see Drążkowska et al. (2023) for a recent review). The formation of the first gravitationally bound building blocks of planets, the planetesimals, which used to be a major bottleneck of the planet formation theory, was addressed by the streaming instability (Johansen et al. 2007). The properties of planetesimals formed in the streaming instability broadly match comets and the Kuiper belt objects (Blum et al. 2017; Nesvorný et al. 2019). The growth timescale of giant planet cores, which was prohibitively long in the classical planetesimal-driven paradigm (Lissauer 1987; Kokubo & Ida 2000, 2002), has been addressed by introducing pebble accretion (Ormel & Klahr 2010; Lambrechts et al. 2014). Despite these new developments, a consistent model covering all the stages of planet formation does not exist yet. Formation of the cores of giant planets remains a challenge, particularly at large orbital distances (Voelkel et al. 2020; Coleman 2021; Eriksson et al. 2023).

Current planet formation models in general fail to meet both the physical and cosmochemical constraints required to explain the Solar System (e.g. Matsumura et al. 2017; Liu et al. 2019; Bitsch et al. 2019; Lau et al. 2024). The main challenge remains to be the ‘migration problem’, where rapid migration occurs for planetary cores of 1 to 10 M⊕, resulting in the formation of super-Earths and hot Jupiters. While the above models generally assume a smooth planetary disc, multiple works (e.g. Coleman & Nelson 2016; Morbidelli 2020; Guilera et al. 2020; Chambers 2021; Andama et al. 2022; Lau et al. 2022) have modelled the formation and evolution of planetary cores retained at the migration trap near a pressure bump in the disc. Nonetheless, the origin of such pressure bumps remains uncertain, the proposed non-planetary causes include late-stage infall of material (Gupta et al. 2023), sublimation (Saito & Sirono 2011), instabilities (Takahashi & Inutsuka 2014; Flock et al. 2015; Dullemond & Penzlin 2018), and the edge of the dead zone, where the magneto-rotational instability (MRI) is suppressed (Pinilla et al. 2016).

More recently, the high-resolution interferometry observations by the Atacama Large Millimeter/submillimeter Array (ALMA) have shown that substructures are typical in protoplanetary discs. Multiple surveys (e.g. Andrews et al. 2018; Long et al. 2018; Cieza et al. 2021), have shown that most of the substructures are presented in the form of axisymmetric rings. While these observations are limited to large and bright discs, disc population synthesis and theoretical models (Toci et al. 2021; Zormpas et al. 2022) have demonstrated that disc substructures may be common in unresolved discs as well.

Chatterjee & Tan (2013) presented an analytic model demonstrating the ‘inside-out’ planet formation scenario, where planet formation starts at the outer edge of the MRI active zone around the star. Although this work focuses on explaining the tightly packed chains of planets commonly seen in exoplanetary systems, the model suggests the possibility of planet formation being triggered by the planet formed in the previous generation. A similar idea was also proposed for the formation of Saturn after the completion of Jupiter (Kobayashi et al. 2012), in which the core of Saturn grows rapidly without significant inward drift in the pressure bump induced by the planetary gap of Jupiter, although they did not consider planetesimal formation from dust and pebble accretion.

Motivated by the current models and observations, a substructure in a protoplanetary disc has recently emerged as an ideal location for planet formation (Chambers 2021; Lau et al. 2022; Jiang & Ormel 2023). Lau et al. (2022) modelled the formation and evolution of planetesimals in an initial axisymmetric disc substructure by coupling the dust and gas evolution code DustPy (Stammler & Birnstiel 2022) with the parallelised symplectic N-body integrator SyMBA parallelised (SyMBAp; Lau & Lee 2023). As the disc evolves and satisfies the condition for the streaming instability and the subsequent gravitational collapse, a fraction of the dust is converted into planetesimals as N-body particles. On top of the full N-body gravitational interactions, additional subroutines are gas drag, planet-disc interactions, and pebble accretion. Lau et al. (2022) showed the rapid formation of planetary cores thanks to the concentrated dust, which are also retained due to the migration trap near the pressure bump. However, this work did not attempt to form giant planets, as the authors focused on the formation of planetary cores, where gas accretion and planetary gap opening are missing in the model.

In this work, we further developed the model in Lau et al. (2022) for the formation of giant planets initiated by a disc substructure and found a scenario of sequential planet formation. In the following, Sec. 2 summarises the methods adopted in Lau et al. (2022) and the new components implemented in this work. The results are presented in Sec. 3, which are followed by the discussions in Sec. 4. The findings of this work are summarised in Sec. 5.

2 Method

We employed the dust and gas evolution code DustPy v1.0.3 (Stammler & Birnstiel 2022) and the symplectic N-body integrator SyMBAp v1.6 (Lau & Lee 2023), a parallelised version of the Symplectic Massive Body Algorithm (SyMBA; Duncan et al. 1998). The coupling of the two codes to construct a consistent planet formation model was first presented in Lau et al. (2022), which only modelled the formation of planetary cores. In this work, we added gas accretion and planetary gap opening to model the subsequent evolution of the embryos formed at a pressure bump. The following summarises the employed method in Lau et al. (2022) and describes the new components in detail.

2.1 Disc model

DustPy simulates the viscous evolution of the gas, coagulation, fragmentation, advection, and diffusion of the dust in a protoplanetary disc. The different parts of the disc model are described in this section.

2.1.1 Gas component

We considered a protoplanetary disc around a solar-type star, which is axisymmetric and in vertical hydrostatic equilibrium. The initial gas surface density ∑g,init is given by (1)

with the distance from the star r, the initial mass of the disc Mdisc, and the characteristic radius rc. We set Mdisc = 0.0263 M⊙ and rc = 50 au, which imply ∑g,init (r = 5 au) ≈ 134.6 g cm−2.

The gas disc viscously evolves in time t according to the advection-diffusion equation (2)

(Lüst 1952; Lynden-Bell & Pringle 1974), while the backreaction from the dust is neglected. The Shakura & Sunyaev (1973) α-parametrisation was adopted for the kinematic viscosity ν such that (3)

with the speed of sound cs and the disc scale height Hg. The viscosity parameter α = {3, 5} × 10−4 was set in this work. The disc scale height is defined by Hg ≡ cs/ΩK, where the local Keplerian orbital frequency ***(eq4)*** with the gravitational constant G. The isothermal sound speed was used and given by ***(eq5)*** with the Boltzmann constant kB, the midplane temperature T and the mean molecular weight of the gas μ = 2.3 mp. The disc was assumed to be passively irradiated by the solar luminosity at a constant angle of 0.05, which gives a midplane temperature profile (4)

We note that the normalisation is about 0.84 of that in Lau et al. (2022), which is the result of a correction made to DustPy since v1.0.2. This setup yields the dimensionless gas disc scale height (5)

The midplane pressure gradient parameter η is then given by (6)

with the midplane gas pressure P, which describes the degree of ‘sub-Keplerity’ of the gas. A logarithmic radial grid was adopted with 133 cells from 3 to 53 au and with an additional 42 cells from 53 to 1000 au.

2.1.2 Dust component

The initial dust surface density ∑d,init is given by (7)

with the global dust-to-gas ratio Z set at the solar metallicity of 0.01. We followed Mathis et al. (1977), that is, the MRN size distribution of the interstellar medium, for the initial size distribution of the dust grains. The maximum initial size was set at 1 μm with the internal density of 1.67 g cm−3 assumed. A total of 141 dust mass bins logarithmically spaced from 10−12 to 108 g were used. Each dust species was evolved in time through transport with the advection-diffusion equation (Clarke & Pringle 1988) coupled to growth and fragmentation with the Smoluchowski equation. The fragmentation velocity was assumed to be 5 ms−1. At collision velocities above which, the dust aggregates are assumed to fragment. The Stokes number Sti of the dust in each dust species i was calculated by considering the Epstein and the Stokes I regimes. The dust scale height of each dust species Hd,i was calculated according to Dubrulle et al. (1995), (8)

assuming Sti < 1. Further details of the algorithms for the disc model are described in Stammler & Birnstiel (2022).

2.1.3 Initial disc gap

An initial axisymmetric gap was introduced to the disc following the model by Dullemond et al. (2018), which is motivated by the commonly observed substructures in protoplanetary discs. To modify the gas profile, we applied a modified α-parameter with radial dependence α′(r) = α/F(r), where the function (9)

with the gap amplitude A = 1, the location r0 = 5.5 au and the gap width w = 0.5 au. The initial gap is removed when the first planet has the gap opening factor K (Kanagawa et al. 2015) of 250. The K factor is further described in Sec. 2.3 on planetary gap opening. This value translates to a perturbation towards the gas disc of about 1/10 of the unperturbed gas surface density.

Since we do not study the physical cause of the initial disc gap, the modified α′(r)-parameter only serves the purpose to attain a target gas profile and does not change the actual turbulence in the disc. Therefore, the modified α-parameter α′(r) is exclusively experienced by the gas, while dust diffusion, dust scale height and turbulent collision speeds are set by α. Nonetheless, the dust evolves according to the resulting gas profile. We note that this treatment is not consistent with the substructure formation scenarios where the actual turbulence is changed, for example, the edge of the dead-zone, but consistent with the cases where the turbulence is unchanged, for example, infall of material.

2.1.4 Planetesimal formation

Lau et al. (2022) adopted the truncated power-law cumulative mass distribution from the fitting by Abod et al. (2019). However, the upper end of the distribution is not limited and Lau et al. (2024) noted that the largest planetesimal in the actual realisation depends on the total number of planetesimals. Therefore, here we adopted the Toomre-like instability criterion Qp for the gravitational collapse of the dense filament induced by streaming instability and the initial mass function from Gerbig & Li (2023). The criterion for collapse is Qp < 1 with (10)

and the mass-averaged Stokes number of the dust in the cell Stavg. The local dust surface density ∑d,local was assumed to be 10 times of the averaged ∑d of the cell. Motivated by the streaming instability simulations in Schreiber & Klahr (2018), the small-scale diffusion parameter δ was set at 10−5. The model converts dust into planetesimals based on the prescription by Drążkowska et al. (2016) and Schoonenberg et al. (2018). The Qp criterion was combined with the smooth planetesimal formation activation function from Miller et al. (2021), which is given by (11)

and evaluated at each radial grid cell. If any cell also satisfies the criterion of ρd/ρg ≥ 1, the local dust surface density for each dust species i is reduced by (12)

The planetesimal formation efficiency per settling time is ζ = 10−3 and the settling time of dust species i is tset,i ≡ 1/(StiΩK).

Then, the removed dust is summed over all dust species and added to the local planetesimal mass surface density. We first draw the location of a new planetesimal using the radial profile of the planetesimal mass surface density. Then, we draw the planetesimal mass according to the initial mass function given by Gerbig & Li (2023), which is resulting from the stability analysis of the dispersion relation for dust influenced by turbulent diffusion (Klahr & Schreiber 2021). The maximum and minimum masses associated with the unstable modes are given by (13)

and the fastest-growing mode is given by (14) (15)

assuming isotropy in the small-scale diffusion. We refer the readers to Eq. (20) in Gerbig & Li (2023) for the expression of the probability density function. Figure 1 shows an example of the differential mass distribution of the planetesimals formed at about 6.5 au near the initial disc substructure for the set of simulations with α = 5 × 10−4.

As described in Lau et al. (2022), the eccentricity e and the inclination i in radian are randomly drawn from two Rayleigh distributions with the scale parameters 10−6 and 5 × 10−7 respectively. The rest of the angles of the orbital elements in radian are drawn randomly from 0 to 2π. The physical radius Rp is calculated by assuming an internal density ρs of 1.5 g cm−3. The drawn mass is subtracted from the surface density from the nearest radial grid cells and the realisation stops when the total remaining mass is less than the drawn mass. Any residue of the planetesimal mass surface density and the last drawn value of planetesimal mass is retained for the next time step to avoid bias towards the lower mass, that is, the drawn mass will be realised as soon as enough planetesimal surface density is available.

thumbnail Fig. 1

Differential mass distribution of the planetesimals formed near the initial disc substructure for the set of simulations with α = 5 × 10−4 drawn from the adopted initial mass function by Gerbig & Li (2023). There are 10 bins per decade in mass m and ΔN is the number of planetesimals in each bin. Each colour corresponds to one of the five simulations.

2.2 Planetesimal evolution

The realised planetesimals, and a Solar-mass star, were then evolved by SyMBAp with full gravitational interactions as well as additional subroutines to include gas drag, the planet-disc interactions, pebble accretion, gas accretion and planetary gap opening. If two bodies collide, they were assumed to merge completely, that is, collisions are perfectly inelastic. At each communication step, on top of the newly formed planetesimals, the radial profiles of the disc were passed to SyMBAp. This included the gas components:

  • the gas surface density ∑g;

  • the midplane temperature T;

  • the gas disc scale height Hg;

  • the midplane gas density ρg, and;

  • the midplane pressure gradient parameter η, and the dust component, for each dust species i,:

  • the Stokes number Sti;

  • the dust disc scale height Hd,i, and;

  • the dust surface density ∑d,i.

Also, the feedback to the disc was also passed to DustPy including

  • the dust mass loss due to pebble accretion (Sect. 2.2.1);

  • the gas mass loss due to gas accretion (Sect. 2.2.2), and;

  • the change in gas surface density due to planetary gap opening (Sec. 2.3),

which are further described in the respective subsections.

2.2.1 Pebble accretion

The treatment for the pebble accretion rate of each planetesimal is identical to that presented in Lau et al. (2022) and summarised below. First, the pebble mass flux of dust species i at any location is given by (16)

The pebble drift speed of dust species i is υdrift,i = 2Sti|η|rΩK (Weidenschilling 1977). Then, we implemented the pebble accretion efficiency factor ϵ PA,i by Liu & Ormel (2018) and Ormel & Liu (2018) to calculate the fraction of the local pebble mass flux being accreted by each planetesimal or planet for dust species i. The pebble accretion rate by a planetesimal or a planet is then given by summing the contributions from all dust species, which is (17)

The mass of the accreted pebbles is then subtracted from the respective dust species and radial cell of the dust disc at the next immediate communication step.

We did not implement the pebble isolation mass explicitly with a prescription (e.g. Lambrechts et al. 2014) in this work. Instead, as dust and gas evolve consistently in this model, the planetary gap opening by a planet (Sec. 2.3) can interrupt the pebble flux capturing the process of pebble isolation within the model.

2.2.2 Gas accretion

We followed Piso & Youdin (2014) and Bitsch et al. (2015) to prescribe the gas accretion rate with the modification by Chambers (2021) to account for the energy released from pebble accretion. Gas accretion generally begins when the energy released from pebble accretion decreases enabling the cooling of the gas envelope. The gas accretion rate in this phase is (18)

with the opacity of the gas envelope κ and the density of the core ρc = 5.5 g cm−3. We note that the solid mass accreted from planetesimal accretion is negligible compared to pebble accretion in our simulations. We followed Brouwers et al. (2021) for the grain size near the Bondi radius of the gas envelope to calculate its Rosseland mean opacity. Assuming the Epstein regime, the incoming solids converge to the size (19)

with the midplane gas density ρg, the thermal velocity υth, the gravity at the Bondi radius gB and the monomer density ρ• = 1.67 g cm−3. The Rosseland mean opacity is given by (20)

with the midplane dust density ρd. The extinction efficiency is Qeff = min(0.6πRs/λpeak, 2) with the peak wavelength of the emission given by (21)

where the temperature T is assumed to be the local disc midplane temperature.

When the gas envelope menv exceeds the solid core mass mc, gas accretion enters the runway phase following the treatment by Bitsch et al. (2015). The gas accretion rate in this phase is determined by the gas stream flowing towards the planet (Tanigawa & Tanaka 2016), which is (22)

with the planet mass m. The accreted gas is then subtracted from the gas disc assuming the half-width of the accretion zone equals twice the Hill radius of the planet.

2.2.3 Physical radius

As the planetesimals or planets grow by many orders of magnitude in mass, the physical radius Rp is evaluated correspondingly. For mass m less than 0.1 Earth mass, we assumed an internal density ρs of 1.5 g cm−3, which is the same when the planetesimals are formed. For m above 0.1 M⊕ but less than 5 M⊕, we followed Seager et al. (2007) for the mass-radius relationship of rocky planets, which is (23)

with the radius of Earth R⊕. For m above 5 M⊕, we followed the mass-radius relationship applied in Matsumura et al. (2017), which is (24)

2.3 Planetary gap opening

For all planetesimals and planets, the non-dimensional gap opening factor (Kanagawa et al. 2015) was evaluated, which is given by (25)

Since the treatment in Sect. 2.1.3 to attain a target gas profile does not change the actual turbulence experienced by the dust and we only consider the torque exerted on the gas disc, planetary gap opening was applied through the modified α′-parameter. When K > 0.25, we impose a planetary gap to the gas disc by dividing α′ by the ratio of the perturbed surface density to the unperturbed one ∑g/∑g,0. In other words, the planetary gap is imposed when the change caused by the corresponding body is more than about 1% of the unperturbed surface gas density.

The effects on α′ were multiplied when more than one planet can open a gap, that is, for all gap opening planets i (26)

with the function to impose the initial disc gap F (see Eq. (9)).

We note that some small but short-period chaotic movements of massive planets can prevent the disc from converging to a quasi-steady state and the speed of code is significantly reduced. Therefore, we allowed α′(r) to reach the target value exponentially on a relaxation timescale of 250 × (10−3/α) yr, which is below the viscous evolution timescales of the disc in the adopted radial domain of the simulations.

We adopted the empirical formula by Duffell (2020) for the gap profile, which is (27)

with ***(eq30)*** evaluated at the planet’s location. The radial profile function ***(eq31)*** is defined by (28)

with the mass ratio q ≡ m/M⊙, the planet’s radial distance from the star rp and the scaling factor ***(eq33)***. The function δ(x) is given by (29)

with the factor ***(eq35)*** and the factor ***(eq36)***. We note that there is a minor discrepancy in the form of δ(x) for the case of x < qNL between the text and the Python code presented in Duffell (2020). We confirm that the expression given above by Eq. (29) is consistent with the mentioned Python code, which is continuous and more appropriate (Duffell 2022, priv. comm.).

2.4 Gas drag and planet-disc interactions

All bodies experience the combined effects of aerodynamic gas drag and the planet-disc interactions. For low-mass planet without gap opening, the treatment for gas drag and planet-disc interaction are identical to that presented in Lau et al. (2022) and summarised below.

We adopt the aerodynamic gas drag by Adachi et al. (1976), which is (30)

with the drag coefficient CD and, the relative velocity between the planetesimal and the gas υrel. The gas flow is assumed to be laminar and cylindrical, where the magnitude is given by rΩK(1 − |η|). As the planetesimals in this work are well larger than a kilometre in size, the large Reynolds number case is generally applicable, that is, CD = 0.5 (Whipple 1972). The gas density ρ at the planetesimal’s position z above the midplane is given by ***(eq38)***.

For type-I damping and migration, we adopted the prescription based on dynamical friction by Ida et al. (2020). The timescales for the isothermal case and finite i, while ***(eq39)***, (Appendix D of Ida et al. 2020) were implemented. The evolution timescales of semimajor axis, eccentricity and inclination are defined, respectively, by (31)

These timescales are given by, with ***(eq41)***, (32) (33) (34)

The characteristic time twav (Tanaka et al. 2002) is given by (35)

where ∑g and ***(eq46)*** are retrieved from the local radial cell of the disc model. The normalised torques CM and CT are given by (36) (37)

with p∑ ≡ −d ln ∑g/d ln r and qT ≡ −d ln T/d ln r. The three timescales were then applied to the equation of motion (38)

in the cylindrical coordinates (r, θ, z) that the velocity of the embryo υ = (υr, υθ, υz). And, the local Keplerian velocity υK was evaluated at the instantaneous r of the particle.

As the planet grows and opens a gap in the disc, Kanagawa et al. (2018) suggested that the magnitude of the torque scales linearly with the local surface density providing a smooth transition to the high-mass (type-II) regime of planet migration. Since the dependence of twav on ∑g as in Eq. (35) is retrieved from the local grid cell, the above treatment combined with gap opening (Sec. 2.3) can already capture this transition.

We note that planetary gap opening (Sec. 2.3) and planet migration are the results of planet-disc interactions, which are physically coupled by the action-reaction pair. However, we have adopted two independent prescriptions for each of them since a prescription for a torque profile which is suitable for a one-dimensional model and a general planet mass is still missing. Further discussions on the adopted treatments are in Sect. 4.4.2.

2.5 Numerical setup

The time step for SyMBAp τ = 0.2 yr was used and particles were removed if the heliocentric distance is less than 4 au or greater than 100 au. The additional subroutines for the evolution of the N-body particles were added to SyMBAp as (39)

The operator handles the effect of pebble accretion, gas accretion and gap opening, handles the effect of gas drag and planet-disc interactions, and is the second-order symplectic integrator in the original SyMBAp. The operators and operate in the heliocentric coordinates and operates in the democratic heliocentric coordinates so coordinate transformation is required at each step.

Since disc dissipation is not included in this work, all simulations stop at 2 Myr, which is the typical timescale that internal photoevaporation becomes significant to the disc (e.g. Owen et al. 2010, 2011; Picogna et al. 2019; Gárate et al. 2021). We note a numerical difficulty when multiple giant planets are produced, which causes multiple deep planetary gaps in the disc and a small integration time step is required for the disc. Each simulation requires a wall-clock time of 2 to 4 weeks. We tested two values of α = {3, 5} × 10−4 and five simulations were conducted for each to evaluate the statistical effect.

3 Results

3.1 The case of α = 5 × 10−4

3.1.1 Formation and evolution of massive bodies

Figure 2 presents one of the simulations with α = 5 × 10−4 with the panels showing the six key timestamps. The solid and dashed lines show the profiles of the gas surface density ∑g and the dust surface density ∑d respectively. The dots show the mass m and the semimajor axis r of the massive bodies, with the error bar indicating the extent of the apoapsis and periapsis for those above 10 M⊕. The final panel (vi), which presents the end results at 2 Myr, also includes the massive bodies from the rest of the simulation set with each colour showing one of the five simulations.

Figure 3 presents the evolution of the semimajor axis of the massive bodies with each colour showing one of the five simulations corresponding to those in the final panel of Fig. 2. The solid lines show the ones reaching more than 10 M⊕ by the end of the simulations. The shaded areas indicate the extents of the lower and upper quartiles of semimajor axes, which is only shown when the total number of bodies is above 50 for a meaningful representation.

At 0.05 Myr (i), the imposed initial substructure has reached the target shape and dust has started to accumulate. At 0.17 Myr (ii), the midplane volumetric dust-to-gas ratio of the disc at about 6.5 au reaches the criteria and planetesimal formation starts. Since these bodies are naturally born in a dust-rich environment, where the dust surface density is more than an order of magnitude higher than the unperturbed case, they can grow rapidly by pebble accretion. The core has also migrated towards the migration trap, which is slightly interior to the peak of the pressure bump, but not further inside as shown by the track in Fig. 3.

At 0.34 Myr (iii), the first massive core has entered the runaway gas accretion phase and opened a significant gap in the disc as it becomes a gas giant (> 100 M⊕). The less massive core and planetesimals, which are also formed from the initial pressure bump, are being scattered mainly to wider orbits as shown by the tracks in Fig. 3. While most planetesimals have been scattered out of the system, the orbit of the second-most massive planetary core is circularised near 8–9 au, and it continues to grow.

At 0.86 Myr (iv), the second core also starts runaway gas accretion but at a much later time relative to the first one. When the second gas giant opens a gap in the disc, the dust near its location follows the sudden change in the gas profile and is pushed away from the forming gas giant to both inner and outer part of the disc. This corresponds to the formation of a small batch of planetesimals near 11 au shown in Fig. 3.

At 1.5 Myr (v), the second gas giant also reached approximately one Jupiter mass with another planetary gap fully opened. A new pressure bump is steadily formed at the outer edge of this gap and dust re-accumulates at about 13–14 au, which contains a part of the leftover dust from the initial dust trap and the dust drifted from the outer disc. Another generation of planetesimals is formed at this location. A minor instability occurred between two newly formed massive cores at around 1.75 Myr that widens their radial separation.

Due to the late formation of the second generation of planetary cores, which ultimately form a pair of ice giants (10–100 M⊕), they remain in the thermal contraction phase of gas accretion at the end of the simulation at 2 Myr (vi). A compact chain of giant planets is produced spanning from 5 to 15 au, with a pair of gas giants formed from the initial pressure bump and a pair of ice giants formed over 1 Myr later from the edge of the planetary gap opened by the outer gas giant. The orbital periods of the inner pair are in near 2:1 commensurability and those of the outer pair are in near 4:3 commensurability.

Across the five random simulations, the final panel of Figs. 2 and 3 show very similar results for the gas giant pair formed in the first generation. For the next generation of planet formation, further stochasticity presents. Two (blue and green) simulations produce a pair of ice giants and another two (red and purple) produce only one ice giant. One simulation (orange) shows no ice giant at 2 Myr but a swarm of still-growing planetesimals and planet embryos.

thumbnail Fig. 2

Six key timestamps demonstrating sequential planet formation in one of the simulations with α = 5 × 10−4. Each panel shows the radial profiles of the gas surface density ∑g (solid line), the dust surface density ∑d (dashed line) and, the mass m and the semimajor axis r of the massive bodies (dot) at the noted time. The final panel also shows the massive bodies from the rest of the simulation set with each colour showing one of the five simulations.
thumbnail Fig. 3

Tracks of the massive bodies in the simulation with α = 5 × 10−4. Each colour shows one of the five simulations corresponding to those in the final panel of Fig. 2, where the key timestamps are also denoted along the time axis on the right. The solid lines show the semimajor axis of the bodies reached 10 M⊕ by the end of the simulations. The shaded areas indicate the extents of the lower and upper quartiles of the semimajor axes of all bodies when the total number is above 50.

3.1.2 Dust mass budget

Figure 4a shows the solid mass budget throughout the simulation presented by the colour blue in Figs. 2 and 3. The solid mass is divided into four categories, which are the dust mass Md inside and outside of 5 au, solids bound in massive bodies, and massive bodies ejected out of the simulation domain. The six timestamps in Fig. 2 are also denoted here on the time axis.

After the initial substructure is imposed, the dust mass between 3 au and 5 au (inner disc) decreases sharply due to the inward drift of dust, while the dust supply from 5 au and beyond (outer disc) is stopped at the pressure bump. At 0.05 Myr (i), the pressure bump is saturated with dust and the dust mass in the inner disc increases due to leakage by turbulent diffusion. Planetesimals start to form and grow by pebble accretion at 0.17 Myr (ii), which start converting dust into massive bodies. The conversion is paused when the first core starts to perturb the disc and stops pebble accretion.

Shortly after 0.2 Myr, there is a spike in the dust mass inside of 5 au. This is caused by the gap opening in a dust rich location as the first planetary core reaches the mass of 10 M⊕. Similarly, another spike occurs at about 0.34 Myr (iii) when the first core enters the runaway gas accretion phase but it is much smaller as dust is already depleted around the planet. Pebble accretion resumes for the second core at about 0.4 Myr after its orbit has been circularised, which also causes the change in the dust mass.

Planetesimal formation occurs again at 1.5 Myr (v) and pebble accretion continues to convert the remaining dust to massive bodies. The final (vi) masses of the four categories show that the majority of solids, or 85% of the initial dust mass beyond 5 au, are eventually incorporated into massive bodies.

Figure 4b shows the cumulative inflow of dust that crossed 5 au after 0.5 Myr, which is the time when the first gas giant has reached approximately one Jupiter mass. The subsequent dust inflow to the inner Solar System is, on average, 0.5–1 M⊕Myr−1 or, in total, about 1.6 M⊕ including two significant episodic inflows to the inner disc.

The first one occurs shortly at about 0.86 Myr (iv) as the second planetary core enters the runaway gas accretion phase and opens a gap in the disc. The dust near its location follows the sudden change in the gas disc and is pushed by the forming gas giant to both inner and outer part of the disc, which is also shown in the profiles of the surface densities (Fig. 2iv). The second one corresponds to a small instability occurs between the two newly formed massive cores at around 1.75 Myr and perturbs the disc (Fig. 3).

thumbnail Fig. 4

Evolution of the solid mass in the simulation with α = 5 × 10−4 that corresponds to the one presented by the colour blue in Figs. 2 and 3. The six timestamps in Fig. 2 are also denoted on the time axis. a) Solid mass budget. The solid mass is divided into four categories: the dust mass Md inside and outside of 5 au, solids bound in massive bodies and ejected massive bodies. It shows a high planet formation efficiency that the majority of solid mass (85% of the initial dust mass beyond 5 au) are converted into massive bodies. b) Cumulative inflow of dust entered 5 au after 0.5 Myr. A total of about 1.6 M⊕ of inflow to the inner disc is recorded over the next 1.5 Myr up to the end of the simulation.

3.2 The case of α = 3 × 10−4

Figure 5 presents the tracks of the massive bodies for the set of simulations with α = 3 × 10−4 in the manner of Fig. 3. Figure 6 presents the end results to the final panel of Fig. 2. The radial profiles of the surface densities are also shown for the whole set of simulations. Compared to the case of α = 5 × 10−4 in Sect. 3.1, a larger variation across the simulations is shown. For all simulations, planetesimal formation occurs at about 0.25 Myr, which is about 0.1 Myr later than the case of α = 5 × 10−4.

In the simulations denoted by the colour red in Fig. 5, a massive core is scattered through the migration trap by another core and is lost to the inner simulation boundary. Later at about 0.6 Myr, the next generation of planetesimals are formed resulting in two massive cores. Similarly, in the simulations denoted by the colour green, only one core is formed from the initial pressure bump and two is formed from the subsequent generation.

In the simulations denoted by the colour blue and purple, two cores are formed from the initial bump. The second generation of planet formation occurs at about 1.75 Myr for the former one, while only a concentrated dust ring presents at the outer edge of the planetary gap for the latter one at the end of the simulation. And, the simulation denoted by the colour orange forms three gas giants from the initial pressure bump and no further planet formation occurs before the simulation ends.

Figure 6 summarises the final architecture of the planets where three out of the five simulations form three gas giants and the remaining two (blue and purple) form two. Also, a significant dust ring remains respectively for all simulations external to their outermost gas giant.

thumbnail Fig. 5

Tracks of the massive bodies in the simulations with α = 3 × 10−4 presented in the same manner as in Fig. 3.
thumbnail Fig. 6

Results of the set of five simulations with α = 3 × 10−4. Similar to the final panel of Fig. 2 except the surface density profiles are shown for all simulations.

4 Discussions

4.1 Sequential planet formation

The results presented in Sec. 3 demonstrate a complete scenario of sequential planet formation. In the case of α = 5 × 10−4 (Sec. 3.1), two gas giants are formed from the initial disc substructure. As the outer gas giant has reached its final mass and a steady planetary gap is opened, dust re-accumulates at a later time near the new pressure bump triggering the next generation of planet formation. The case of α = 3 × 10−4 shows similar trend despite of the greater degree of stochasticity.

Comparing the results with the ‘inside-out’ planet formation scenario by Chatterjee & Tan (2013) for Kepler systems, we note that gas giant that has reached its final mass is more likely to trigger the next generation of planet formation. For low-mass planets that cannot open a significant gap in the disc (≲ 100 M⊕), dust leakage from the pressure bump at the outer edge of the gap is significant and requires a large supply from the outer disc to reach the conditions for planetesimal formation. Even if planetesimals may form, they are under greater gravitational influence of the planets as the width of the gap scales with m1/2 (Kanagawa et al. 2016; Duffell 2020) while the Hill radius scales with m1/3. In this case, the perturbation from the planet is more likely to prevent the growth of the planetesimals as pebble accretion, particularly when the planet continues to grow, is sensitive to the relative velocity as noted in Lau et al. (2022). These planetesimals will likely be scattered out of the system as well, if the planet enters the runaway gas accretion phase to become a gas giant. Therefore, the outer edge of the planetary gap opened by a steady gas giant is a much more favourable environment for the next generation of planet formation.

4.1.1 Architecture of the resulting systems

In the case of α = 5 × 10−4 (Sec. 3.1), the delay in the formation of the second generation of planets directly shortens the time available for their growth. Therefore, they remain at about 10 M⊕ by the end of the simulations. Although disc dissipation is not included in the current model, the sequential planet formation scenario provides the delay in formation time of the ice giants required by the models explaining their masses (e.g. Lee et al. 2014; Ogihara et al. 2020; Raorane et al. 2024). At the end of the simulations, the second generation planetesimals and embryos formed still remain in the system as the ice giants are not able to scatter them. This results in a system with diversity that consists of gas giants, ice giant(s) and small massive bodies.

Since the radial separation of the two generations of planets is determined by the gap width, the resulting system is compact with the four giants planets from 5 to 15 au. The gas giants formed from the initial substructure are also commonly in the 2:1 mean-motion resonance. While a further test with a larger sample size is required, this may justify the compact chain of giant planets adopted in the initial conditions of the Nice model (e.g. Tsiganis et al. 2005; Morbidelli et al. 2005) and the early instability model (e.g. Clement & Kaib 2017; Deienno et al. 2018) from the formation point of view.

Due to the computational cost, only two values of α are tested in this work where a larger degree of stochasticity is presented for the case of α = 3 × 10−4 (Sec. 3.2). Other than α, which determines the evolution timescale of the disc, the resulting planetary system should also be sensitive to the initial disc mass Mdisc, the characteristic radius rc and the location of the initial disc substructure. Upon the availability of TriPoD (Pfeil et al. submitted), which is a simplified three-parameter dust coagulation model, a more extensive parameter study shall become computationally feasible to study the diversity of planetary systems.

4.1.2 Dust mass budget

The resulting solid mass budget (Fig. 4a) shows a high planet formation efficiency for the presented simulation with α = 5 × 10−4. The common gap opened by the two gas giants is very effective in preventing dust from drifting through. As a result, the remaining dust is retained in the disc for a prolonged period of time and preserved solids for the subsequent planet formation.

After the formation of the first gas giant, the quasi-steady inflow and the episodic inflows result in a total of about 1.6 M⊕ of outer disc dust flowing into the inner disc over a time period of 1.5 Myr (Fig. 4b). The leftover planetesimals formed from the initial disc substructure are also generally scattered outward by the rapid formation of the first gas giant (Fig. 3). This indicates the chemical division caused by the first gas giant is robust. While further tests with the correct masses of the giant planets and their time of formation are required, this formation scenario may provide the required rapid formation of Jupiter’s core to prevent significant exchange of the non-carbonaceous and carbonaceous reservoirs in the early Solar System (Kruijer et al. 2017).

Stammler et al. (2023) studied the efficiency of dust trapping by planetary gaps corresponding to different planetary masses and values of α. Their results show a dust leakage rate of about 1 M⊕ Myr−1 for a gap created by a Saturn-mass planet with α = 10−4. This is broadly consistent with the presented leakage rate of about 0.5 M⊕ Myr−1 with the episodic inflows due to the dynamical instabilities excluded, since the planetary gap prescription used here is about 30% deeper than the one given by Kanagawa et al. (2016) (see Fig. 8 of Duffell 2020) and the planet is about three times more massive in our case.

The small inflow of dust from the outer disc is likely to be preferentially accreted by a proto-Earth and proto-Venus as pebble accretion is more efficient for bodies with higher mass and dynamically colder orbits. This may explain Earth’s chemical abundances relative to Mars and Vesta (Kleine et al. 2023) without causing a significant growth to form super-Earth. This result also suggests against the scenario of significant growth by pebble accretion in the inner Solar System (e.g. Johansen et al. 2021).

4.2 Comparison with the Solar System

We note that the results of α = 5 × 10−4 (Sec. 3.1) show a pair of Jupiter-mass giants, instead of one Jupiter- and one Saturn-mass giants. Although the second gas giant entered the runway gas accretion about 1 Myr later than the first one, there is no significant mass difference in the end. This is likely due to the absence of disc dissipation, which results in a high gas surface density throughout the simulation. Further developments of the model to include photoevaporation are required to provide a complete scenario of the formation of the early outer Solar System and planet formation in general. We anticipate that the remaining dust in the protoplanetary disc will form a belt of planetesimals resembling the scenario proposed by Carrera et al. (2017). In this case, these planetesimals will remain dynamically cold while growth by pebble accretion is prohibited due to the depletion of gas and dust. After disc dissipation, the N-body part (SyMBAp) can continue to model the long-term evolution of the system and show if a Nice model-like instability can occur among the giant planets.

The formation of the gas giants in the results is likely too quick compared to the meteoritic record (Kruijer et al. 2017, 2020). We note that the composition and the opacity of the envelope of gas giant, which are critical to the gas accretion rate, are still an active field of research (e.g. Szulágyi et al. 2016; Lambrechts et al. 2019; Schulik et al. 2019; Ormel et al. 2021). And, different gas accretion prescriptions are adopted among recent planet formation models (e.g. Liu et al. 2019; Bitsch et al. 2019; Matsumura et al. 2021; Chambers 2021; Lau et al. 2024). Further investigations on the different recipes and their consequences are required to match the formation history of the Solar System’s giant planets.

While the source of the initial disc substructure is not investigated in this work, the water ice line in the early Solar System has been proposed as a key feature in reproducing the Solar System by multiple works (e.g. Morbidelli et al. 2016, 2022; Brasser & Mojzsis 2020; Charnoz et al. 2021; Chambers 2023). Further investigations are required to determine the criteria at the water ice line in the early Solar System to trigger planet formation, particularly, the change in the surface density required.

4.3 Other recent works

4.3.1 Predictions of planet formation at pressure bumps

Xu & Wang (2024) made theoretical predictions on the architecture of the planetary systems assuming efficient planet formation at pressure bumps. They concluded three main pathways: slow core formation, fast core formation but slow gas accretion and, fast core formation and gas accretion.

While Lau et al. (2022) and Jiang & Ormel (2023) show that the high dust concentration at a pressure bump will likely favour rapid growth of core, the slow core formation can be possible if the planetesimals formed are dynamically heated to an extent that pebble accretion is halted but they still remain close to their birthplace.

The case of fast core formation but slow gas accretion is suggested to form a chain of super-Earths or potentially Saturn-mass planets over a prolonged period of time. This case is similar to the scenario proposed by Chatterjee & Tan (2013). As discussed in Sec. 4.1, the planetary gap formed by super-Earths may not be able to trap a significant amount of dust and trigger planetesimal formation. And, even if planetesimals could form, they are likely much closer to the core and cannot grow efficiently by pebble accretion.

The case of fast core formation and fast gas accretion is suggested to form a chain of gas giants. This prediction is the closest to the results shown in Sec. 3 while two to three gas giants can form from the pressure bump in the presented work. We note that the number of cores that can form from each pressure bump likely depends on the amplitude and the width of the dust trap, while this has not been tested in the presented work. This requires an extensive parameter study and more random simulations per set of parameters.

Although dust trap favours planetesimal formation, we emphasise that in situ planet formation at pressure bumps is unlikely a general solution to the diversity of exoplanetary systems. In particular, for the observed compact planet chains in resonance, a more probable formation scenario is that the cores are formed at a temporary pressure bump which migrate subsequently through the disc. Multiple works have studied the scenario that the inner most planet can be trapped at the disc edge (e.g. Terquem & Papaloizou 2007; Cossou et al. 2013; Brasser et al. 2018; Huang & Ormel 2023) and the subsequent inward-migrating planets can form a resonant chain through convergent inward migration (e.g. Tamayo et al. 2017; Delisle 2017; MacDonald & Dawson 2018; Wong & Lee 2024).

4.3.2 Sandwiched planet formation

With gas and dust hydrodynamics simulations, Pritchard et al. (2024) proposed the ‘sandwiched planet formation’ scenario where planet formation can occur with the dust trapped between two massive planets that each creates a pressure maximum. The authors already noted that formation of the planets and dust fragmentation are not modelled, which may have critical effect to the dust concentration between the planets. Furthermore, from the results presented above (Sec. 3), we also note that if planetesimals could form between the planets, they are likely under a great gravitational influence from the massive planets. This will prevent them from growing efficiently by accreting pebbles, or, more likely, scatter them. Nonetheless, their work confirms dust rings can be created by massive planets with hydrodynamics simulations and these are preferred locations of planet formation.

4.4 Caveats

4.4.1 Initial disc substructure

In this work, we studied the consequence of a substructure in the disc that can trigger planetesimal formation, with the location motivated by that of Jupiter. Although multiple non-planetary mechanisms are proposed in the field as discussed in Sec. 1, the parameter space, including the location, amplitude, width and lifetime, is not explored in this work and requires future investigations. For instance, from some test runs, we note that the amplitude of the substructure needs to be large enough to trap dust effectively and trigger planetesimal formation, although non-axisymmetric features that may aid dust concentration such as vortices (e.g. Barge & Sommeria 1995; Tanga et al. 1996) are not considered. Sequential planet formation also cannot occur if the dust mass remaining is not enough to trigger the next generation of planet formation. We emphasise that the criteria to form planetesimals are not trivial to satisfy in typical disc conditions. This work only focuses on a case where planetesimal formation is possible. However, by combining the unknowns in the shape and location of the initial disc substructure with different disc parameters, we expect the model can produce a variety of planetary systems through a parameter study, where sequential planet formation may not always occur.

4.4.2 Planet migration and gap opening

At the end of Sec. 2.4, we note that planetary gap opening and planet migration are treated by independent prescriptions while they are both the results of planet-disc interactions and are physically coupled. Also, multiple works (e.g. Lin & Papaloizou 1986; Armitage & Bonnell 2002; D’Angelo & Lubow 2010) have studied both effects consistently and provided formulae for the torque density profile exerted by a planet on the disc. However, upon applying the formula given by D’Angelo & Lubow (2010) in our 1-D model, we note that the gap opened by a Jupitermass planet is much narrower than that described in Duffell (2020), which is also given by Eq. (27). Since this gap profile is tested against a set of 2-D hydrodynamical simulations in a general parameter space and is consistent with other works (e.g. Kanagawa et al. 2015), we have opted to prescribe the gap profile and planet migration separately. We also note that a fixed temperature profile is adopted in this work, which implies that the effect of shock heating (e.g. Zhu et al. 2015; Rafikov 2016) is neglected while its effect is likely more significant in the outer disc where less irradiation from the star is received. Nonetheless, upon the availability of a general torque formula applicable to a 1-D model, this part of the model shall be modified for consistency, especially in the case of having multiple gap-opening planets in the disc.

5 Conclusions

This work demonstrates a scenario of sequential giant planet formation that is triggered by an initial disc substructure. We further extended the model in Lau et al. (2022) by including the effects of planetary gas accretion and gap opening. We employed DustPy to model a protoplanetary disc initially with micronsized dust, and SyMBAp was employed to model the evolution of the planetesimals upon formation.

Consistent with the previous results, planetary cores are formed rapidly from the initial disc substructure, which can then be retained at the migration trap and start gas accretion. The results show multiple (up to three) cores can form and grow into giant planets in each generation. As the first generation of gas giants has formed and opened a steady gap, the new pressure bump at the outer edge of the planetary gap becomes the next location of planet formation.

In the case of the higher value of α = 5 × 10−4, the second generation of planet formation occurs about 1 Myr after the first one, and only ice giants were formed instead of gas giants. This case also shows a high planet formation efficiency where more than 85% of the dust beyond 5 au is converted into massive bodies. As the first generation of gas giants effectively prevent dust from flowing through to reach the inner disc, the retained dust is then available for the next generation of planet formation. In the case of a lower value of α = 3 × 10−4, a larger degree of stochasticity was shown, while the general scenario of sequential giant planet formation remains. In both cases, a compact chain of giant planets are formed at the end of the simulations. While the simulations were stopped at 2 Myr, a natural continuation to the model would be to include the effect of photoevaporation to physically dissipate the disc and stop gas accretion.

Although the formation mechanisms of disc substructure are beyond the scope of this work, further investigations are required to study the possible shape and location produced by physical processes. It is unlikely that any disc substructure can trivially provide the conditions required for planetesimal or planet formation. Also, the parameter space and the number of random simulations in this work are limited by the computational costs. Further code optimisation is required to study the statistical effects and to model the diversity of planetary systems. And, planetary gas accretion is still an active field of research. Further investigations specifically on gas accretion are required to model the formation time of the Solar System’s giant planets.

Acknowledgements

We thank T. Kleine, C. Ormel, R. Teague and A. Vazan for the insightful discussions. We acknowledge funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme under grant agreement no. 714769, under the European Union’s Horizon Europe Research and Innovation Programme 101124282 (EARLYBIRD) and funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under grant 325594231 and Germany’s Excellence Strategy – EXC-2094 -– 390783311. J.D. was funded by the European Union under the European Union’s Horizon Europe Research & Innovation Programme 101040037 (PLANETOIDS). Views and opinions expressed are however those of the authors only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them.

References

  1. Abod, C. P., Simon, J. B., Li, R., et al. 2019, ApJ, 883, 192 [Google Scholar]
  2. Adachi, I., Hayashi, C., & Nakazawa, K. 1976, PThPh, 56, 1756 [NASA ADS] [Google Scholar]
  3. Andama, G., Ndugu, N., Anguma, S. K., & Jurua, E. 2022, MNRAS, 512, 5278 [Google Scholar]
  4. Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJ, 869, L41 [Google Scholar]
  5. Armitage, P. J., & Bonnell, I. A. 2002, MNRAS, 330, L11 [Google Scholar]
  6. Barge, P., & Sommeria, J. 1995, A&A, 295, L1 [NASA ADS] [Google Scholar]
  7. Bitsch, B., Lambrechts, M., & Johansen, A. 2015, A&A, 582 [Google Scholar]
  8. Bitsch, B., Izidoro, A., Johansen, A., et al. 2019, A&A, 623, A88 [Google Scholar]
  9. Blum, J., Gundlach, B., Krause, M., et al. 2017, MNRAS, 469, S755 [Google Scholar]
  10. Brasser, R., & Mojzsis, S. J. 2020, NatAs, 4, 492 [NASA ADS] [Google Scholar]
  11. Brasser, R., Matsumura, S., Muto, T., & Ida, S. 2018, ApJ, 864, L8 [Google Scholar]
  12. Brouwers, M. G., Ormel, C. W., Bonsor, A., & Vazan, A. 2021, A&A, 653 [Google Scholar]
  13. Carrera, D., Gorti, U., Johansen, A., & Davies, M. B. 2017, ApJ, 839, 16 [Google Scholar]
  14. Chambers, J. 2021, ApJ, 914, 102 [Google Scholar]
  15. Chambers, J. 2023, ApJ, 944, 127 [NASA ADS] [CrossRef] [Google Scholar]
  16. Charnoz, S., Avice, G., Hyodo, R., Pignatale, F. C., & Chaussidon, M. 2021, A&A, 652, A35 [Google Scholar]
  17. Chatterjee, S., & Tan, J. C. 2013, ApJ, 780, 53 [Google Scholar]
  18. Cieza, L. A., González-Ruilova, C., Hales, A. S., et al. 2021, MNRAS, 501, 2934 [Google Scholar]
  19. Clarke, C. J., & Pringle, J. E. 1988, MNRAS, 235, 365 [Google Scholar]
  20. Clement, M. S., & Kaib, N. A. 2017, Icarus, 288, 88 [NASA ADS] [CrossRef] [Google Scholar]
  21. Coleman, G. A. L. 2021, MNRAS, 506, 3596 [Google Scholar]
  22. Coleman, G. A. L., & Nelson, R. P. 2016, MNRAS, 460, 2779 [Google Scholar]
  23. Cossou, C., Raymond, S. N., & Pierens, A. 2013, A&A, 553, A2 [Google Scholar]
  24. D’Angelo, G., & Lubow, S. H. 2010, ApJ, 724, 730 [Google Scholar]
  25. Deienno, R., Izidoro, A., Morbidelli, A., et al. 2018, ApJ, 864, 50 [NASA ADS] [CrossRef] [Google Scholar]
  26. Delisle, J.-B. 2017, A&A, 605, A96 [Google Scholar]
  27. Drążkowska, J., Alibert, Y., & Moore, B. 2016, A&A, 594, A105 [Google Scholar]
  28. Drążkowska, J., Bitsch, B., Lambrechts, M., et al. 2023, Protostars Planets VII, eds. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, ASP Conf. Ser., 534, 717 [NASA ADS] [Google Scholar]
  29. Dubrulle, B., Morfill, G., & Sterzik, M. 1995, Icarus, 114, 237 [Google Scholar]
  30. Duffell, P. C. 2020, ApJ, 889, 16 [Google Scholar]
  31. Dullemond, C. P., & Penzlin, A. B. T. 2018, A&A, 609, A50 [Google Scholar]
  32. Dullemond, C. P., Birnstiel, T., Huang, J., et al. 2018, ApJ, 869, L46 [Google Scholar]
  33. Duncan, M. J., Levison, H. F., & Lee, M. H. 1998, AJ, 116, 2067 [Google Scholar]
  34. Eriksson, L. E. J., Mol Lous, M. A. S., Shibata, S., & Helled, R. 2023, MNRAS, 526, 4860 [NASA ADS] [CrossRef] [Google Scholar]
  35. Flock, M., Ruge, J. P., Dzyurkevich, N., et al. 2015, A&A, 574, A68 [Google Scholar]
  36. Gárate, M., Delage, T. N., Stadler, J., et al. 2021, A&A, 655 [Google Scholar]
  37. Gerbig, K., & Li, R. 2023, ApJ, 949, 81 [Google Scholar]
  38. Guilera, O. M., Sándor, Z., Ronco, M. P., Venturini, J., & Bertolami, M. M. M. 2020, A&A, 642, A140 [Google Scholar]
  39. Gupta, A., Miotello, A., Manara, C. F., et al. 2023, A&A, 670, A8 [Google Scholar]
  40. Huang, S., & Ormel, C. W. 2023, MNRAS, 522, 828 [NASA ADS] [CrossRef] [Google Scholar]
  41. Ida, S., Muto, T., Matsumura, S., & Brasser, R. 2020, MNRAS, 494, 5666 [Google Scholar]
  42. Jiang, H., & Ormel, C. W. 2023, MNRAS, 518, 3877 [Google Scholar]
  43. Johansen, A., Oishi, J. S., Low, M.-M. M., et al. 2007, Nature, 448, 1022 [Google Scholar]
  44. Johansen, A., Ronnet, T., Bizzarro, M., et al. 2021, Sci. Adv., 7, eabc0444 [Google Scholar]
  45. Kanagawa, K. D., Tanaka, H., Muto, T., Tanigawa, T., & Takeuchi, T. 2015, MNRAS, 448, 994 [Google Scholar]
  46. Kanagawa, K. D., Muto, T., Tanaka, H., et al. 2016, PASJ, 68, 43 [NASA ADS] [Google Scholar]
  47. Kanagawa, K. D., Tanaka, H., & Szuszkiewicz, E. 2018, ApJ, 861, 140 [Google Scholar]
  48. Klahr, H., & Schreiber, A. 2021, ApJ, 911, 9 [Google Scholar]
  49. Kleine, T., Steller, T., Burkhardt, C., & Nimmo, F. 2023, Icarus, 397, 115519 [NASA ADS] [CrossRef] [Google Scholar]
  50. Kobayashi, H., Ormel, C. W., & Ida, S. 2012, ApJ, 756, 70 [NASA ADS] [CrossRef] [Google Scholar]
  51. Kokubo, E., & Ida, S. 2000, Icarus, 143, 15 [Google Scholar]
  52. Kokubo, E., & Ida, S. 2002, ApJ, 581, 666 [Google Scholar]
  53. Kruijer, T. S., Burkhardt, C., Budde, G., & Kleine, T. 2017, Proc. Natl. Acad. Sci. U.S.A., 114, 6712 [Google Scholar]
  54. Kruijer, T. S., Kleine, T., & Borg, L. E. 2020, Nat. Astron., 4, 32 [Google Scholar]
  55. Lambrechts, M., Johansen, A., & Morbidelli, A. 2014, A&A, 572, A35 [Google Scholar]
  56. Lambrechts, M., Lega, E., Nelson, R. P., Crida, A., & Morbidelli, A. 2019, A&A, 630, A82 [Google Scholar]
  57. Lau, T. C. H., & Lee, M. H. 2023, RNAAS, 7, 74 [Google Scholar]
  58. Lau, T. C. H., Drążkowska, J., Stammler, S. M., Birnstiel, T., & Dullemond, C. P. 2022, A&A, 668, A170 [Google Scholar]
  59. Lau, T. C. H., Lee, M. H., Brasser, R., & Matsumura, S. 2024, A&A, 683, A204 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  60. Lee, E. J., Chiang, E., & Ormel, C. W. 2014, ApJ, 797, 95 [Google Scholar]
  61. Lin, D. N. C., & Papaloizou, J. 1986, ApJ, 309, 846 [Google Scholar]
  62. Lissauer, J. J. 1987, Icarus, 69, 249 [NASA ADS] [CrossRef] [Google Scholar]
  63. Liu, B., & Ormel, C. W. 2018, A&A, 615, A138 [Google Scholar]
  64. Liu, B., Ormel, C. W., & Johansen, A. 2019, A&A, 624, A114 [Google Scholar]
  65. Long, F., Pinilla, P., Herczeg, G. J., et al. 2018, ApJ, 869, 17 [Google Scholar]
  66. Lüst, R. 1952, Zeitschrift Naturforschung Teil A, 7, 87 [Google Scholar]
  67. Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603 [Google Scholar]
  68. MacDonald, M. G., & Dawson, R. I. 2018, AJ, 156, 228 [Google Scholar]
  69. Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425 [Google Scholar]
  70. Matsumura, S., Brasser, R., & Ida, S. 2017, A&A, 607, A67 [Google Scholar]
  71. Matsumura, S., Brasser, R., & Ida, S. 2021, A&A, 650, A116 [Google Scholar]
  72. Miller, E., Marino, S., Stammler, S. M., et al. 2021, MNRAS, 508, 5638 [Google Scholar]
  73. Morbidelli, A. 2020, A&A, 638, A1 [Google Scholar]
  74. Morbidelli, A., Levison, H. F., Tsiganis, K., & Gomes, R. 2005, Nature, 435, 462 [Google Scholar]
  75. Morbidelli, A., Bitsch, B., Crida, A., et al. 2016, Icarus, 267, 368 [Google Scholar]
  76. Morbidelli, A., Baillié, K., Batygin, K., et al. 2022, Nat. Astron., 6, 72 [Google Scholar]
  77. Nesvorný, D., Li, R., Youdin, A. N., Simon, J. B., & Grundy, W. M. 2019, Nat. Astron., 3, 808 [Google Scholar]
  78. Ogihara, M., Kunitomo, M., & Hori, Y. 2020, ApJ, 899, 91 [Google Scholar]
  79. Ormel, C. W., & Klahr, H. H. 2010, A&A, 520, A43 [Google Scholar]
  80. Ormel, C. W., & Liu, B. 2018, A&A, 615, A178 [Google Scholar]
  81. Ormel, C. W., Vazan, A., & Brouwers, M. G. 2021, A&A, 647, A175 [Google Scholar]
  82. Owen, J. E., Ercolano, B., Clarke, C. J., & Alexander, R. D. 2010, MNRAS, 401, 1415 [Google Scholar]
  83. Owen, J. E., Ercolano, B., & Clarke, C. J. 2011, MNRAS, 412, 13 [Google Scholar]
  84. Picogna, G., Ercolano, B., Owen, J. E., & Weber, M. L. 2019, MNRAS, 487, 691 [Google Scholar]
  85. Pinilla, P., Flock, M., Ovelar, M. d. J., & Birnstiel, T. 2016, A&A, 596, A81 [Google Scholar]
  86. Piso, A.-M. A., & Youdin, A. N. 2014, ApJ, 786, 21 [Google Scholar]
  87. Pritchard, M., Meru, F., Rowther, S., Armstrong, D., & Randall, K. 2024, MNRAS, 528, 6538 [NASA ADS] [CrossRef] [Google Scholar]
  88. Rafikov, R. R. 2016, ApJ, 831, 122 [Google Scholar]
  89. Raorane, A., Brasser, R., Matsumura, S., et al. 2024, Icarus, submitted [arXiv:2403.17771] [Google Scholar]
  90. Saito, E., & Sirono, S.-i. 2011, ApJ, 728, 20 [NASA ADS] [CrossRef] [Google Scholar]
  91. Schoonenberg, D., Ormel, C. W., & Krijt, S. 2018, A&A, 620, A134 [Google Scholar]
  92. Schreiber, A., & Klahr, H. 2018, ApJ, 861, 47 [Google Scholar]
  93. Schulik, M., Johansen, A., Bitsch, B., & Lega, E. 2019, A&A, 632, A118 [Google Scholar]
  94. Seager, S., Kuchner, M., Hier-Majumder, C. A., & Militzer, B. 2007, ApJ, 669, 1279 [Google Scholar]
  95. Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
  96. Stammler, S. M., & Birnstiel, T. 2022, ApJ, 935, 35 [Google Scholar]
  97. Stammler, S. M., Lichtenberg, T., Drążkowska, J., & Birnstiel, T. 2023, A&A, 670, L5 [Google Scholar]
  98. Szulágyi, J., Masset, F., Lega, E., et al. 2016, MNRAS, 460, 2853 [Google Scholar]
  99. Takahashi, S. Z., & Inutsuka, S.-i. 2014, ApJ, 794, 55 [NASA ADS] [CrossRef] [Google Scholar]
  100. Tamayo, D., Rein, H., Petrovich, C., & Murray, N. 2017, ApJ, 840, L19 [Google Scholar]
  101. Tanaka, H., Takeuchi, T., & Ward, W. R. 2002, ApJ, 565, 1257 [Google Scholar]
  102. Tanga, P., Babiano, A., Dubrulle, B., & Provenzale, A. 1996, Icarus, 121, 158 [Google Scholar]
  103. Tanigawa, T., & Tanaka, H. 2016, ApJ, 823, 48 [Google Scholar]
  104. Terquem, C., & Papaloizou, J. C. B. 2007, ApJ, 654, 1110 [Google Scholar]
  105. Toci, C., Rosotti, G., Lodato, G., Testi, L., & Trapman, L. 2021, MNRAS, 507, 818 [Google Scholar]
  106. Tsiganis, K., Gomes, R., Morbidelli, A., & Levison, H. F. 2005, Nature, 435, 459 [Google Scholar]
  107. Voelkel, O., Klahr, H., Mordasini, C., Emsenhuber, A., & Lenz, C. 2020, A&A, 642, A75 [Google Scholar]
  108. Weidenschilling, S. J. 1977, MNRAS, 180, 57 [Google Scholar]
  109. Whipple, F. L. 1972, in Plasma Planet, ed. A. Elvius, 211 [Google Scholar]
  110. Wong, K. H., & Lee, M. H. 2024, AJ, 167, 112 [Google Scholar]
  111. Xu, W., & Wang, S. 2024, ApJ, 962, L4 [NASA ADS] [CrossRef] [Google Scholar]
  112. Zhu, Z., Dong, R., Stone, J. M., & Rafikov, R. R. 2015, ApJ, 813, 88 [Google Scholar]
  113. Zormpas, A., Birnstiel, T., Rosotti, G. P., & Andrews, S. M. 2022, A&A, 661, A66 [Google Scholar]

All Figures

thumbnail Fig. 1

Differential mass distribution of the planetesimals formed near the initial disc substructure for the set of simulations with α = 5 × 10−4 drawn from the adopted initial mass function by Gerbig & Li (2023). There are 10 bins per decade in mass m and ΔN is the number of planetesimals in each bin. Each colour corresponds to one of the five simulations.
In the text
thumbnail Fig. 2

Six key timestamps demonstrating sequential planet formation in one of the simulations with α = 5 × 10−4. Each panel shows the radial profiles of the gas surface density ∑g (solid line), the dust surface density ∑d (dashed line) and, the mass m and the semimajor axis r of the massive bodies (dot) at the noted time. The final panel also shows the massive bodies from the rest of the simulation set with each colour showing one of the five simulations.
In the text
thumbnail Fig. 3

Tracks of the massive bodies in the simulation with α = 5 × 10−4. Each colour shows one of the five simulations corresponding to those in the final panel of Fig. 2, where the key timestamps are also denoted along the time axis on the right. The solid lines show the semimajor axis of the bodies reached 10 M⊕ by the end of the simulations. The shaded areas indicate the extents of the lower and upper quartiles of the semimajor axes of all bodies when the total number is above 50.
In the text
thumbnail Fig. 4

Evolution of the solid mass in the simulation with α = 5 × 10−4 that corresponds to the one presented by the colour blue in Figs. 2 and 3. The six timestamps in Fig. 2 are also denoted on the time axis. a) Solid mass budget. The solid mass is divided into four categories: the dust mass Md inside and outside of 5 au, solids bound in massive bodies and ejected massive bodies. It shows a high planet formation efficiency that the majority of solid mass (85% of the initial dust mass beyond 5 au) are converted into massive bodies. b) Cumulative inflow of dust entered 5 au after 0.5 Myr. A total of about 1.6 M⊕ of inflow to the inner disc is recorded over the next 1.5 Myr up to the end of the simulation.
In the text
thumbnail Fig. 5

Tracks of the massive bodies in the simulations with α = 3 × 10−4 presented in the same manner as in Fig. 3.
In the text
thumbnail Fig. 6

Results of the set of five simulations with α = 3 × 10−4. Similar to the final panel of Fig. 2 except the surface density profiles are shown for all simulations.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.