© K. Gerbig et al. 2019

1 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. 2009, 2013; Ansdell 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⁷ yr to be in line with typical lifetimes of protoplanetary disks (Hernández et al. 2007; Mamajek 2009; Fedele et al. 2010; Pfalzner et al. 2014).

The model principles are described in detail in Sect. 2, in particular the three processes of dust growth, planetesimal formation, and planetesimal collisions. They lead to a set of differential equations, which we present for two scenarios in Sect. 3. There we additionally discuss our numerical setup, that is our parameters, initial conditions, and simulation run time. We present our results obtained from a parameter study in Sect. 4, where we further estimate how the total mass in the disk is distributed among the three species dust, pebbles, and planetesimals. In Sect. 5, we discuss model limitations. Finally, we summarize and conclude in Sect. 6.

2 Model principles

2.1 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 (1)

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 (2)

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 (3)

where we introduced the initial dust-to-gas ratio (4)

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), (5)

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 molecularhydrogen is given by (6)

with Boltzmann constant k_B = 1.3807 × 10⁻¹⁶erg K⁻¹, and proton mass m_p = 1.673 × 10⁻²⁴ g. Further, the vertical pressure scale height of the gas disk h_g is defined as (7)

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 (8)

where G ≈ 6.674 × 10⁻⁸ cm³ g⁻¹ s⁻². 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): (9)

to which the column density is inversely proportional in the inner disk. We assume a temperature profile of (Chiang & Goldreich 1997) (10)

where we choose stellar properties of R_star = 1.25 × R_⊙≈ 1.25 × 7 × 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⁻⁴ ≤ α ≤ 10⁻² 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) and Flaherty et al. (2017).

2.2 Dust growth

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; 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 (11)

where a₀ 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.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 (12)

where ρ is the mass density of dust. We choose ρ = 1.2 g cm⁻³, 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-to-pebbles growth timescale in the midplane with (13)

The approximation in Eq. (13) holds if relative dust velocities are set by disk turbulence with St_dst ≪ 1, i.e. for sufficiently large α-values.

The dust to pebble growth rate is calculated via (14)

The differential equation (15)

illustrates the decrease of the local dust supply due to particle growth to pebble size. It is solved analytically by (16)

Here, is the initial dust column density from the profile in Eq. (5) and (17)

is constant in time for a fixed St_pbb.

2.3 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 (18)

where we introduced the exponent of the gas pressure power law (19)

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⁻¹ quickly driftinward toward the central star on very short timescales. Such a particle at 10 AU can be expected to reach the star in under 10⁴ 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 Hill-density, it can collapse and form a planetesimal as discussed in Johansen et al. (2006, 2007). 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 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 (20)

We further express the planetesimal formation rate by (21)

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 (22)

where τ_ff is the free-fall timescale (Shariff & Cuzzi 2015; Klahr & Schreiber 2015). Hence, we expect pebbles to dominate the mass contributionof 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): (23)

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 (24)

If Eq. (24) is not fulfilled, 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 objects – relics 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 (Johansen et al. 2007; Cuzzi et al. 2008; Klahr & Schreiber 2015). Assuming homogeneous and spherical planetesimals with mass density ρ = 1.2 g cm⁻³, the planetesimal mass is approximated with (25)

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.

2.4 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 (26)

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 . 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₀ ≤ m < m₁, for pebbles m₁ ≤ m < m₂ and for planetesimals m₂ ≤ m ≤ M. The derivations of m₁ and m₂ are shown in Appendix C.2.

With these requisites, one can solve the integral 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 andplanetesimals, respectively. For typical disk parameters (see Table 1), m₁, m₂ ≪ 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 collisionsbetween particles i and j is σ_ij, then (Birnstiel et al. 2016) (27)

The collision timescale τ_col,i is defined as the average time during which particle i experiences one encounter (28)

where Δv_ij 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: (29)

where the escape velocity v_esc of the planetesimals is (30)

By approximating the vertical planetesimal distribution with a Gaussian with a root mean square width of h_pls = v_hillΩ⁻¹ (Goldreich et al. 2004), the number density of planetesimals can be written as (31)

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) (32)

Plugging the above expressions into Eq. (28) we get for the collision timescale of equal-sized planetesimals (33)

We follow Morbidelli et al. (2009) and insert the Hill velocity for the relative velocity of planetesimals Δv_ij. 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 must also depend on the available planetesimal surface density. It is given by (34)

3 Numerical model setup

3.1 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 and each appear once as a sink and a sourceterm and p_dst + p_pbb + p_pls = 1 per definition, the total column density in the system is conserved in time (36)

In this configuration, we expect a column density steady state to occur.

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 equations

where

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. Contraryto 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.

3.2 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.

Fig. 1 Flowchart visualizing different column density transfer processes in scenario 2 configuration, where collisional dust can not grow. Scenario 1 would simply add an arrow connecting collisional dust to the dust growth triangle, rendering the primordial populations equivalent to their collisional counterpart.

4 Results

4.1 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⁶ 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³ to approximately 2 × 10⁴ 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 planetesimalpopulation comprises over 98% of the total available column density, i.e. the mass available to the system starting around 5 × 10⁴ 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 oncethey 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 asource term.

We notethat until ~10⁶ 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⁶ yr.

Fig. 2 Local evolution of the different species with column densities normalized by the initial dust column density versus time. Left panel: scenario 1 configuration where collisional dust can grow back to (primordial) pebbles, and reach an equilibrium state. Right panel: scenario 2 configuration, where compact collisional dust can not grow pebble sizes. This simulation was done using the default parameter set R = 10 AU, Mdisk = 0.01 M⊙ and ɛ = 0.1.

4.2 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.

4.2.1 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³ 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⁵ 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 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.

Fig. 3 Column density of different populations is displayed against R for several snapshots. Top left panel: primordial dust, top right, bottom left and bottom right: primordial pebbles, planetesimals and collisional dust, respectively. We depict a combined result of multiple 0-d simulations executed at 20 different radii to allow an insight in a possible global evolution of the dust profile. All simulations shown in this figure used ɛ = 0.01.

4.2.2 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⁴ to 10⁶ 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⁴, 5 × 10⁴, 6 × 10⁵ and 10⁷ 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) . 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).

Fig. 4 Local evolution at R = 10 AU for Mdisk = 0.01M⊙ of the normalized column density for different values of the trap efficiency parameter ɛ (indicated by different colors). Top panel: evolution of the planetesimal column density. Middle panel: both primordial (solid lines) and collisional pebble populations (dashed lines). Bottom panel: primordial dust (solid lines) and its collisional counterpart (dashed lines).

4.2.3 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 star-forming 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⁻⁴ M_⊙ while gas has a total mass of 1.3 × 10⁻²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 (40)

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. (2009, 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 planetesimalcolumn 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⁵ 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⁵ 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⁶ 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⁴ yr for the most massive disk compared to 10⁵ 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 and r_pls since their influence on the column density evolution is not the main focus of this work.

4.3 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₁ and R_N at a certain time t is given by integrating the column density at t over r (41)

We canapproximate this integral by performing multiple simulation runs k at different radii R_k and then summing over them, i.e., using the trapezoidal rule, (42)

We set N = 20, R₁ = 1 AU, and R₂₀ = 100 AU. For radii larger than 100 AU, the column density profile in Eq. (5) declines exponentially and prevails the r² -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 non-negligible 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 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 thatthe 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 expectedfrom 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⁵ 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⁶ 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 ɛ.

Fig. 5 Local evolution at R = 10 AU, ɛ = 0.01 of normalized (left panels) and absolute (right panels) column density of the species. Panels are arranged in analogy toFig. 4: top panels: evolution of planetesimal column densities; middle panels: both primordial (solid lines) and collisional (dashed lines) pebble population. Bottom panels: species of primoridal dust (solid lines) and collisional dust (dashed lines). Different color indicate different values for the disk mass μ = Mdisk∕M⊙∈ [0.01, 0.03, 0.05, 0.07]. Bright green lines approximately correspond to the MMSN model first discussed in Weidenschilling (1977).

Fig. 6 Evolution of the normalized total mass of solid particles in the disk for different values of ɛ ∈ {0.001, 0.01, 0.1, 0.5} and disk mass μ = Mdisk∕M⊙∈{0.01, 0.02, 0.05}. Different colors correspond to the different species (red: planetesimals, yellow: dust, blue: pebbles, black dotted lines: dust and pebbles combined), the solid lines to primordial and the dashed lines to collisional populations. Left panels: disk mass increases from top to bottom while the trap efficiency ɛ remains constant. Likewise, in the right panels: trap efficiency increases downward, however the middle left panel would also fit within this sequence between the top right and the middle right panel.

Fig. 7 Mass fraction in planetesimals (top panels), pebbles (middle panels) and dust (bottom panels) vs. trap efficiency ɛ (left panels) and disk mass μ = Mdisk∕M⊙ (right panels). 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.

5 Limitations and advantages of the model

Our model is limited by various factors. Most strikingly, its locality and the resulting absence of spatial transport ofmaterial 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 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₂ 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 (Klahr et al. 2018). 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⁶ 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.

6 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 inthe second this was forbidden. In the former scenario, this leads to a steady state with constant column densities as visualized inthe 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⁵ 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⁶ 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⁷ yr. Beyond 10⁷ 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 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. (2010, 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 proto- planetary 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⁷ yr ranging from the birth of the protoplanetary disk until its dissipation.

Acknowledgements

K.G. thanks Ruth A. Murray-Clay and U.C. Santa Cruz for hosting him as a graduate visitor. C.L. acknowledges funding from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) as part of the Schwerpunktprogramm (SPP, Priority Program) SPP 1833 “Building a Habitable Earth”. Finally, the authors are grateful to the anonymous referee for an insightful and detailed report which greatly helped to improve this paper.

Appendix A Comparison to the planetesimal synthesis model in Lenz et al. (2019)

In Fig. A.1, we compare the integrated mass of planetesimals and smaller particles in our model to Lenz et al. (2019), where planetesimal collisions are not implemented. Moreover, Lenz et al. (2019) solve the Smoluchowski equation forgrowth and fragmentation of particles (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).

Fig. A.1 Evolution of normalized total mass of solid particles in the disk for our model (solid lines) and the model in Lenz et al. (2019) (dashed lines) for α = 10−3. Both models used ɛ = 0.1, μ = 0.05, RC = 35 AU, and the same temperature profile (Lenz et al. 2019, Eq. (25)). Times that are not accounted for in Lenz et al. (2019) are shaded in grey.

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 (B.1)

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⁻³ (see Lenz et al. 2019, Fig. 3). In more turbulent disks (α ≥ 10⁻²) 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.

Fig. B.1 Local evolution at R for ɛ = 0.01 and μ = 0.01 of normalized column density of the species for different pebble Stokes numbers (dotted lines: Stpbb = 1, solid lines (see Fig. 2): Stpbb = 0.1, dashed lines:Stpbb = 0.01) and different radii (colors). Each panel corresponds to a different particle species: primordial dust and pebbles, planetesimals, collisional dust and pebbles (from top to bottom).

For overestimates of the pebble Stokes number, the dust-to-pebbles growth timescale will also be overestimated, since after Eq. (13). Because of this, the pebble population in Fig. B.1 grows faster for smaller Stokes numbers. 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.

Appendix C Details of the planetesimal collision model

C.1 Normalizing the specific number column density power law

In our model, colliding planetesimals are equally sized and therefore equally massive. The surface density must be conserved during a collision, i.e. (C.1)

where σ_pls is a parameter specifying the surface density of a single planetesimal, and m₀ and M the masses of the smallest and largest possible fragment respectively. We set and M = m_pls for destructive collisions. We solve the integral in Eq. (C.1), giving (C.2)

leading to (C.3)

We note that the choice of σ_pls is irrelevant, because it cancels when calculating the relative fractions p_dst, p_pbb and p_pls.

C.2 Calculation of the transition masses m₁ and m₂

The drift velocity of a particle such that it can just be considered a pebble and not dust, is equal to the required velocity for it to cover one trap distance d during exactly one trap lifetime (C.4)

If this condition is met, a particle with mass m₁ has just enough time to participate in planetesimal formation via graviational instability and can therefore be considered a pebble. Plugging Eq. (12) into the drift velocity in Eq. (18) yields a second order polynomial with the smaller solution being (C.5)

where we used .

Following Birnstiel et al. (2012), we compare a particle’s drift timescale to its growth timescale, which yields the maximum Stokes number, that can be reached before it is removed by radial drift, i.e. (C.6)

Particles more massive than m₂ grow faster than they are removed from their position, potentially forming planetesimals by growth and not by gravitational collapse, which is not the focus of this paper. Hence, we assign all particles with mass m > m₂ to the planetesimal population. The condition in Eq. (C.6) also results in a quadratic equation after using Eqs. (12), (13) and (18). We then define the larger solution as the transition of pebbles to planetesimals in the fragment distribution: (C.7)

References

All Tables

All Figures

	Fig. 1 Flowchart visualizing different column density transfer processes in scenario 2 configuration, where collisional dust can not grow. Scenario 1 would simply add an arrow connecting collisional dust to the dust growth triangle, rendering the primordial populations equivalent to their collisional counterpart.
In the text

Fig. 2 Local evolution of the different species with column densities normalized by the initial dust column density versus time. Left panel: scenario 1 configuration where collisional dust can grow back to (primordial) pebbles, and reach an equilibrium state. Right panel: scenario 2 configuration, where compact collisional dust can not grow pebble sizes. This simulation was done using the default parameter set R = 10 AU, Mdisk = 0.01 M⊙ and ɛ = 0.1.

In the text

Fig. 3 Column density of different populations is displayed against R for several snapshots. Top left panel: primordial dust, top right, bottom left and bottom right: primordial pebbles, planetesimals and collisional dust, respectively. We depict a combined result of multiple 0-d simulations executed at 20 different radii to allow an insight in a possible global evolution of the dust profile. All simulations shown in this figure used ɛ = 0.01.

In the text

Fig. 4 Local evolution at R = 10 AU for Mdisk = 0.01M⊙ of the normalized column density for different values of the trap efficiency parameter ɛ (indicated by different colors). Top panel: evolution of the planetesimal column density. Middle panel: both primordial (solid lines) and collisional pebble populations (dashed lines). Bottom panel: primordial dust (solid lines) and its collisional counterpart (dashed lines).

In the text

Fig. 5 Local evolution at R = 10 AU, ɛ = 0.01 of normalized (left panels) and absolute (right panels) column density of the species. Panels are arranged in analogy toFig. 4: top panels: evolution of planetesimal column densities; middle panels: both primordial (solid lines) and collisional (dashed lines) pebble population. Bottom panels: species of primoridal dust (solid lines) and collisional dust (dashed lines). Different color indicate different values for the disk mass μ = Mdisk∕M⊙∈ [0.01, 0.03, 0.05, 0.07]. Bright green lines approximately correspond to the MMSN model first discussed in Weidenschilling (1977).

In the text

Fig. 6 Evolution of the normalized total mass of solid particles in the disk for different values of ɛ ∈ {0.001, 0.01, 0.1, 0.5} and disk mass μ = Mdisk∕M⊙∈{0.01, 0.02, 0.05}. Different colors correspond to the different species (red: planetesimals, yellow: dust, blue: pebbles, black dotted lines: dust and pebbles combined), the solid lines to primordial and the dashed lines to collisional populations. Left panels: disk mass increases from top to bottom while the trap efficiency ɛ remains constant. Likewise, in the right panels: trap efficiency increases downward, however the middle left panel would also fit within this sequence between the top right and the middle right panel.

In the text

Fig. 7 Mass fraction in planetesimals (top panels), pebbles (middle panels) and dust (bottom panels) vs. trap efficiency ɛ (left panels) and disk mass μ = Mdisk∕M⊙ (right panels). 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.

In the text

	Fig. A.1 Evolution of normalized total mass of solid particles in the disk for our model (solid lines) and the model in Lenz et al. (2019) (dashed lines) for α = 10−3. Both models used ɛ = 0.1, μ = 0.05, RC = 35 AU, and the same temperature profile (Lenz et al. 2019, Eq. (25)). Times that are not accounted for in Lenz et al. (2019) are shaded in grey.
In the text

Fig. B.1 Local evolution at R for ɛ = 0.01 and μ = 0.01 of normalized column density of the species for different pebble Stokes numbers (dotted lines: Stpbb = 1, solid lines (see Fig. 2): Stpbb = 0.1, dashed lines:Stpbb = 0.01) and different radii (colors). Each panel corresponds to a different particle species: primordial dust and pebbles, planetesimals, collisional dust and pebbles (from top to bottom).

In the text