Constraining the parameter space for the Solar Nebula

If we want to understand planetesimal formation, the only data set we have is our own Solar System. It is particularly interesting as it is so far the only planetary system we know of that developed life. Understanding the conditions under which the Solar Nebula evolved is crucial in order to understand the different processes in the disk and the subsequent dynamical interaction between (proto-)planets, once the gas disk is gone. Protoplanetary disks provide a plethora of different parameters to explore. The question is whether this parameter space can be constrained, allowing simulations to reproduce the Solar System. Models and observations of planet formation provide constraints on the initial planetesimal mass in certain regions of the Solar Nebula. By making use of pebble flux-regulated planetesimal formation, we perform a parameter study with nine different disk parameters like the initial disk mass, initial disk size, initial dust-to-gas ratio, turbulence level, and more. We find that the distribution of mass in planetesimals in the disk depends on the planetesimal formation timescale and the pebbles' drift timescale. Multiple disk parameters can influence pebble properties and thus planetesimal formation. However, it is still possible to draw some conclusions on potential parameter ranges. Pebble flux-regulated planetesimal formation seems to be very robust, allowing simulations with a wide range of parameters to meet the initial planetesimal constraints for the Solar Nebula. I.e., it does not require a lot of fine tuning.


Introduction
In order to form planets, tiny micron sized dust grains have to grow to hundreds or thousands of kilometers. First, grains grow by collisions with other grains. But at some point they cannot continue to grow because either (1) the relative velocities become so high that a collision leads to fragmentation (Blum & Münch 1993;Blum & Wurm 2008;Gundlach & Blum 2014), or (2) they start to drift faster toward the star, than they could potentially grow (Klahr & Bodenheimer 2006;Birnstiel et al. 2012). There is also the bouncing barrier (Zsom et al. 2010), but charging effects might lead to growth to to sizes an order of magnitude above this barrier's limit (Steinpilz et al. 2019). Laboratory experiments point in the direction of low fragmentation speeds for icy particles (Musiolik & Wurm 2019) of around 1 m s −1 . This could cause particles to hit the fragmentation barrier first, which is why we are not considering the bouncing in this paper.
It is believed that planets are formed by so-called planetesimals, of a few to hundreds of kilometers in size. These planetesimals are the building blocks of planets. Once they have formed, accretion of pebbles (Ormel & Klahr 2010) may become important too (for a review see, e.g., Ormel 2017). But since grains stop growing at some point, continuous growth leaves a Member of the International Max Planck Research School for Astronomy and Cosmic Physics at the Heidelberg University missing link between roughly mm-dm to planetesimal diameters (∼ 100 km) via continuous growth.
If there is a pressure bump in the disk, e. g., caused by a vortex or zonal flow, particles can get trapped around the center of the bump (Whipple 1972;Barge & Sommeria 1995). After the accumulation of enough pebbles, the streaming instability (Youdin & Goodman 2005) can potentially be the dominant turbulent process to trigger fragmentation in the laminar case (α t = 0), or in the turbulent case (α t > 0) gravoturbulent planetesimal formation can occure (Johansen et al. 2006;Johansen et al. 2007).
To summarize, we follow the idea that pebbles form planetesimals in a gravoturbulent process leading to a Gaussian-like size distribution of planetesimals that peaks around 100 km in diameter (Klahr & Schreiber 2015;Schreiber 2018). These planetesimals can then build planetary embryos, which can grow via further accretion of both planetesimals and pebbles to form (proto-)planets.
The initial distribution of planetesimals is one of the biggest unknowns in planet formation models. The radial distribution in the disk is important for embryo formation and subsequent accretion of planetesimals onto embryos. While we can observe protoplanetary disks and debris disks around other stars, we cannot observe planetesimal populations. There is only one system where we have relatively good knowledge of the present day small body population -the Solar System. This is also

Mass Constraints for the Initial Planetesimal Population
In this section we will review literature studies to infer constraints on the initial planetesimal mass in various regions of the Solar Nebula. A summary of this literature research is depicted in Fig. 1.

Mercury Region: Interior to 0.7 au
There is no observed stable population of asteroids within the orbit of Mercury (Steffl et al. 2013), even though e.g. Campins et al. (1996) found a dynamically stable region in the Solar System's inner region. On top of that, no models for terrestrial planet formation require any planetesimals inside of Mercury's orbit (∼ 0.4 au). Additionally, terrestrial planet models fit observations better if the planetesimal disk is truncated around 0.7 AU (Hansen 2009;Walsh et al. 2011;Levison et al. 2015;Morbidelli et al. 2016). Even though there are some models suggesting how to clear out this region (e.g. Ida & Lin 2008;Batygin & Laughlin 2015;Volk & Gladman 2015), none of these are demonstrated in a comprehensive way. However, this implies that we can't define an upper limit on the mass of initial planetesimals in this region. Low mass short-period planets around other stars may indicate that in other planetary systems there were planetesimals in short period orbits that formed planets via in situ planetesimal accretion (e.g. Chiang & Laughlin 2013;Hansen & Murray 2012) or pebble accretion (e.g. Chatterjee & Tan 2013), but the presence of migration implies that there is no evidence that a population of initial planetesimals was present within Mercury's orbit in our own Solar System. Bottom line: 0 to an unknown upper limit 2.2. Earth/Venus Region 0.7 − 1 au As it was shown by Hansen (2009), placing ∼ 2 M ⊕ of oligarchs within 0.7 − 1 au can lead to good matches to the sizes and spacing of the terrestrial planets of the Solar System. The results from the parameter study of Kokubo et al. (2006), in terms of the number of Earth-like planets within 0.5 to 1.5 au, seem to be very robust with respect to mass and radial profile of oligarchs. They have shown that one can still obtain reasonable results with a total initial mass of 2.77 M ⊕ in that region. For ∼ 23 M ⊕ , Super-Earths form.
Pebble flux can allow a significant amount of mass to be transported into this region. If there are enough pebbles and appropriate disk structure, it is possible to produce reasonable Solar System analogs beginning with < 10 −2 M ⊕ planetesimal masses (Levison et al. 2015). If Jupiter migrated inwards and then out of the Asteroid Belt (known as the grand tack; Walsh et al. 2011), to leave about the needed mass in the current Asteroid Belt it would have implanted ∼ 1 M ⊕ of material into the Earth/Venus forming region. This suggests that either primordially or after early pebble accretion there was initially 1 M ⊕ of planetesimals/embryos in this region (Walsh et al. 2011). Since the grand tack removes most objects in the Asteroid Belt in simulations, the inital planetesimal population needs to be massive enough early on. This probably implies that pebble accretion only grew the mass of the Asteroid Belt by a factor of a few at most, and thus the terrestrial planet region probably only grew by a factor of a few as well. As a rough estimate we can assume the inital planetesimal mass to be around 0.1 M ⊕ .
Bottom Line: 0.1 (if there are enough pebbles that can be accreted) to 2.77 M ⊕ (if there aren't)

Asteroid Belt: 2 − 3 au
The Asteroid Belt currently has a mass of about 5 · 10 −4 M ⊕ (e.g. Kresak 1977), where roughly 50% of the mass is in the 4 largest objects, with 1/3 in Ceres. Over the history of the Solar System, it has potentially been depleted by the following effects. (1) Over the last 4 Gyr dynamical chaos in the current structure of the Solar System has removed about ∼ 50% of the mass of the Asteroid Belt (Minton & Malhotra 2010). Vesta's crust indicates that the Asteroid Belt population was only modestly larger than it is today at the time the mean collision velocities were pumped up to 5 km/s (i.e. the current mean impact velocity in the main belt region ;Bottke Jr et al. (1994)). If the Asteroid Belt had substantially more 30 km-sized planetesimals in it over the last 4 Gyr than it has today, Vesta would be expected to have more than 1 large basin (Bottke Jr et al. 2005a,b;O'Brien & Greenberg 2005). If planetesimals were "born big" (Morbidelli et al. 2009), i.e. 80 km, this suggests that collisional evolution should not be particularly important in removing material; at least not more than a factor of a few. The late Jupiter-Saturn interaction in terms of reshuffling of the Giant planets likely depleted the Asteroid Belt by a factor of ∼ 2 to ∼ 10 (Minton & Malhotra 2010). If the Grand Tack happened (Walsh et al. 2011), then only few times 10 −3 to a few times 10 −4 of the population would have survived (ignoring newly implanted plantesimals from other regions; see e.g. Fig. 7 in Morbidelli et al. 2015). One should also note that in this model the C-complex asteroids, which comprise about 75% of the asteroids (Gradie et al. 1989) and include Ceres, are implanted from the outer Solar System. But these modifications are swamped by the uncertainty in the clearing rate on the migration efficiency. If pebble accretion plays a crucial role for the growth of large asteroids (> 200 km in diameter), then this also would reduce the "initial mass" of planetesimals that is needed. This effect probably would be a factor of ∼ 2 (Johansen et al. 2015).
Bottom line: ∼ 2 · 10 −3 M ⊕ (4 times current mass) to ∼ 5M ⊕   Fig. 1, but converted into column densities, assuming a power-law shape ∝ r −2.25 (see Eq. (37) of Lenz et al. (2019) for a motivation). The blue (orange) line shows (three times) the minimum mass Solar Nebula profile for solids (Weidenschilling 1977b;Hayashi 1981). core has at least 7M ⊕ (Wahl et al. 2017). So far, we do not know how much of this mass was originally in planetesimals from the appropriate region. Hence, for the lower limit, we will ignore Saturn and take 50% (assuming that the other half stems from pebbles, Bordukat 2019) of the lower mass estimate from Wahl et al. (2017) which gives 3.5M ⊕ .
In order to reach critical masses for strong gas accretion, around 5 times the mass of the MMSN seems to be needed (Thommes & Duncan 2006, these simulations used planetesimals of 20 km in diameter). This leads to 132 M ⊕ within 4 to 15 au. Assuming that 50% of the mass was contributed by pebbles (Bordukat 2019), the lower limit would be around 66 M ⊕ . This lower limit is still high, but desublimation effects just outside the ice line can lead to a pile-up in planetesimals by a factor of ∼ 5 (Drążkowska & Alibert 2017;Schoonenberg et al. 2018). This effect is not included in this work. Raymond & Izidoro (2017) found that ∼ 10% of the asteroids around the giant planets were scattered into the asteroid belt. So in order to explain the mass of C-type asteroids in the Asteroid Belt, there probably had to have been around few times 10 −2 M ⊕ of asteroids in the giant planet forming region. This value would be the absolute minimum for this region.
Bottom line: ∼ 66 M ⊕ to unknown high mass 2.5. The Nice Disk ∼ 15 − 30 au In order to match the observed structure of the Kuiper Belt, where many objects are in resonance with Neptune, outward migration of the giant planets is the preferred explanation. For this type of outward migration, planetesimal driven migration is the leading explanation (Fernandez & Ip 1984;Malhotra 1995;Thommes et al. 1999). the Nice Model Morbidelli et al. 2005;Gomes et al. 2005) is a comprehensive model explaining how this outward planetesimal driven migration could have occurred. In this model the giant planets initially formed closer to the sun than their current locations, and migrated outwards due to interactions with a planetesimal disk known as "the Nice disk". Even though the details of the model have been changed (e.g. Morbidelli et al. 2007;Levison et al. 2011 detailed characteristics of the objects in the 3:2 resonance (Nesvornỳ & Vokrouhlickỳ 2016). -Comets: Gas drag prevents km-sized planetesimals from being scattered into the Oort Cloud while the gas disk is still around (Brasser et al. 2007). This suggests that the long period comets were scattered into the Oort cloud after the gas disk went away, a natural outcome of something like the Nice Model. An initial population of around 30 M ⊕ is needed to populate the Oort Cloud (Dones et al. 2004), however the existence of more massive planets in the inner Oort Cloud (e.g. a planet 9, Batygin & Brown 2016) could decrease the required reservoir size. However, since comets could have potentially been shared between stars in the birth cluster under favorable conditions, it is possible comets are not a reliable constraint (Levison et al. 2010). -Ice giant ejection: In the models in which an ice giant is ejected from the Solar System, the best overall structure of the Solar System, for example surviving terrestrial planets, needs ∼ 20 M ⊕ in planetesimals (Nesvornỳ & Morbidelli 2012).
Additionally, it has been found that the column density profile of planetesimals has a minimal effect on the outcomes for a relatively broad range of power-laws. This was tested by Batygin & Brown (2010) for Σ p ∝ r −k , where k = 1 . . . 2, with inner edge r in ∼ 12 au and outer radius r out = 30 au. The Nice scattering occurred after the disk went away so that pebble accretion could have increased the total mass of the Nice disk. However, because the disk was very likely flaring in the outer regions, it is unlikely that pebble accretion was efficient and increased the total mass in these regions more than a factor of ∼ 2 (Lambrechts & Johansen 2012). Bottom line: ∼ 10 M ⊕ seems to be needed 2.6. The Cold Classical Kuiper Belt ∼ 30 au − 50 au There is a population of objects in the Kuiper belt with low eccentricities and inclinations which look as though they are not transplanted, though likely to be primordial. Observations indicate that the mass of the current classical population is 8·10 −3 M ⊕ (Fuentes & Holman 2008). If the larger KBOs were formed by coagulation from small planetesimals, there must have been significantly more mass in this region in the early stages. For instance, Pan & Sari (2005) suggested that the high end size distribution could be matched by collisional evolution. However, if one combines more modern description laws with the need to preserve wide binaries, one cannot match the observed population in this type of collisional environment. Therefore, this suggests that planetesimals formed as large bodies and that the total mass of the cold classcial Kuiper Belt objects (CKBO), dominated by bodies larger than diameters of ∼ 100 km has not evolved significantly ).
The Nice migration may have dynamically depleted the Kuiper Belt by up to an order of magnitude (Morbidelli et al. 2008). Singer et al. (2019) found a lack of small craters on Pluto and Charon, indicating that planetesimals in the Kuiper Belt are not a collisionally evolved population, or that collisions destroyed small planetesimals.
Bottom line: 0.008 to ∼ 0.1 M ⊕ 2.7. Beyond 50 au A radial distance of roughly 50 au appears to be a real edge to the cold classical Kuiper belt (Jewitt et al. 1998;Trujillo & Brown 2001;Fuentes & Holman 2008). If there is a population with similar size and albedo's to the observed KBO at 60 au, it's mass cannot be more than 8% of the observed KBOs, as otherwise it would have been detected (Fuentes & Holman 2008). There are small bodies with semi-major axes greater than 50 au in the Solar System, but most of them are dynamically coupled to the giant planets, suggesting that they have been scattered into their large orbits. Hence, they do not represent primordial orbits. A possible exception to the objects coupled to giant planets are the Sedna type objects (Brown et al. 2004), but these objects are on highly eccentric orbits, suggesting that they were scattered to their current locations and only decoupled from the rest of the Solar System after being scattered outward, e.g. by the tidal influence of the Sun's birth cluster (Brasser et al. 2006(Brasser et al. , 2007Kaib & Quinn 2008;Brasser et al. 2012).
Bottom line: No evidence that anything formed at these distances initially

The Model
We use a new python based version of the Birnstiel et al. (2010) code called DustPy (Stammler & Birnstiel in prep.). This code allows us to compute the radial motion and growth of particles, as well as gas evolution. DustPy is a 1-d (radial) code with analytical vertical integration, solving the Smoluchowski equation (von Smoluchowski 1916) for particle growth. For more details see Birnstiel et al. (2010). In the following we are describing basics of the dust model, a simple accretion heating model. The sink term we have chosen for the gas due to photoevaporative winds is shown in Appendix A.

Basics
For simplicity we assume spherical compact particles with mass m = (4/3)πρ s a 3 , where ρ s is the material density and a the particle radius. Epstein (1924) derived a friction force under the condition that a λ g and v rel v th for spherical particles The drag force particles feel while moving through a fluid (a λ g ) is ρ g being the gas mass density and v rel the relative velocity to the gas. C D is called the drag coefficient. This drag law was already expressed -in the same form but with a constant drag coefficient -by Newton (1729) in section 2 and 7 of his second book, for the impact of air on the falling motion of hollow glass spheres, where inertia is dominant over viscous forces. The first formulation of this drag formula in the form F D = a 2 ρ g v rel 2 · f (Re) -here C D is given by some function depending on the Reynolds number with molecular viscosity ν mol -was given by Rayleigh (1892). The drag coefficient for Re ≤ 2 · 10 5 is given by (Cheng 2009) The molecular viscosity for hard spheres, neither attracting nor repulsing, is roughly given by (Chapman 1916, see his Eq. (249)), where the gas mean free path is and we further assume that the geometrical cross section σ g is that of molecular hydrogen σ H2 = 2 · 10 −15 cm 2 .
Massey & Mohr (1933) pointed out that this classical approximation is good enough for helium and hydrogen over a large range of temperatures (see their table III on p. 450), i. e., quantum mechanics is not required. For cold temperatures (∼ 10 K) quantum mechanics is important, but these temperatures are typical for the outer disk, where the gas density is so low that particles are in the Epstein drag regime anyway.
We define the stopping time as following Whipple (1972). According to Newton's second law it iṡ i. e., τ s is the time it takes for the velocity of the particle relative to the gas to be reduced from v rel to v rel /e. How well particles are coupled to the gas is described by their Stokes number, which we define by where Ω = GM r 3 (10) is the Keplerian frequency. Since the Stokes number is the ratio of the stopping time, over which particles couple to the gas, and the dynamical gas timescale, small values (St 1) mean that particles are coupled to the gas. Large values (St 1) indicate that particles are decoupled from the gas. I.e., particles are coupled to the gas motion in less than an orbit for St 1, and large Stokes numbers (St 1) would need many orbits to synchronize to the gas motion. If the mean free path of gas molecules λ g is large enough, particles are in the Epstein drag regime. If λ g is small compared to the particle radius a, they are in the fluid regime. The transition between the two regimes occurs around 1 λ g = 4a/9. The Stokes number is thus If the fluid regime is reached, we follow Birnstiel et al. (2010) and assume the Stokes drag law regime, i. e., 1 Following Weidenschilling (1977a), this condition can be obtained by setting either the stopping time in the Stokes drag law and the Epstein drag regime or the two drag forces equal and making use of Eq. (5).
for Re ≤ 1 (Stokes 1851, his Eq. (126)). If λ g < 4a/9, this leads to This way the velocity of particles, their relative velocity to the gas, and the Stokes number don't have to be solved iteratively together.

Column Densities
We define the column density as the mass per 3-dim volume, density ρ, vertically integrated over height z of the disk, where i = {d, g, p} can be dust (d), gas (g), and planetesimals (p). We use Σ d as the column density, including all solid particles without planetesimals. If it has St as an argument, it is the column density of particles with this Stokes number. Following Birnstiel et al. (2010), we define the dust column density distribution per logarithmic bin of grain radius a as where n a is the number density per grain size bin. This way, knowledge of the used size grid is not needed and the total dust column density is given by As initial condition for the gas we use the self-similar profile (Lynden-Bell & Pringle 1974) where Z 0 = Σ d /Σ g is the initial solid-to-gas ratio in terms of column densities. The initial dust column density is then given by Z 0 Σ g (t 0 ).

Drift velocities
Particles, from tiny dust grains up to boulders, are embedded in the gas disk. With the force of gravity from the central star balanced by the centrifugal force, particles move on Keplerian Orbits. The action of the gas pressure gradient on the particles can be neglected because the internal density of the particles is so much larger than the gas density. The gas does feel gravity, centrifugal force, and the pressure gradient force. If these forces balance each other, the gas moves on slightly sub-Keplerian orbits. Particles with St 1 are coupled to the gas, thus they feel a centrifugal deficiency due to sub-Keplerian gas motion and drift radially inward. As long as St < 1 this leads to a stronger radial drift for increasing St. If the particle Stokes number is larger than unity, they decouple from the gas and feel a headwind from the surrounding gas. The mass-to-surface ratio increases with size and this effect becomes weaker for increasing St (i. e. increasing stopping time τ s ). The steady-state solution for radial drift reads (Nakagawa et al. 1986) for low dust-to-gas ratios (Weidenschilling 1977a). We use the latter expression for this paper to save computation time.

Planetesimal Formation Rate
For the planetesimal formation rate we follow Lenz et al. (2019). The model is based on the idea that pebble traps appear and disappear on a given timescale. In those pebble traps, pebble clouds can then collapse to planetesimals. In this model, the pebble flux (in mass per time) is transformed into planetesimals over a conversion length : We assume that this conversion length is proportional to the gas pressure scale height h g . Mass conversion from pebbles to planetesimals according to this recipe is only allowed if the condition is fulfilled, where m p is the mass of a single 100 km planetesimal and τ l is the lifetime of traps. In this paper, we assume that τ l = 100 t orb for all simulations. ε is the efficiency with which pebbles are transformed into planetesimals. For more details we refer to Lenz et al. (2019). With help of the mean radial trap separation d, one can relate this parameter to the conversion length, = d/ε. The jump from pebble-size to objects 100 km in diameter is a direct result of the particle diffusion timescale within the particle cloud and the collapse timescale (Klahr & Schreiber 2015;Schreiber 2018;Gerbig et al. 2020).

Comparison to other Planetesimal Formation Rate Models
The model for the planetesimal formation rate from Lenz et al. (2019) differs from other models. In Lenz et al. (2019) planetesimal formation is regulated by a conversion length scale, over which drifting particles are converted into planetesimals. The conversion length scale depends on the radial density of pebble traps and the efficiency of concentrating particles and converting pebble clouds into bound objects. Drążkowska et al. (2016) and Schoonenberg et al. (2018) suggest models for which planetesimal formation occurs with a certain efficiency per orbit from the local particle density. These models do assume that particles are not trapped while drifting. Potentially, an equivalent situation could be reached for explicit traps that build up and vanish on a given timescale everywhere in the disk with some average radial distance between each other. Lenz et al. (2019) parameterized this via the conversion length , see Eq.
(3) in their paper. Adding a gas gap to the simulation, Stammler et al. (2019) used the model of Schoonenberg et al. (2018) to produce planetesimals just outside this gap and were able to reproduce the observed optical depth of HD 163296. Eriksson et al. (2020) used the criterion from Yang et al. (2017) and assumed that all the available local mass is transformed into planetesimals once the condition is met for which particles in the midplane can concentrate to a particle-to-gas mass ratios of more than 10.
For further discussion of other planetesimal formation models see e.g. section 5.2 of Lenz et al. (2019).

Advection-Diffusion Equation
The particle diffusion coefficient D d for species i can be estimated with help of the gas diffusion coefficient as (Youdin & Lithwick 2007) This means small particles are diffused with the gas and larger particles are less influenced by gas diffusion. As first described by Fick (1855, p. 66), reviewed in more modern notation by Tyrrell (1964, his Eq. (1)), and derived from fundamental principles by Reeks (1983, his Eq. (25)) 2 the diffusive flux is given by (see also Cuzzi et al. (1993)) which gives the z-integrated version in radial direction In the last step we used the fact that, due to the Gaussian shape of ρ i d in z-direction, the highest contribution of the integral comes from the region within , within which the gas density does not change by much (especially because h i d < h g ). If the gas density is roughly constant, the particles Stokes number also stays roughly constant. If at z = h i d the temperature is similar to the one of the midplane, D i d can be considered to be zindependent. As long as these conditions are met, the right hand side of Eq. (26) gives a good approximation. Since particles exhibit diffusive mixing due to turbulent gas motion, they are not Table 1: Parameters that are checked in this study. Standard values are marked in gray. Those for the second fiducial set are shown in bold. Disk mass M disk , characteristic radius r c , viscosity power-law index γ, which is also the power-law index of the column density of our initial condition for r r c , breakup speed of grains v f , initial dust-to-gas ratio Z 0 , trap formation time τ f , turbulent viscosity parameter α t , planetesimal formation efficiency ε, and X-ray luminosity L X . For comparison, the disk mass of the MMSN is roughly 0.013M . In a separate row, we highlight our most appealing case, which includes a simple model for accretion heating that is not used in any other simulation.
We would like to point out that in the expressions in e. g. Desch et al. (2017), which are based on Morfill & Völk (1984), the diffusive flux is not given by This expression is strictly speaking only valid for small particles, -which couple to the gas motion on timescales shorter than the correlation time of the fastest turbulent eddy τ Kolmogorov , i.e. the smallest eddy at the dissipation scale of turbulence (Kolmogorov scale) τ s < τ Kolmogorov -or for constant gas densities. Otherwise, Eq. (26) should be used. Unfortunately, we do not know the value of τ Kolmogorov . For further details on the different turbulence regimes see e.g. Ormel & Cuzzi (2007). The difference between the two diffusion terms can be significant if the gas density drops quickly, as is the case for gap opening due to photoevaporation. Despite Dubrulle et al. (1995) using this diffusive flux for the vertical direction (z derivative instead of r For v f 10 m s −1 γ ∼ 0.5 − 1. For v f 1 m s −1 and M disk 0.1M , and gas density power-law index γ ∼ 0.5 could work but γ ∼ 1 seems more likely v f frag. speed 1 m s −1 to allow pebbles with St 10 −2 to form Z 0 initial dust-to-gas ratio 0.01 Z 0 0.03 works more or less equally well, whereas Z 0 0.003 fails τ f trap formation time Traps needed at least 300 t orb to form outside of 50 au or never formed there α t turbulence parameter α t ∼ 10 −5 − 10 −3 (or only up to a few 10 −4 if v f ∼ 1 m s −1 ) ε planetesimal formation efficiency 0.002 < εh g /d 0.06 (if ε and d/h g are constant) L X X-ray luminosity For r c 20 au, photoevaporation does not affect the final planetesimal profile significantly derivative), their result for the particle scale height is still valid since the gas density does not change much in the vertical direction within one particle scale height. By making use of Eq. (26), the advection-diffusion equation reads Here, θ(·) is the Heaviside function anḋ is the critical pebble flux to allow planetesimal formation ). We introduce a smoothing function for the Stokes number dependency of the efficiency parameter ε This pre factor is displayed in Fig. 3. The idea is to smooth out the strong dependence on the fragmentation speed-which is similar to the idea presented in Windmark et al. (2012), where particles have a velocity distribution. The evolution of the gas is given by (Pringle 1981) whereΣ w is a loss term due to photoevaporative winds. The photoevaporation model is based on Picogna et al. (2019) and described in Appendix A. For the viscosity ν we choose the turbulent viscosity according to Shakura & Sunyaev (1973) which is the same expression as Eq. (23).

Temperature Model
In order to calculate the midplane gas temperature, one needs the contribution from radiation (internal and external) as well as from accretion heating. From pure radiation heating (e.g., Armitage 2010, his section 2.4.2) one obtains where θ ≈ tan θ ≈ h g /r ≈ 0.04 (e.g. Chiang & Goldreich 1997;Pfeil & Klahr 2019). We set the background temperature due to external sources to T bg = 10 K. Gough (1981) gives a luminosity evolution of the sun, As the age of the sun is roughly t ≈ 4.6 · 10 9 yr and our simulations run for a few 10 6 yr, we can make the approximation L ≈ 5L /7 ≈ 2.73 · 10 33 erg s −1 .
For pure accretion heating (i.e. ignoring radiation heating for the moment), and without taking optical depth effects into account and assuming that T acc ∝ c 2 s , the local midplane temperature can be calculated as (Nakamoto & Nakagawa 1994;Pringle 1981) Following Ostriker (1963) and Armitage (2010, his Eq. (3.37)), we approximate the midplane temperature due to accretion and radiation heating as The Rosseland optical depth τ is approximated by where κ R is the size and wavelength averaged Rosseland opacity (Birnstiel et al. 2018) that is calculated in every time step based on the local size distribution, and Σ d is the column density of all particles except planetesimals. We use this accretion heating model only for some further test cases. For the majority of presented simulations in the main text, we stick to radiation heating only, see Eq. (32).

Analyzed Parameters
For the total disk masses we used values between the MMSN (0.013 M ) and roughly the critical value at which disk fragmentation can occur, ∼ 0.1 M (Toomre 1964;Goldreich & Lynden-Bell 1965). However, for collapse due to the disks own gravity, the cooling time is also an important criterion (Baehr et al. 2017).
For the viscosity power law index γ we also allowed extreme cases, i.e. 0.5 ≤ γ ≤ 1.5, and made the turbulence parameter disk radius dependent for the cases γ 1: where g = γ + q − 3/2 and T ∝ r −q . For the fragmentation speed, recent work by Musiolik & Wurm (2019) indicated that the value should be around 1 m s −1 . We still analyze values up to the former default of 10 m s −1 .
Values for the solar metallicity span from Z = 1.34% (Asplund et al. 2009) to Z = 2% (for a review see Vagnozzi 2019). We will use Z = 0.0134 as our fiducial initial dust-to-gas ratio.
For the trap formation time τ f , we took typical timescales for the significant evolution of disk instabilities -such as the convective overstability, vertical convective instability, subcritical baroclinic instability, or vertical shear instability (Pfeil & Klahr 2019), as well as Hall MHD (Bai & Stone 2014;Béthune et al. 2016). The values can span from a few hundred to thousands of orbits, depending on how fast the instability evolves and how fast it can then create pressure bumps. Therefore, we should also consider the viscous timescale (e.g. Armitage 2010) on which structures could form. For example, for α t ∼ 10 −4 this would give ∼ 1600t orb . We look at turbulence levels that represent an almost laminar case (α t = 10 −5 ) up to a very turbulent state (α t = 10 −2 ).
For the planetesimal formation efficiency, here defined as ε = 5h g / , we rely on numerical experiments in order to judge whether values are high or low. We found that ε = 0.3 is already high, with almost all the mass that was originially in dust being in planetesimals by the end of the simulations, see Fig. 4. If this ratio is around 0.1, we consider the efficiency to be rather low, which is the case for ε ≈ 0.01. X-ray luminosities are found in the range (Güdel et al. 1997;Vidotto et al. 2014) 10 28 L X /(erg/s) 10 31 .

Effect of Planetesimal Formation Efficiency
Before analyzing the results of all nine different disk parameters, we would first like to concentrate on the planetesimal formation efficiency. Since many processes that we do not fully understand yet are hidden in this parameter, we first have to clarify which values are low or high. In the left panel of Fig. 4 we show, for the first fiducial parameter set, the final mass in planetesimals within a given disk region over the total initial dust mass as, a function of planetesimal formation efficiency. All simulations were stopped at 10 Myr, or when essentially all the gas was drained. One can identify a linear regime for small ε, which makes sense since the formation rate scales linearly with this parameter (Eq. (21)). For very high efficiency values, the mass in all the regions should reach a plateau, as the planetesimal profile should be very close to the initial condition of the dust. For an extreme case, the conversion length is infinitaly small and pebbles are all instantly transformed into planetesimals. Once a critical large value of ε/d is reached, the initial structure is basically reproduced, which leads to a plateau in this plot. Given these extreme cases, we expect there to be a sweet spot, i.e. for a given efficiency the final planetesimal mass reaches a local maximum. This behavior can be described by looking at timescales. How fast planetesimals can be built locally from pebbles is determined by the planetesimal formation timescale Article number, page 9 of 21 A&A proofs: manuscript no. aanda If τ ptes < τ dr , pebbles are transformed into planetesimals faster than particles are removed from their location by radial drift, i.e. the planetesimal profile becomes closer to the initial dust profile. Setting those timescales equal and making use of = ε/d, one obtains which means that, for a fixed location in the disk, only the gas temperature matters, since d ∝ h g ∝ √ T . Starting from very low ε, the mass increases linearly until the sweet spot where planetesimal formation and drift occur on similar timescales. If ε is increased even more, pebbles are transformed into planetesimals before they can significantly drift, leading to a profile that is closer to the initial dust profile the higher ε is.
However, considering the simplicity of the estimate for the maximum, the prediction works surprisingly well for both parameter sets (compare the left and right panels of Figure 4). The difference in the shape of the curves may be related to whether the disk is mostly limited by drift or fragmentation, as well as the mass budget in pebbles. For the fragmentation limited case, material stays within a certain region in the disk longer, as fragmentation events force it to start growing again from tiny, very slowly drifting dust grains, or small dust being swept up by larger grains. As can be seen in Figure 6, the high fragmentation speed in the first parameter set allows the disk to be mostly fragmentation limited in the inner disk and drift limited in the outer disk, over the typical time span of planetesimal formation. However, for the second parameter set the fragmentation speed is so low that during the time of planetesimal formation basically the entire disk is fragmentation limited. Also compare to the findings of Birnstiel et al. (2012), which did not include planetesimal formation.
The vertical purple lines in Figure 4 mark the point beyond which all the initial dust mass has ended up in planetesimals, since there is not much planetesimal mass outside of 50 au and since this is a prediction for the maximum in mass at 50 au. For the radial positions further in, the required planetesimal formation efficiency to reach a maximum in mass is even smaller.

Deeper analysis of special cases
In this section we will focus on the two fiducial runs and the most appealing simulation, where only the latter contains a simple model for accretion heating, that we described in section 3.7. Fig. 5 shows the time evolution of the gas, dust, and planetesimal profiles (top panels) as well as the local values of these at three different disk radii (bottom panels). The kinks in the planetesimal profiles (solid lines) and dust profiles (dotted lines) inidcate the position of the water ice line in that simulation. Interior to the water ice line, loss of water ice due to sublimation is assumed.  the first fiducial run we used a turbulence strength parameter of α t = 10 −3 and a fragmentation speed of v f = 10 m s −1 , whereas the second fiducial run we have set α t = 10 −4 and v f = 1 m s −1 . Despite the fact that for the latter case both values are one order or magnitude smaller than in the former, the smaller v f in the second fiducial run leads to much smaller maximum particle sizes, since the fragmentation barrier scales quadratically with v f but only inversely linearly with α t . As a result, in the second parameter set almost the entire disk is limited by fragmentation over the major time of planetesimal formation, while the first becomes drift limited much faster (see Figure 6). In Fig. 5 one can also see that higher α t lets the disk spread faster, with material being removed in the very outer regions due to the constant external FUV sink term that we used.
Article number, page 11 of 21 A&A proofs: manuscript no. aanda  Table 1). The fiducial values for this set of simulations are shown in the header of the plot. In each panel, the simulation with those parameters is shown as dashed gray line, and the solid lines with colors show simulation results where only one parameter of the set was changed. The gray area shows values more massive than the solids of the MMSN (Weidenschilling 1977b;Hayashi 1981). Figure 6 shows the particle flux for particles of different sizes in the Stokes number space as a function of disk radius at three different snapshots. The fragmentation (red lines) and drift limit (purple lines) set the maximum size of the flux dominating particle species. The horizontal gray line marks the Stokes number beyond which particles are assumed to contribute to particle trapping and planetesimal formation, see also Fig. 3. Within ∼ 3 au, some particles have higher Stokes numbers, forming a kink feature because they enter the Stokes drag regime (see Figure 6 and B.3). However, it looks much less extreme in the grain size space.
More details on the special case including accretion heating can be found in Appendix B, where we show the time evolution of the pebble flux and of the planetesimal, dust, and gas profiles.

Mass Evolution
It might be valuable to know if and when the total mass in planetesimals saturates. This saturated mass can be compared with the MMSN solid mass. If this mass is not reached, we consider the parameter set of that simulation to be unable to reproduce the Solar System, as usually several times the MMSN is needed in order to get results that are comparable to the Solar System. Figures 7 and 8 show the time evolution in mass of the first and second sample, respectively. In each panel only one parameter was changed compared to the ficucial parameter set shown in the title of both figures. The final value of the total disk mass that is in planetesimals can be compared to the minimum mass for the solar nebula based on Weidenschilling (1977b) and Hayashi (1981). The required mass for initial planetesimals is marked by gray regions. If the final mass is below the gray region for a given parameter set, this set can be excluded for the solar nebula. In both figures one can see that initial dust-to-gas ratios of 0.003 and lower are not able to lead to the MMSN mass in planetesimals. A fragmentation speed of v f = 10 2 cm s −1 seems to be a critical value, under which the MMSN mass cannot be reached unless the initial dust-to-gas ratio Z 0 > 0.0134, the turbulence strength α t ≤ 10 −4 (here this parameter is also used for vertical  Table 1). The fiducial run produces just enough planetesimal mass in the disk to reach the MMSN mass in solids, even though the disk mass is at the high end already. From this point of view, one needs either a higher initial dust-to-gas ratio, a higher planetesimal formation efficiency, or a larger fragmentation speed. For the bottom left panel, no line is visible for α t = 0.01 as no planetesimals are formed in that case. and radial particle diffusion, as well as relative turbulent velocities), or ε > 0.03. For a disk starting with M disk < 0.1 M , ε or Z 0 might need even higher values to reach a final planetesimal disk mass more massive than the MMSN.
The yellow line in the middle lower panel of Fig. 8 shows a case of low planetesimal formation efficiencies in a mostly fragmentation limited disk. In this case, after about 10 5 yr photoevaporation allows another phase of planetesimal formation after the planetesimal mass in the disk reached a plateau. The same effect shows up for high initial dust-to-gas ratios, see the case of Z = 0.03 in the centered panel of Fig. 8. In both cases the reason for the second planetesimal formation phase is the higher mass budget. For low ε, particles survive longer in a fragmentation limited disk -since their average radial drift velocity is much lower, due to disruptive collisions that replenish slowly drifting dust grains -and less mass is transformed into planetesimals. For Z 0.03, the initial particle mass budget is already so high that there remains enough mass for planetesimal formation at later times, when photoevaoration has removed a significant amount of gas mass. However, this effect of a second planetesmial formation phase only occurs if the disk is mostly fragmentation limited, which applies for the results shown in Fig. 8 but not for those shown in Fig. 7, in which case the disks are mostly drift limited. Additionally, the second planetesimal formation phase induced by photoevaporation is not sufficient to reach the mass of the MMSN for low planetesimal formation efficiencies, i.e. for ε 0.01. L X /erg s −1 = 0,3 · 10 28 , 10 29 ,3 · 10 29 , 10 31 Fig. 9: Final (at 10 Myr) planetesimal column density as a function of disk radius for the first sample (gray values in Table 1). The fiducial values for this set of simulations are shown in the header of the plot. In each panel, the simulation with those parameters is shown as dashed gray line, and the solid lines with colors show simulation results where only one parameter of the set was changed. The gray areas represent the constraints that we described in sec. 2, where the mass in the given region was translated into a column density, assuming a planetesimal profile ∝ r −2.25 (see Eq. (37) of Lenz et al. (2019)). For the outermost region, we also overplotted a box ∝ r −8 .
Pinilla et al. 2012), even if these traps appear and disappear on a given timescale. This could lead to longer planetesimal formation and greater impact of photoevapotation. However, this might not change the results significantly, leaving the presented conclusions untouched.

Deep Parameter Analysis
By looking at the final planetesimal profiles for all nine parameters, we find a huge variety of possible parameters for the Solar Nebula. It is reassuring that the model works not only for a very finely tuned subset of parameter choices. Though the different parameters can influence each other, it is still possible to draw some conclusions. Table 2 shows which parameters fail to fulfill the outer Solar System constraints or the MMSN mass. In Table 3 we present disk parameter ranges that could potentially reproduce the Solar System. Those conclusions are based on Ta-ble 2. How much mass the initial disk should contain depends on the fragmentation speed, since the latter determines how much mass is in particles with St 0.01.
Our parameter analysis is based on Figs. 7, 8, 9, and 10. The last two of these figures shows the column density profiles of the final planetesimal population. In each panel only one parameter is varied compared to the fiducial parameter set (dashed lines). In the background, the gray boxes represent the mass constraints discussed in Sec. 2.
The initial characteristic radius r c has two major effects. One is the radial position beyond which the dust and gas density drops exponentially. The second is that for smaller (larger) r c there is more mass in the inner (outer) disk region. If the disk is too large, there is simply too much mass available around 40 au and beyond to form planetesimals. Too much mass in the outer disk then leads to violation of the upper CCKBO constraint. Whether this constraint is indeed violated depends also on the  Table 1). If ε is high enough, photoevaporation does not change the result significantly. But it seems to have an effect for low ε, see yellow line in the bottom middle panel of Fig. 8. In the top left panel we also plot the most appealing case as dotted line (last row in Tab. 1). initial disk mass, the planetesimal formation efficiency, the initial dust-to-gas ratio, and the viscosity power-law index γ. However, for a narrow set of these four parameters, finding constraints for r c is possible.
The fragmentation speed v f changes the outcome a lot since the fragmentation limit depends quadratically on this parameter. For v f 10 m s −1 γ ∼ 0.5 is actually beneficial due to the stronger density drop in the outer disk, see orange line in the top right panel of Fig. 9. However, for v f 1 m s −1 one would need more than M disk 0.1M to create enough mass to build all planets, specifically in the inner disk part ( 15 au), see orange line in the top right panel of Fig. 10.
If v f < 1 m s −1 , too few or no particles with St 0.01 would be formed, which is the necessary Stokes number to make trapping, and the collapse of pebble clouds into planetesimals work. Already, with v f = 1 m s −1 it is difficult to meet all the constraints, especially a total mass in planetesimals larger than the minimum mass Solar Nebula (Weidenschilling 1977b;Hayashi 1981) -see the lowest line in the middle left panel of Fig. 7 and the dashed lines in Fig. 8.
The initial dust-to-gas ratio Z 0 determines how much mass is initially in particles, and also the dust dynamics, as a low dust-togas ratio leads to a drift limited disk that loses particles quickly due to drift. Large Z 0 of up to roughly 0.03 seem to allow fulfilling the mass constraints on initial planetesimals. However, values of 0.003 lead to too little mass in the final planetesimal population as can be seen in Figs. 7,8,9,and 10.
The more time particles have to drift from the region outside of ∼ 50 au to the inner parts before planetesimal formation, the better the CCKB constraints can be met. Alternatively, traps might not occur outside of 50 au at all (Pfeil & Klahr 2019). From the simulations we conclude that the trap formation time must be τ f > 300 t orb or even 1000 t orb . Other values lead to masses between ∼ 30 au and 50 au, which are orders of magnitudes higher than the upper limit.
Constraining values for the turbulence parameter α t is also linked to the fragmentation speed, because the fragmentation limit scales inversely linear with α t but quadratically with v f . For this limit, only the relative velocity matters. However, we assumed that vertical and radial diffusion, as well as the viscos-ity parameter for the gas, have the same value α t that we used for the turbulent velocities. For smaller α t , particles can settle closer to the midplane. If this value is low enough, growth is not limited by relative turbulent velocties but by relative settling speeds or relative radial drift. When α t is high (∼ 10 −2 ), relative turbulent velocities are too high to allow St > 0.01 particles. At the same time, the radial viscous gas motion drags dust along to the outer regions of the disk, leading to too much mass in planetesimals outside of 30 au. Hence, for v f ∼ 10 m s −1 we suggest α t ∼ 10 −5 − 10 −3 , while the upper end should be up to a few 10 −4 if v f ∼ 1 m s −1 .
Values of ε/d ≤ 0.002/h g can be excluded for the solar nebula (if ε and d/h g are constant), because the total final planetesimal mass in the disk is below the MMSN, and the mass required for the Nice disk can not be reached. Unless the solar nebula was not very small (r c 10 au), planetesimal formation should not have been too efficient, i.e., ε/d 0.06/h g . Otherwise too many planetesimals are formed outside of 30 au. Thus, the range of possible values is 0.002 < εh g /d 0.06.
Photoevaporation did not influence the final planetesimal profile significantly for small enough disks (r c 20 au). However, for large disks (r c ∼ 100 au) it can make a difference. But this case is not interesting for finding similar conditions to the solar nebula, as in this case too much mass ends up in planetesimals in the outer disk regions anyway.
In the left panel of Fig. 4, the plateau is reached for smaller values of ε/d than in the right panel, which is linked to the smaller pebble mass that is available for planetesimal formation, since the smaller fragmentation velocity leads to smaller maximum Stokes numbers.

With Accretion Heating
Accretion heating leads to a hotter inner disk, which is good since not a single planetesimal was detected inside the orbit of Mercury. I.e., the planetesimal profile has to drop drastically in the inner region before reaching Mercury's current radial position. This can be satisfied due to the higher gas temperatures, as these are forcing the fragmentation limit to be at lower Stokes numbers, and to a larger conversion length ( ∝ h g ∝ √ T ). In this simulation the ice line is moving radially over time, which is why there is no distinct kink feature in the profile. The constraints in the outer disk, that is the Nice disk and the CCKB, are the strongest ones we have. The other constraints may be a bit more flexible. These constraints are roughly met by this simulation.
Since the temperature model presented in this paper was not tested in a comprehensive way, e.g. by comparing with Hubeny (1990) or Nakamoto & Nakagawa (1994), we only use it to show one special case. In addition, planetesimal-planetesimal collisions would replenish the small dust population (Gerbig et al. 2019). This effect is not taken into account in this paper but could change the gas midplane temperature via the mean dust opacity.

Summary
We used an extended version of the Lenz et al. (2019) model, including Stokes drag for particles, and allowing the gas to evolve viscously while photoevaporation is removing gas over time. The analyzed parameter space was largely increased. While this paper provides a parameter study for pebble flux-regulated planetesimal formation, we focused on meeting Solar System constraints for initial planetesimals. Therefore, we used two different default parameter sets and varied one out of nine parameters per simulation. Overall, while some parameters can be excluded, the model seems to be very robust, thus it does not require parameter fine tuning in order to fulfill the constraints.
The computation times of the presented simulations were between roughly a week and six months, while running on ten cores each. Using more than ten cores for a simulation would not decrease the computation time significantly since the code cannot make use of further parallelization. The runs of the second sample, in particular, were running for months. To shorten the computation time, a simple model must be used such as the two population model presented by Birnstiel et al. (2012). However, this simplified model was only tested for a narrow set of parameters and causes deviations from DustPy simulations for certain parameters. Additionally, the two population model was not yet tested in detail with the inclusion of the planetesimal formation model that was used in this study. A simple model reproducing the results shown in this study is likely possible, however, we preferred to use DustPy in order to rely on fundamental physics principles and a sophisticated growth and fragmentation model rather than simplified and untested models.
In sec. 2 we suggested mass constraints for initial planetesimals in different regions of the disk: The fragmentation speed v f that leads to breakup in particle collisions and the turbulence parameter of relative velocities α t determine how large particles can grow. If the combination of both leads to a fragmentation limit that is close to Stokes numbers of 0.01, the available mass for planetesimal formation will be affected by these parameters. This is why a constraint in total initial disk mass has to be linked to (mostly) v f . To fulfill the constraints we suggested, we need M disk 0.1M for v f ∼ 1 m s −1 and M disk 0.02M for v f 10 m s −1 . In addition, the solar nebular was not larger than r c 50 au (r c is the initial transition radius between a power-law and a dropping exponential profile). The power-law index of that inner region was likely around γ ∼ 1, but for large fragmentation speeds γ ∼ 0.5 can be beneficial for the outer region due to the density drop (if traps can be formed outside of 50 au). To allow pebbles with St 10 −2 to form, which is roughly the needed Stokes number for trapping and subsequent planetesimal formation, one needs v f 1 m s −1 . For the initial dust-to-gas ratio, many values 0.01 Z 0 0.03 could work, but Z 0 0.003 leads to too little mass in planetesimals. Outside of 50 au traps needed at least 300 t orb , or never formed there. For the turbulence parameter we find a wide range of possible values α t ∼ 10 −5 − 10 −3 (or only up to a few 10 −4 if v f ∼ 1 m s −1 ). Since disk parameters can affect each other, we also find a wide range for the radial pebble to planetesimal conversion length: 0.002 < h g / 0.06. If the disk is sufficiently small (r c 20 au), photoevaporation does not change the final planetesimal profile by much.
We estimated the maximum mass in planetesimals by equating the planetesimal formation and drift timescale. This approach leads to = r, see Eq. (42), which seems to fit our simulation results well (see Figure 4). If the planetesimal formation timescale is much shorter than the drift timescale, the planetesimal profile reproduces the initial dust profile. Planetesimal formation efficiencies smaller than the value corresponding to this sweet-spot lead to planetesimal profiles steeper than the initial dust profile or even steeper than the minimum mass Solar Nebula profile. This effect was already observed in Lenz et al. (2019), and this study provides an estimate for the transition to more local planetesimal formation, which is linked to slopes closer to the initial dust profile.
Within the model, further limitations are that no pebble accretion was included, which could especially affect the planetesimal profile in the inner disk as pebbles get accreted before reaching that zone. In addition, our simulations did not consider planetesimal-planetesimal collisions, which would lead to multiple generations in planetesimals, pebbles, and dust.

Conclusions
The MMSN is not consistent with viscous disk evolution models and does not provide enough mass in the giant planet forming region to allow strong gas accretion (see Fig. 2). While typically the MMSN distribution is assumed to be present from the beginning, the timing of substantial planetesimal formation could also matter for further embryo formation and evolution. We have shown that pebble flux-regulated planetesimal formation produces beneficial planetesimal distributions for a wide range of parameters, both with respect to planetesimal formation and initial conditions of the disk. Even though the impact of disk parameters on the evolution of initial planetesimals influence each other, some constraints on these parameters were found in this study. Having only a narrow set of parameters that could reproduce the Solar System would have indicated model fine tuning. This stresses the applicability of our parameterization to models of planet formation, e.g. population synthesis models.  due to external FUV radiation. Also for the case of an inner hole, we follow Picogna et al. (2019). The hole radius r h is implicitly defined via the radially integrated midplane gas number density I.e., once the critical flux for planetesimal formation is reached the flux is orders of magnitude larger than the critical valueṀ cr (above the shaded areas in both panels). Photoevaporation (an effect not included in Lenz et al. (2019)) leads to a small increase of the pebble flux at late times, see the evolution after ∼ 8 · 10 6 yr in the lower panel of Fig. B.1. However, this increase has only a negligible effect on the final planetesimal population since the pebble flux has dropped by many orders of magnitude compared to its maximum value. In Fig. B.2, since in this simulation the gas temperature depends on dust evolution, the ice line radially moves over time. Hence, the kink feature in planetesimals that is clearly visible at early times (∼ 10 4 yr) is smeared out at late times (∼ 10 6 yr). At all three locations shown in the bottom panel, planetesimal formation is going on for around ∼ 10 6 yr with significant mass contributions. Note that a higher X-ray luminosity (up to ∼ 10 30 erg s −1 ) would not change the results by much as the disk would not vanish before ∼ 2 Myr. At this time the planetesimal population has saturated already.
For our most appealing case which includes accretion heating (Fig. B.3), the higher gas midplane temperatures in the inner disk region are leading to smaller maximum Stokes numbers compared to a situation with pure radiation heating. At late times (∼ 1 Myr), enough dust was converted into planetesimals causing the opacity to drop and thus gas temperatures are much lower than in the initial phase of disk evolution. As a result, higher Stokes numbers can be reached and significantly more planetesimals are formed within ∼ 1 au. This effect is also visible in Figures B.1