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

TL;DR

This study uses a local toy model to explore how planetesimals form and evolve in protoplanetary disks, revealing the timing, mass distribution, and dependence on physical parameters, crucial for understanding planet formation.

Contribution

The paper introduces a simplified 0-D model to analyze the evolution of dust and planetesimal mass ratios, highlighting the rapid formation and significant mass of planetesimals in disks.

Findings

Planetesimals form quickly and are widespread in disks.

A significant portion of mass is in planetesimals after 10^4-10^6 years.

The timing and mass depend on distance, initial mass, and formation efficiency.

Abstract

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.

Tables

Table 1. Table 1: Default parameter set

Parameter	Value	Unit	Parameter name	References
$R$	$10$	AU	Distance to star	Arbitrarily set to Saturn’s position in Solar System
$M_{\mathrm{disk}}$	$0.02$	$M_{\mathrm{star}}$	Disk mass (gas + solids)	MMSN (Weidenschilling, 1977)
$ϵ$	$0.01$		Trap efficiency	Chosen arbitrarily
${\mathrm{St}}_{\mathrm{pbb}}$	$0.1$		Pebble Stokes number	Maximum Stokes number in Birnstiel et al. (2012)
$M_{\mathrm{star}}$	$1$	$M_{\odot }$	Stellar mass	Set to mass of Sun
$T_{\mathrm{star}}$	$4000$	K	Stellar effective temperature	Beckwith et al. (1990); Chiang & Goldreich (1997)
$R_{\mathrm{star}}$	$1.25$	$R_{\odot }$	Stellar radius	Fiducial value for pre-main sequence stars
${\alpha }_{\mathrm{irr}}$	$0.1$		Irradiation angle	E.g. Pfeil & Klahr (2019)
${ϵ}_{\mathrm{dg}}^0$	$0.01$		Initial dust-to-gas ratio	Savage & Jenkins (1972); Draine et al. (2007)
$R_C$	$40$	AU	Characteristic radius of initial dust profile	Set to Kuiper Belt’s location in Solar System
$n$	$1$		Power-law exponent of initial dust profile	After radial viscosity profile in Shakura & Sunyaev (1973)
$a_0$	${10}^{-4}$	cm	Dust grain radius	Mathis et al. (1977)
$r_{\mathrm{pls}}$	$50$	km	Planetesimal radius	Morbidelli et al. (2009); Klahr & Schreiber (2015)
$\rho$	$1.2$	$g{\mathrm{cm}}^{-3}$	Volume material density	Carry (2012)
$d$	$5$	$h_g$	Trap distance	Dittrich et al. (2013)
${\tau }_{\mathrm{trap}}$	$100$	Orbits	Trap lifetime	Dittrich et al. (2013); Manger & Klahr (2018)

Equations

$$\Sigma (R,\varphi )={\int }_{-\infty }^{\infty }\rho (R,\varphi ,z)\mathrm{d}z.$$

$$\frac{\partial {\Sigma }_{\mathrm{solids}}(R,t)}{\partial t}=0.$$

$${\Sigma }_{\mathrm{total}}={\Sigma }_{\mathrm{g}}+{\Sigma }_{\mathrm{solids}}={\Sigma }_{\mathrm{g}}+{\Sigma }_{\mathrm{dst}}^0=\left(1+{ϵ}_{\mathrm{dg}}^0\right){\Sigma }_{\mathrm{g}},$$

$${ϵ}_{\mathrm{dg}}^0=\frac{{\Sigma }_{\mathrm{dst}}^0}{{\Sigma }_{\mathrm{g}}}=10^{-2},$$

$${\Sigma }_{\mathrm{total}}(R)=C_{\Sigma }{\left(\frac{R}{R_{\mathrm{C}}}\right)}^{-n}\exp \left[-{\left(\frac{R}{R_{\mathrm{C}}}\right)}^{2-n}\right],$$

$$c_{\mathrm{s}}=\sqrt{\frac{k_{\mathrm{B}}T}{2m_{\mathrm{p}}}},$$

$$h_{\mathrm{g}}=\frac{c_{\mathrm{s}}}{\Omega },$$

$$\Omega =\sqrt{\frac{G{M}_{\mathrm{star}}}{R^3}},$$

$$\nu =\alpha c_{\mathrm{s}}{h}_{\mathrm{g}}\propto \frac{c_{\mathrm{s}}^2}{\Omega }\propto T\cdot R^{3/2}\propto R^{3/2-q}=R^1,$$

$$T=T_{\mathrm{star}}{\left(\frac{R_{\mathrm{star}}}{R}\right)}^q{\alpha }_{\mathrm{irr}}^{0.25},$$

$$a=a_0\exp \left(\frac{t}{{\tau }_{\mathrm{growth}}}\right),$$

$${\mathrm{St}}_{\mathrm{dst}}=\frac{\pi }{2}\frac{a_0\rho }{{\Sigma }_{\mathrm{g}}},$$

$${\tau }_{\mathrm{growth}}\simeq \ln \left(\frac{{\mathrm{St}}_{\mathrm{pbb}}}{{\mathrm{St}}_{\mathrm{dst}}}\right)\frac{1}{{ϵ}_{\mathrm{dg}}\Omega }.$$

$${\dot{\Sigma }}_{\mathrm{growth}}=\frac{{\Sigma }_{\mathrm{dst}}}{{\tau }_{\mathrm{growth}}}=\ln \left(\frac{{\mathrm{St}}_{\mathrm{dst}}}{{\mathrm{St}}_{\mathrm{pbb}}}\right){ϵ}_{\mathrm{dg}}\Omega {\Sigma }_{\mathrm{dst}}\propto {\Sigma }_{\mathrm{dst}}^2.$$

$${\dot{\Sigma }}_{\mathrm{dst}}=-{\dot{\Sigma }}_{\mathrm{growth}}$$

$${\Sigma }_{\mathrm{dst}}(t)={\left[C\cdot t+\frac{1}{{\Sigma }_{\mathrm{dst}}^0}\right]}^{-1}.$$

$$C\text{mathchar}58=\frac{\Omega }{{\Sigma }_{\mathrm{g}}\ln \left(\frac{{\mathrm{St}}_{\mathrm{pbb}}}{{\mathrm{St}}_{\mathrm{dst}}}\right)}$$

$$v_{\mathrm{drift}}=\frac{{\mathrm{St}}_{\mathrm{pbb}}}{{\mathrm{St}}_{\mathrm{pbb}}^2+{\left(1+{ϵ}_{\mathrm{dg}}\right)}^2}\frac{h_{\mathrm{g}}}{R}\gamma c_{\mathrm{s}},$$

$$\gamma =\frac{\partial \ln P}{\partial \ln R},$$

$$l\text{mathchar}58=\frac{d}{ϵ}.$$

$${\dot{\Sigma }}_{\mathrm{form}}=\frac{|v_{\mathrm{drift}}|}{l}{\Sigma }_{\mathrm{pbb}}.$$

$${\tau }_{\mathrm{sed}}\approx \frac{{\tau }_{\mathrm{ff}}}{\mathrm{St}},$$

$$M_{\mathrm{trapped}}=2\pi Rϵ{\tau }_{\mathrm{trap}}{\int }_{{\mathrm{St}}_{\min }}^{{\mathrm{St}}_{\max }}|v_{\mathrm{drift}}|{\Sigma }_{\mathrm{pbb},\mathrm{St}}\mathrm{dSt}\ge m_{\mathrm{pls}},$$

$$M_{\mathrm{trapped}}=2\pi Rϵ{\tau }_{\mathrm{trap}}|v_{\mathrm{drift}}|{\Sigma }_{\mathrm{pbb}}\ge m_{\mathrm{pls}}.$$

$$m_{\mathrm{pls}}=\frac{4\pi }{3}\rho r_{\mathrm{pls}}^3,$$

$$n_m(m)=C_n\cdot m^{-\xi },$$

$${\lambda }_{\mathrm{mfp}}=\frac{1}{n_j{\sigma }_{ij}}.$$

$${\tau }_{\mathrm{col},i}=\frac{{\lambda }_{\mathrm{mfp}}}{\Delta v_{ij}}=\frac{1}{n_j{\sigma }_{ij}\Delta v_{ij}},$$

$${\sigma }_{\mathrm{grav}}=4\pi r_{\mathrm{pls}}^2\left[1+{\left(\frac{v_{\mathrm{esc}}}{\Delta v_{ij}}\right)}^2\right],$$

$$v_{\mathrm{esc}}=\sqrt{\frac{2G{m}_{\mathrm{pls}}}{r_{\mathrm{pls}}}}.$$

Full text

*Aims. *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.

*Methods. *We perform a parameter study based on a model considering dust growth, planetesimal formation and collisional fragmentation of planetesimals, while neglecting radial transport processes.

*Results. *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$. Planetesimal collisions are more significant in more massive disks, leading to lower relative planetesimal fractions compared to less massive disks. After approximately $10^{6}$ yr, our model predicts planetesimal collisions to dominate, which resupplies small particles.

*Conclusions. *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.

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 for 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 to determine the amount of mm- and $\mathrm{\SIUnitSymbolMicro}$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 are too small to be be detected individually. Models of planetesimal formation (e.g. Dr\każ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 perform 0-d local simulations to study how the planetesimal to dust and pebble ratio is evolving in 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 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. In addition, we there discuss our numerical setup, i.e. 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_{\mathrm{star}}$. We use the concept of column densities, which is commonly used in astrophysics when describing accretion disks. In cylindrical coordinates ($R,\phi,z$) it is defined as the integral of the volumetric mass density $\rho$ of gas or solids integrated along the vertical path $z$ that goes through the disk

[TABLE]

By assuming a cylindrical symmetric disk, one can eliminate the $\phi$-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 $\Sigma_{\mathrm{solids}}$ must be conserved

[TABLE]

Initially, the total column density of solids is made up purely by dust. We relate the dust column density $\Sigma_{\mathrm{dst}}$ to the gas column density $\Sigma_{\mathrm{g}}$ with the dust-to-gas ratio $\epsilon_{\mathrm{dg}}$, which we assume to be $R$-independent. The total column density of all material $\Sigma_{\mathrm{total}}$, including solids and gas, is then given by

[TABLE]

where we introduced the initial dust-to-gas ratio

[TABLE]

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 $\Sigma_{\mathrm{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),

[TABLE]

where $R_{\mathrm{C}}$ is the so-called characteristic radius, marking the transition of the power law to the exponential. $C_{\Sigma}$ is a normalization constant chosen such that the disk contains a total mass of solids and gas of $M_{\mathrm{disk}}$. For our simulations we choose $R_{\mathrm{C}}=40\mathrm{\ AU}$, roughly corresponding to Kuiper Belt’s location in the Solar System.

For a vertically isothermal disk with temperature $T$, the sound speed $c_{\mathrm{s}}$ of molecular hydrogen is given by

[TABLE]

with Boltzmann constant $k_{\mathrm{B}}=1.3807\cdot 10^{-16}\mathrm{erg}\,\mathrm{K}^{-1}$, and proton mass $m_{\mathrm{p}}=1.673\cdot 10^{-24}\mathrm{\ g}$. Further, the vertical pressure scale height of the gas disk $h_{\mathrm{g}}$ is defined as

[TABLE]

which originates from considering vertical hydrostatic balance (v. Weizsäcker, 1948). Here, we introduced the Keplerian angular velocity $\Omega$, which is derived from balance of centrifugal and gravitational acceleration

[TABLE]

where $G\approx 6.674\cdot 10^{-8}\ \mathrm{cm}^{3}\,\mathrm{g}^{-1}\,\mathrm{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):

[TABLE]

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

[TABLE]

where we choose stellar properties of $R_{\mathrm{star}}=1.25\cdot R_{\odot}\approx 1.25\cdot 7\cdot 10^{10}\mathrm{\ cm}$ and $T_{\mathrm{star}}=4000\mathrm{\ K}$ (Beckwith et al., 1990), $q=0.5$ after Chiang & Goldreich (1997) and an irradiation angle of $\alpha_{\mathrm{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 $\alpha$ quantifies disk turbulence. Typical values are in the range $10^{-4}\leq\alpha\leq 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); 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 $\mathrm{St}\gg 1$, they are decoupled from the gas, while for $\mathrm{St}\ll 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

[TABLE]

where $a_{0}$ is an initial reference size and $\tau_{\mathrm{growth}}$ a certain growth timescale. Equation (11) is only valid for constant $\epsilon_{\mathrm{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\text{\,}\mathrm{\SIUnitSymbolMicro m}$ after Mathis et al. (1977) and the pebble species to $a=a(\mathrm{St}_{\mathrm{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 $\mathrm{St}_{\mathrm{dst}}$ via

[TABLE]

where $\rho$ is the mass density of dust. We choose $\rho=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 $\mathrm{St}_{\mathrm{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

[TABLE]

The approximation in Eq. (13) holds if relative dust velocities are set by disk turbulence $\Delta v_{\mathrm{dst}}\approx c_{s}\sqrt{3\alpha\mathrm{St}_{\mathrm{dst}}}$ with $\mathrm{St}_{\mathrm{dst}}\ll 1$, i.e. for sufficiently large $\alpha$-values.

The dust to pebble growth rate is calculated via

[TABLE]

The differential equation

[TABLE]

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

[TABLE]

Here, $\Sigma_{\mathrm{dst}}^{0}$ is the initial dust column density from the profile in Eq. (5) and

[TABLE]

is constant in time for a fixed $\mathrm{St}_{\mathrm{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 towards 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

[TABLE]

where we introduced the exponent of the gas pressure power law

[TABLE]

which is for $R\ll R_{\mathrm{C}}$ and $T\propto R^{-1/2}$ given by $\gamma=-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 $\epsilon_{\mathrm{dg}}$ changes in time.

One finds that pebbles with $\mathrm{St}_{\mathrm{pbb}}\approx 10^{-1}$ quickly drift inward towards 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 towards 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) and Johansen et al. (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 $\epsilon$. We will also refer to $\epsilon$ 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 $\tau_{\mathrm{trap}}$ a trap’s typical lifetime. We set $d=5h_{\mathrm{g}}$ and $\tau_{\mathrm{trap}}$ to 100 orbits, as found numerically by Dittrich et al. (2013); Manger & Klahr (2018).

Following Lenz et al. (2019), pebbles are transformed into planetesimals over the conversion length

[TABLE]

We further express the planetesimal formation rate by

[TABLE]

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 $\tau_{\mathrm{sed}}$ for $\mathrm{St}<1$ particles increases due to aerodynamic friction on the order of

[TABLE]

where $\tau_{\mathrm{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_{\mathrm{pls}}$ (Lenz et al., 2019):

[TABLE]

where $\Sigma_{\mathrm{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

[TABLE]

If Eq. (24) is not fulfilled, $\dot{\Sigma}_{\mathrm{form}}$ is set to 0. We deduce the fixed parameter $m_{\mathrm{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 $\rho=1.2\mathrm{\ g}\,\mathrm{cm}^{-3}$, the planetesimal mass is approximated with

[TABLE]

where $r_{\mathrm{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

[TABLE]

where $C_{n}$ is a normalization constant and $\xi=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, $\xi=1.83$ is robust against variations of physical parameters like the impact velocity (Dohnanyi, 1969), and our collision model is robust against variations of $\xi$, 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 $\int_{m_{i}}^{m_{j}}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}\leq m<m_{1}$, for pebbles $m_{1}\leq m<m_{2}$ and for planetesimals $m_{2}\leq m\leq M$. The derivations of $m_{1}$ and $m_{2}$ are shown in Appendix C.2.

With these requisites, one can solve the integral $\int_{m_{i}}^{m_{j}}n_{m}mdm$ for each population, which leads to three fractions $p_{\mathrm{dst}},p_{\mathrm{pbb}}$ and $p_{\mathrm{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 tab. 1), $m_{1},m_{2}\ll m_{\mathrm{pls}}$. With $\xi=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 $\lambda_{\mathrm{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 $\sigma_{ij}$, then (Birnstiel et al., 2016)

[TABLE]

The collision timescale $\tau_{\mathrm{col},i}$ is defined as the average time during which particle $i$ experiences one encounter

[TABLE]

where $\Delta 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_{\mathrm{pls}}$, is given by the gravitational cross section $\sigma_{\mathrm{grav}}$ (Safronov, 1969), which considers both the geometry of the system but also gravitational focusing:

[TABLE]

where the escape velocity $v_{\mathrm{esc}}$ of the planetesimals is

[TABLE]

By approximating the vertical planetesimal distribution with a Gaussian with a root mean square width of $h_{\mathrm{pls}}=v_{\mathrm{hill}}\Omega^{-1}$ (Goldreich et al., 2004), the number density of planetesimals can be written as

[TABLE]

Here, we introduced the Hill velocity $v_{\mathrm{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_{\mathrm{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)

[TABLE]

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

[TABLE]

We follow Morbidelli et al. (2009) and insert the Hill velocity for the relative velocity of planetesimals $\Delta v_{ij}$. We note that the relative velocity appears once directly in (28) and again as its inverse squared in (29), leading to the perhaps counter intuitive implication of higher relative velocities causing longer collision timescales. The rate of planetesimal collisions $\dot{\Sigma}_{\mathrm{col}}$ must also depend on the available planetesimal surface density. It is given by

[TABLE]

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:

[TABLE]

Since $\dot{\Sigma}_{\mathrm{growth}}$ and $\dot{\Sigma}_{\mathrm{form}}$ each appear once as a sink and a source term and $p_{\mathrm{dst}}+p_{\mathrm{pbb}}+p_{\mathrm{pls}}=1$ per definition, the total column density in the system is conserved in time

[TABLE]

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

[TABLE]

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 (35a) - (35c), collisional dust is too compact and dense to grow further, whereas collisional pebbles may participate in planetesimal formation again.

We will denote the configuration depicted in equations (35a), (35b) and (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 (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

[TABLE]

according to the radial column density profile in Eq. (5). All other populations are set to zero, i.e.

[TABLE]

Our simulation utilizes parameters introduced in Sect. 2 and summarized in Table 1. In our parametric study, we vary $R$, $M_{\mathrm{disk}}$ and $\epsilon$. 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_{\mathrm{disk}}=0.02\cdot M_{\odot}$ 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 $\epsilon=0.01$. A brief discussion of the influence of $\mathrm{St}_{\mathrm{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.

4 Results

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.

4.1 Local evolution

Figure 2 displays the local evolution of the normalized column densities $\Sigma_{\mathrm{dst,prim}},\Sigma_{\mathrm{dst,col}},\Sigma_{\mathrm{pbb,prim}},\Sigma_{\mathrm{pbb,col}}$ and $\Sigma_{\mathrm{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\cdot 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 equation (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\cdot 10^{3}$ to approximately $2\cdot 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 $\epsilon$. The planetesimal population comprises over 98 % of the total available column density, i.e. the mass available to the system starting around $5\cdot 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 equation (37e) only consists of a source term.

We note that until $\sim 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.

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\mathrm{AU}\leq R\leq 100\mathrm{AU}$ all using $\epsilon=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\cdot 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 towards 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.

4.2.2 Variation of the efficiency of planetesimal formation

Further, we investigate the influence of the trap efficiency $\epsilon$ on the column density evolution. All other parameters are kept constant according to the default parameter set. The trap efficiency is per definition $0\leq\epsilon\leq 1$, with $\epsilon=1$ implying that the entire incoming pebble flux is trapped and converted into planetesimals and $\epsilon=0$ meaning that the trap is not capturing and forming planetesimals.

A comparison of the local column density evolution for different $\epsilon$-values ($\epsilon$ = 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 $\epsilon$ 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 $\epsilon$. During the initial phase of planetesimal formation (until $10^{4}$ to $10^{6}$ yr depending on the value of $\epsilon$) the planetesimal column density is linear in $\epsilon$ 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 $\epsilon$-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 $\epsilon$. Here, planetesimal formation is insignificant, because there is only few pebbles available. We point out, that for all $\epsilon$-values almost all the mass ends up within the planetesimal population — for $\epsilon=0.5,0.1,0.01,0.001$ starting at approximately $2\cdot 10^{4},5\cdot 10^{4},6\cdot 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 $\epsilon$. In the middle panel, one can see an equivalent plot for primordial pebbles (solid lines) and collisional pebbles (dashed lines). The influence of $\epsilon$ on the evolution of primordial dust is negligible, since the dust growth timescale does not directly depend on $\epsilon$, 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 $\epsilon$ 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) $\dot{\Sigma}_{\mathrm{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 $\epsilon$, because $\epsilon$ is what determines when enough planetesimals can be formed for collisions to produce noteworthy amounts of fragments (see top panel of Fig. 4).

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_{\mathrm{disk}}=0.004-0.143M_{\odot}$ in the $\sim$ 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\cdot 10^{-4}M_{\odot}$ while gas has a total mass of $1.3\cdot 10^{-2}M_{\odot}$. 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

[TABLE]

and investigate values ranging from $\mu=0.01$, approximately corresponding to the MMSN after Weidenschilling (1977) up to $\mu=0.07$. More massive disks are expected to be gravitationally unstable and are additionally uncommon when taking recent observations by Andrews et al. (2009) and Andrews 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 $\sim 5\cdot 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\cdot 10^{5}$ yr). For MMSN-like disks, the peak height is not only lower ($\sim$ 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\cdot 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 $\alpha$, 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 $\mathrm{St}_{\mathrm{pbb}}=0.1$ for our parameter study. We discuss implications in B. A more complex model would introduce a dependence on $\alpha$ 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 $\epsilon_{\mathrm{dg}}^{0}$ and $r_{\mathrm{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_{1}$ and $R_{N}$ at a certain time $t$ is given by integrating the column density at $t$ over $r$

[TABLE]

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,

[TABLE]

We set $N=20$, $R_{1}=1\mathrm{\ AU}$, and $R_{20}=100\mathrm{\ 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 ($\epsilon\in\{0.001,0.01,0.1,0.5\}$) and disk mass $\mu=M_{\mathrm{disk}}/M_{\odot}$ ($\mu\in\{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 $\mu$ 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 ( $\epsilon_{\mathrm{dg}}M_{\mathrm{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 $\epsilon$-values for a constant disk mass of $\mu=0.02$, one discovers, that $\epsilon$ both influences the time where the maximal planetesimal mass is reached and its peak value. Similar to Fig. 4, a higher $\epsilon$-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 $\epsilon=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_{\mathrm{Earth}}$. A tenfold increase of $\epsilon$ 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 $\epsilon$ remains constant at $\epsilon=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_{\mathrm{Earth}}$ similar to the top right panel, implying a doubled normalized disk mass $\mu$ would be approximately balanced by a tenfold decrease of $\epsilon$. 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 $\epsilon=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 $\epsilon$ (left panels) and the dimensionless disk mass $\mu$ (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 $\epsilon$ and $\mu$ for times $\gtrsim 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 $\epsilon$ and $\mu$ as seen in the top panels of Fig. 7. Here, while the mass in planetesimals always increases with $\epsilon$, the relation is clearly non-linear. This is explained due to the fact, that the pebble supply is limited and therefore the maximum planetesimal column density is independent of $\epsilon$, 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 $\mu$, because the initial growth phase from dust to pebbles only depends on the dust-to-gas ratio, which is picked the same for all $\mu$. The small differences originate from the fact that the significance of planetesimal collisions depends on total mass. However, the pebble fraction correlates strongly with $\epsilon$ (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\cdot 10^{6}$ yr, the mass in pebbles is rather independent of $\epsilon$. This is due to planetesimal formation being insignifcant in comparison to planetesimal collisions, which do not depend on $\epsilon$.

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 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. $\mathrm{St}\gg\alpha$. 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 $\mathrm{St}_{\mathrm{pbb}}=0.1$ during the entire disk evolution. The influence $\mathrm{St}_{\mathrm{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 modelled 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 et al. (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 e.g. 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 H2O would correspond to approximately $T_{\mathrm{g}}\gtrsim$ 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 $\sim$ 100 km, the efficiency of pebble accretion is minimal. Therefore, our assumption of a fixed planetesimal size of $\sim$ 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 $\epsilon$. 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 $\epsilon$ 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 ($\sim 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.

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 Eqns. (37a) to (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 $\mathrm{St}_{\mathrm{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. 2. 2.

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 $\sim 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. 3. 3.

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 $\epsilon$ and disk mass $M_{\mathrm{disk}}$. Note again, that $\epsilon$ 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 $\epsilon$-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 in 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 Fig. 6 and Fig. 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 $\epsilon$ and $M_{\mathrm{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 $\epsilon$ as seen in Fig. 6 and Fig. 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 $\epsilon$ and disk mass $M_{\mathrm{disk}}$.

Despite many simplifications made in our model, it generated results connecting important model parameters like trap efficiency $\epsilon$ and disk mass $M_{\mathrm{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.

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. 8, 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 for growth 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. 8. 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

[TABLE]

with $\mathrm{St}_{\mathrm{frag}}$ being the maximum Stokes number in a fragmentation dominated size distribution, and $\mathrm{St}_{\mathrm{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 $\mathrm{St}_{\mathrm{pbb}}=0.1$. This is a very good approximation for $\alpha=10^{-3}$ (see Lenz et al., 2019, Fig. 3). In more turbulent disks ($\alpha\geq 10^{-2}$) the higher relative velocities of dust particles lead to a decrease of the fragmentation size and our approximation of $\mathrm{St}_{\mathrm{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. 9.

For overestimates of the pebble Stokes number, the dust-to-pebbles growth timescale will also be overestimated, since $\tau_{\mathrm{growth}}\propto\ln{\left(\mathrm{St}_{\mathrm{pbb}}\right)}$ after Eq. (13). Because of this, the pebble population in Fig. 9 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.

[TABLE]

where $\sigma_{\mathrm{pls}}$ is a parameter specifying the surface density of a single planetesimal, and $m_{0}$ and $M$ the masses of the smallest and largest possible fragment respectively. We set $m_{0}=4/3\pi a_{\mathrm{0}}^{3}$ and $M=m_{\mathrm{pls}}$ for destructive collisions. We solve the integral in Eq. (44), giving

[TABLE]

leading to

[TABLE]

We note that the choice of $\sigma_{\mathrm{pls}}$ is irrelevant, because it cancels when calculating the relative fractions $p_{\mathrm{dst}},p_{\mathrm{pbb}}$ and $p_{\mathrm{pls}}$.

C.2 Calculation of the transition masses $m_{1}$ and $m_{2}$

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

[TABLE]

If this condition is met, a particle with mass $m_{1}$ 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

[TABLE]

where we used $(\epsilon_{\mathrm{dg}}+1)^{2}\approx 1$.

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.

[TABLE]

Particles more massive than $m_{2}$ 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_{2}$ to the planetesimal population. The condition in (49) also results in a quadratic equation after using Eqns. (12), (13) and (18). We then define the larger solution as the transition of pebbles to planetesimals in the fragment distribution:

[TABLE]