The impact of pebble flux regulated planetesimal formation on giant planet formation

Forming gas giant planets by the accretion of 100 km diameter planetesimals, a typical size that results from self-gravity assisted planetesimal formation, is often thought to be inefficient. Many models therefore use small km-sized planetesimals, or invoke the accretion of pebbles. Furthermore, models based on planetesimal accretion often use the ad hoc assumption of planetesimals distributed radially in a minimum mass solar nebula fashion. We wish to investigate the impact of various initial radial density distributions in planetesimals with a dynamical model for the formation of planetesimals on the resulting population of planets. In doing so, we highlight the directive role of the early stages of dust evolution into pebbles and planetesimals in the circumstellar disk on the following planetary formation. We have implemented a two population model for solid evolution and a pebble flux regulated model for planetesimal formation into our global model for planet population synthesis. This framework is used to study the global effect of planetesimal formation on planet formation. As reference, we compare our dynamically formed planetesimal surface densities with ad-hoc set distributions of different radial density slopes of planetesimals. Even though required, it is not solely the total planetesimal disk mass, but the planetesimal surface density slope and subsequently the formation mechanism of planetesimals, that enables planetary growth via planetesimal accretion. Highly condensed regions of only 100 km sized planetesimals in the inner regions of circumstellar disks can lead to gas giant growth. Pebble flux regulated planetesimal formation strongly boosts planet formation, because it is a highly effective mechanism to create a steep planetesimal density profile. We find this to lead to the formation of giant planets inside 1 au by 100 km already by pure planetesimal accretion.


Physical background
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 . The accretion of smaller objects with a higher effective cross section, like either km-sized planetesimals (Ida & Lin 2004) or 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 restrain to make a statement on the efficiency of pebble accretion, the scenario of a planetary core accreting 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), thus a pebble accreting embryo could well have formed from planetesimal collisions. This crucial step adds room to discuss 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 their size is an active field of research. While some studies infer that the current size of asteroid belt objects is well constrained and found to be in the order of magnitude of 100 km in diameter (Bottke Jr et al. 2005;Walsh et al. 2017;Delbo' et al. 2017), other studies find that the size distribution found today merely reflects which sizes are most resilient to clearing, and therefore suggest a smaller primordial size (Zheng et al. 2017). The observed size distribution could also arise from the growth of planetesimals of originally 100 m in size (Weidenschilling 2011). In the Kuiper belt, the size distribution has a similar shape as predicted by simulations including the streaming instability between 10 and 100 km (Schäfer et al. 2017), indicating large initial sizes. On the other hand, recent discoveries of Kuiper belt objects via stellar occultations rather indicate a size of 1-2 km (Arimatsu et al. 2019). Small initial sizes of 0.4 -4 km are also inferred theoretically by Schlichting et al. (2013). Also, the surface density profile of planetesimals for extrasolar systems is unknown. Studies of our own solar system motivated the minimum mass solar nebula (mmsn) hypothesis (Weidenschilling 1977;Hayashi 1981) that results in a power law drop of the planetesimal surface density with a decay of Article number, page 1 of 14 arXiv:2004.03492v1 [astro-ph.EP] 7 Apr 2020 A&A proofs: manuscript no. Consistent_planetesimal_formation Σ 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 however is not planetesimals, but the dust in the circumstellar disk, as the distribution of planetesimals in protoplanetary disks is currently unobservable. Lenz et al. (2019) model 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 find 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 orders of magnitude without increasing the total mass in planetesimals. Since the accretion rate of planetesimals is proportional to the local planetesimal surface density, these highly condensed planetesimal zones are promising to have a drastic effect on planetary growth.

Previous models
Before discussing some of the previous work, we would like to 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 yield vice versa. Key to the population synthesis approach is that the model is complex enough to take into account the physical effects that are deemed crucial for planet formation, yet its single system computational cost is low enough so that it can 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 like Johansen & Bitsch (2019), Mordasini (2018) or Ida & Lin (2004) all use initial distributions of planetesimals and initially placed planetary embryos, while neglecting the presence of pebbles. Other formation models like Bitsch et al. (2015) or Brügger et al. (2018) model planetary growth by the accretion of pebbles and initially set planetary embryos, while neglecting the formation, or accretion of planetesimals. Yet, 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 have chosen to thus improve our planet population synthesis model by a "disk consistent" 1 model for solid evolution (Birnstiel et al. 2012) and planetesimal formation (Lenz et al. 2019) to take the early stages of the disks evolution into account. This early phase determines the planetesimal surface density distribution, the radial pebble flux evolution, the formation of planetary cores and therefore planet formation as a whole. For our study, we will focus on the formation and accretion of planetesimals. We will display the impact of the planetesimal surface density and its formation on the population of planets. We will show that the accretion by 100 km sized planetesimals is in fact a highly efficient growth mechanism for planets, due to the highly condensed planetesimal regions in the disk. Furthermore we will give an overview over the future possibilities that arise from our newly implemented modules. This paper is outlined as followed: in Sect. 2 the planetesimal formation model is explained, as well as the newly implemented solid evolution model on which it is based on. Sect. 3 will give insight into the population synthesis framework and how it was modified for our purpose. The changes in Σ P in the population synthesis code, as well as the newly computed synthetic populations are presented in Sect. 4. Sect. 5 will discuss the results followed by a brief summary and an outlook on our new possibilities and future work in Sect. 6.

The 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). This chapter gives a brief outline of the assumptions and displays the most important reasons why we have chosen to use it in our framework. Our goal is to implement a fast computing, one dimensional, parameterized algorithm for solid evolution that is well tested and in good agreement with more sophisticated models. Key of the performance of the two population approximation is a parameterized mass ratio f m (r) as a function of orbital distance r between two populations of solids, that depends on whether the growth of the particles is limited by drift or by fragmentation. Each time step, the model solves one advection-diffusion equation given by with Σ s as the total solid surface density without planetesimals, Σ g as the gas surface density and D g as the gas diffusion coefficient and t and r as time and radial distance.ū describes the weighted velocity of the total solid density and is defined as where f m is the fit parameter for the mass ratio between the two populations. u 0 and u 1 describe their velocities, while the surface densities of the two populations are given as 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. Σ 0 describes the smaller population, that can be seen as dust, subject to diffusion and transport with the gas, while Σ 1 describes the larger population, that can be seen as pebbles , which on top of being diffused by the gas are also sedimenting towards the midplane and drifting towards pressure maxima, for instance towards the star. The fit parameter f m has been 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 = 0.97, drift-limited case 0.75, fragmentation-limited case .
These are also the ones that we used in our simulations. One can see the effect of this implementation in Fig. 3, where the ratio between dust and pebbles varies with space and time, visible in the two blue curves. The full model and its results are described in Lenz et al. (2019) in greater detail, we thus outline here only the basic physical assumptions behind 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 drifting by into planetesimals. Thus in principle it acknowledges that pebbles want to drift inward and that one can form more planetesimals if more material comes by. Many different planetesimal formation prescriptions can therefore be parameterised in such a fashion. Be it in the frame work of turbulent clustering (Cuzzi et al. 2010;Hartlep & Cuzzi 2020) or streaming instabilities (Johansen et al. 2009;Schäfer et al. 2017) or local trapping in zonal flows (Johansen et al. 2007(Johansen et al. , 2011Dittrich et al. 2013; or in vortices (Raettig et al. 2015;Lyra et al. 2018), the formation is always limited by how much a region receives in fresh pebbles, after consuming the locally available ones. Thus our parameterisation is per se 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 are slowly evolving coherent flow structures in protoplanetary disks like vortices and zonal flows , which can form everywhere, which live only for a limited time and thus only trap a fraction of drifting pebbles. In these traps, pebbles get sufficiently concentrated, that regulated by streaming and Kelvin Helmholtz instabilities, planetesimal formation will be triggered. The planetesimal formation rate is generally proportional to the radial pebble fluẋ with v drift as the drifting velocity of the particles and St min and St max as the minimal and maximal Stokes number for which a particle is considered a pebble. v drift is given as with P(r) as gas pressure, h g (r) as gas pressure scale height (h g (r) = c s (r)/Ω(r)) and c s (r) as sound speed. Ω(r) is given as the orbital frequency at the radial distance r. The source term for planetesimals, i.e. for Σ P is then given as (Lenz et al. 2019) with d(r) as the radial separation of the pebble traps and as the efficiency parameter, that describes how much of the pebble flux is transformed into planetesimals after drifting over a distance of d. We choose a constant value of = 0.05 as a good value to form a sufficient number of planetesimals as 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, if the formation of planetesimals might follow a different underlying mechanism, like, e.g., around the water iceline as described by  or Schoonenberg & Ormel (2017). This flexibility allows us to study a broad range of planetesimal formation scenarios, using the same implementation. Right now our two-poppy implementation has no proper treatment of the processes of evaporation and possible recondensation. The only effect of the existing iceline is incorporated into the parameter f ice (T ): in effect to reduce the pebble flux inside the iceline to compensate for the evaporation of water ice. Therefore the ice line is visible in the distribution of planetesimals, even so it is not visible in the pebbles themselves (See Fig. 3). We also use a fixed planetesimal size of 100 km in diameter as in (Lenz et al. 2019). As a consequence, there is a threshold of transformed mass to be reached to build at least one planetesimal. From this we can derive a critical pebble flux, that is necessary for Σ P to change. It is given as (Lenz et al. 2019) where τ t describes the average lifetime of a trap, which is given as 100 local orbits and m P the mass of a single planetesimal. For simplicity we assume spherical planetesimals with a uniform density of ρ s = 1.0g/cm 3 . 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 where the sink term L is defined as and where θ (·) is the heavyside function. This combines the above mentioned conditions for planetesimal formation. The surface density can only change while a critical mass is transformed (θ (Ṁ peb −Ṁ cr )) and if the Stokes numbers of the particles are within St min and St max (θ (St 0 − St min ) · θ (St mas − St 0 )).

The planet formation and evolution model
The most up to date version of our planet population synthesis model can be found in Emsenhuber et al. (in prep.), which corresponds to an update of the model presented in Mordasini (2018). This model combines planet formation (Alibert et al. 2005 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 Emsenhuber et al. (in prep.). We will provide here only an overview of the physical processes that are tracked in the model, while focusing on the solid components of the protoplanetary disk model in Sect 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 α parametrization (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 away 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 planet's mass and once 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 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 life time of the protoplanetary disk, but at least 10 Myr. Once this is passed, the model switches to the evolution stage (Mordasini et al. 2012) where 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 use a method similar to Mordasini et al. (2009), with several adaptations. 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 is fixed to 0.03 au. The initial embryo mass is 0.0123 M ⊕ and the location is random with a uniform distribution in the logarithm of the distance between 0.06 and 40 au. Embryos are placed directly at the beginning of the simulations.

The solid component
A schematic overview of the different modules can be seen in Fig. 1. Previous generations of the model including the upcoming Emsenhuber et al. (in prep.) use 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 mmsnlike distribution (Emsenhuber et al. in prep.). 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 is chosen to be uniform and with a radius of r P = 300 m. Important to point out is that Σ 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 go beyond the standard implementation on Emsenhuber et al. (in prep.). We now included two additional solid quantities (dust & pebbles) that are evolving along with the gas evolution of the disk model. The initial mass in dust and pebbles is given as the metalicity 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 also do planetesimals now form over time, which opens a completely new level of dynamical interaction with the disk. The size of planetesimals that we assume in the following simulations is given as 100 km in diameter.
Thus the main differences between Emsenhuber et al. (in prep.) and this paper is 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. (in prep.). Emsenhuber et al. (in prep.) on the other side includes an N-body integration for multiple simultaneously evolving cores, an option we have not used in the present paper, as we wanted to focus on the effect of dynamical planetesimal formation.

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 , generally it depends on the individual evolution of the disk. Due to the steeper slope, we find a remarkable increase of Σ P in the inner regions of a protoplanetary disk and a corresponding decrease further 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, yet only a fraction of that is transformed into planetesimals. We will 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 Fig. 3: Exemplary disk evolution including our dynamical model for planetesimal formation after 0.1 Myrs, 1 Myrs and 2 Myrs. We show the 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 /year. The planetesimal and solid evolution parameters can be found in table 1. Note the effect of the iceline 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.  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). The bottom right panel shows the population in which planetesimals form over time using the model described in Sec. 2. The circles given around the datapoints show the mass fraction of envelope mass over core mass. The number of systems are 1999 (Σ P ∝ r −0.9 ), 1990 (Σ P ∝ r −1.5 ), 1961 (Σ P ∝ r −2.1 ) and 1945 (Σ P -dyn). as a fraction of the gas surface density. Considering /d > 1, this would lead to local pebble-to-planetesimal conversion and the outer material could not drift into the inner regions of the disk, which would have changed the density profile. For a more detailed treatment of this behaviour 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 can be seen in Fig. 7, as well as the masses within 10 au and 1 au.
We find that the mean total disk masses are lower for the steeper density profiles by a factor of M −2.1 tot /M −0.9 tot ≈ 0.62 or M −1.5 tot /M −0.9 tot ≈ 0.87. This is to be expected as more material is inside the iceline, which is taken care of in these models. Still the masses within 1 au of the steeper models are by orders of magnitude higher (M −2.1 1au /M −0.9 1au ≈ 21, 58 or M −1.5 1au /M −0.9 1au ≈ 6, 01). M tot and M 1au refer to the median masses from Fig. 7. The lowest total median mass ratios of planetesimals can be found in the dynamically formed simulation with M dyn tot /M −0.9 tot = 0.504, the mass ratio within 1 au however is the second highest with M dyn 1au /M −0.9 1au = 8.27. The smaller total masses for the steeper planetesimal surface density can be explained by the smaller amount 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 further out regions to the inner disk. This would evaporate the icy planetesimals within the iceline, leaving only the silicate planetesimals, therefore effectively loosing mass. The mass loss here is therefore only due to icy planetesimals within the iceline, whereas the amount of silicate planetesimals stays the same. Regardless of this mass loss, we find that the mass in the inner disk (r < 1 au) is signigicantly higher for the steeper density slopes. The lifetimes of the gas disks studied in our case can be seen in Fig. 2 The global effect on planet formation of these changes in Σ P is presented in Sect. 4.2.

Synthetic populations
In the following we will present several synthetic populations that have been 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 Article number, page 6 of 14 growth of planetary embryos by the accretion of solids is only given by the accretion of planetesimals in these simulations. To ensure the right 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. in prep.), we chose to stay with 1 embryo per run for our study, as it makes no sense to mix our study with effects of multiple planets in that forthcoming paper. Therefore, we focus on the general distribution of masses and semimajor axes and the overall mass occurrences of planets. Figure 4 shows the mass and semi major axis distribution of four synthetic populations around a solar type star using 1 lunar mass (0.0123 M ⊕ ) planetary embryo for each system. We simulate 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 for the dynamic planetesimal formation run. The initial conditions of the four populations are the same, except for the initial Σ P and the formation of planetesimals respectively. 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 of 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 heavy 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 further out like e.g. beyond the water iceline as it can be observed in the population of exoplanets (Winn & Fabrycky 2015), which will probably change once we will incorporate recondensation of water vapor, effectively boosting the pebble flux regulated birth of planetesimals.

Mass occurences
For a more quantitative analysis we study the mass occurrences for the different planetesimal density slopes. Here we focus on the planetary mass and the core mass. Fig. 8 and Fig. 9 show histograms with the occurrences of the different masses for the various populations from Fig. 4. As Fig. 4 shows, most of the high masses are found in the inner parts of the protoplanetary disk, whereas the outer placed embryos fail to grow. We therefore also focus our study on the inner region within 1 au. Fig. 8 takes the complete population into account whereas Fig. 9 only contains planets with a semimajor axis below 1 au. We also give the median masses for the planets and their cores. A cumulative function of the planetary masses is shown in Fig. 5. 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  Fig. 4. The yaxis shows how many planets for each density profile are above the current planetary mass, normalised 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.

Influence of the starting location
In Fig. 6 we see 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 most heavy planets 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 au to 10 au for the Σ P ∝ r −2.1 run and mostly from around 4 au to 10 au for the dynamical formation model. Embryos that are placed at a distance from 2 au to 4 au appear to have a lower probability to become gas giant planets in both cases, but as the probability of their formation at this location is also low, due to 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 model in the future also the generation of embryos in a consistent fashion. Ultimately the above mentioned effect of recondensation beyond the iceline can further change this picture.

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. Fig 10 shows planetary growth tracks, the mass growth over time and the corresponding semimajor axis evolution. The embryo in these sys-Article number, page 7 of 14 A&A proofs: manuscript no. Consistent_planetesimal_formation Fig. 6: Semimajor axis distribution and starting location of planets that have grown to gas giant masses (M P > 100M ⊕ ) 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 100M ⊕ in the Σ P ∝ r −2.1 runs is given as 41 out of 1961. This however is heavily biased by the placement of the planetary seeds which, also occurs at far out regions with low planetesimal surface densities.
tems was placed initially at 8.2 au which seems to be a preferential starting location for giant planets, see Fig. 6. We can see that the higher planetesimal surface density has a drastic impact 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 ⊕ .

Discussion
In our models with a fixed initial density slope for the planetesimals we find that we can not form gas giant planets from planetesimal accretion with 100 km sized planetesimals, if we assume that the surface density distribution of the planetesimals is shallow, varying as r −0.9 . This is in agreement with studies from Johansen & Bitsch (2019), in which planetesimal accre-tion of large planetesimals is an inefficient accretion mechanism for low mass planetary embryos. Yet 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, using only 100 km sized planetesimals, while also forming a large amount 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. But more importantly we find that pebble flux regulated planetesimal formation leads automatically from a shallow distribution of dust to a steep planetesimal distribution, leading 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 ⊕ . 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), if we choose the formation of planetesimals to be consistent with the disks evolution.
One has to keep in mind that the total mass in planetesimals is the lowest for the dynamical planetesimal formation model, since 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, since pebble accretion on protoplanets is currently neglected.
We also find that our current models assuming 100 km sized planetesimals do not form cold giants around the water iceline in any scenario due to orbital migration, although giant planets migrate trough that area, as a study of the initial embryo location in Fig. 6 shows. This might also indicate that the formation of planetesimals could be enhanced by the mechanisms around the iceline, as already predicted by  and Schoonenberg & Ormel (2017). They suggest, that sublimation and recondensation of icy pebbles at the iceline can have a drastic effect on the formation of planetesimals. This effect on planetesimal formation can be incorporated with our implemenatation by locally adapting the formation efficiency and promises to have a significant impact on the formation of heavy planets around the iceline. Finally, we find that the placement of planetary embryos appears to be a strong component for giant planet formation, see Fig. 6. 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 iceline and pebble accretion, we believe these features to have a drastic impact on our synthetic planet populations. We expect this to explain the abundance of cold/hot giants and terrestrial planet diversity.

Summary & Outlook
Using the two population solid evolution model by Birnstiel et al. (2012)  effect of planetesimal formation using our model for planetary population synthesis. 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 cirscumstellar 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 impact of the planetesimal surface density distribution and formation on the population of planets. The main results of planetesial formation for single embryo planet population synthesis are: -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, due to highly condensed planetesimal areas in the inner regions of circumstellar disks.
-Pebble flux regulated planetesimal formation fails to form cold giant planets outside the iceline. The reason however is not a too long core accretion timescale 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 like the mmsn model for the distribution of planetesimals in protoplanetary disks, but can start with much shallower mass distributions in agreement with observations of disks around young stars.
-Dynamic planetesimal formation increases the amount 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 biggest technical advantages that the newly implemented solid evolution model brings are: -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 us the opportunity to study increased planetesimal formation around the iceline, or other dynamically evolving planetesimal surface density profiles, like e.g. 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. Fig. 7: 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. 4. 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. 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 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 show the median planet and median core masses. The histograms show clear shifts towards the higher mass ranges for steeper planetesimal surface densities and for the dynamically formed planetesimals, compared to the Σ P ∝ r −0.9 or even the Σ P ∝ r −1.5 distribution.
Article number, page 11 of 14 A&A proofs: manuscript no. Consistent_planetesimal_formation Fig. 9: Planetary mass and core mass occurrences of the four different populations within 1 au. From Fig. 4 we see that the heavier planets are located in the inner disk since 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. We choose to focus on the planetary masses within 1 au to neglect the "failed" planetary cores that are randomly placed at far out regions of the disk and do not grow substantially from the initial lunar mass.  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 . The upper panel shows the mass and semimajor axis change during the evolution of the system, while the middle panel shows the growth of the embryo over time. The lower panel shows the semimajor axis evolution over time. Over Gyr timescales 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.