© O. Voelkel et al. 2022

1 Introduction

The latest observational constraints on grain growth in young protostellar disks via thermal dust emission (Harsono et al. 2018) imply that the formation of planets may begin in the earliest embedded phases in the life of young stars. The idea that previous generations of gas giant planets may have been accreted by the host star as a result of inward migration has already been introduced by Lin & Papaloizou (1986). This hypothesis states that the final system of planets that can be observed around a star may only reflect a small subset of the planets that initially formed. The possibility of a previous protogiant planet in the solar system is mentioned as well, however, the topic of the possible existence of previous super-Earths or other terrestrial mass planets has not been considered thus far. The number of planets that form and survive during the lifetime of a circumstellar disk is unknown. It is only the minimum number of survived planets per system that stands as a lower constraint, as it is given as the number of exoplanet detections. While this number often lacks completeness due to the low detectability of low-mass planets in certain systems, it completely lacks information on the previous history of the system. Another recent topic of discussion is the notion that planet-bearing stars might be polluted and show higher metalicities (Gonzalez 1997; Murray & Chaboyer 2002). It is thus possible that the currently observed population of planets merely reflects a small fraction of the planets that had initially been formed. While the research conducted on how individual planets grow and evolve that is based on an initially placed embryo continues to flourish, the preliminary formation of the used embryo is typically an initial assumption. As recently shown in Schlecker et al. (2021), this initial embryo location is the initial condition with the highest predictive power with regard to the resulting formation of a planet. The question on how many embryos form and how many of those survive cannot be neglected, however, when we are looking to study the formation of planets in a consistent fashion. An essential foundation of this approach is is a model that predicts the number of planetary embryos from the previously evolved system and tracks their combined evolution until the dispersal of the circumstellar disk. This study will present a self-consistent global model of planet formation and disk evolution that allows for such a study. The presented model enables us to investigate the number of planets that form and evolve within a circumstellar disk over its entire lifetime. We find that the formation history of planetary systems is far more complex than an initial distribution of planetary embryos could reflect.

To clarify the terminology, here we give a brief overview on current models of planet formation, focusing on their differences, similarities, and limitations. One similarity that all planet formation models bring into the picture is that the approach aims to combine multiple physical processes in a common framework. This results from the complexity of the problem, as the formation of planets cannot be described in an isolated way. Not only does the process range from a dust grain to a gas giant over numerous orders of magnitude in mass, it also needs to be embedded in the global evolution of a circumstellar disk. In addition to the evolution of the disk itself, interactions with other simultaneously growing planets can influence planet formation. Planet-planet interaction, as well as planet―disk interactions, can decide the fate of a planet during its formation and its later evolutionary stages.

Current global models of planet formation focus either on a specific time in the course of planetary formation, a specific accretion mechanism (i.e., pebble accretion or planetesimal accretion), or a specific location in the disk. Models introduced by Ida & Lin (2004), Alibert et al. (2005), Mordasini et al. (2012b) and Emsenhuber et al. (2021a) focus on the accretion of planetesimals on initially placed planetary cores. Even though the size of these planetesimals has its own ongoing field of research, here we refer to planetesimals as objects in the size range from 600 m in diameter (Emsenhuber et al. 2021a) to 100km in diameter Voelkel et al. (2020). Their size plays a major role in the accretion mechanism, as more massive objects are less likely to be accreted and the stirring by a protoplanet increases their eccentricities and inclinations. Even at sizes of several hundred meters to km however, the gas disk can significantly damp the planetesimal dynamical states and make them a highly efficient mechanism for protoplanetary growth (Emsenhuber et al. 2021a). The accretion of smaller particles for which gas drag can cause the object to even spiral onto the accreted protoplanet is called pebble accretion (Ormel & Klahr 2010). Planet formation models that are built around the accretion of pebbles onto initially placed planetary embryos have been introduced in Bitsch et al. (2015), Ndugu et al. (2017) and Jungo et al. (2020). The aforementioned planet formation models either focus on the accretion of planetesimals or the accretion of pebbles. Hybrid accretion models were recently introduced by, for instance, Alibert et al. (2018) or Guilera et al. (2020). A major drawback of these models, however (which also pertains to models aimed at the study of pebble or planetesimal accretion in an isolated fashion), is the initial placement of planetary embryos. This initial assumptions skips the earliest phase of circumstellar disk evolution whereby planetesimals form and later accumulate to compose the planetary embryos. As the results from Bitsch et al. (2015) show, the location and the time when an embryo is placed plays a dominant role in the subsequent evolution. The recent work presented in Voelkel et al. (2021a) and Voelkel et al. (2021b) studies the formation of planetary embryos from plan-etesimals that form from an evolving pebble disk. They find that more distant embryos (>2−3 au) form after the pebble flux has largely vanished. While the accretion of pebbles on planetesi-mals and planetary embryos is included in Voelkel et al. (2021b), only the innermost planetary embryos can benefit from pebble accretion, as the outer embryos fail to form during the lifetime of the pebble flux.

A recent work by Guilera et al. (2020) studied the formation of giant planets around pressure bumps. They combined a global disk evolution model containing gas, dust, and pebble dynamics with the formation of planetesimals due to streaming instability behind the pressure bumps. The embryo was included once the mass in planetesimals became equivalent to the mass of a planetary embryo. They also used a hybrid accretion model that combines pebble and planetesimal accretion, as well as a global disk evolution model. While in their extensive model, they did not use the embryo as an initial assumption, the embryo was placed in a specific location and subject to an initial assumption as well. The formation time of the embryo in their work is given by the time it takes until a lunar mass of 100 km planetes-imals has formed around their pressure bump. While this is a first constraint on the initial placement time, it does not account for planetesimal growth by planetesimal collisions, which can take significantly longer than it takes to form the 100 km planetesimals themselves. Additionally, their placed embryo was the only embryo in the system, which neglects planet–planet interactions. While this model attempts to model all stages of the single planet in the system (beginning from dust and pebbles to a final gas giant) using a global disk evolution model, the planet formation studied remains local, as the location and the number of planets in the system remains fixed.

A planet formation model based on planetesimal accretion that forms planetesimals consistent with the radial evolution of dust and pebbles has also been presented in Voelkel et al. (2020). This approach connects a two-population model for dust and pebble dynamics (Birnstiel et al. 2012) with the pebble flux regulated model for planetesimal formation from Lenz et al. (2019). The evolution of dust, pebbles, and planetesimals was merged with the planet formation model from Emsenhuber et al. (2021a) to study the impact of different planetesimal distributions on planet formation. While planetesimals were formed consistently with the disk evolution, pebble accretion was neglected and planetary embryos remained an initial assumption. The formation model of Emsenhuber et al. (2021a), however, is capable of also tracking the growth and N-body dynamics and interactions of up to 100 planetary embryos. A model that forms planetary embryos based on the local planetesimal surface density evolution is presented and discussed in Voelkel et al. (2021a) and Voelkel et al. (2021b). To bridge the gap between disk evolution and the accretion of pebbles and planetesimals onto planetary embryos, we decided to implement the embryo formation model from Voelkel et al. (2021a) into the planet formation model described in Voelkel et al. (2020) and model the accretion of both pebbles and planetesimals.

Here, we highlight that the planet formation model presented in this work combines the currently known accretion mechanisms (pebble accretion, planetesimal accretion, gas accretion) during the entire lifetime of the circumstellar disk, which includes the formation of planetary embryos, as well as late stage gas accretion. The number, as well as the formation time and location of planetary embryos are no longer an assumption, but the result of the analytic embryo formation model from Voelkel et al. (2021a). This stands in contrast to the models of Emsenhuber et al. (2021b) where a fixed number of embryos (1, 20, 50, 100) is inserted at t = 0 throughout the disk uniformly in log. In the following, we discuss the individual stages of the planet formation model used in this study. A detailed description of the different stages of planet formation that are covered in the model presented in this paper are described in Sect. 2.

2 Our global model of planet formation

In the following, we discuss the different stages of planet formation that are covered in our global formation model. As a reminder, the computation of the disk evolution, the accretion of material, the formation of embryos, and planetary migration are computed simultaneously. The existence of a planet changes the pebble flux due to pebble accretion. The formation of planetesimals (as it is regulated by the pebble flux) changes accordingly and, thus, the formation of other planetary embryos is affected as well.

2.1 Disk evolution model

A detailed description on the implementation of the gas and solid disk evolution model used in our planet formation framework can be found in Voelkel et al. (2020). Here, we focus on the discussion of the underlying fundamentals.

We used a one-dimensional disk evolution model that tracks the evolution of gas, dust, pebbles, and planetesimals. The viscously evolving gas disk (Lüst 1952; Lynden-Bell & Pringle 1974) uses an α-prescription for turbulence (Shakura & Sunyaev 1973) and includes internal and external photoevaporation (Picogna et al. 2019).

The inner edge of the disk is modeled as in Alibert et al. (2005), setting the surface density at the boundary to the floor value of the gas and the mass flux through the inner boundary equal to the equilibrium flux. Coupled to the evolution of the gas disk is the two-population solid evolution setup from Birnstiel et al. (2012). This model solves an advection diffusion equation of a combined solid density. Depending on whether the particles at a given radial location can be considered in the drift or in the fragmentation limit of growth, a fixed mass relation is applied. This relation splits the solid density into two populations. Depending on the individual Stokes number of the particles, these populations can be considered as dust (St ≪ 1) or pebbles (St ≥ 1), respectively.

The formation of planetesimals in our framework is regulated by the local radial pebble flux. The model we use was introduced by Lenz et al. (2019) and it does not specify which physical mechanism (e.g., Kelvin Helmholtz instability or streaming instability) drives the formation of planetesimals. Its underlying assumption is that planetesimals form in trapping zones that can appear at any location of the disk. These trapping zones appear for a given lifetime and with a radial separation of d(r). In Lenz et al. (2019), planetesimals form proportional to the radial pebble flux and the formation rate of planetesimals is given as: $${\dot{\Sigma }}_{\text{p}}\left(r\right)=\frac{\in }{d\left(r\right)}\frac{{\dot{M}}_{\text{peb}}\left(r\right)}{2\pi r},$$(1)

with Σ_p(r) as the local planetesimal surface density at a heliocentric distance, r, and, M_peb(r), as the local radial pebble flux; ϵ describes the amount of the pebble flux that is transformed into planetesimals over the trap distance d(r). The distance of pebble traps is given as five gas pressure scale heights in our approach, as in Lenz et al. (2019). The reason for this assumption stems from the typical separation of zonal flows found in Dittrich et al. (2013). Planetesimals are assumed to be in the oligarchic regime (Ida & Makino 1993; Thommes et al. 2003; Chambers 2006).

They are described in a fluid type fashion using a surface density, Σ_P, with eccentricity and inclination. For their dynamical state, we use the approach from Fortier et al. (2013). Planetesimals are stirred by the embryos, as well as by each other and damped by the gas disk. This stirring by the protoplanet follows Guilera et al. (2010) while the planetesimal planetesimal stirring follows Ohtsuki et al. (2002). The damping of plan-etesimals follows Inaba et al. (2001) in the quadratic regime, Adachi et al. (1976) and Rafikov (2004) in the Stokes and Epstein regime.

The size at which planetesimals form is given as 100 km in diameter. While the size of planetesimals is an ongoing field of research, we choose a size of 100 km in diameter, as observational constraints from the solar systems infer (Bottke Jr et al. 2005; Walsh et al. 2017; Delbo' et al. 2017) and what numerical simulations suggest (Schäfer et al. 2017; Klahr & Schreiber 2020). Other works suggest smaller sizes, in the range of several 100m to kilometres in diameter (Arimatsu et al. 2019; Schlichting et al. 2013; Weidenschilling 2011; Zheng et al. 2017).

2.2 Planetesimals to planetary embryos

Planetesimals are described as a one-dimensional (1D) surface density (Σ_p), analogous to gas, dust, and pebbles. While we did not track the N-body evolution of this large number of planetesimals, we tracked the dynamical N-body evolution of up to 100 planetary embryos. These embryos are introduced over time into the simulation, consistent with the evolution of the planetesimal surface density and its dynamical state. The embryo formation model that we use was introduced by Voelkel et al. (2021a) and we briefly describe it in the following.

Once planetesimals begin to form, we track their growth by integrating the local mass growth rate of a planetesimal in the oligarchic regime within a swarm of planetesimals (Lissauer 1993): $$\frac{\text{d}{M}_{\text{p}}\left(r,t\right)}{\text{d}t}=\frac{\sqrt{3}}{2}{\Sigma }_{\text{P}}\left(r,t\right)\cdot \text{Ω}\left(r\right)\pi r_{\text{b}}{}^2\left(1+\frac{{\upsilon }_{\text{esc}}^2\left(M_{\text{p}},R_{\text{b}}\right)}{{\upsilon }_{\infty }^2\left(r,t\right)}\right),$$(2)

with M_p(r, t) as the mass of the largest object at a heliocentric distance, r, at a time, t; Σ_P(r, t) is given as the local planetesimal surface density; Ω(r) as the orbital Kepler frequency; R_b the radius of the largest object; υ_esc as the escape velocity of M_p at its surface; and υ_∞ as the dispersion velocity of planetesimals, which we give as υ_∞ = e(r) Ω(r) and with e(r) as the planetesimals eccentricity; M_p is initially set to the mass of a 100 km planetesimal with a solid density of ρ_s = 1.0 g cm⁻³: $$M_{\text{p}}\left(r,t_0\right)=M_{100\text{\hspace{0.17em}km}}.$$(3)

Once Mp locally surpasses the mass of a planetary embryo, which in our case, is given as a lunar mass (M_emb = 0.0123 M_⊕), a new N-body object is introduced into the simulation and Mp is reset to M_p(r, t) = M_100km, within 15 Hill radii of the placed embryo. An additional constraint from Voelkel et al. (2021a) is that an embryo cannot form within 15R_Hill of any other embryo in the system. The orbital separation of planetary embryos in the oligarchic growth regime has already been found and confirmed in Kokubo & Ida (1998); Kobayashi et al. (2011) and Walsh & Levison (2019). It serves as a good constraint on the number of planetary embryos within a given spatial distribution. While the model has been derived without including the effect of pebble accretion, the effect of pebble accretion on embryo formation was studied in Voelkel et al. (2021b) in a similar framework, including pebble accretion. It has been shown that the accretion of pebbles largely affects the mass growth rate of planetary embryos with masses >0.01 M_⊕ and, thus, their physical spacing to each other. When expressed in the embryo Hill radii, however, the orbital separations remain similar to what has been found in Voelkel et al. (2021a); Kokubo & Ida (1998); Kobayashi et al. (2011); Walsh & Levison (2019). Pebble accretion only begins to be an effective accretion mechanism at masses much larger than that of a 100 km planetesimal. The initial planetesimal growth from 100 km to a lunar mass object therefore remains dominated by planetesimal collisions. The initial formation time of a lunar mass object is therefore only weakly influenced by pebble accretion (Voelkel et al. 2021b). Its subsequent growth however begins to be dominated by pebble accretion. Including the embryo formation model from Voelkel et al. (2021a) into our planet formation framework therefore ensures that the total number of embryos in the system, their formation time and their spatial distribution is consistent with the evolution of the planetesimal surface density and its dynamical state. The initial eccentricity of an embryo is chosen randomly between e = 10⁻³ and e = 10⁻⁵. The initial inclination is given as i = e/2.

2.3 Embryos and beyond

Once a planetary embryo has formed (as described in Sect. 2.2), it is subject to several simultaneously occurring processes. Its mass growth is given by the accretion of pebbles, planetesimals, and gas (Pollack et al. 1996). To avoid confusion, there is no physical difference in our model between a planetary embryo and a planet. The terminology of an embryo only refers to the object at its initial lunar mass of 0.0123 M_⊕. It is then treated as a single N-body object and is hereafter referred to as planet. Every planet that formed in our model was initially introduced to the systems as a lunar mass embryo based on the model described in Sect. 2.2. During its growth, the planet is subject to planetary migration and the dynamical interaction with other planets in the system. Planets can be scattered out of the system if they reach an orbital distance to the host star of above 1000 au. Mergers with other simultaneously forming planets are included as well, if the orbital distance between two planets is less than the sum of their radii. The density of the core is computed as in Mordasini et al. (2012a) by using a modified polytropic EOS from Seager et al. (2007).

As mentioned, planetesimals are considered to be in the oligarchic regime. They are accreted by embryos and evolve their eccentricity and inclination by self-stirring, by the interaction with the embryo, and by the damping of the gas disk (see Sect. 2.1). Next to the accretion of planetesimals, we included the accretion of pebbles from the disk, based on the prescription of Ormel (2017). The accretion of pebbles is considered until the planet reaches its local pebble isolation mass (Lambrechts et al. 2014). Once planets are large enough to accrete gas from their surrounding, the 1D structure of the gas envelope is retrieved by solving the internal structure equations (Bodenheimer & Pollack 1986). The accretion of solids affects the accretion of gas during the initial phase via accretional heating (Pollack et al. 1996; Lee & Chiang 2015; Alibert et al. 2018). As seen in Voelkel et al. (2020), this effect can suppress runaway gas accretion for high planetesimal surface densities, as runaway gas accretion can only occur if the accretion rate of gas surpasses the accretion rate of planetesimals (Pollack et al. 1996).

As these planets that are embedded in the gas disk grow, they can undergo type I and type II migration. Here, type I migration is treated as described in Coleman & Nelson (2014), type II migration as in Dittkrist et al. (2014) – and to distinguish between them we use the prescription by Crida et al. (2006) for gap opening. As discussed in Coleman & Nelson (2014), the formation of strong corotation torques can also lead to type I outward migration within our framework.

The total type I torque exerted onto the planet is given as (Paardekooper et al. 2011; Coleman & Nelson 2014): $${\Gamma }_1=F_{\text{L}}{\Gamma }_{\text{L}}+F_{\text{e}}{F}_i\left({\Gamma }_{\text{c,baro}}+{\Gamma }_{\text{c,ent}}\right)\text{,}$$(4)

with Γ_L as the Lindblad torque and Γ_c,baro, and Γ_c,ent as the barotropic and entropic corotation torque. The reduction of the Lindblad torque caused by the eccentricity and inclination of the planets orbit follows Cresswell & Nelson (2008) and is given as: $$F_{\text{L}}^{-1}=P_e+\left(\frac{P_e}{\left|P_e\right|}\right)\left(0.07\widehat{i}+0.085{\widehat{i}}^4-0.08\widehat{e}{\widehat{i}}^2\right)\text{,}$$(5)

with $$P_e=\frac{1+{\left(\frac{\widehat{e}}{2.25}\right)}^{1/2}+{\left(\frac{\widehat{e}}{2.84}\right)}^6}{1-{\left(\frac{\widehat{e}}{2.02}\right)}^4}\text{,}$$(6)

using ê = e/h_asp and î = i/h_asp with h_asp = h/r as the aspect ratio of the disk; F_e and F_i in Eq. (4) denote the reduction of the corotation torque as caused by the planets eccentricity and inclination (Bitsch & Kley 2010). They are given as: $$F_e=\exp \left(-\frac{e}{h_{\text{asp}}/2+0.01}\right)\text{,}$$(7) $$F_i=1-\tanh \widehat{i}\text{,}$$(8)

see Fendyke & Nelson (2014) for Eq. (7) and Coleman & Nelson (2014) for Eq. (8). The damping of the planets eccentricity and inclination is included as in Cresswell & Nelson (2008).

The vertical structure of gas disk is computed at each time step of the disks evolution using the approach of Nakamoto & Nakagawa (1994). The vertical scale height is then given as h = c_s/Ω with $c_{\text{s}}=\sqrt{k_{\text{b}}{T}_{\text{mid}}/\left(\left(\mu m_{\text{H}}\right)\right)}$ (using μ = 2.24 as the molecular weight and m_H as the hydrogen atom mass). The midplane temperature, T_mid, is given as: $${\sigma }_{\text{SB}}{T}_{\text{mid}}^4=\frac{1}{2}\left(\frac{3}{8}{\tau }_{\text{R}}+\frac{1}{2{\tau }_{\text{P}}}\right)\dot{E}+{\sigma }_{SB}{T}_{\text{s}}^4\text{,}$$(9)

with T_s as the irradiation temperature, σ_SB as the Stefan-Boltzmann constant, τ_R and τ_P as the Rosseland and Planck mean optical depths, and Ė as the viscous dissipation rate. The entire framework is described in great detail in Emsenhuber et al. (2021a); here, we only give a brief overview of the related physics. Strictly speaking, Emsenhuber et al. (2021a) refer to a model for planet population synthesis that contains a planet formation model at its heart. In this sense, Emsenhuber et al. (2021a) is an updated version of Mordasini (2018), in which a model for planet formation (Alibert et al. 2005, 2013) and planet evolution (Mordasini et al. 2012b) are combined to carry out a planet population synthesis approach, as in Mordasini et al. (2009). As we are not carrying out a population synthesis in this study, focusing, rather, on the investigation of a single system, we do not make use of the population synthesis capabilities of our used framework. However, we do wish to mention here that the entire framework presented in this paper is capable of conducting the same population synthesis studies as presented, for instance, in Emsenhuber et al. (2021b) or Schlecker et al. (2021), since our additional physical models (pebble and dust dynamics, pebble accretion, planetesimal formation, and planetary embryo formation) do not require high computational costs.

3 Numerical setup and initial conditions

We focus on a set of disk parameters introduced in Lenz et al. (2020) (see Table 1). Their goal was to constrain the parameters of the solar nebula by conducting a broad parameter study. The resulting parameters led to what Lenz et al. (2020) refer to as the “most appealing solar nebula” (MASN). The distribution of planetesimals that resulted from the aforementioned parameters performed optimally in constraining the solar nebula, based on the distribution of planets and asteroids in the solar system today.

4 Simulation results

We investigate the effect of dynamic embryo formation on the formation of a planetary system. While this extensive model allows for a multitude of effects to be investigated in greater detail, our study focuses on the composition, mass, and semimajor axis evolution of the resulting planetary system (Sect. 4.1), the number of active and formed planets over time (Sect. 4.2), the evolution of the disk surface densities and masses (Sect. 4.3), the evolution of the solid mass components (Sect. 4.4), and the final system of planets (Sect. 4.5).

4.1 Planetary system evolution

Figure 1 shows the mass, pebble mass fraction, and semi-major axis over time of all planets during the lifetime of the gas disk. Figure 2 shows the mass and semimajor axis distribution of the system, including their growth tracks and the initial formation time of the corresponding planetary embryo. Figure 3 shows the evolution of the corresponding planet masses. In Fig. 1, within the first 0.35 Myr, we find that several embryos form within 1 au and grow primarily through pebble accretion, as shown by the color of the dots. As seen in Fig. 3, the most massive planet reaches up to 20 M_⊕. These early formed planets experience strong type I migration and eventually end at the inner edge of the gas disk. As more planets migrate to the inner edge, the innermost planets are accreted by the host star, as shown by the triangular markers in Fig. 2. The same evolution is undertaken by the next set of planets that forms after 0.35 Myr. Within the first 1 Myr, we see that several new embryos form in the terrestrial planet zone. Initially dominated by pebble accretion, those super-Earth-mass planets first experience outward migration due to positive corotation torques (Coleman & Nelson 2014), followed by inward migration before halting at the inner edge of the gas disk again. Those planets form within 1.2 au but temporarily reach a semimajor axis of 2–3 au. Within the first 3.5 Myr, we see that one planet does not migrate to the inner edge of the disk within the first 1 Myr, but over a significantly longer timescale. As it can be seen by the pebble mass fraction evolution of the outermost planet at 1 Myr, it was initially dominated by pebble accretion. Over the course of the next 2.5 Myr, its pebble mass fraction strongly decreases due to ongoing planetesimal accretion. During that phase, a set of sub-Earth-mass planets has formed in the area around 1 au. After the outer super-Earth goes into another phase of type I inward migration, it pushes the sub-Earth-mass planets to the inner edge of the gas disk as well, eventually clearing the terrestrial planet zone of planets. The migration of planets beyond the inner edge of the gas disk is due to N-body interactions with other planets. This involves resonant chains and scattering up to the point of reaching the surface of the host star. After the last super-Earth has migrated to the inner edge of the gas disk at 3.5 Myr, a large set of sub-Earth-mass planets emerges over the next 13 Myr. Those planets stay at smaller masses than their super-Earth predecessors for the lifetime of the gas disk and consequentially experience significantly slower type I migration. As seen in Fig. 3, the two most massive planets of the remaining system have formed at 0.35 Myr as part of the second generation of embryo formation and no planet from the first generation of embryo formation survived beyond the dispersal of the gas disk.

Fig. 1 Semimajor axis over time evolution for the planetary system. The evolution of a planet is linked via the grey line and the size of the dots indicate the mass of the planets every 10 ky. The color of the planets indicates their pebble mass fraction Mpeb/MP. The horizontal lines are drawn at 0.35 Myr, 1 Myr, and 3.5 Myr and show the moments at which most currently active planets are accreted by the host star or were subject to mergers. Thus we find four distinct generations of planet formation.

4.2 Number of planets over time

In Fig. 4, we show the total number of active planets and the cumulative number of planets that formed during the lifetime of the gas disk. As the number of formed planets continuously increases, the number of active planets shows three significant moments of decrease. The first decrease can be found at 0.35 Myr, the second at 1 Myr and the third at 3.5 Myr. The local minima of active planets in the disk are indicated by the vertical lines. The total number of planets that formed during the lifetime of the disk is given as 78, whereas the total number of active planets after the lifetime of the disk is given as 16. The largest number of planets formed in the last generation after 3.5 Myr. While the number of formed planets keeps increasing after 5 Myr, the number of active planets after that time remains almost constant at ~30 due to giant impacts and mergers. In the latest stages after 16 Myr, the number of active planets drops to 16, as planets continue to collide and merge, but no more embryos are forming.

Fig. 2 Mass and semimajor axis distribution of the entire planetary system until the dispersal of the gas disk. The large circle markers indicate planets that remain active until the end of the simulation. Planets that were accreted by the host star are shown as triangles and planets that merged via collisions with other planets are indicated as squares. The track of the planets is shown as the solid grey line for active planets and dotted grey lines for accreted or merged planets. The color of the final marker indicates the initial formation time of the corresponding planetary embryo.

Fig. 3 Planet mass over time. The vertical lines are drawn at 0.35 Myr, 1 Myr, and 3.5 Myr and show the moments at which most currently active planets are accreted by the host star or were subject to mergers. The vertical steps are caused by giant impacts that result in planet mergers. The large circle markers indicate planets that remain active until the end of the simulation. Planets that were accreted by the host star are shown as triangles and planets that merged via collisions with other planets are indicated as squares. The track of the planets is shown as the solid grey line for active planets and dotted grey lines for accreted or merged planets. The color of the final marker indicates the pebble mass fraction of the corresponding planet.

Fig. 4 Number of planets over time. The blue line shows the number of currently active planets in the system, while the orange line shows the total number of planets that have formed. The vertical lines are drawn at 0.35 Myr, 1 Myr, and 3.5 Myr and show the moments at which most currently active planets are accreted by the host star or were subject to mergers.

4.3 Disk evolution

Figure 5 shows the surface densities and planetesimal eccentricity of the circumstellar disk at 0.35 Myr, 1 Myr, 3.5 Myr, and 16.8 Myr. The semimajor axes of the active planets in the system are displayed as dashed vertical lines. Figure 5 also shows the disk component masses at the various snapshots. The very long gas disk lifetime of this setup is due to the large initial gas disk mass (0.1 M_⊙) in combination with the small a = 3 × 10⁻⁴. As discussed in Lenz et al. (2019), a higher photoevaporation rate would not greatly influence the formation of planetesimals, as most planetesimals form within the first Myr of the system’s evolution. In order to stay consistent with Lenz et al. (2019), we chose to use the same parameters. However, a higher photoevaporation rate, meant to induce a shorter disk lifetime, would not affect the initial Myr of the system’s evolution. Within the first 1 Myr, the mass of the pebble and dust disk drops from an initial value of 450.93 M_⊙ to only 5.92 M_⊙. The mass of the planetesimal disk after 1 Myr is given as 73.61 M_⊙. The largest fraction of the dust and pebble disk is accreted by the host star due to continuous inward drift. The inner region of the circumstellar disk is largely depleted of planetesimals when the gas disk has vanished. After 16.8 Myr, we still find 67.59 M_⊙ of planetesimals in the entire disk, most of which are between 5 au and 10 au. As can be seen in any snapshot in which planets are present, the eccentricity of planetesimals greatly increases at the location of active planets. Once the planets have migrated, however, the planetesimals eccentricity is again reduced via damping from the gas disk.

Figure 6 shows the fraction of the planetesimal disk and the combined dust and pebble disk mass over the gas disk mass within 1 au, 2.5 au, and 5 au during the lifetime of the gas disk. As we know from Lenz et al. (2019), the planetesimal surface density profile is steeper than the surface density profile of the gas disk, due to the influx of pebbles from distant regions of the disk. As a consequence, we find that in the early phase of the evolution, the disk mass fraction of planetesimal mass over gas mass within a smaller region shows the highest value. After 200 ky, we find that the fraction of planetesimal mass over gas mass within 1 au is >1.5%, whereas the fraction of planetesimal mass over gas mass within 5 au only reaches 0.5% in that time. If it were not for the formation of planetary embryos, we would expect this trend to continue throughout the lifetime of the disk. After the formation of planetary embryos, however, the planetesimal disk mass within 1 au strongly varies as a consequence of planetesimal accretion and continuous planetesimal formation. As can be seen in Fig. 5, the planetesimal disk within 1 au experiences the greatest depletion due to planets, as most embryos within 1 Myr form within 1.5 au. The sharp increase in the mass fraction at later times is due to the depletion of the gas disk as a consequence of photoevaporation and accretion to the host star. The mass of the combined dust and pebble disk vastly exceeds the mass of the planetesimal disk for the first 0.35 Myr within 5 au, 2.5 au, and 1 au. As the planets that formed within the first 0.35 Myr also happened to form within 2.5 au, their mass growth is dominated by pebble accretion. Between 0.35 Myr and 1 Myr, the planetesimal disk mass exceeds the combined dust and pebble disk mass for every shown radius. As the total solid disk mass starts to become dominated by planetesimals, the accretion of pebbles is no longer the dominant mechanism of growth. Planets that form in the second generation begin to reduce their pebble mass fraction a a result of the depletion of the pebble reservoir.

Fig. 5 Surface density and planetesimal eccentricity evolution of the circumstellar disk. We show the gas surface density (blue), planetesimal surface density (black), and the combined dust and pebble surface density (orange) at t = 0.35 Myr, t = 1 Myr, t = 3.5 Myr, and t= 16.8 Myr. The initial corresponding density is shown as the dashed line. The initial mass of the gas disk is given as M(Σg) = 34102.64 M⊕, the initial solid mass in dust and pebbles is given as Μ(Σtp) = 450.93 M⊕, and the initial planetesimal disk mass is given as M(Σtp) = 0 M⊕. The location of an active planet is indicated via dashed vertical lines. The planetesimal eccentricity is given as the dotted line, respectively.

Fig. 6 Fraction of the planetesimal disk mass (MD(ΣP)) and combined dust and pebble disk mass (MD(Σtp)) over the gas disk mass (MD(Σg)) within 1 au, 2.5 au, and 5 au over time in percent. The grey dotted line indicates the global initial dust to gas ratio dg = 1.34%. The dashed vertical lines are drawn at 0.35 Myr, 1 Myr, and 3.5 Myr and show the moments at which most currently active planets are accreted by the host star or were subject to mergers. The total mass in the planetesimal disk after the dispersal of the gas disk is given as 67.59 M⊕. The final planetesimal disk mass within 5 au is given as 10.5 M⊕, within 2.5 au as 0.349 M⊕, and within 1 au as 2.35 × 10−3 M⊕.

4.4 Solid mass evolution

Figure 7 shows the evolution of the different solid mass components of the system in Fig. 1. This includes the total pebble and dust disk mass, the total planetesimal disk mass, the mass in all active planetary cores, the mass in all formed planetary cores, and the mass of all planetary cores that were accreted by the host star. We find that the mass of cores accreted by the host star is larger than the remaining mass in active cores at the end of the gas disks lifetime. The mass of the planetesimal disk surpasses the mass of the dust and pebble disk after ~0.35 Μyr. The mass of active planetary cores never surpasses the mass of the planetesimal disk; however, the mass of all cores that have formed does indeed surpass the mass of the planetesimal disk within 1 Myr. The planetesimal disk mass reaches its highest value after ~0.45 Μyr and then decreases. This is due to low planetesimal formation as a result of the largely depleted pebble and dust disk and planetesimal accretion onto planetary cores.

Fig. 7 Solid mass evolution during the lifetime of the gas disk. We show the total pebble and dust disk mass (MD(Σtp)), the total planetesimal disk mass (MD(ΣP)), the mass in active planetary cores (Mc(active)), the mass of all formed planetary cores (MC(formed)), and and the mass of planetary cores that have been accreted by the host star (MC(accreted)). The dashed vertical lines are drawn at 0.35 Myr, 1 Myr, and 3.5 Myr and show the moments at which most currently active planets are accreted by the host star or were subject to mergers.

Fig. 8 Mass and semimajor axis distribution of the planetary system after the gas disk has dispersed. The perihel and aphel of the planet as caused by its eccentricity is displayed via the error bars. The color shows the planet’s pebble mass fraction.

4.5 Final planetary system

Figure 8 shows the mass, semimajor axis, and eccentricity of the final system at 16.8 Myr. The color map indicates the pebble mass fraction of the individual planet. We find that the sub-Earth-mass planets in the terrestrial planet zone are dominantly composed of planetesimals and the inner super-Earths show a composition that stems from both pebbles and planetesimals. There is a clear dichotomy to be found between the close in super-Earths and the planets outside of 0.1 au. This dichotomy is related to orbital properties, such as the semimajor axis and individual eccentricities but also to planetary properties, such as their total mass, pebble content, and embryo formation time. The close-in planets show significantly higher eccentricities than the planetesimal composed planets in the terrestrial planet zone. Most planets outside 0.1 au share very low eccentricities.

Figure 9 shows the pebble mass fraction over the planet mass for the final system after 16.8 Myr. The color map indicates the formation time of the corresponding embryo. The highest mass planets also contain the highest pebble mass fraction and the earliest formation time. The planets that formed at a later stage of the disk evolution remain at small masses and their pebble mass fraction remains below 10%.

Fig. 9 Mass and pebble mass fraction of the planetary system after the gas disk has dispersed. The color shows the formation time of the planets.

5 Discussion

5.1 Formation of multiple generations of planets

As a result of the embryo formation model, the first embryos form only in the inner region of the disk, were they are subject to effective pebble accretion. Embryos then rapidly grow in mass and migrate to the inner edge of the disk. Since the formation of planetary embryos depends both on the planetesimal surface density and the heliocentric distance, embryos at larger distances (>10 au) do not form within the lifetime of the gas disk.

As the formation of embryos in our applied model does not occur within a given orbital separation to other embryos or planets, we find no embryo formation once the terrestrial planet region is populated by simultaneously growing planets. Further out, embryo formation cannot take place within that time, as planetesimals could not grow to a lunar mass as a consequence of larger growth time scales with orbital distances, as well as lower corresponding planetesimal surface densities. Once the super-Earth-mass planets have migrated to the inner edge of the gas disk and have been accreted by the star, the inner region of the disk is free from planets and embryo formation from the remaining planetesimals occurs. Since there is still a greater amount of gas by the time the first super Earth planets have migrated inwards, the eccentricity damping of the remaining planetesimals occurs, which had been excited by the super-Earth-mass planets that rapidly grew to >10 M_⊕. The embryo growth rate depends on the dispersion velocity (see Eq. (2)), which again is given as υ_∞ = e(r) Ω(r) in our framework. An increase in the eccentricity thus reduces the growth of planetary embryos.

Eccentricity damping by the gas, however, leads again to shorter embryo growth time scales, since eccentricity damping reduces the dispersion velocity.

In our model, we thus find multiple generations of planets. The first generation of planets that forms in the terrestrial planet zone grows rapidly by pebble accretion, followed by rapid inward migration and subsequent accretion onto the star. This process is followed by a second generation with a similar fate as the first generation. The second generation of embryos also forms within the lifetime of the pebble flux. This generation grows to super-Earth masses via pebble accretion as well and is then subject to migration. As the pebble flux vanishes over time, the next generation cannot grow as massive as the previous ones and its migration speed is thus largely reduced. We find a set of sub-Earth-mass planets growing by planetesimal accretion from 1 Myr to 3.5 Myr. Those planets, however, are also pushed toward the host star eventually, as a more massive planet from the previous generation migrates inward as well. After the last pebble-based super-Earth has migrated to the inner edge of the disk after 3.5 Myr, the terrestrial planet zone is free from planets once again and eccentricity damping enables the last generation to form. Until the end of the gas disk lifetime at about 16.8 Myr, planetary embryos can form up to a distance of 7 au. The last generation grows dominantly by planetsimal accretion and remains largely at sub-Earth masses. The last generation of planets does not experience strong type I migration because of their low masses; therefore, the clearing of planets in the terrestrial planet zone as occurs with previous generations, does not take place in this phase. We find that for as long as we have an active pebble flux, pebble accretion on embryos and fast type I migration clear the terrestrial zone from planets and the eccentricity damping of planetesimals enables the next generation of embryos to form.

The planetary composition, in terms of whether the mass stems from pebble or planetesimal accretion, reflects this picture. The early generations of planets are mostly composed of pebbles, whereas the lower mass, later generations of planets are predominantly formed by planetesimal accretion. The close-in super-Earth-mass planets in the final system are composed both of pebbles and planetesimals. They stem from the second generation of planetary formation, when the pebble flux had largely vanished. While the most massive remaining planet initially grew mostly by pebble accretion as well, a large fraction of its mass stems from planetesimals as it continues to accrete planetesimals after the pebble flux has vanished. We also find that the highest number of planets forms in the latest generation, since type I migration no longer forces planets to the inner edge of the disk and their small masses allow for a smaller orbital spacing.

5.2 Relevance of pebble accretion

The studied system raises several questions on the role of pebble accretion for planet formation. Pebble accretion itself is a very efficient mechanism for planetary growth. However, it is most efficient in the early phase of the disk evolution because it requires a pebble reservoir. During this phase, planetary migration is also highly efficient, due to the presence of a massive gas disk. The earliest planets that formed in our simulation thus experience both rapid growth and rapid migration. This combination is found to be detrimental to the probability of their survival in our study. We refrain from making a general statement here on the effectiveness of pebble accretion for planetary formation because our study only involves a single set of disk parameters. Nonetheless, we wish to highlight that efficient growth via pebble accretion might be a destructive mechanism in planetary formation for embryos forming early on in the inner disk.

In our study, most final planets are dominantly shaped via planetesimal accretion. These planets maintain their lower masses because the planetesimal reservoir in the terrestrial planet zone is rather small at the time of their formation. The reason for this small planetesimal reservoir lies in its depletion via planetesimal accretion onto the previously formed super-Earths. While pebble accretion does not play a direct role in the formation of the last generation of planets, it does so implicitly. The planetesimal reservoir in the terrestrial planet zone would contain much more mass if it was not for the previous generations of super-Earths that depleted it. The late, low-mass planetesimal-based generation thus required the earlier super-Earth, pebble-based generation to remain at such small masses and not to be carried away via migration as well. In this scenario, pebble accretion plays a crucial role in the global process of planet formation, even though most final planets are largely composed of planetesimals. Identifying the role of pebble accretion for a wider range of disk parameters will be subject to future work.

5.3 Embryo formation and migration

The model for embryo formation that is used in this study has been derived using N-body simulations including planetesimal formation (Voelkel et al. 2021a) and has been compared to simulations that also include the effect of pebble accretion (Voelkel et al. 2021b). Both of these studies did not show multiple generations of planetary embryo formation within 1 Myr. Even though several setups in Voelkel et al. (2021b) showed the rapid growth of super-Earth-mass planets in the terrestrial planet zone, these super-Earth-mass planets did not migrate inwards, since planetary migration due to planet disk interaction was not included. The formation of a next generation was therefore suppressed due to the presence of the super Earth mass planets in the terrestrial region. The inclusion of planetary migration in the more sophisticated N-body simulations, which contain both the formation of planetesimals and the accretion of pebbles, should form multiple generations of embryos within 1 Myr as well. A study along these lines to support and underline the findings of this paper will be conducted in future work.

5.4 Long-term evolution

We chose to end our simulation after the dispersal of the gas disk because our focus lies on the dynamic embryo formation of the first couple of Myr. However we still wish to briefly discuss the long-term evolution of the system. After the gas disk has vanished we find 16 active planets in the system and 68.59 M_⊕ in planetesimals. Three of those planets are very close in, with masses of 20.7 M_⊕ at 0.014 au, 4.5 M_⊕ at 0.029 au, and 0.93 M_⊕ at 0.022 au and with eccentric orbits. Those (super) Earth-mass planets are possibly to be accreted by the host star due to tidal interactions over timescales of Gyrs. The remaining system would then consist of the 13 planetesimal composed planets and the remaining 68.59 M_⊕ of planetesimals. The remaining planetesimal disk mass greatly exceeds that of the 13 planetesimal composed planets, a long-term integration of the system should, thus, also include the remaining planetesimals and their potential embryo formation to formulate a concise statement on the final system after several hundred Myr. Without the damping effect of the gas disk, however, the higher eccentricities of the planetesimals would reduce planetesimal accretion and embryo formation due to higher dispersion velocities.

5.5 Considering the architecture of the Solar System

We wish to discuss our simulations in consideration to the initial setup of the Grand Tack model (Walsh et al. 2011). In the Solar System, we find two gas giant planets at distances of 5–10 au, followed by two ice giants at 19–30 au. The inner region is populated with four smaller terrestrial planets. As our model suggests, the last generation resembles a large set of Earth-mass and sub-Earth-mass terrestrial embryos that are believed to form the four terrestrial planets in the Grand Tack model. The profound difference between the Grand Tack and our formation model is that the sub-Earth-mass terrestrial embryos did not form as a second or third generation in the Grand Tack scenario. Instead, they were merely the first generation of embryos and the reason why they did not grow to super Earth mass planets due to pebble accretion is due to Jupiter shielding the pebble flux. In contrary to the Grand Tack, our simulation suggests that a first generation of super Earth mass planets (eventually accreted by the host star) may have populated the terrestrial planet region during the first stages of the Solar System’s evolution.

Since the Solar System contains gas giants, which we do not include within our model framework, this hypothesis is subject to further investigation. The non-formation of Jupiter in the model presented here raises several profound challenges. The formation of planetary embryos is the result of planetesimals growth. This results in the very late formation of an embryo at larger distances (>5 Myr at 5 au). By that time, the flux of pebbles has vanished and the growth of a Jupiter core would extend the lifetime of the gas disk. The early formation of a core at a larger distance to form Jupiter would either require a single, initially much larger planetesimal to form or a local overdensity in planetesimals to reduce the formation time. Such an overdensity of planetesimals has been described as the result of a pressure bump in the gas disk (Guilera et al. 2020). This pressure bump, however, would also have major implications with regard to the evolution of the inner system. The question of whether or not we also find multiple generations of terrestrial planets when a pressure bump is included will be addressed in a future work.

Another possible way for giant planets to form would require the treatment of additional physics, namely, via further outward migration due to orbital resonances. If two planets are captured in mean motion resonance, they may form a gap in the gas disk (Walsh et al. 2011). In cases where the inner planet is the more massive one, this gap can cause the outward migration of both planets. Effectively, this might cause the formation of giant planets outside the initial orbit of their embryo formation. As gap opening in the disk is currently not included in our model, we cannot observe this process within our simulation.

6 Summary and outlook

In this paper, we investigate the effect of dynamic planetary embryo formation during the lifetime of the gaseous disk. To pave the way for our approach, we start with a self-consistent global model of planet formation that begins with an initial circumstellar disk of gas, dust, and pebbles. The model presents a viscous evolution of the gas disk and uses a two-population approach to model the evolution of dust and pebbles. Planetesimals form on the basis of the radial pebble flux and planetary embryos are introduced based on the evolution of the planetesimal surface density and their dynamical state. The eccentricities and inclinations of planetesimals are increased thanks to nearby planetary embryos and self-stirring. Concurrently, the eccentricity and inclination damping by the evolving gas disk is also considered here. Once planetary embryos have formed, they can grow via pebble, planetesimal, and gas accretion. Planets follow N-body dynamics with other planets and are subject to planet disk interactions, such as planetary migration. The number of embryos in the system, their initial location, and formation time, are no longer part of an initial assumption, but the results of the disk’s evolution. Our main findings can be summarized as followed:

• We find distinct generations of planet formation in the terrestrial planet region within the lifetime of the gas disk. Earlier generations grow dominantly by pebble accretion and are largely accreted by the host star due to migration. Later generations are composed largely of planetesimals, as those planets form after the pebble flux has mostly vanished;

• We find close-in super-Earth-mass planets composed of both pebbles and planetesimals and mostly planetesimally composed sub-Earth-mass planets in the terrestrial region. A first generation of embryos mostly formed of pebbles could not survive the gaseous disk, as they were accreted by the host star. The formation of close-in super-Earths and mini-Neptunes is a likely outcome for the early generations of planet formation;

• The majority of planetary embryos that form do not outlive the gas disk. Out of the 78 embryos that formed in total, only 16 remained after the disk vanished. The rest fell victim to either accretion to the star or mergers.

These findings mark the onset of a wide range of possibilities for the planetary formation and disk evolution model presented here. While the parameter space that we cover in this paper is focused on a single set of disk parameters, our model can be used within the framework of planet population synthesis as well. Additionally, it can be used to study individual features of single systems in a more detailed fashion as, for instance, the formation of planets in primordial rings due to pressure bumps. Next to our presented planet formation model, we will study the possibility of multiple generations of embryo formation using large scale N-body simulations in a future work. As multiple generations appear to already form within the first 1 Myr, a sophisticated N-body study, similar to that of Voelkel et al. (2021b), is computationally feasible and should confirm our findings. The underlying hypothesis of most planet formation models states that the final planets are the end-product from the initially placed bodies. This hypothesis is heavily challenged by our results. Dynamic planetary embryo formation shows the possibility for multiple distinct generations of planet formation. This finding is certain to have a fundamental effect on the formation history and composition of planets both in the solar system and exoplanet systems. It therefore needs to be accounted for in future studies.

Even though we claim to start with a nebula that was designed to create a solar system (Lenz et al. 2020), our simulations did not lead to a planetary system that resembles our own Solar System. Three effects can be responsible for this: a statistical effect of the N-body solver, our systematic initial condition, missing physics (or a combination of these effects). The statistical effect could be tested by performing numerous similar simulations and check whether this leads to a more Solar System like state for a number of outcomes. Then, Markov chain Monte Carlo (MCMC) simulations can be used to further constrain the potential initial conditions that formed the solar system, such as disk mass, size, profile, etc. However, the culprit is most likely the missing physics – even without improving the turbulence model, the viscous evolution of the disk, or the dust growth physics – which may be crucial for the non-formation of a complete solar system that would resemble our own. A major drawback in our framework are the tidal forces acting on the disk, which have not yet been implemented. Gap formation with pebble trapping or resonant outward migration of planet pairs as in the Grand Tack model, therefore cannot occur in this case; however, we note that this missing process could have feasibly led to a more Solar System-like outcome based on the chosen initial conditions, rather than stand as an impediment.

Acknowledgements

We wish to thank Remo Burn and Rogerio Deienno for many fruitful discussions. This research has been supported by the Deutsche Forschungsgemeinschaft via priority program (DFG SPP) SPP 1992 “Exploring the diversity of extrasolar planets” under contract: KL 1469/17-1 Consistent Planetesimal Formation from Pebbles for Synthetic Population Syntheses of Exo-Planets and SPP 1833 “Building a Habitable Earth” under contract: KL 1469/13-1 and KL 1469/13-2. This research was supported by the Deutsche Forschungs-gemeinschaft through the Major Research Instrumentation Programme and Research Unit FOR2544 “Blue Planets around Red Stars” for H.K. under contract DFG KL1469/15-1.

References

All Tables

All Figures

Fig. 1 Semimajor axis over time evolution for the planetary system. The evolution of a planet is linked via the grey line and the size of the dots indicate the mass of the planets every 10 ky. The color of the planets indicates their pebble mass fraction Mpeb/MP. The horizontal lines are drawn at 0.35 Myr, 1 Myr, and 3.5 Myr and show the moments at which most currently active planets are accreted by the host star or were subject to mergers. Thus we find four distinct generations of planet formation.

In the text

Fig. 2 Mass and semimajor axis distribution of the entire planetary system until the dispersal of the gas disk. The large circle markers indicate planets that remain active until the end of the simulation. Planets that were accreted by the host star are shown as triangles and planets that merged via collisions with other planets are indicated as squares. The track of the planets is shown as the solid grey line for active planets and dotted grey lines for accreted or merged planets. The color of the final marker indicates the initial formation time of the corresponding planetary embryo.

In the text

Fig. 3 Planet mass over time. The vertical lines are drawn at 0.35 Myr, 1 Myr, and 3.5 Myr and show the moments at which most currently active planets are accreted by the host star or were subject to mergers. The vertical steps are caused by giant impacts that result in planet mergers. The large circle markers indicate planets that remain active until the end of the simulation. Planets that were accreted by the host star are shown as triangles and planets that merged via collisions with other planets are indicated as squares. The track of the planets is shown as the solid grey line for active planets and dotted grey lines for accreted or merged planets. The color of the final marker indicates the pebble mass fraction of the corresponding planet.

In the text

	Fig. 4 Number of planets over time. The blue line shows the number of currently active planets in the system, while the orange line shows the total number of planets that have formed. The vertical lines are drawn at 0.35 Myr, 1 Myr, and 3.5 Myr and show the moments at which most currently active planets are accreted by the host star or were subject to mergers.
In the text

Fig. 5 Surface density and planetesimal eccentricity evolution of the circumstellar disk. We show the gas surface density (blue), planetesimal surface density (black), and the combined dust and pebble surface density (orange) at t = 0.35 Myr, t = 1 Myr, t = 3.5 Myr, and t= 16.8 Myr. The initial corresponding density is shown as the dashed line. The initial mass of the gas disk is given as M(Σg) = 34102.64 M⊕, the initial solid mass in dust and pebbles is given as Μ(Σtp) = 450.93 M⊕, and the initial planetesimal disk mass is given as M(Σtp) = 0 M⊕. The location of an active planet is indicated via dashed vertical lines. The planetesimal eccentricity is given as the dotted line, respectively.

In the text

Fig. 6 Fraction of the planetesimal disk mass (MD(ΣP)) and combined dust and pebble disk mass (MD(Σtp)) over the gas disk mass (MD(Σg)) within 1 au, 2.5 au, and 5 au over time in percent. The grey dotted line indicates the global initial dust to gas ratio dg = 1.34%. The dashed vertical lines are drawn at 0.35 Myr, 1 Myr, and 3.5 Myr and show the moments at which most currently active planets are accreted by the host star or were subject to mergers. The total mass in the planetesimal disk after the dispersal of the gas disk is given as 67.59 M⊕. The final planetesimal disk mass within 5 au is given as 10.5 M⊕, within 2.5 au as 0.349 M⊕, and within 1 au as 2.35 × 10−3 M⊕.

In the text

Fig. 7 Solid mass evolution during the lifetime of the gas disk. We show the total pebble and dust disk mass (MD(Σtp)), the total planetesimal disk mass (MD(ΣP)), the mass in active planetary cores (Mc(active)), the mass of all formed planetary cores (MC(formed)), and and the mass of planetary cores that have been accreted by the host star (MC(accreted)). The dashed vertical lines are drawn at 0.35 Myr, 1 Myr, and 3.5 Myr and show the moments at which most currently active planets are accreted by the host star or were subject to mergers.

In the text

	Fig. 8 Mass and semimajor axis distribution of the planetary system after the gas disk has dispersed. The perihel and aphel of the planet as caused by its eccentricity is displayed via the error bars. The color shows the planet’s pebble mass fraction.
In the text

	Fig. 9 Mass and pebble mass fraction of the planetary system after the gas disk has dispersed. The color shows the formation time of the planets.
In the text