Linking planetesimal and dust content in protoplanetary disks via a local toy model

If planetesimal formation is an efficient process, as suggested by several models involving gravitational collapse of pebble clouds, then, before long, a significant part of the primordial dust mass should be absorbed in many km sized objects. A good understanding of the total amount of solids in the disk around a young star is crucial for planet formation theory. But as the mass of particles above the mm size cannot be assessed observationally, one must ask how much mass is hidden in bigger objects. We perform 0-d local simulations to study how the planetesimal to dust and pebble ratio is evolving in time and to develop an understanding of the potentially existing mass in planetesimals for a certain amount of dust and pebbles at a given disk age. We perform a parameter study based on a model considering dust growth, planetesimal formation and collisional fragmentation of planetesimals, while neglecting radial transport processes. While at early times, dust is the dominant solid particle species, there is a phase during which planetesimals make up a significant portion of the total mass starting at approximately $10^4 - 10^6$ yr. The time of this phase and the maximal total planetesimal mass strongly depend on the distance to the star $R$, the initial disk mass, and the efficiency of planetesimal formation $\epsilon$. After approximately $10^6$ yr, our model predicts planetesimal collisions to dominate, which resupplies small particles. In our model, planetesimals form fast and everywhere in the disk. For a given $\epsilon$, we were able to relate the dust content and mass of a given disk to its planetesimal content, providing us with some helpful basic intuition about mass distribution of solids and its dependence on underlying physical processes.


Introduction
Planetesimals are the smallest solid objects that are bound by their own gravitational attraction in lieu of material binding forces. Since planetesimals are believed to constitute the building blocks of planets (e.g., Kokubo & Ida 2012), studying the physical processes involved in their formation and evolution within a protoplanetary disk around a young star is of great importance for planet-formation theory and thus also to the question of the origin of our Solar System's planets. Understanding the solid material available in a protoplanetary disk is a key aspect of this issue. However, observations are limited in what particle sizes can be detected. While for instance the Atacama Large Millimeter/submillimeter Array allows us to determine the amount of mm-and µm-sized particles via observations in the scattered light (see e.g., Andrews et al. 2009Andrews et al. , 2013Ansdell et al. 2016), it is unknown how much solid material is in the size of planetesimals since they do not show characteristic absorption features and they are too small to be be detected individually. Models of planetesimal formation (e.g., Drążkowska et al. 2016;Drażkowska & Alibert 2017) based on the dust evolution model presented by Birnstiel et al. (2012) show rapid transformation of solid material from dust via pebbles into planetesimals. The ideal scenario of a model combining dust growth, pebble formation, and planetesimal formation with evolution models, such as Morbidelli et al. (2009) or Kobayashi et al. (2016), would be numerically very expensive because size ranges would be significantly larger than what is currently covered in evolution models that include planetesimal interactions (e.g., Kokubo & Ida 2002;Leinhardt & Stewart 2009;Levison et al. 2012;San Sebastián et al. 2019). Therefore, before starting this endeavor, we performed 0-d local simulations to study how the planetesimal to dust and pebble ratio evolves over time. Hereby, we consider dust evolution based on the two-population model by Birnstiel et al. (2012), pebble flux-regulated planetesimal formation after Lenz et al. (2019), and collisional fragmentation of planetesimals to resupply small particles. Our aim is to develop an understanding of the potentially existing mass in planetesimals for a certain amount of dust and pebbles, by exploring two scenarios and various parameters.
In our model we ignore spatial transport of material, which is only important if inward drift of pebbles occurs on shorter timescales than their transformation into planetesimals. If the latter process is sufficiently efficient, spatial transport becomes unimportant. Further, we assume a constant gas column density by ignoring viscous gas evolution (Lüst 1952;Pringle 1981;Hueso & Guillot 2005;Birnstiel et al. 2010) or photoevaporation (see e.g., Ercolano et al. 2009;Owen et al. 2011;Nakatani et al. 2018;Picogna et al. 2019). This would be of particular importance during the late stages of the disk evolution and in the inner disk regions. We, therefore, limit the run time of our simulations to 10 7 yr to be in line with typical lifetimes of protoplanetary

Protoplanetary disk setup
A protoplanetary disk initially contains gas and dust circulating around a central star with mass M star . We use the concept of column densities, which is commonly used in astrophysics when describing accretion disks. In cylindrical coordinates (R, φ, z) it is defined as the integral of the volumetric mass density ρ of gas or solids integrated along the vertical path z that goes through the disk By assuming a cylindrical symmetric disk, one can eliminate the φ-dependence of the column density. If spatial transport of material is neglected, conservation of mass implies that also the total column density of all solid material Σ solids must be conserved ∂Σ solids (R, t) ∂t = 0.
Initially, the total column density of solids is made up purely by dust. We relate the dust column density Σ dst to the gas column density Σ g with the dust-to-gas ratio dg , which we assume to be R-independent. The total column density of all material Σ total , including solids and gas, is then given by where we introduced the initial dust-to-gas ratio as found by Savage & Jenkins (1972) or Draine et al. (2007) in the interstellar medium (ISM). The gas profile is non-evolving in our model, hence also the total column density Σ total in Eq.
(3) must be constant in time. The initial column density depends on R according to a self-similar initial profile derived by Lynden-Bell & Pringle (1974), where R C is the so-called characteristic radius, marking the transition of the power law to the exponential. C Σ is a normalization constant chosen such that the disk contains a total mass of solids and gas of M disk . For our simulations we choose R C = 40 AU, roughly corresponding to Kuiper Belt's location in the Solar System. For a vertically isothermal disk with temperature T , the sound speed c s of molecular hydrogen is given by with Boltzmann constant k B = 1.3807 × 10 −16 erg K −1 , and proton mass m p = 1.673 × 10 −24 g. Further, the vertical pressure scale height of the gas disk h g is defined as which originates from considering vertical hydrostatic balance (Weizsäcker 1948). Here, we introduced the Keplerian angular velocity Ω, which is derived from balance of centrifugal and gravitational acceleration where G ≈ 6.674 × 10 −8 cm 3 g −1 s −2 . We set the exponent of the polynomial decline in Eq. (5) to n = 1 after the radial viscosity profile of turbulent disks presented by Shakura & Sunyaev (1973): to which the column density is inversely proportional in the inner disk. We assume a temperature profile of (Chiang & Goldreich 1997) where we choose stellar properties of R star = 1.25 × R ≈ 1.25 × 7 × 10 10 cm and T star = 4000 K (Beckwith et al. 1990), q = 0.5 after Chiang & Goldreich (1997) and an irradiation angle of α irr = 0.1 (e.g., Pfeil & Klahr 2019). This profile is valid for disks where heating is dominated by radiation from the star, in lieu of viscous heating in accretion disks (Pringle 1981), which we assume to be insignificant in comparison. Lastly, the dimensionless parameter α quantifies disk turbulence. Typical values are in the range 10 −4 ≤ α ≤ 10 −2 as concluded from simulations (Johansen & Klahr 2005;Dzyurkevich et al. 2010;Nelson et al. 2013;Flock et al. 2017) or observations such as Andrews et al. (2009) andFlaherty et al. (2017). Birnstiel et al. (2012) showed that growth of dust particles in protoplanetary disks can be sufficiently described by representing the entire dust population by only two sizes. The small size (henceforth referenced to as "dust") represents those particles which are tightly coupled to the gas via aerodynamic friction. The larger size (henceforth referenced to as "pebbles") represents particles subject to a substantial radial drift. These species are characterized by their respective Stokes numbers, a crucial quantity for describing the aerodynamic behavior of particles. For St 1, they are decoupled from the gas, while for St 1, they are strongly affected by gas dynamics. Assuming monodisperse coagulation (Stepinski & Valageas 1997;Kornet et al. 2001; A116, page 2 of 17 K. Gerbig et al.: Linking planetesimal and dust content in protoplanetary disks Brauer et al. 2008;Birnstiel et al. 2012) the size a of a dust particle at a certain time t can be described by an expression of exponential growth

Dust growth
where a 0 is an initial reference size and τ growth a certain growth timescale. Equation (11) is only valid for constant dg , which we assume to be approximately true, as the resulting error is insignificant in comparison to other assumptions in model. We simplify the two-population model by assigning the dust species to a 0 = 0.1 µm after Mathis et al. (1977) and the pebble species to a = a(St pbb ).
As the size of dust grains in our model remains smaller than the mean free path of the gas, the Epstein regime of gas drag is valid and the molecular nature of the gas has to be considered (Epstein 1924). Near the midplane we can then follow Birnstiel et al. (2012) and calculate St dst via where ρ is the mass density of dust. We choose ρ = 1.2 g cm −3 , which was found by Carry (2012) for asteroids. Fragmentation (Blum & Münch 1993;Blum & Wurm 2008) and radial drift (Klahr & Bodenheimer 2006;Birnstiel et al. 2012) may limit the maximum pebble size. In the inner disk, the former is typically more significant, whereas in the outer disk, the latter typically is more important (Birnstiel et al. 2012;Lenz et al. 2019). As these limits and their evolution are not the focus of this paper, we will simplify and treat the Stokes number of pebbles as a fixed parameter of our model with St pbb = 0.1 approximately corresponding to the maximum Stokes number of the larger particle species as found by Birnstiel et al. (2012). Implications are discussed in Appendix B.
Following Birnstiel et al. (2012), we estimate the dust-topebbles growth timescale in the midplane with The approximation in Eq. (13) holds if relative dust velocities are set by disk turbulence ∆v dst ≈ c s √ 3αSt dst with St dst 1, i.e. for sufficiently large α-values.
The dust to pebble growth rate is calculated viȧ The differential equatioṅ illustrates the decrease of the local dust supply due to particle growth to pebble size. It is solved analytically by Here, Σ 0 dst is the initial dust column density from the profile in Eq. (5) and is constant in time for a fixed St pbb .

Pebble flux-regulated planetesimal formation
A gas parcel at a certain distance from the central star is in a force balance between gravitational, centrifugal and thermal pressure forces. A solid particle does not feel this outward oriented pressure force and should therefore move on a Keplerian orbit, whereas the gas moves on a sub-Keplerian orbit. Because of aerodynamic friction, solid particles lose angular momentum to the gas and spiral inward toward the central star. After Nakagawa et al. (1986), the radial drift velocity of pebbles can be derived by considering the equations of motions of particles and gas v drift = St pbb where we introduced the exponent of the gas pressure power law which is for R R C and T ∝ R −1/2 given by γ = −2.75. Due to the assumption of a constant pebble Stokes number, the drift velocity of pebbles in Eq. (18) remains almost constant as only dg changes in time.
One finds that pebbles with St pbb ≈ 10 −1 quickly drift inward toward the central star on very short timescales. Such a particle at 10 AU can be expected to reach the star in under 10 4 yr. If pebbles can not grow an order of magnitude in a shorter time than this, they fall victim to evaporation in the inner disk. Particle growth also leads to higher relative velocities, so that upon collision the bodies tend to break up rather than stick together (Blum & Münch 1993;Chokshi et al. 1993;Blum & Wurm 2008). Therefore, further coagulation from pebbles to planetesimals seems unlikely (Homma & Nakamoto 2018), though not impossible as for example shown in Kataoka et al. (2013).
Our work focuses on a different way to overcome the fragmentation and drift barrier, as first suggested by Safronov (1969) and later Goldreich & Ward (1973). Here, when particles settle toward the midplane, they undergo a gravitational instability and form planetesimals in a spontaneous event. Balbus & Hawley (1998) found that magneto-rotational instability in the disk forms short-lived turbulent eddies, vortices and pressure bumps which may trap pebble-sized particles (Johansen & Klahr 2005). Further, streaming instability (Youdin & Goodman 2005;Squire & Hopkins 2018;Umurhan et al. 2019), where the feedback of particles onto gas is crucial, may be able to create significant particle over-densities (Johansen & Youdin 2007;Carrera et al. 2015;Simon et al. 2016;Nesvorny et al. 2019). If a particle cloud is not ripped apart by the star's tidal forces, i.e. has reached Hilldensity, it can collapse and form a planetesimal as discussed in Johansen et al. (2006. The initial size of the resulting planetesimal may then be given by the balance of contraction and particle diffusion timescale (Klahr & Schreiber 2015). Instead of considering different types of traps and examining physical properties of the disk in detail, we will follow Lenz et al. (2019) and describe the trapping mechanism with the help of parameters, most important of which is the planetesimal formation efficiency . We will also refer to as the trap (or trapping) efficiency, but it is important to keep in mind, that it quantifies not only how efficient the trap can accumulate pebbles, but also how efficient trapped pebbles can be converted into planetesimals. The trap distance d quantifies typical radial separation of the trap and τ trap a trap's typical lifetime. We set d = 5h g and τ trap to 100 orbits, as A116, page 3 of 17 A&A 629, A116 (2019) found numerically by Dittrich et al. (2013) and Manger & Klahr (2018).
Following Lenz et al. (2019), pebbles are transformed into planetesimals over the conversion length We further express the planetesimal formation rate bẏ Here, we ignore the contribution of any dust that might also be trapped in the collapsing cloud. Dust is slowed down more significantly than pebbles, as the sedimentation (or contraction) timescale τ sed for St < 1 particles increases due to aerodynamic friction on the order of where τ ff is the free-fall timescale (Shariff & Cuzzi 2015;Klahr & Schreiber 2015). Hence, we expect pebbles to dominate the mass contribution of the collapsing cloud.
Of course, the recipe in Eq. (21) is only valid once the flux of pebbles through the trap has reached a critical value, in other words once the trap has accumulated enough mass to form at least a single planetesimal with mass m pls (Lenz et al. 2019): where Σ pbb,St is the pebble column density per unit Stokes number. Since in this model, we assume all pebbles to have the same constant Stokes number, the integral vanishes and the condition simplifies to If Eq. (24) is not fulfilled,Σ form is set to 0. We deduce the fixed parameter m pls from the initial planetesimal size, which can be estimated by studying Asteroid Belt and Kuiper Belt objectsrelics of the planet formation process in the Solar System. Observations of the cumulative size distribution of bodies in these regions show a decrease of slope for increasing sizes at roughly 100 km (Bottke et al. 2005;Nesvorný et al. 2011;Fraser et al. 2014;Delbo' et al. 2017). It was concluded that this kink cannot be obtained by collisional evolution alone and instead originates from the primordial size distribution. Morbidelli et al. (2009) proposed that planetesimals were born with a minimum diameter of 80 km, which is proven to be possible Cuzzi et al. 2008;Klahr & Schreiber 2015). Assuming homogeneous and spherical planetesimals with mass density ρ = 1.2 g cm −3 , the planetesimal mass is approximated with where r pls = 50 km is the planetesimal radius. We note that this recipe for planetesimal formation also works for smaller masses, as the pebble flux is typically super-critical (Lenz et al. 2019). Our model does not include a mass grid and therefore has to also exclude growth of planetesimals to larger sizes, meaning the planetesimal population will exclusively contain 100 km sized objects.

Planetesimal collision model
Collisions between two planetesimals are complex physical processes with many potential outcomes. The simplified picture used in this paper's model imagines the column density of the two colliding planetesimals to be redistributed among fragments of different sizes upon a collisional event. Collisions are thus assumed to be destructive, as the total planetesimal column density in the system decreases upon planetesimal collision. The most convenient approach is to apply an inverse power law for the function of the specific number column density n m (m), which describes the fragment distribution by assuming a collisional cascade, as seen in e.g., Zvyagina et al. (1974). The specific number column density reads where C n is a normalization constant and ξ = 1.83 as found by Dohnanyi (1969) by analytically investigating the collective dynamical interactions of asteroids for inelastic collisions and fragmentation. This value was also found experimentally via collision experiments with basalt rocks (Fujiwara et al. 1977). The slope of the power law incorporates specifics of a collisional event, such as impact velocities, which Dohnanyi (1969) assumes to be on the order of kilometers per second. We point out, that this assumption is overestimating typical planetesimal impact velocities in protoplanetary disks by about 2 orders of magnitude (Wetherill & Stewart 1993;Morbidelli et al. 2009). Although, ξ = 1.83 is robust against variations of physical parameters like the impact velocity (Dohnanyi 1969), and our collision model is robust against variations of ξ, we note that we therefore also overestimate the steepness of the power law. We derive the normalization C n in Appendix C.1.
As we are only interested in an average size distribution of fragments, we express the column density that is redistributed into objects with masses between m i and m j via the integral m j m i n m mdm. Now, the column density redistribution ratio upon a collisional event can be estimated by defining three mass intervals corresponding to the three species dust, pebbles and planetesimals, i.e. for dust m 0 ≤ m < m 1 , for pebbles m 1 ≤ m < m 2 and for planetesimals m 2 ≤ m ≤ M. The derivations of m 1 and m 2 are shown in Appendix C.2.
With these requisites, one can solve the integral m j m i n m mdm for each population, which leads to three fractions p dst , p pbb and p pls corresponding to the column density participating in a single collision that is redistributed to dust, pebbles and planetesimals, respectively. For typical disk parameters (see Table 1), m 1 , m 2 m pls . With ξ = 1.83, the mass in fragments is therefore significantly dominated by objects which we still consider to be planetesimals.
The effect of planetesimal collisions on the surface density evolution is affected not only by the outcome of a single collision, but also on the frequency of planetesimal encounters in a protoplanetary disk, quantified by planetesimal collision timescale and rate. We consider the mean free path λ mfp between collisions. If a single particle i is moving through a cluster of other particles j with number density n j and cross section for collisions between particles i and j is σ i j , then (Birnstiel et al. 2016)   (2018) The collision timescale τ col,i is defined as the average time during which particle i experiences one encounter where ∆v i j is the relative velocity of the particles i and j to each other. Since all planetesimals in our model are equal-sized, the cross section two planetesimals i and j with radius r pls , is given by the gravitational cross section σ grav (Safronov 1969), which considers both the geometry of the system but also gravitational focusing: where the escape velocity v esc of the planetesimals is v esc = 2Gm pls r pls .
By approximating the vertical planetesimal distribution with a Gaussian with a root mean square width of h pls = v hill Ω −1 (Goldreich et al. 2004), the number density of planetesimals can be written as Here, we introduced the Hill velocity v hill , which is the velocity of particle i relative to particle j when they are just close enough to each other that the total gravitational force acting on them is not dominated by the mass of the central star M star , but instead their own masses m i and m j , i.e. when entering their respective Hill spheres. For equally massive planetesimals the hill velocity reads (Hill 1878) v hill = RΩ 2m pls 3M star Plugging the above expressions into Eq. (28) we get for the collision timescale of equal-sized planetesimals We follow Morbidelli et al. (2009) and insert the Hill velocity for the relative velocity of planetesimals ∆v i j . We note that the relative velocity appears once directly in Eq. (28) and again as its inverse squared in Eq. (29), leading to the perhaps counter intuitive implication of higher relative velocities causing longer collision timescales. The rate of planetesimal collisionsΣ col must also depend on the available planetesimal surface density. It is given bẏ

Balance equations for two scenarios
The evolution of the column densities due to the dust growth to pebbles, planetesimal formation and planetesimal collisions can be assembled to two systems of coupled differential equations corresponding to two different scenarios. The rates for the model processes are represented by sink and source terms: SinceΣ growth andΣ form each appear once as a sink and a source term and p dst + p pbb + p pls = 1 per definition, the total column density in the system is conserved in time In this configuration, we expect a column density steady state to occur. A116, page 5 of 17 A&A 629, A116 (2019) Another view arises when one considers that planetesimal collisions may produce very compact fragments, in lieu of fluffy aggregates. In particular, if planetesimal formation compactifies the material involved in the gravitational collapse, then collisional dust would be much more dense than primordial dust aggregates. Material properties of dust grains are important for the microphysics of growth (Ormel et al. 2007;Paszun & Dominik 2009). For very compact dust particles, growth via sticking may be more difficult than for fluffy aggregates, because they may not be able to absorb all of the collisional energy and consequently break apart or restructure rather than grow. This is especially the case if they are not encased by an ice mantle which could be destroyed in collisions, provided there is no recondensation. Further, even if two compact grains collide gently enough, they may be very loosely packed such that future collisions are more likely to be destructive.
Therefore, we differentiate dust and pebbles further into the subspecies "primordial" and "collisional" dust and pebbles respectively. We note that primordial pebbles do not exist at the beginning of the simulation. Rather, they grow directly from primordial dust grains. This is in contrast to collisional pebbles, which originate from material that was previously comprised in planetesimals. We introduce two additional differential equationṡ where f pbb,prim = Σ pbb,prim Σ pbb,prim + Σ pbb,col and (37f) correspond to the fraction of pebbles that is considered primordial and collisional, respectively. This configuration is visualized per flowchart in Fig. 1. In contrast to Eqs. (35a)-(35c), collisional dust is too compact and dense to grow further, whereas collisional pebbles may participate in planetesimal formation again.
We denote the configuration depicted in Eqs. (35a)-(35c), where no differentiation between primordial and collisional dust and pebbles is performed, as scenario 1. Scenario 2 will describe the configuration in Fig. 1. Contrary to scenario 1, scenario 2 will not reach an equilibrium state as Eq. (37c) only contains a loss term, while Eq. (37e) only has a source term. These two scenarios represent limiting cases. Intermediate scenarios where a fraction of collisional dust can grow, for example by coagulating with primordial dust, are also conceivable but not subject of this work.

Initial conditions and default parameters
Due to locality of the model, it is very easy to implement and computationally inexpensive. The code written specifically for this local model is based on the 2nd order Runge Kutta algorithm, a one-step procedure for solving differential equations with boundary conditions. We set the initial primordial dust column density to according to the radial column density profile in Eq. (5). All other populations are set to zero, i.e.
Our simulation utilizes parameters introduced in Sect. 2 and summarized in Table 1. In our parametric study, we vary R, M disk and . The default value for the distance to the star R = 10 AU was chosen arbitrarily and corresponds roughly to the position of the planet Saturn in the solar system (semi-major axis of Saturn: 9.537 AU). M disk = 0.02 × M is roughly twice the mass of the Minimum Mass Solar Nebula (MMSN) as studied in Weidenschilling (1977) and Hayashi (1981). The MMSN will be discussed in more detail later on, when looking at the influence of the disk mass on the simulation results. For the trap efficiency we will assume a default value of = 0.01. A brief discussion of the influence of St pbb is shown in Appendix B. The influence of the other parameters in Table 1 is not investigated in this work and they remain constant.

Local evolution
We present our simulation results by first describing the column density evolution of a fixed parameter set and then later discussing the change resulting from parameter variations. Figure 2 displays the local evolution of normalized column densities Σ dst,prim , Σ dst,col , Σ pbb,prim , Σ pbb,col and Σ pls . The left panel depicts scenario 1, where the dust-sized fragments of planetesimal collisions grow further to pebble-sizes. This implies the existence of an equilibrium state, which begins after 2 × 10 6 yr with this particular set of parameters. It is not surprising that the population of primordial dust follows a monotone decline, since the corresponding differential Eq. (37c) contains only a sink term. There is a short period during the local evolution where pebble population makes up a significant portion of the column density -in Fig. 2 this is from 4 × 10 3 to approximately 2 × 10 4 yr after the start of the simulation. As we will discuss later, both the duration and the peak height of the pebble column density of this period will vary strongly with the choice of parameters, in particular the trapping efficiency . The planetesimal population comprises over 98% of the total available column density, i.e. the mass available to the system starting around 5 × 10 4 yr. Collisions between planetesimals naturally will be more common with more planetesimals available leading to the emergence of the two fragment species: collisional pebbles and collisional dust. However, both only make up less than 2% of the total column density once they reach their peak during the equilibrium state. The right panel in Fig. 2 shows scenario 2. We see that per definition the collisional dust column density monotonically increases, since its differential Eq. (37e) only consists of a source term. We note that until ∼10 6 yr, the left and right panel are almost congruent with each other. This is expected. The growth of collisional dust, which is the fundamental distinction of scenario 1 and 2, only becomes relevant once planetesimals collisions produce a noteworthy amount of fragments, which for this set of parameters takes about 10 6 yr.

Parameter study
To develop an understanding of the column density evolution in a system it is crucial to explore various parameters. For this purpose, we chose to focus on the scenario 2 configuration (see Fig. 1), and base all simulation runs presented in this section on it.

Variation of the distance to the star
An insight in a possible global evolution can be enabled by executing several simulations for multiple distances to the star R.
The result of such a simulation run (20 different local simulations at various radii 1 AU ≤ R ≤ 100 AU all using = 0.01) is presented in Fig. 3, where the column density of the populations is displayed versus R for several snapshots during the simulation. The top left panel shows the primordial dust population, where the decline in time after Eq. (16) is clearly visible. Moreover, the dust column density declines more rapidly closer to the star, resulting in an inside-out growth of material (also found in e.g., Birnstiel et al. 2012). This behavior is expected when one considers the R-dependence of the growth timescale in Eq. (13), which is much shorter for small R than in the outer disk.
The same effect can be seen upon inspection of the top right panel in Fig. 3, which displays the primordial pebble column density. In the beginning at 10 3 yr, the pebble population increases in the inner disk much faster than in the outer disk resulting in a much steeper profile compared to the primordial dust population. Further, we can identify that removal of primordial pebbles via planetesimal formation likewise is an inside-out process, as a direct implication of inside-out dust growth.
The planetesimal population naturally mirrors this behavior too, as pictured in the bottom left panel of Fig. 3: until 2.5 × 10 5 yr, the planetesimal column density increases more rapidly at smaller radii. However, once the local pebble supply is depleted, planetesimal formation slows down drastically and collisions start to decrease the planetesimal local column density again. This also happens first at smaller radii, because here, the pebble supply is depleted first. Therefore, the peak of the planetesimal column density also moves toward the outer disk with increasing time, resulting in an inside-out formation of planetesimals as also shown by Lenz et al. (2019). The same cannot be argued for further growth to planets. Here, in addition to the radial planetesimal mass distribution, dynamical stirring and the available pebble supply become relevant, introducing complex R-dependencies.
Finally, the bottom right panel of Fig. 3 displays the column density of collisional dust particles, which will eventually make A116, page 7 of 17  up the bulk of the column density as already seen in the right panel of Fig. 2. For a fixed R, it increases all the time, again first in the inner disk and later in the outer disk, where it never surpasses the planetesimal column density.

Variation of the efficiency of planetesimal formation
Further, we investigate the influence of the trap efficiency on the column density evolution. All other parameters are kept constant according to the default parameter set. The trap efficiency is per definition 0 ≤ ≤ 1, with = 1 implying that the entire incoming pebble flux is trapped and converted into planetesimals and = 0 meaning that the trap is not capturing and forming planetesimals. A comparison of the local column density evolution for different -values ( = 0.001, 0.01, 0.1, 0.5) is shown in Fig. 4. The top panel depicts the evolution of the planetesimal column density. For all values of the general trend already seen in the right panel of Fig. 2 is replicated. The column density of planetesimals increases until they make up the vast majority of the column density available in the system. However, the rate of planetesimal formation strongly depends on . During the initial phase of planetesimal formation (until 10 4 to 10 6 yr depending on the value of ) the planetesimal column density is linear in after Eq. (21). However, this linearity ceases to exist once the normalized column density approaches unity.
On the top panel of Fig. 4, the curves for the different -values seem to roughly coincide once the decline of the planetesimal population sets in. This is not surprising, since the collision rate in Eq. (34) does not depend directly on . Here, planetesimal formation is insignificant, because there is only few pebbles available. We point out, that for all -values almost all the mass ends up within the planetesimal population -for = 0.5, 0.1, 0.01, 0.001 starting at approximately 2 × 10 4 , 5 × 10 4 , 6 × 10 5 and 10 7 yr, respectively.
The bottom panel in Fig. 4 displays the evolution of primordial dust (solid lines) as well as collisional dust (dashed lines) for the same four values of the trap efficiency parameter . In the middle panel, one can see an equivalent plot for primordial pebbles (solid lines) and collisional pebbles (dashed lines). The influence of on the evolution of primordial dust is negligible, since the dust growth timescale does not directly depend on , as seen in Eq. (13). However, both the peak height and width, as well as the rate of decrease of the primordial pebble population strongly depend on as shown in the plot. Inefficient traps leave behind a significant population of primordial pebbles, while more efficient traps can drain the pebble supply more rapidly. The collisional pebble population peaks once production of collisional pebbles can not keep up with drainage of pebbles through planetesimal formation, i.e. in Eq. (37d)Σ pbb,col = 0. This is strongly influenced by the column density available in primordial pebbles, since they outnumber their collisional counterparts at early times. The intersection with the horizontal axis of both the collisional dust and the collisional pebbles curve depends on , because is what determines when enough planetesimals can be formed for collisions to produce noteworthy amounts of fragments (see top panel of Fig. 4).

Variation of the initial disk mass
Protoplanetary disks are observed with a number of different masses. Andrews et al. (2010) found disk masses in the range of M disk = 0.004-0.143M in the ∼1 Myr old Ophiuchus starforming region. The MMSN model discussed in Weidenschilling (1977) describes a disk of solar composition containing the minimum amount of solids necessary to form the eight planets and the asteroid belts of today's solar system, by only considering their masses and today's positions. In the MMSN model, rocky and icy objects have a total mass of roughly 2 × 10 −4 M while gas has a total mass of 1.3 × 10 −2 M . Although effects such mass loss via photoevaporation are not considered, it is still expedient to execute parameter variations of the disk mass.
Such a parameter study is shown in Fig. 5. We define a dimensionless disk mass and investigate values ranging from µ = 0.01, approximately corresponding to the MMSN after Weidenschilling (1977) up to µ = 0.07. More massive disks are expected to be gravitationally unstable and are additionally uncommon when taking recent observations by Andrews et al. (2009Andrews et al. ( , 2010 into account. The other model parameters are kept constant. The normalized column density is depicted in the left panels whereas the right panels show its absolute value. Panels are arranged in analogy to Fig. 4: the top panels show the evolution of the planetesimal column densities. The maximum absolute planetesimal column density naturally increases with disk mass, since a higher disk mass implies more available material to form planetesimals. Interestingly this trend is reversed when looking at the normalized planetesimal column density evolution in the top left panel of Fig. 5. Less massive disks have a lower relative planetesimal fraction than more massive disks , because the effect of planetesimal collisions becomes more important with increasing disk mass. The middle panels of Fig. 5 display the primordial (solid lines) and the collisional pebble population (dashed lines). In the bottom panels, one can see the primordial (solid lines) and collisional (dashed lines) dust populations. The absolute primordial dust and primordial pebble content increases with increasing disk mass. However, there seems to be no strong correlation between normalized primordial column density and disk mass (middle left and bottom left panels of Fig. 5). Only starting at ∼5 × 10 5 yr, the panel shows that the primordial dust supply is drained faster for more massive disks. The peak height of the absolute collisional pebble column density in the middle right panel correlates strongly with disk mass. This trend is replicated for the normalized column density. For the most massive disks pebbles make up around 30% of the total available column density at the peak of the pebble population (roughly 4 × 10 5 yr). For MMSN-like disks, the peak height is not only lower (∼2% of available column density), it also occurs slightly later during the evolution (10 6 yr). Similarly, both absolute and normalized collisional dust column density evolution correlate with disk mass. For massive disks the fraction of collisional dust is generally larger than for less massive disks and also starts to increase somewhat sooner (5 × 10 4 yr for the most massive disk compared to 10 5 yr for the MMSN disk) and with a steeper slope.
In our model, evolution of dust, pebbles and planetesimals is independent of the turbulence parameter α, but only valid for disks where turbulence is strong enough to set relative dust velocities as discussed in Sect. 2.2. Further and more importantly, the pebble Stokes number is for simplicity kept constant at St pbb = 0.1 for our parameter study. We discuss implications in Appendix B. A more complex model would introduce a dependence on α as seen in Birnstiel et al. (2009). Additionally, turbulence influences radial transport of disk material strongly, which is also neglected in this simple model. Finally, in this paper, we will refrain from performing parameter studies with 0 dg and r pls since their influence on the column density evolution is not the main focus of this work.

Estimating the evolution of the mass distribution
The evolution of the mass that is contained within planetesimals is of particular interest since planetesimals are too cold to be observed in the infrared and they also provide too little column density to be detected via scattered light observations. Hence, this section presents an approximation of the evolution of the mass distribution. The total mass M i of a given population i that is contained within R 1 and R N at a certain time t is given by integrating the column density at t over r We can approximate this integral by performing multiple simulation runs k at different radii R k and then summing over them, i.e., using the trapezoidal rule, We set N = 20, R 1 = 1 AU, and R 20 = 100 AU. For radii larger than 100 AU, the column density profile in Eq. (5) declines exponentially and prevails the r 2 -dependence in Eq. (41), rendering the contribution to the total mass beyond this point insignificant. We point out that this procedure is only expected to deliver a rough estimate of the total mass evolution in the disk.
For a more precise picture, spatial mass transport certainly has to be considered. Figure 6 displays the evolution of the normalized mass of solid particles in the disk, as obtained via the algorithm explained above. Different colors correspond to different species (yellow for dust, blue for pebbles and red for the planetesimal population). Black dashed lines represent the sum of all dust and pebble particles, since we can not expect to be able to distinguish these species with observations. Different panels originate from simulation runs with different values for trap efficiency ( ∈ {0.001, 0.01, 0.1, 0.5}) and disk mass µ = M disk /M (µ ∈ {0.01, 0.02, 0.05}, where the lower boundary roughly approximates the MMSN-mass from Weidenschilling (1977), which is indicated in the top right corner of each panel. We note that while µ represents the total disk mass containing both solids and gas, the vertical axis in each panel of Fig. 6 only displays the mass of each of the solid species. The data is normalized with respect to the total mass of solid particles ( dg M disk ). The mass evolution for each panel is comparable to the column density evolution in Fig. 2. In the beginning, all mass in the system is contained within the combined dust population. After a period during which the combined pebble population makes up a nonnegligible part of the total mass, the planetesimal population starts to dominate. Once the collisions take over, the mass is being transferred back to the combined dust population. When Left panels: the disk mass is fixed at µ = 0.02 while the trap efficiency is fixed to = 0.01 in the right panels. Colors correspond to different times in the disk evolution. Data is taken from the same simulation runs that were depicted in Fig. 6. However, Fig. 7 does not differentiate between primordial and collisional populations and depicts their sum.
comparing the three panels on the right side as well as the middle left panel, each corresponding to different -values for a constant disk mass of µ = 0.02, one discovers, that both influences the time where the maximal planetesimal mass is reached and its peak value. Similar to Fig. 4, a higher -value implies that the planetesimal population peaks at earlier times and with a greater maximal peak height. In the top right plot, where the trap efficiency parameter was set to the low value of = 0.001, the planetesimal population (red line) never becomes more massive than combined dust-and-pebble population (green line). It's peak is roughly equivalent to a mass of 26.8 M Earth . A tenfold increase of implies approximately a 20% increase of the maximal planetesimal mass and this peak occurring roughly ten times earlier.
In the left panels, disk mass increases from the top to the bottom panel while trap efficiency remains constant at = 0.01. As expected from Fig. 5, more massive disks tend to have a smaller relative planetesimal fraction than less massive disks, leading to a lower ratio of planetesimal mass to dust-and-pebble combined mass for more massive disks. An interesting observation is that in the top left panel, the maximal planetesimal mass is also approximately 26.8 M Earth similar to the top right panel, implying a doubled normalized disk mass µ would be approximately balanced by a tenfold decrease of . When observing these model disks at the times of the respective planetesimal mass peak, one would conclude significantly different dust masses, even though the planetesimal mass is comparable. The evolution of the dust mass in the top left panel would for example compare better to the bottom right panel with = 0.5, where planetesimals reach the mass fraction peak earlier than in the top left panel.
In Fig. 7, we further visualize how the mass in solids is distributed among the populations and how this depends on the efficiency (left panels) and the dimensionless disk mass µ (right panels). Data are taken from the same simulation runs as Fig. 6, while colors indicate the time the snapshot was taken, i.e. the disk age. The evolution of the dust content (bottom panels) does only depend on and µ for times 10 5 yr. This is because only then planetesimal collisions start to produce a significant amount of dust. The amount of these fragments is determined by the planetesimal column density, which depends on and µ as seen in the top panels of Fig. 7. Here, while the mass in planetesimals always increases with , the relation is clearly nonlinear. This is explained due to the fact, that the pebble supply is limited and therefore the maximum planetesimal column density is independent of , as seen Fig. 4. In the top right panel of Fig. 7, we can again observe that the planetesimal fraction is indeed decreasing with increasing disk mass, i.e. less massive disks are more efficient in forming planetesimals, which is in accordance to Fig. 5. The behavior of the pebble fraction is depicted in the middle panels of Fig. 7. It is approximately constant in µ, because the initial growth phase from dust to pebbles only depends on the dust-to-gas ratio, which is picked the same for all µ. The small differences originate from the fact that the significance of planetesimal collisions depends on total mass. However, the pebble fraction correlates strongly with (middle left panel), because more efficient planetesimal formation implies that the pebble supply is drained faster. We also observe that during late times at t = 9.2 × 10 6 yr, the mass in pebbles is rather independent of . This is due to planetesimal formation being insignifcant in comparison to planetesimal collisions, which do not depend on .

Limitations and advantages of the model
Our model is limited by various factors. Most strikingly, its locality and the resulting absence of spatial transport of material constitutes a strong limitation, as our model can not currently cover scenarios where the pebble drift timescale is shorter than the conversion timescale for pebbles into planetesimals. The 0-dimensional model can also not consider any effects proceeding perpendicular to the disk plane, such as dust settling and vertical turbulent stirring. The dependence of the model on disk turbulence was eliminated by assuming that particle growth does not scale with disk turbulence, i.e. St α. Additionally, viscous gas evolution and photoevaporation (see e.g., Ercolano et al. 2009;Owen et al. 2011;Nakatani et al. 2018;Picogna et al. 2019) were also not considered. In comparison to Birnstiel et al. (2012), who derive the maximum pebble Stokes number by taking the dust growth limiting effects drift and fragmentation into account, we fix the pebble Stokes number to a guiding value at St pbb = 0.1 during the entire disk evolution. The influence St pbb A116, page 12 of 17 is discussed in Appendix B. This could be improved by implementing these growth barriers, which depend on gas and particle density. This would allow for a time-and space-dependent pebble Stokes number as in Birnstiel et al. (2012). Alternatively, dust coagulation could be modeled in more detail by implementing a particle size grid and solving the Smoluchowski equation (Smoluchowski 1916), as was done by Lenz et al. (2019), who we further compare our results to in Appendix A.
Many observed circumstellar disks contain substructures, mainly rings as first seen by ALMA Partnership (2015) in HL Tauri or recently in Andrews et al. (2018) presenting results from the "Disk Substructure at High Angular Resolution Project". As the underlying physical process is still in debate, as recently discussed in van der Marel et al. (2018), disk substructures are not considered in this work's model for gas column density evolution. Furthermore, ice lines, which mark the border where condensation of volatiles such as for example water is possible, were not considered. The Clausius-Clapeyron (Clausius 1850) relation describes the slope of the tangent dividing the two phases in a pressure-temperature diagram, which for typical concentrations of H 2 O would correspond to approximately T g 150 K-200 K. The addition of the water ice line would imply a significant kink in the dust profile as seen for example in the model by Drażkowska & Alibert (2017). In their work, ice lines mark a favorable location for planetesimal formation. Further not considered is the accretion of solids onto the central star, removing material in the inner disk and decelerating or even preventing planetesimal formation in these regions.
Chemical composition of solids and how it depends on R is also not part of this model. One may expect the composition of a dust particle to influence its growth growth rate. Additionally, the composition of a planetesimal affects the size distribution of collisional fragments (Johnson et al. 2012).
The fact that the local model does not consider different planetesimal sizes, is also a restriction, as the outcome of a collision is dependant on the mass and sizes of the colliding planetesimals. The temporal evolution of the planetesimal size distribution is not considered. Further, we neglect any planet-forming processes such as pebble accretion (see e.g., Ormel & Klahr 2010;Ormel 2017;Rosenthal et al. 2018;Lambrechts et al. 2019), which may drain the pebble supply and therefore decelerating the birth of new planetesimals. At planetesimal sizes of ∼100 km, the efficiency of pebble accretion is minimal. Therefore, our assumption of a fixed planetesimal size of ∼100 km does not allow for pebble accretion to be efficient. Likewise, planetesimal accretion (e.g., Kokubo & Ida 2012) is also neglected.
Lastly, various limiting assumptions go into this work's model for planetesimal formation, as we condense unknown information about the physical nature of particle traps in the trapping efficiency parameter . Besides the streaming instability (Youdin & Goodman 2005), other hydrodynamical processes such as subcritical baroclinic instability (Klahr & Bodenheimer 2003), convective overstability (Klahr & Hubbard 2014), or vertical shear instability (Urpin & Brandenburg 1998) may also affect the pebble trapping mechanism . While Pfeil & Klahr (2019) show that hydrodynamic instabilities can act throughout the entire disk, the one that dominates the formation of traps may depend on the location in the disk. Therefore, it is plausible that may also depend on R, as well as on the pebble Stokes number.
Strength of our local model is, that it considers three defining processes in planetesimal formation and evolution theory and integrates them into one simple mode, which is both easy to understand and due to its locality easy to implement as well as computational inexpensive. Our model is able to form significant column densities of planetesimals everywhere in the disk and fast (∼10 6 yr), i.e. well before the onset of gas dissipation. Therefore, our toy model can be used to inform initial conditions of late-stage protoplanetary disk evolution models.

Summary and conclusions
We investigated a local model of pebble flux-regulated planetesimal formation (Lenz et al. 2019), dust growth (Birnstiel et al. 2012) and planetesimal collisions to study the evolution of the planetesimal to dust and pebble ratio using an enclosed set of coupled differential equations. Aim of the performed parameter study was to develop an understanding of the potentially existing mass in the planetesimal population for a certain amount of dust and pebbles and how it relates to our model parameters.
In Sect. 2, we discussed the principles on which our local model is based on. To each of the three processes dust growth, planetesimal formation, and planetesimal collisions a timescale is attached to, which allows the formulation of rate equations (see Eqs. (37a)-(37e)) describing the evolution of the enclosed system. To solve this set of differential equations a code following a one-step procedure was written. To keep the code simple and numerically inexpensive, radial transport of material and gas evolution was neglected, besides planetesimals, only two small particle populations -dust and pebbles -were considered, and the Stokes number of the latter was assumed to stay constant at St pbb = 0.1, which is roughly the maximum Stokes number in Birnstiel et al. (2012). Implications are discussed in Appendix B. The small particle populations were further differentiated into primordial/ pristine particles and fragments from planetesimal collisions. We discussed two scenarios. In the first, collisional dust was allowed to grow back to pebbles and form planetesimals again and in the second this was forbidden. In the former scenario, this leads to a steady state with constant column densities as visualized in the left panel of Fig. 2. Here, column densities remain constant. However, for our parameter study we used the latter scenario, which is the more realistic approach, because fragments of planetesimal collisions are assumed to be too compact to undergo growth via sticking. Here, the local evolution can be broken down into three phases.
During the first stage of the disk evolution, primordial dust grows to pebble sizes, the dust population shrinks to one hundredth of its initial value in roughly 10 5 yr. For the behavior of primordial dust the simple analytical solution in Eq. (16) can be found. The phase in which planetesimals are the dominating species naturally follows the increase of pebbles that are available for planetesimal formation. This occurs fast, i.e. within ∼10 6 yr, and everywhere in the disk. As a feedback, the pebble population is drained again, limiting planetesimal formation. At the same time collisions between planetesimals produce collisional fragments in both pebble and dust sizes. While the former species can form new planetesimals and prolong the planetesimal dominated phase, the collisional dust species starts to steadily increase. Once the pebble supply is depleted, the planetesimal population shrinks and dusty fragments start to become the dominating species in the system. This occurs at approximately 10 7 yr. Beyond 10 7 yr gas dissipation cannot be neglected (Hernández et al. 2007;Mamajek 2009;Fedele et al. 2010;Pfalzner et al. 2014) and our model assumptions fail.
In our parameter study, we focused on the effects of distance to the star R, trap efficiency and disk mass M disk . We note again, that incorporates both the efficiency of the pebble trapping mechanism as well as the efficiency of planetesimal A116, page 13 of 17 A&A 629, A116 (2019) formation itself. We summarize our parameter study below. The dust growth timescale is much shorter close to the star than in outer regions, which has the direct implication that the entire local evolution is occuring on shorter timescales for small R than for larger R as seen in Fig. 3. This inside-out growth of dust is in line with e.g., Brauer et al. (2008); Birnstiel et al. (2010Birnstiel et al. ( , 2012 or Krijt et al. (2016). Inside-out formation of planetesimals agrees with Lenz et al. (2019). High -values generally imply an acceleration of the planetesimal formation process, leading to a longer planetesimal dominated phase (see Fig. 4). When studying different disk masses as shown in Fig. 5, we found that while more massive disks have a higher absolute planetesimal column density than less massive disks, their normalized planetesimal column density is lower, because planetesimal collisions are more significant in massive disks compared to less massive disks.
Lastly, we showed an estimate for the mass evolution of the disk in Figs. 6 and 7. Here, we displayed how the distribution of material among the three populations evolves in time and depends on planetesimal formation efficiency and disk mass. Therefore, our local toy model can be a potential tool for constraining the material in planetesimals in observed protoplanetary disks.
By investigating the influence of and M disk on the mass evolution in Fig. 6, we are able to make statements on the potentially existing mass in planetesimals for a given amount of pebbles and dust. Provided the total mass of a certain disk is known and measured independent from dust mass by observations such as e.g., Pascucci et al. (2016), or modeling, the more dust is observed the less mass can be expected to be compromised in planetesimals relative to the total disk mass. If one can determine the age of the protoplanetary disk, e.g., via the position of the central star (or its neighbors) in the Hertzsprung-Russel diagram, and provided both total disk mass and dust portion relative to it are known from empirical data, the toy model may constrain the value of the trap efficiency as seen in Figs. 6 and 7. If one, however, is unsure about the total mass of a certain disk, which work by Andrews et al. (2010) suggest is typically the case, it becomes challenging to predict the planetesimal mass for a given observed dust mass, since the dust and pebble column density evolution can look almost identical for two different sets of trap efficiency and disk mass M disk . Despite many simplifications made in our model, it generated results connecting important model parameters like trap efficiency and disk mass M disk to the local column density evolution of three different species dust, pebbles and planetesimals. Additonally, it helped us to further our understanding for the evolution of the mass distribution over the course of over 10 7 yr ranging from the birth of the protoplanetary disk until its dissipation.  (Smoluchowski 1916), instead of approximating dust coagulation via the two population model, which tends to underestimate growth timescales and thus overestimates particle growth (Birnstiel et al. 2012). For this reason, dust growth and therefore also planetesimal formation sets in sooner in our model compared to Lenz et al. (2019). Further, Lenz et al. (2019) use the same prescription for planetesimal formation, which is reflected in the slope of the planetesimal mass evolution of the two models in Fig. A.1. However, once planetesimal collisions become relevant, i.e. when about 30% of the total mass available in solids is comprised in planetesimals, our model predicts a smaller planetesimal fraction than Lenz et al. (2019).

Appendix B: On treating the pebble Stokes number as a fixed parameter
Following Birnstiel et al. (2012), the most abundant pebble size a resides at the fragmentation limit (typically in the inner disk) or at the drift limit (typically in the outer disk), i.e. in terms of Stokes numbers with St frag being the maximum Stokes number in a fragmentation dominated size distribution, and St drift the maximum Stokes number that can be reached before drift removes the particle. In our model, we treat the pebble Stokes number as a fixed parameter with St pbb = 0.1. This is a very good approximation for α = 10 −3 (see Lenz et al. 2019, Fig. 3). In more turbulent disks (α ≥ 10 −2 ) the higher relative velocities of dust particles lead to a decrease of the fragmentation size and our approximation of St pbb = 0.1 is too high, especially in the inner disk. Hence, we investigate the influence of the pebble Stokes number on the local evolution in Fig. B.1. For overestimates of the pebble Stokes number, the dustto-pebbles growth timescale will also be overestimated, since τ growth ∝ ln St pbb after Eq. (13). Because of this, the pebble population in Fig Conversely, planetesimals form slower for smaller Stokes numbers, as we expect a decrease in drift velocity as portrayed in Eq. (18). As this decrease is dominating over the increase of available pebbles, the effect of the pebble Stokes number is similar to the effect of the trapping efficiency depicted in Fig. 4.
In summary, for smaller pebble Stokes numbers, we expect more pebbles, though they drift significantly slower resulting in a smaller pebble flux and less planetesimals. In our parameter study, this error gains importance for disks with moderate or high turbulence, and typically is in the inner disk more significant than in the outer disk.