© O. Voelkel et al. 2020

1 Introduction

A current conundrum of planetesimal accretion in the core-accretion scenario of planet formation is that for 100 km planetesimals it appears to require an unreasonably high disk mass to be an effective mechanism for giant planet formation within the lifetime of a circumstellar disk (Fortier et al. 2013). The accretion of smaller objects with a higher effective cross section, such as either kilometer (km)-sized planetesimals (Ida & Lin 2004) or centimeter (cm)-sized bodies known as pebble accretion (Ormel & Klahr 2010) is often described as the solution for giant planet formation and has been studied widely by Klahr & Bodenheimer (2006), Lambrechts & Johansen (2012), Levison et al. (2015), and Bitsch et al. (2015), to name just a few. While we refrain from making a statement on the efficiency of pebble accretion, the scenario of a planetary core that accretes inward-drifting pebbles also lacks an explanation on how, where, and when this planetary core first forms. Planetesimals are typically too small for efficient pebble accretion (Ormel & Klahr 2010), therefore a pebble-accreting embryo might well have formed from planetesimal collisions. This crucial step adds room for discussing the formation of planetesimals and subsequently their role in planetary core and planet formation. From Tanaka & Ida (1999) we know that the accretion rate of planetesimals depends on the planetesimal size and linearly on the planetesimal surface density. Constraining the size of planetesimals is an active field of research. While some studies infer that the current size of asteroid belt objects is well constrained and that they are about 100 km in diameter (Bottke Jr et al. 2005; Walsh et al. 2017; Delbo’ et al. 2017), other studies reproted that the size distribution found today merely reflects the sizes that are most resilient to clearing and therefore suggest a smaller primordial size (Zheng et al. 2017). The observed size distribution might also arise from the growth of planetesimals that originally measured 100 m (Weidenschilling 2011). In the Kuiper belt, the size distribution has a similar shape as predicted by simulations that include the streaming instability between 10 and 100 km (Schäfer et al. 2017), which indicates large initial sizes. On the other hand, recent discoveries of Kuiper belt objects through stellar occultations rather indicate a size of 1–2 km (Arimatsu et al. 2019). Small initial sizes of 0.4–4 km have also been inferred theoretically by Schlichting et al. (2013). The surface density profile of planetesimals for extrasolar systems is likewise unknown. Studies of our own Solar System motivated the minimum-mass solar-nebula (mmsn) hypothesis (Weidenschilling 1977; Hayashi 1981), which results in a power-law drop of the planetesimal surface density with a decay of Σ_P ∝ r^−1.5. Observations of solid material in disks (Andrews et al. 2010) and the widely used α-disk model for the viscous evolution of an accretion disk (Shakura & Sunyaev 1973) suggest a shallower density distribution of Σ_P ∝ r^−0.9 for radially constant α. The observed solid material is not planetesimals, however, but the dust in the circumstellar disk, as the distribution of planetesimals in protoplanetary disks is currently unobservable. Lenz et al. (2019) modeled the formation of planetesimals based on the solid evolution of a viscously evolving disk, assuming that planetesimals form proportional to the time-dependent local radial pebble flux. They reported that the profile of the planetesimal surface density becomes significantly steeper (Σ_P ∝ r^−2.1) than the initial dust, pebble, and gas density (Σ ∝ r^−0.9). This mass transfer results in an increase in the planetesimal surface density in the inner circumstellar disk by several orders of magnitude without increasing the total mass in planetesimals. Because the accretion rate of planetesimals is proportional to the local planetesimal surface density, these highly condensed planetesimal zones are promising candidates to exert a drastic effect on planetary growth.

Before we discuss some of the previous work, we distinguish between a global planet formation model and a model for planet population synthesis. While a model for planet population synthesis contains (or should contain) a global formation model, this does not apply vice versa. Key to the population synthesis approach is that the model is complex enough to take the physical effects into account that are deemed crucial for planet formation, but its single-system computational cost is low enough for it to be used to study a wide range of parameters. Only this will enable a statistical comparison with observational data. For this purpose, it is vital to find ways to simplify complex physical processes and merge them to a more complex framework without loosing the essence of their nature. The formation of planetesimals is such a process, and the one-dimensional formation model by Lenz et al. (2019) is such an attempt. Previous work on the accretion of planetesimals for planetary growth such as Johansen & Bitsch (2019), Mordasini (2018), or Ida & Lin (2004) used initial distributions of planetesimals and initially placed planetary embryos, while neglecting the presence of pebbles. Other formation models such as those of Bitsch et al. (2015) or Brügger et al. (2018) modeled planetary growth by the accretion of pebbles and initially set planetary embryos, while neglecting the formation, or accretion, of planetesimals. A model that contains both pebble and planetesimal accretion while also taking the formation of planetesimals and planetary embryos into account is still pending.

We therefore chose to improve our planet population synthesis model by a “disk-consistent”¹ model for solid evolution (Birnstiel et al. 2012) and planetesimal formation (Lenz et al. 2019) to take the early stages of the disk evolution into account. This early phase determines the planetesimal surface density distribution, the radial pebble-flux evolution, and the formation of planetary cores and therefore planet formation as a whole. For our study, we focus on the formation and accretion of planetesimals. We display the effect of the planetesimal surface density and its formation on the population of planets. We show that the accretion by 100 km sized planetesimals is in fact a highly efficient growth mechanism for planets because of the highly condensed planetesimal regions in the disk. Furthermore, we give an overview of future possibilities that arise from our newly implemented modules.

This paper is structured as followed: in Sect. 2 we explain the planetesimal formation model together with the newly implemented solid-evolution model on which it is based. Section 3 provides insight into the population synthesis framework and how it was modified for our purpose. The changes in Σ_P in the population synthesis code and the newly computed synthetic populations, are presented in Sect. 4. Section 5 discusses the results, followed by a brief summary and an outlook on our new possibilities and future work in Sect. 6.

2 Planetesimal formation model

2.1 Two-population solid-evolution model

The two-population model for solid evolution by Birnstiel et al. (2012) is a parameterized approach to model the evolution and growth of dust and cm sized bodies in circumstellar disks. A detailed description of the model can be found in Birnstiel et al. (2012) and Lenz et al. (2019). We briefly outline the assumptions here and list the most important reasons for choosing it in our framework. Our goal is to implement a fast-computing, one-dimensional, parameterized algorithm for solid evolution that is well tested and agrees well with more sophisticated models. Key for the performance of the two-population approximation is a parameterized mass ratio f_m (r) as a function of orbital distance r between twopopulations of solids, which depends on whether the growth of the particles is limited by drift or by fragmentation. In each time step, the model solves one advection-diffusion equation that is given by $$\frac{\partial {\Sigma }_{\text{s}}}{\partial t}+\frac{1}{r}\frac{\partial }{\partial r}\left[r\left({\Sigma }_{\text{s}}\overline{u}-D_{\text{g}}{\Sigma }_{\text{g}}\frac{\partial }{\partial r}\left(\frac{{\Sigma }_{\text{s}}}{{\Sigma }_{\text{g}}}\right)\right)\right]=0,$$(1)

where Σ_s is the total solid-surface density without planetesimals, Σ_g is the gas surface density, D_g is the gas diffusion coefficient, and t and r are time and radial distance. ū describes the weighted velocity of the total solid density and is defined as $$\overline{u}=(1-f_m(r))\cdot u_0+f_m(r)\cdot u_1,$$(2)

where f_m is the fit parameter for the mass ratio between the two populations. u₀ and u₁ describe their velocities, while the surface densities of the two populations are given as $$\begin{array}{ll}{\Sigma }_0(r)\hfill & ={\Sigma }_{\text{s}}(r)\cdot (1-f_m(r))\hfill \\ {\Sigma }_1(r)\hfill & ={\Sigma }_{\text{s}}(r)\cdot f_m(r).\hfill \end{array}$$

The two populations are defined by their Stokes number. Particles with a small Stokes number of St ≪1 are strictly coupled to the evolution of the gas, whereas particles with St ≥ 1 are not. Σ₀ describes the smaller population, which can be seen as dust, subject to diffusion and transport with the gas, while Σ₁ describes the larger population, which can be seen as pebbles, which in addition to being diffused by the gas are also sedimenting toward the midplane and drifting toward pressure maxima, for instance, toward the star. The fit parameter f_m was derived by comparing the two-population model to the more sophisticated dust model from Birnstiel et al. (2010). The values for f_m that were the best fit are given as $$f_m=\{\begin{array}{ll}0.97,\hfill & \text{drift-limited case}\hfill \\ 0.75,\hfill & \text{fragmentation-limited case}\text{.}\hfill \end{array}$$(5)

These are also the values that we used in our simulations. The effect of this implementation is shown in Fig. 1, where the ratio between dust and pebbles varies with space and time. This is shown by the two blue curves.

Fig. 1 Exemplary disk evolution including our dynamical model for planetesimal formation after 0.1, 1, and 2 Myr. We showthe surface density for the dust, pebbles, planetesimals, gas, and their individual disk masses. The dashed lines refer to the initial profile of the corresponding density. This run does not contain a planetary embryo, it only evolves the disk dynamically. The total disk gas mass is given as 0.012 M⊙ with a dust-to-gas ration of 1.5% and α = 10−3. The exponential cutoff radius of the disk is at 137 au, the inner radius at 0.03 au, and the evaporation rate is given as 2.87 × 10−5 M⊙ yr−1. The planetesimal and solid-evolution parameters can be found in Table 1. The effect of the ice line is visible in the kink in the planetesimal distribution around 1 au and the effect of drift vs. fragmentation-limited pebble size in the radially varying dust-to-pebble ratio.

2.2 Pebble flux-regulated planetesimal formation

The full model and its results are described in Lenz et al. (2019) in greater detail. We therefore outline here only the basic physical assumptions of this one-dimensional approach and summarize the most important equations and results. The principle behind this parameterized model is that planetesimals form by a local continuous mechanism that converts a certain fraction of the pebbles that drift by into planetesimals. In principle, it thus acknowledges that pebbles will drift inward and that more planetesimals can be formed when more material comes by. Many different planetesimal formation prescriptions can therefore be parameterized in this way. Wheter in the framework of turbulent clustering (Cuzzi et al. 2010; Hartlep & Cuzzi 2020), streaming instabilities (Johansen et al. 2009; Schäfer et al. 2017), local trapping in zonal flows (Johansen et al. 2007, 2011; Dittrich et al. 2013; Drążkowska & Alibert 2017), or in vortices (Raettig et al. 2015; Lyra et al. 2018), the formation is always limited by the number of fresh pebbles that a region receives after consuming the locally available pebbles. Our parameterization is thus by definition model independent. Different scenarios might lead to the same conversion rates for the pebble flux. The parameters we need is the fraction ɛ of pebbles that is converted into planetesimals after having drifted over a distance of d within the disk. We can motivate these parameters easily in our paradigm of trapping zones that slowly evolve coherent flow structures in protoplanetary disks, such as vortices and zonal flows (Klahr et al. 2018), which can form everywhere, live only for a limited time, and thus only trap a fraction of drifting pebbles. In these traps, pebbles become sufficiently concentrated for planetesimal formation to be triggered, regulated by streaming and Kelvin Helmholtz instabilities. The planetesimal formation rate is generally proportional to the radial pebble flux, $${\dot{M}}_{\text{peb}}:=2\pi r{\sum }_{{\text{St}}_{\text{min}}\le \text{St}\le {\text{St}}_{\text{max}}}|v_{\text{drift}}(r,\text{St})|{\Sigma }_{\text{s}}(r,\text{St}),$$(6)

where v_drift is the drifting velocity of the particles and St_min and St_max are the minimum and maximum Stokes number for which a particle is considered a pebble. v_drift is given as $$v_{\text{drift}}(r,\text{St})=\frac{\text{St}}{{\text{St}}^2+1}\frac{h_{\text{g}}(r)}{r}\frac{\partial \text{ln}P(r)}{\partial \text{ln}r}{c}_s(r),$$(7)

where P(r) is the gas pressure, h_g(r) is the gas pressure scale height (h_g(r) = c_s(r)∕Ω(r)), and c_s(r) is the sound speed. Ω(r) is given as the orbital frequency at the radialdistance r. The source term for planetesimals, that is, for Σ_P, is then given as (Lenz et al. 2019) $${\dot{\Sigma }}_{\text{p}}(r)=f_{\text{ice}}(T)\frac{\varepsilon }{d(r)}\frac{{\dot{M}}_{\text{peb}}}{2\pi r,}$$(8)

where d(r) is the radial separation of the pebble traps and ɛ is the efficiency parameter that describes how much of the pebble flux is transformed into planetesimals after drifting over a distance of d. We chose a constant value of ɛ = 0.05 as a good value to form a sufficient number of planetesimals, as was found in Lenz et al. (2019) for d(r) = 5.0 pressure scale heights, motivated by our findings in the detailed numerical simulations of zonal flows (Dittrich et al. 2013). Generally, we can change ɛ locally when the formation of planetesimals might follow a different underlying mechanism, for instance, around the water-ice line, as described by Drążkowska & Alibert (2017) or Schoonenberg & Ormel (2017). This flexibility allows us to study a broad range of planetesimal formation scenarios, using the same implementation. Our two-population implementation currently has no proper treatment of the processes of evaporation and possible recondensation. The only effect of the existing ice line is incorporated into the parameter f_ice(T): $$f_{\text{ice}}(T)=\{\begin{array}{ll}1\hfill & \text{for }T<=170\text{K}\hfill \\ \frac{1}{3}\hfill & \text{for }T>170\text{K}\hfill \end{array}\text{}\text{}\text{}\text{},$$(9)

in effect to reduce the pebble flux inside the ice line to compensate for the evaporation of water ice. Therefore the ice line is visible in the distribution of planetesimals, even though it is not visible in the pebbles themselves (see Fig. 1). We also used a fixed planetesimal size of 100 km in diameter as in Lenz et al. (2019). As a consequence, a threshold of transformed mass must be reached to build at least one planetesimal. From this we can derive a critical pebble flux that is required for Σ_P to change. It is given as (Lenz et al. 2019) $${\dot{M}}_{\text{cr}}:=\frac{m_{\text{p}}}{\varepsilon {\tau }_{\text{t}},}$$(10)

where τ_t describes the average lifetime of a trap, which is given as 100 local orbits, and m_P is the mass of a single planetesimal. For simplicity, we assumed spherical planetesimals with a uniform density of ρ_s = 1.0 g cm⁻³. The mass that is transformed into planetesimals arises as a sink term in the advection-diffusion equation (Eq. (1)). The new advection-diffusion equation is then given as $$\frac{\partial {\Sigma }_{\text{s}}}{\partial t}+\frac{1}{r}\frac{\partial }{\partial r}\left[r\left({\Sigma }_{\text{s}}\overline{u}-D_{\text{g}}{\Sigma }_{\text{g}}\frac{\partial }{\partial r}\left(\frac{{\Sigma }_{\text{s}}}{{\Sigma }_{\text{g}}}\right)\right)\right]=L,$$(11)

where the sink term L is defined as $$L=(1-f_m(r))\cdot L_0+f_m(r)\cdot L_1,$$(12)

with $$\begin{array}{ll}{L}_0=\hfill & \frac{\varepsilon }{d(r)}\cdot v_{\text{drift 0}}{\Sigma }_0\cdot \theta ({\dot{M}}_{\text{peb}}-{\dot{M}}_{\text{cr}})\hfill \\ \hfill & \times \theta ({\text{St}}_{\text{0}}-{\text{St}}_{\text{min}})\cdot \theta ({\text{St}}_{\text{max}}-{\text{St}}_{\text{0}})\hfill \end{array}$$(13)

and $$\begin{array}{ll}{L}_1=\hfill & \frac{\varepsilon }{d(r)}\cdot v_{\text{drift 1}}{\Sigma }_1\cdot \theta ({\dot{M}}_{\text{peb}}-{\dot{M}}_{\text{cr}})\hfill \\ \hfill & \times \theta ({\text{St}}_{\text{1}}-{\text{St}}_{\text{min}})\cdot \theta ({\text{St}}_{\text{max}}-{\text{St}}_{\text{1}}),\hfill \end{array}$$(14)

where θ (⋅) is the Heaviside function. This combines the conditions described above for planetesimal formation. The surface density can only change while a critical mass is transformed (θ (Ṁ_peb −Ṁ_cr)) and when the Stokes numbers of the particles are within St_min and St_max (θ (St₀ −St_min) ⋅ θ (St_mas −St₀)).

3 Planet formation and evolution model

3.1 Model components

The current version of our planet population synthesis model can be found in Emsenhuber et al. (2020). It is an update of the model presented in Mordasini (2018). This model combines planet formation (Alibert et al. 2005, 2013) and evolution (Mordasini et al. 2012). Descriptions of the model can be found in Benz et al. (2014), Mordasini et al. (2015), Mordasini (2018), and in upcoming work by Emsenhuber et al. (2020). We provide here an overview of the physical processes that are tracked in the model and focus on the solid components of the protoplanetary disk model in Sect. 3.2.

The formation part of the model follows the core-accretion scenario of planetary embryos in viscously evolving circumstellar disks (Lüst 1952; Lynden-Bell & Pringle 1974). The macroscopic viscosity is given by the α parameterization (Shakura & Sunyaev 1973). Planetesimals are assumed to be in the oligarchic regime (Ida & Makino 1993; Thommes et al. 2003; Chambers 2006; Fortier et al. 2013). The structure of the envelope is retrieved by solving the internal structure equations (Bodenheimer & Pollack 1986). During the initial phase, gas accretion is governed by the ability to radiate the potential energy gained by the accretion of both solids and gas (Pollack et al. 1996; Lee & Chiang 2015). The efficiency of cooling increases with the planetary mass, and when the gas-accretion rate is limited by the supply of the gas disk, the planet contracts (Bodenheimer et al. 2000).

Planets embedded in a gas disk will undergo migration (e.g., Baruteau et al. 2014). The model uses the prescription of Dittkrist et al. (2014). For type I migration, it is based on the work by Paardekooper et al. (2010), while for type II migration, planets move in equilibrium with the gas disk. The switch between the two follows the criterion of Crida et al. (2006).

The formation stage lasts for the entire lifetime of the protoplanetary disk, but for 10 Myr at least. When this is completed, the model switches to the evolution stage (Mordasini et al. 2012), in which the planets are followed until 10 Gyr. This stage follows the thermodynamical evolution of the planets, with atmospheric escape (Jin et al. 2014) and tidal migration.

To perform population synthesis, we used a method similar to that of Mordasini et al. (2009), with severaladaptations. The distribution of disk gas masses and the relationship between the mass and the exponential cutoff radius follow Andrews et al. (2010). The inner radius was fixed to 0.03 au. The initial embryo mass was 0.0123 M_⊕ and the location was random, with a uniform distribution in the logarithm of the distance between 0.06 and 40 au. Embryos were placed directly at the beginning of the simulations.

Fig. 2 Schematic display of the different stages in the formation model for solid-evolution development. Upper panel: previously published model from Mordasini (2018). Middle panel: currently improved version in this work, including the two-population solid evolution for dust and pebbles, as well as the formation of planetesimals (see Sect. 2). Lower panel: outlook on possible future development stages. The new modules and functions are highlighted in red, and future possibilities are highlighted in blue.

Fig. 3 Planetesimal disk masses within 1 au, 10 au, and the complete disk for three different analytic density slopes and the dynamically formed planetesimal mass. The analytic masses are given at the start of the simulation, while the dynamically formed disk masses are shown after one million years, after most planetesimals have already formed. The dynamic runs do not contain a planetary embryo, they only simulate the disk evolution. The disk parameters, however, are the same as in the population in Fig. 5. We show the mass in planetesimals in the whole disk in the upper panel, the planetesimal mass within 10 au in the middle panel, the planetesimal mass within 1 au in the lower panel, and the corresponding median masses for every setup.

3.2 Solid component

A schematic overview of the different modules is shown in Fig. 2. Previous generations of the model, including the upcoming version of Emsenhuber et al. (2020), used an initial planetesimal surface density slope that was set either to be equal to the initial gas-density slope (Mordasini et al. 2009) or used a Σ_P ∝ r^−1.5 mmsn-like distribution (Emsenhuber et al. 2020). For the first case, this gave a planetesimal surface density distribution of Σ_P ∝ r^−0.9 up to an exponential cutoff radius, which depends on the given disk size. The total mass in planetesimals was chosen to be the metalicity (in the following, dust-to-gas ratio d_g) of the host star times the total gas disk mass, modulo the effect of condensation fronts. The size of the planetesimals was chosen to be uniform and with a radius of r_P = 300 m. Importantly, Σ_P only evolved while being accreted or ejected by embryos. Planetesimal formation or drift were not included, which left us with a static distribution of planetesimals and a complete lack of a physical description of the early phases of planet formation.

With our newly implemented model for planetesimal formation we proceed beyond the standard implementation in Emsenhuber et al. (2020). We now include two additional solid quantities (dust and pebbles) that evolve along with the gas evolution of the disk model. The initial mass in dust and pebbles21 is given as the metallicity of the host star times the gas-disk mass. Their density slope is set to be equal to that of the gas disk, giving an initial solid density profile of Σ_s ∝ r^−0.9. There are no initially placed planetesimals. Planetesimals only form based on the evolution of dust and pebbles. This ensures that planetesimals form consistently with the disk evolution. Not only is the final distribution of planetesimals highly different than the static assumption of the previous disk model (see Sect. 4.1), but planetesimals now also form over time, which opens a completely new level of dynamical interaction with the disk. The size of planetesimals that we assumed in the following simulations is given as 100 km in diameter.

The main differences between Emsenhuber et al. (2020) and this paper therefore are the size of planetesimals (r_P = 300 m vs. r_P = 50 km) and the option for dynamic planetesimal formation, which is not yet implemented in Emsenhuber et al. (2020). Emsenhuber et al. (2020), on the other hand, include an N-body integration for multiple simultaneously evolving cores. We did not use this option here because we wished to focus on the effect of dynamical planetesimal formation.

Fig. 4 Cumulative distribution of gas-disk lifetimes for our synthetic population. We used a constant value of α = 10−3 in our runs. We find that 50% of the lifetimes are shorter than 4.1 Myr.

4 Results

4.1 Disk evolution

Previous simulations with our model used an initial Σ_P of d_g ⋅ Σ_g, where d_g is the dust-to-gas ratio. The slope in Σ_P was therefore given as the slope of the initial gas surface density.

The density slope that arises from the pebble flux-regulated model for planetesimal formation can have a slope as steep as Σ_P ∝ r^−2.1, and it generally depends on the individual evolution of the disk. Because of the steeper slope, we find a remarkable increase in Σ_P in the inner regions of a protoplanetary disk and a corresponding decrease farther out. Another profound difference to the previous implementation of our model is the total mass in planetesimals. The initial mass in dust in the planetesimal formation runs is equal to the initial mass in planetesimals with the analytically given planetesimal surface density, but only a fraction of this is transformed into planetesimals. We therefore always undershoot the total mass in planetesimals for our dynamically formed simulations compared to the previous implementation. Choosing higher values for the planetesimal formation efficiency can result in a shallower density profile, similar to that of the initial gas distribution. The initial dust density is given as a fraction of the gas surface density. If we were to consider ɛ∕d > 1, this would lead to local pebble-to-planetesimal conversion and the outer material would be unable to drift into the inner regions of the disk, which would have changed the density profile. For a more detailed treatment of this behavior, we refer to Lenz et al. (2019). To find similar densities to the analytic Σ_P ∝ r^−2.1 runs, we would have to increase our disk masses to match the final mass in planetesimals.

The total disk masses for the different density distributions are shown in Fig. 3, as well as the masses within 10 and 1 au. We find that the mean total disk masses are lower for the steeper density profiles by a factor of $M_{\text{tot}}^{-2.1}/M_{\text{tot}}^{-0.9}\approx 0.62$ or $M_{\text{tot}}^{-1.5}/M_{\text{tot}}^{-0.9}\approx 0.87$. This is to be expected as more material is inside the ice line, which is taken care of in these models. The masses within 1 au of the steeper models are still higher by several orders of magnitude ($M_{\text{1au}}^{-2.1}/M_{\text{1au}}^{-0.9}\approx 21,58$ or $M_{\text{1au}}^{-1.5}/M_{\text{1au}}^{-0.9}\approx 6,01$). M_tot and M_1au refer to the median masses from Fig. 3. The lowest total median mass ratios of planetesimals can be found in the dynamically formed simulation with $M_{\text{tot}}^{\text{dyn}}/M_{\text{tot}}^{-0.9}=0.504$, the mass ratio within 1 au, however, is the second highest with $M_{\text{1au}}^{\text{dyn}}/M_{\text{1au}}^{-0.9}=8.27$. The lower total masses for the steeper planetesimal surface density can be explained by the smaller number of icy planetesimals in these setups. Choosing a steeper density slope for the same mass as in the Σ_P ∝ r^−0.9 shifts material (icy planetesimals and silicate planetesimals) from farther out regions to the inner disk. This would evaporate the icy planetesimals within the ice line, leaving only the silicate planetesimals, therefore effectively loosing mass. The mass loss here is therefore only due to icy planetesimals within the ice line, whereas the amount of silicate planetesimals remains the same. Regardless of this mass loss, we find that the mass in the inner disk (r < 1 au) is significantly higher for the steeper density slopes. The lifetimes of the gas disks studied in our case are shown in Fig. 4 The global effect on planet formation of these changes in Σ_P is presented in Sect. 4.2.

Fig. 5 Mass vs. semimajor axis of synthetic planet populations for different ΣP distributions after 100 million years. Each setup contains one single planetary embryo. The initially set distributions for ΣP are ΣP = Σ0 ⋅ r−0.9 (initial gas density slope), ΣP = Σ0 ⋅ r−1.5 (mmsn) and ΣP = Σ0 ⋅ r−2.1 (Lenz et al. 2019). Bottom right panel: population in which planetesimals form over time using the model described in Sect. 2. The circles given around the data points show the mass fraction of envelope mass over core mass. The numbers of systems are 1999 (ΣP ∝ r−0.9), 1990 (ΣP ∝ r−1.5), 1961 (ΣP ∝ r−2.1), and 1945 (ΣP- dyn).

4.2 Synthetic populations

In the following we present several synthetic populations that were computed with different initial planetesimal surface density profiles and the dynamic planetesimal formation model from Lenz et al. (2019). It is important to mention that the growth of planetary embryos by the accretion of solids is only given by the accretion of planetesimals in these simulations. To ensure the correct comparison of planetesimal accretion with different slopes of Σ_P, we neglect the accretion of pebbles for this part of our study. We also consider systems with one embryo each because our focus lies on the changes to the previous implementation. Although populations with a much higher number of embryos are possible in the new version of the model (Emsenhuber et al. 2020), we chose to stay with one embryo per run for our study because mixing our study with effects of multiple planets in that forthcoming paper avails us nothing. We therefore focus on the general distribution of masses and semimajor axes and on the overall mass occurrences of planets.

Fig. 6 Planetary and core-mass occurrences of the four different populations. Every planet in each systems starts with a core mass of 0.0123 M⊕ and no envelope. The quantities that arise from the three analytical planetesimal surface density profiles are shown in blue (ΣP ∝ r−2.1), orange (ΣP ∝ r−1.5), and green (ΣP ∝ r−0.9), whereas the properties of the planetesimal formation population are shown in black. The dashed lines in the plots showthe median planet and median core masses. The histograms show clear shifts toward the higher mass ranges for steeper planetesimal surface densities and for the dynamically formed planetesimals than the ΣP ∝ r−0.9 or even the ΣP ∝ r−1.5 distribution.

4.2.1 Mass semimajor axis distributions

Figure 5 shows the mass and semimajor axis distribution of four synthetic populations around a solar-type star using one planetary embryo of lunar mass (0.0123 M_⊕) for each system. We simulated a total number of 1999 systems for the Σ_P ∝ r^−0.9 distribution, 1990 for Σ_P ∝ r^−1.5, 1961 for Σ_P ∝ r^−2.1, and 1945 forthe dynamic planetesimal formation run. The initial conditions of the four populations are the same, except for the initial Σ_P and the formation of planetesimals.

The upper left green panel refers to an initial Σ_P of Σ_P ∝ r^−0.9, the upper right orange to Σ_P ∝ r^−1.5 and the lower left blue to Σ_P ∝ r^−2.1. The lower right panel in black refers to the final planets that formed using the pebble flux-regulated model for planetesimal formation. We find a large number of planets that exceed a mass of ten Earth masses (necessary for runaway gas accretion, see Pollack et al. 1996) and sometimes even reach several hundreds of Earth masses when the slope of Σ_P is given with a slope of r^−2.1. The simulation in which the slope is given with the r^−0.9 does not even produce one single planet with a mass higher than that of ten Earth masses. Overall, this plot shows an immense increase in planetary masses for steeper planetesimal density profiles. It is important to mention here that the high-mass gas giant planets all end up within 1 au, which is due to the high masses in planetesimals in the inner disk and planetary migration. In our synthetic runs we do not see gas giants farther out, for example, beyond the water-ice line, as can be observed in the population of exoplanets (Winn & Fabrycky 2015). This will probably change when recondensation of water vapor is included, which effectively boosts the birth of planetesimals that is regulated by pebble flux.

Fig. 7 Planetary mass and core-mass occurrences of the four different populations within 1 au. Figure 5 shows that the higher mass planets are located in the inner disk because these are the regions with the highest planetesimal surface density. The dashed lines in the plots show the median planet and median core masses within 1 au.

4.2.2 Mass occurrences

For a more quantitative analysis, we studied the mass occurrences for the different planetesimal density slopes. Here we focus on the planetary mass and the core mass. Figures 6 and 7 show histograms with the occurrences of the different masses for the various populations from Fig. 5. As Fig. 5 shows, most of the high masses are found in the inner parts of the protoplanetary disk, whereas embryos placed farther out fail to grow. We therefore also focused our study on the inner region within 1 au. Figure 6 takes the complete population into account, whereas Fig. 7 only contains planets with a semimajor axis smaller than1 au. We also give the median masses for the planets and their cores. A cumulative function of the planetary masses is shown in Fig. 8. We find that the number of planets above 10 M_⊕ is given as 0 (Σ_P ∝ r^−0.9), 159 (Σ_P ∝ r^−1.5), 565 (Σ_P ∝ r^−2.1), and 301 (Σ_P - dyn). The number of planets above 20 M_⊕ is given as 31 (Σ_P ∝ r^−1.5), 383 (Σ_P ∝ r^−2.1), and 138 (Σ_P- dyn).

4.3 Effect of the starting location

In Fig. 9 we show the semimajor axis distribution and the initial starting location distribution of high-mass planets in the Σ_P ∝ r^−2.1 run and the dynamic formation model. We find that planets with the highest masses end up at the inner edge of the disk due to migration. There is no in situ giant formation, but rather a preferential zone in which planetary embryos need to be placed in order to grow to giant planets. This preferential area appears to be around 1 au and from around 4 to 10 au for the Σ_P ∝ r^−2.1 run and mostly from around 4 to 10 au for the dynamical formation model. Embryos that are placed at a distance from 2 to 4 au appear to have a lower probability of becoming gas giant planets in both cases, but because the probability of their formation at this locationis also low because of the local deficiency in planetesimals, they should not have been placed there in the first place. Now that we have a distribution of planetesimals, we can use this information to also model the generation of embryos in a consistent fashion. Ultimately, the effect of recondensation beyond the ice line described above can further change this picture.

4.4 Gas giant growth

Here we focus on a system that forms a 997.6 M_⊕ mass planet for the Σ_P ∝ r^−2.1 density distribution and a 281.7 M_⊕ planet for the dynamical planetesimal formation run. The initial disk parameters for the setup are given in Table 1. Figure 10 shows planetary growth tracks, the mass growth over time, and the corresponding semimajor axis evolution. The embryo in these systems was placed initially at 8.2 au, which seems to be a preferential starting location for giant planets, see Fig. 9. The higher planetesimal surface density has a drastic effect on the early stages of planetary growth. The planets in the Σ_P ∝ r^−2.1 setup and the dynamical planetesimal formation run can grow fast enough to undergo runaway gas accretion, whereas the planet in system Σ_P ∝ r^−1.5 fails to do so, even though its core reaches a core mass of 43 M_⊕. The planet in system Σ_P ∝ r^−0.9 fails to build a large enough core for significant gas accretion and ends up at 4.64 M_⊕. The 43 M_⊕ planet does not end up as a gas giant because of planetary migration and the continuous accretion of planetesimals during the migration phase. Continuous planetesimal accretion reduces the gas-accretion rate due to accretional heating, and while the accretion rate of solids surpasses the accretion rate of gas, the planet cannot proceed into runaway gas accretion (Pollack et al. 1996). We find that it is not the total mass available or the accretion timescales for 100 km planetesimals that prevent giant planet formation, but fast type I migration. The upper right panel of Fig. 5 shows that most planets above 10 M_⊕ end up at the inner edge of the gas disk. Their evolution is similar to that of the planet in the Σ_P ∝ r^−1.5 system in Fig. 10 in the sense that they grow well above 10 M_⊕, but then migrate too quickly to the inner edge of the disk. Figure 9 also shows that most massive planets in the Σ_P ∝ r^−2.1 run end up at the inner edge of the disk due to migration as well. The change in migration type allows the planets in the Σ_P ∝ r^−2.1 and the dynamical run in Fig. 10 to undergo runaway gas accretion. They switch from type I to type II migration before they end up too close to the star. The planet in the Σ_P ∝ r^−1.5 run also transitions from type I to type II migration at 0.26 au, but this is already too close to the star, and type II migration continues to transport the planet to the disk edge before entering runaway gas accretion. The planet in the Σ_P ∝ r^−2.1 run switches to type II migration at 0.65 au, giving it enough time to enter runaway gas accretion before migrating too closely to the star. When the migration rate from type I to type II is slowed down, a decrease in the accretion rate of planetesimals is caused that reduces the accretional heating and thus enables runaway gas accretion if the planet has not already migrated too closely to the star as in the Σ_P ∝ r^−1.5 case of Fig. 10. The limitation for the formation of gas giant planets by 100 km planetesimal accretion is therefore mostly given by planetary migration and not by growing 10 M_⊕ planets within the lifetime of a gaseous disk.

Fig. 8 Cumulative function of planetary masses within 1 au for the four different synthetic populations from Fig. 5. The y-axis shows how many planets for each density profile are above the current planetary mass, normalized by the total number for planets in each population. The planets that have formed in the ΣP ∝ r−0.9 run are shown in green, the ΣP ∝ r−1.5 population is shown in orange, and the ΣP ∝ r−2.1 population is shown in blue. The dynamic planetesimal formation population is shown in black.

Fig. 9 Semimajor axis distribution and starting location of planets that have grown to gas giant masses (MP >100 M⊕) in the ΣP ∝ r−2.1 and the planetesimal formation runs. We find that most massive planets end up at the inner edge of the disk. The total number of planets that have reached over 100 M⊕ in the ΣP ∝ r−2.1 runs is given as 41 out of 1961. This is heavily biased by the placement of the planetary seeds, however, which also occurs in far out regions with low planetesimal surface densities.

5 Discussion

In our models with a fixed initial density slope for the planetesimals we find that we cannot form gas giant planets from planetesimal accretion with 100 km sized planetesimals when we assume that the surface density distribution of the planetesimals is shallow and varies as r^−0.9. This agrees with studies from Johansen & Bitsch (2019), in which planetesimal accretion of large planetesimals is an inefficient accretion mechanism for low-mass planetary embryos. On the other hand, we can clearly show that a change in the planetesimal surface density slope has a drastic effect on the global evolution of planetary systems. A steeper profile in the initial planetesimal surface density distribution can lead to gas giant growth in the inner region of protoplanetary disks when only 100 km sized planetesimals are used, while also forming a large number of terrestrial planets and super-Earths. This result indicates that planetesimal accretion alone can be a very effective mechanism for planetary growth in the inner regions of circumstellar disks and can explain large diversities in the population of planets.

More importantly, however, we find that pebble flux-regulated planetesimal formation leads automatically from a shallow distribution of dust to a steep planetesimal distribution. This in turn leads to much higher planetary masses than in the Σ_P ∝ r^−1.5 density profile.

The largest planets still can be formed using the Σ_P ∝ r^−2.1 density slope and reach 1062.8 M_⊕. The most massive planet in the Σ_P ∝ r^−1.5 run reaches only 166.2 M_⊕ and 7.8 M_⊕ for Σ_P ∝ r^−0.9. The maximum planetary mass for our dynamic simulation peaks at 317.1 M_⊕. By comparing the mmsn (Σ_P ∝ r^−1.5) profile with the dynamic formation model, we find that we increase the number of planets above 10 M_⊕ by 89% (from 159 to 301) and the number of planets above 20 M_⊕ by 345% (from 31 to 138) when we choose the formation of planetesimals to be consistent with the disk evolution.

We recall that the total mass in planetesimals is lowest for the dynamical planetesimal formation model because only a fraction of the dust and pebbles is transformed into planetesimals. The slope of the planetesimals that form over time, however, is steeper than the r^−1.5 slope. The total mass that is available for accretion is therefore lower in the planetesimal formation run because pebble accretion onto protoplanets is currently neglected.

We also find that our current models assuming 100 km sized planetesimals do not form cold giants around the water-ice line in any scenario due to orbital migration, although giant planets migrate through this area, as a study of the initial embryo location in Fig. 9 shows. This might also indicate that the formation of planetesimals could be enhanced by the mechanisms around the ice line, as has been predicted by Drążkowska & Alibert (2017) and Schoonenberg & Ormel (2017).

These authors suggested that sublimation and recondensation of icy pebbles at the ice line can have a drastic effect on the formation of planetesimals. This effect on planetesimal formation can be incorporated with our implementation by locally adapting the formation efficiency ɛ and promises to have a significant effect on the formation of high-mass planets around the ice line. Finally, we find that the placement of planetary embryos appears to be a strong component for giant planet formation, see Fig. 9. The effect of the starting location of planetary embryos in combination with the formation of planetesimals can be studied in future work in greater detail, including the dynamical placement of planetary seeds during the evolution of the disk. In combination with the increased planetesimal formation around the ice line and pebble accretion, these features have a drastic effect on our synthetic planet populations. We expect this to explain the abundance of cold or hot giants and the diversity of terrestrial planets.

Fig. 10 Planetary growth tracks, mass over time, and semimajor axis evolution for a giant planet system. The system that is studied leads to a gas giant planet of 997.6 M⊕ for the ΣP ∝ r−2.1 density distribution and a 281.7 M⊕ planet for the dynamic model. The other systems lead to 51.1 M⊕ for ΣP ∝ r−1.5 and 4.64 M⊕ for ΣP ∝ r−0.9. Upper panel: mass and semimajor axis change during the evolution of the system, while middle panel: growth of the embryo over time. Lower panel: semimajor axis evolution over time. On timescales of million years, the giant planet in the dynamical planetesimal formation and the ΣP ∝ r−1.5 simulation falls into the star due to tidal forces, which is no longer shown in the lower right panel.

6 Summary and outlook

Using the two-population solid-evolution model by Birnstiel et al. (2012) and the model for pebble flux-regulated planetesimal formation by Lenz et al. (2019), we have studied the effect of planetesimal formation using our model for planetary population synthesis. By comparing the dynamical planetesimal formation with different ad hoc planetesimal surface density distributions, we find strong differences for the formation of planets in the inner parts of circumstellar disks for a planetesimal size of 100 km. This can be linked directly to the steeper slope in Σ_P, as reference simulations with shallower surface density profiles show. We hereby show the effect of the planetesimal surface density distribution and formation on the population of planets. The main results of planetesimal formation for single embryo planet population synthesis are listed below

• Planetesimal accretion with 100 km sized planetesimals can be a very efficient planetary growth mechanism in the inner regions of circumstellar disks and creates a large variety of planets.

• Pebble flux-regulated planetesimal formation enables gas giant formation by accreting only 100 km sized planetesimals because of highly condensed planetesimal areas in the inner regions of circumstellar disks.

• Pebble flux-regulated planetesimal formation fails to form cold giant planets beyond the ice line. The reason is not a core-accretion timescale that is too long compared to the disk lifetimes, but orbital migration that removes the cores faster than they can grow.

• We no longer rely on an ad hoc assumption such as the mmsn model for the distribution of planetesimals in protoplanetary disks, but can start with much shallower mass distributions that agree with observations of disks around young stars.

• Dynamicplanetesimal formation increases the number of planets above 10 M_E by 89% and the number of planets above 20 M_E by 345% compared to the mmsn hypothesis.

The greatest technical advantages that the newly implemented solid evolution model brings are listed below.

• Pebble accretion can be included next to planetesimal accretion into our population synthesis framework to study their individual contributions to planetary growth.

• Locally adapting the planetesimal formation efficiency ɛ gives usthe opportunity of studying increased planetesimal formation around the ice line, or other dynamically evolving planetesimal surface density profiles such as rings in disks.

• Planetary embryo formation based on the local planetesimal surface density evolution can be incorporated.

These improvements will enable us to consistently study the full size range of planet formation in a globally coupled framework, beginning from a disk of gas and dust.

Acknowledgements

This research is supported by the German Science Foundation (DFG) under the priority program SPP 1992: “Exoplanet Diversity” under contract KL 1469/17-1, by the priority program SPP 1385 “The first ten million years of the Solar System” under contract KL 1469/4-(1-3) “Gravoturbulent planetesimal formation in the early solar system” and SPP 1833 “Building a Habitable Earth" under contract KL 1469/13-(1-2) “Der Ursprung des Baumaterials der Erde: Woher stammen die Planetesimale und die Pebbles? Numerische Modellierung der Akkretionsphase der Erde.” and DFG Research Unit FOR2544 “Blue Planets around Red Stars” under contract KL 1469/15-1. This research was also supported by the Munich Institute for Astro- and Particle Physics (MIAPP) of the DFG cluster of excellence “Origin and Structure of the Universe” and was performed in part at KITP Santa Barbara by the National Science Foundation under Grant No. NSF PHY11-25915. C.M. and A.E. acknowledge the support from the Swiss National Science Foundation under grant BSSGIO_155816 “PlanetsInTime”. C.L. was supported by the DFG as part of the priority program SPP 1833 “Building a Habitable Earth”. The authors thank Martin Schlecker and Remo Burn for many fruitful discussions.

References

All Tables

All Figures

Fig. 1 Exemplary disk evolution including our dynamical model for planetesimal formation after 0.1, 1, and 2 Myr. We showthe surface density for the dust, pebbles, planetesimals, gas, and their individual disk masses. The dashed lines refer to the initial profile of the corresponding density. This run does not contain a planetary embryo, it only evolves the disk dynamically. The total disk gas mass is given as 0.012 M⊙ with a dust-to-gas ration of 1.5% and α = 10−3. The exponential cutoff radius of the disk is at 137 au, the inner radius at 0.03 au, and the evaporation rate is given as 2.87 × 10−5 M⊙ yr−1. The planetesimal and solid-evolution parameters can be found in Table 1. The effect of the ice line is visible in the kink in the planetesimal distribution around 1 au and the effect of drift vs. fragmentation-limited pebble size in the radially varying dust-to-pebble ratio.

In the text

Fig. 2 Schematic display of the different stages in the formation model for solid-evolution development. Upper panel: previously published model from Mordasini (2018). Middle panel: currently improved version in this work, including the two-population solid evolution for dust and pebbles, as well as the formation of planetesimals (see Sect. 2). Lower panel: outlook on possible future development stages. The new modules and functions are highlighted in red, and future possibilities are highlighted in blue.

In the text

Fig. 3 Planetesimal disk masses within 1 au, 10 au, and the complete disk for three different analytic density slopes and the dynamically formed planetesimal mass. The analytic masses are given at the start of the simulation, while the dynamically formed disk masses are shown after one million years, after most planetesimals have already formed. The dynamic runs do not contain a planetary embryo, they only simulate the disk evolution. The disk parameters, however, are the same as in the population in Fig. 5. We show the mass in planetesimals in the whole disk in the upper panel, the planetesimal mass within 10 au in the middle panel, the planetesimal mass within 1 au in the lower panel, and the corresponding median masses for every setup.

In the text

	Fig. 4 Cumulative distribution of gas-disk lifetimes for our synthetic population. We used a constant value of α = 10−3 in our runs. We find that 50% of the lifetimes are shorter than 4.1 Myr.
In the text

Fig. 5 Mass vs. semimajor axis of synthetic planet populations for different ΣP distributions after 100 million years. Each setup contains one single planetary embryo. The initially set distributions for ΣP are ΣP = Σ0 ⋅ r−0.9 (initial gas density slope), ΣP = Σ0 ⋅ r−1.5 (mmsn) and ΣP = Σ0 ⋅ r−2.1 (Lenz et al. 2019). Bottom right panel: population in which planetesimals form over time using the model described in Sect. 2. The circles given around the data points show the mass fraction of envelope mass over core mass. The numbers of systems are 1999 (ΣP ∝ r−0.9), 1990 (ΣP ∝ r−1.5), 1961 (ΣP ∝ r−2.1), and 1945 (ΣP- dyn).

In the text

Fig. 6 Planetary and core-mass occurrences of the four different populations. Every planet in each systems starts with a core mass of 0.0123 M⊕ and no envelope. The quantities that arise from the three analytical planetesimal surface density profiles are shown in blue (ΣP ∝ r−2.1), orange (ΣP ∝ r−1.5), and green (ΣP ∝ r−0.9), whereas the properties of the planetesimal formation population are shown in black. The dashed lines in the plots showthe median planet and median core masses. The histograms show clear shifts toward the higher mass ranges for steeper planetesimal surface densities and for the dynamically formed planetesimals than the ΣP ∝ r−0.9 or even the ΣP ∝ r−1.5 distribution.

In the text

	Fig. 7 Planetary mass and core-mass occurrences of the four different populations within 1 au. Figure 5 shows that the higher mass planets are located in the inner disk because these are the regions with the highest planetesimal surface density. The dashed lines in the plots show the median planet and median core masses within 1 au.
In the text

Fig. 8 Cumulative function of planetary masses within 1 au for the four different synthetic populations from Fig. 5. The y-axis shows how many planets for each density profile are above the current planetary mass, normalized by the total number for planets in each population. The planets that have formed in the ΣP ∝ r−0.9 run are shown in green, the ΣP ∝ r−1.5 population is shown in orange, and the ΣP ∝ r−2.1 population is shown in blue. The dynamic planetesimal formation population is shown in black.

In the text

Fig. 9 Semimajor axis distribution and starting location of planets that have grown to gas giant masses (MP >100 M⊕) in the ΣP ∝ r−2.1 and the planetesimal formation runs. We find that most massive planets end up at the inner edge of the disk. The total number of planets that have reached over 100 M⊕ in the ΣP ∝ r−2.1 runs is given as 41 out of 1961. This is heavily biased by the placement of the planetary seeds, however, which also occurs in far out regions with low planetesimal surface densities.

In the text

Fig. 10 Planetary growth tracks, mass over time, and semimajor axis evolution for a giant planet system. The system that is studied leads to a gas giant planet of 997.6 M⊕ for the ΣP ∝ r−2.1 density distribution and a 281.7 M⊕ planet for the dynamic model. The other systems lead to 51.1 M⊕ for ΣP ∝ r−1.5 and 4.64 M⊕ for ΣP ∝ r−0.9. Upper panel: mass and semimajor axis change during the evolution of the system, while middle panel: growth of the embryo over time. Lower panel: semimajor axis evolution over time. On timescales of million years, the giant planet in the dynamical planetesimal formation and the ΣP ∝ r−1.5 simulation falls into the star due to tidal forces, which is no longer shown in the lower right panel.

In the text