Planetesimal Population Synthesis: Pebble Flux Regulated Planetesimal Formation

TL;DR

This paper introduces a model linking local planetesimal formation rates to pebble flux, exploring how disk turbulence influences planetesimal formation timing and distribution, with implications for planetary embryo development.

Contribution

It presents a new formulation for planetesimal formation driven by pebble flux and investigates the effects of disk turbulence and pebble flux thresholds on planetesimal distribution.

Findings

Lower turbulence (α≈0.001) favors rapid planetesimal formation across the disk.

Planetesimals form after dust grows to pebbles and pebble flux exceeds a critical value.

Final planetesimal distribution is steeper than initial dust and gas profiles.

Abstract

We propose an expression for a local planetesimal formation rate proportional to the instantaneous radial pebble flux. The result --- a radial planetesimal distribution --- can be used as initial condition to study the formation of planetary embryos. We follow the idea that one needs particle traps to locally enhance the dust-to-gas ratio sufficiently such that particle gas interactions can no longer prevent planetesimal formation on small scales. The location of these traps can emerge everywhere in the disk. Their occurrence and lifetime is subject of ongoing research, thus they are implemented via free parameters. This enables us to study the influence of the disk properties on the formation of planetesimals, predicting their time dependent formation rates and location of primary pebble accretion. We show that large $\alpha$-values of $0.01$ (strong turbulence) prevent the formation of planetesimals in the inner part of the disk, arguing for lower values of around $0.001$ (moderate turbulence), at which planetesimals form quickly at all places where they are needed for proto-planets. Planetesimals form as soon as dust has grown to pebbles ($\sim\mathrm{mm}$ to $\mathrm{dm}$) and the pebble flux reaches a critical value, which is after a few thousand years at $2-3\,$AU and after a few hundred thousand years at $20-30\,$AU. Planetesimal formation lasts until the pebble supply has decreased below a critical value. The final spatial planetesimal distribution is steeper compared to the initial dust and gas distribution which helps to explain the discrepancy between the minimum mass solar nebula and viscous accretion disks.

Tables

Table 1. Table 1 Parameters which are kept constant in this paper. Values of the disk mass and of r c subscript 𝑟 c {r_{\mathrm{c}}} can be compared with the statistics in Fig. 3 of Andrews et al. ( 2010 ) . Per size decade we use 27 grid points.

Symbol	Value	Meaning
$M_{\star }\left(t_0\right)$	$1M_{\odot }$	Initial mass of the central star
$M_{\mathrm{disk}}$	$0.05M_{\odot }$	Total initial disk mass (gas+particles)
$r_c$	$35\mathrm{au}$	Transition radius from power-law to exponential drop off for ${\Sigma }_g$, Eq. (20)
$d$	$5h_g$	Pebble trap distance, Eq. (1)
${\tau }_l$	$100t_{\mathrm{orb}}$	Lifetime of pebble traps
$f_{\mathrm{ice}}$	$1/3$ for $T>170K$	Dust ice parameter, Eq. (26)
$Z\left(t_0\right)$	${10}^{-2}$	Initial dust-to-gas column density ratio, ${\Sigma }_d/{\Sigma }_g$
$\gamma$	$1$	Power-law index in gas density profile, Eq. (20)
$q$	$0.4$	Power-law index of temperature profile, Eq. (25)
$\xi$	$1.83$	Power-law index for fragment distribution, Eq. (40)
$v_f$	$1-10m/s$	Fragmentation threshold speed, Eq. (18)
$a_0$	$0.1m$	Smallest grain radius for all times
$a_1$	$1m$	Largest grain radius for initial condition
$n$	$\propto a^{-3.5}$	Initial number density distribution, Eq. (24)
${\rho }_s$	$1.2g/{\mathrm{cm}}^3$	Internal density of all solids (dust, pebbles, planetesimals)
$m_p$	$1.05\cdot {10}^{-7}{M}_{\oplus }$	Planetesimal mass of $100\mathrm{km}$ diameter planetesimals, Eq. (4)
${\mathrm{St}}_{\min }$	${10}^{-2}$	Minimum Stokes number participating in trapping & planetesimal formation
${\mathrm{St}}_{\max }$	$10$	Maximum Stokes number participating in trapping & planetesimal formation

Table 2. Table 2 Overview of different approaches leading to streaming instability (SI) or planetesimals (ptes). With α t = 0 subscript 𝛼 t 0 \alpha_{\mathrm{t}}=0 we mean that there is no turbulence active, or only turbulence stemming from streaming instability itself. The second quantity is the sub-Keplerianness, defined as η : = 1 / 2 ⋅ ( h g / r ) 2 ⋅ d ln ρ g / d ln r \eta\mathrel{\mathop{\mathchar 58\relax}}=1/2\cdot(h_{\mathrm{g}}/r)^{2}\cdot\operatorname{d\!}{\ln{{\rho_{\mathrm{g}}}}}/\operatorname{d\!}{\ln{r}} . In the third column with "feedback" we mean particle-gas coupling such that momentum is also transferred to the gas if particles loose momentum due to the interaction with the gas. The velocity dispersion of particles is δ ​ v 𝛿 𝑣 {\delta v} , Ω Ω \Omega is the Kepler frequency, and G 𝐺 G the gravitational constant. For larger particles δ ​ v 𝛿 𝑣 {\delta v} decreases, allowing gravitational instability for large particles (Weidenschilling, 1995 ) . The question mark in the third row indicates that one cannot be entirely sure that streaming instability is active.

input				result
${\alpha }_t$	$\eta$	feedback	$Z={\Sigma }_d/{\Sigma }_g$	SI	ptes	references
$=0$	$=0$	no	$>\Omega \delta v/\left(\pi G{\Sigma }_g\right)$	no	yes	Safronov (1969); Goldreich & Ward (1973); Youdin & Shu (2002)
$>0$	$>0$	no	0.01	no	yes	Johansen et al. (2006)
$>0$	$>0$	yes	0.01	?	yes	Johansen et al. (2007)
$=0$	$>0$	yes	0.01	yes	no	Johansen et al. (2007)
$=0$	$>0$	yes	$\gtrsim 0.03$	yes	yes	Johansen et al. (2009b); Simon et al. (2016); Abod et al. (2018, $Z=0.1$, $\eta$ varies)

Equations

$${\dot{\Sigma }}_{\mathrm{p}}(r)=\frac{\epsilon }{d(r)}\frac{{\dot{M}}_{\mathrm{peb}}}{2\pi r},$$

$${\dot{M}}_{\mathrm{peb}}\text{mathchar}58=2\pi r\sum _{{\mathrm{St}}_{\min }\le \mathrm{St}\le {\mathrm{St}}_{\max }}\left|v_{\mathrm{drift}}(r,\mathrm{St})\right|{\Sigma }_{\mathrm{d}}(r,\mathrm{St}).$$

$$\ell \text{mathchar}58=d/\epsilon .$$

$$m_{\mathrm{p}}=\frac{4\pi }{3}(50\mathrm{km}{)}^3{\rho }_{\mathrm{s}}\approx 6.28\cdot 10^{20}\mathrm{g}\approx 1.05\cdot 10^{-7}{M}_{\oplus }$$

$$\epsilon {\tau }_{\mathrm{l}}{\dot{M}}_{\mathrm{peb}}>m_{\mathrm{p}}.$$

$${\dot{M}}_{\mathrm{cr}}\text{mathchar}58=\frac{m_{\mathrm{p}}}{\epsilon {\tau }_{\mathrm{l}}}$$

$$\nu ={\alpha }_{\mathrm{t}}{c}_{\mathrm{s}}{h}_{\mathrm{g}}.$$

$${\rho }_i(r,z)={\rho }_{i,0}(r)\cdot \exp \left[-\frac{1}{2}{\left(\frac{z}{h_i}\right)}^2\right],$$

$${\Sigma }_i(r,t)\text{mathchar}58={\int }_{-\infty }^{\infty }{\rho }_i(r,z,t)dz.$$

$${\rho }_{i,0}(r)=\frac{1}{\sqrt{2\pi }}\frac{{\Sigma }_i(r)}{h_i}.$$

$$h_{\mathrm{g}}=c_{\mathrm{s}}/\Omega$$

$$h_{\mathrm{d}}={\left(1+\frac{\mathrm{St}}{{\alpha }_{\mathrm{t}}}\right)}^{-1/2}{h}_{\mathrm{g}}.$$

$$\mathrm{St}=\frac{{\rho }_{\mathrm{s}}}{{\rho }_{\mathrm{g}}}\frac{a}{v_{\mathrm{th}}}\Omega \stackrel{z=0}{\approx }\frac{\pi }{2}\frac{a{\rho }_{\mathrm{s}}}{{\Sigma }_{\mathrm{g}}}.$$

$$u_{\mathrm{r}}=-\frac{3}{{\Sigma }_{\mathrm{g}}\sqrt{r}}\frac{\partial }{\partial r}({\Sigma }_{\mathrm{g}}\nu \sqrt{r}).$$

$$v_{\mathrm{drift}}=\frac{\mathrm{St}}{{\mathrm{St}}^2+1}\frac{h_{\mathrm{g}}}{r}\frac{\partial \ln P}{\partial \ln r}{c}_{\mathrm{s}},$$

$$v_{\mathrm{r}}=\frac{u_{\mathrm{r}}}{1+{\mathrm{St}}^2}+v_{\mathrm{drift}}.$$

$$p_{\mathrm{f}}\text{mathchar}58=\{\begin{array}{ll}0.5\cdot \exp \left[2(\Delta v-v_{\mathrm{f}})/(cm/s)\right]& \Delta v\le v_{\mathrm{f}}\\ 1-0.5\cdot \exp \left[-2(\Delta v-v_{\mathrm{f}})/(cm/s)\right]& \Delta v>v_{\mathrm{f}}\end{array},$$

$$\frac{v_{\mathrm{f}}(r)}{m/s}=\{\begin{array}{ll}10^{0.5\cdot \exp \left[4(r-r_{\mathrm{ice}})/\mathrm{au}\right]}& r\le r_{\mathrm{ice}}\\ 10^{1-0.5\cdot \exp \left[-4(r-r_{\mathrm{ice}})/\mathrm{au}\right]}& r>r_{\mathrm{ice}}\end{array}.$$

$$D_{\mathrm{d}}=\frac{D_{\mathrm{g}}}{1+{\mathrm{St}}^2}$$

$${\Sigma }_{\mathrm{g}}(r)={\Sigma }_{\mathrm{c}}{\left(\frac{r}{r_{\mathrm{c}}}\right)}^{-\gamma }\exp \left[-{\left(\frac{r}{r_{\mathrm{c}}}\right)}^{2-\gamma }\right].$$

$${\Sigma }_{\mathrm{c}}=(2-\gamma )M_{\mathrm{disk}}/[2\pi {r_{\mathrm{c}}}^2\cdot (1+Z(t_0))].$$

$${\Sigma }_{\mathrm{d}}(t_0)=Z(t_0){\Sigma }_{\mathrm{g}}(t_0).$$

$$\begin{array}{rl}\mathrm{St}& \le \frac{6\left|2-q-\gamma -(2-\gamma ){\left(\frac{r}{r_{\mathrm{c}}}\right)}^{2-\gamma }\right|}{2\gamma +3+q+2(2-\gamma ){\left(\frac{r}{r_{\mathrm{c}}}\right)}^{2-\gamma }}{\alpha }_{\mathrm{t}}\\ & \stackrel{\gamma =1}{\underset{q=0.4}{=}}\frac{6\left|0.6-r/r_{\mathrm{c}}\right|}{5.4+2r/r_{\mathrm{c}}}{\alpha }_{\mathrm{t}}\end{array}$$

$$n(a)\propto a^{-3.5}$$

$$T(r)=158.9\mathrm{K}\cdot {\left(\frac{r}{\mathrm{au}}\right)}^{-2/5}.$$

$$f_{\mathrm{ice}}=1/3\cdot \theta (T-170\mathrm{K})+\theta (170\mathrm{K}-T).$$

$$\begin{array}{rl}& \frac{\partial {\Sigma }_{\mathrm{d}}^i}{\partial t}+\frac{1}{r}\frac{\partial }{\partial r}\left\{r\left[{\Sigma }_{\mathrm{d}}^i{v}_{\mathrm{r}}^i-D_{\mathrm{d}}^i\frac{\partial }{\partial r}\left(\frac{{\Sigma }_{\mathrm{d}}^i}{{\Sigma }_{\mathrm{g}}}\right){\Sigma }_{\mathrm{g}}\right]\right\}\\ & =-\frac{\epsilon }{d}\left|v_{\mathrm{drift}}^i\right|{\Sigma }_{\mathrm{d}}^i\cdot \theta ({\dot{M}}_{\mathrm{peb}}-{\dot{M}}_{\mathrm{cr}})\\ & \theta ({\mathrm{St}}_i-{\mathrm{St}}_{\min })\cdot \theta ({\mathrm{St}}_{\max }-{\mathrm{St}}_i).\end{array}$$

$$\begin{array}{rl}{L}_i& =-\frac{\epsilon }{d}\left|v_{\mathrm{drift}}^i\right|\cdot \theta ({\dot{M}}_{\mathrm{peb}}-{\dot{M}}_{\mathrm{cr}})\\ & \theta ({\mathrm{St}}_i-{\mathrm{St}}_{\min })\theta ({\mathrm{St}}_{\max }-{\mathrm{St}}_i)\end{array}$$

$$d{\Sigma }_{\mathrm{p}}={\dot{\Sigma }}_{\mathrm{p}}dt=-\sum _{{\mathrm{St}}_{\min }\le {\mathrm{St}}_i\le {\mathrm{St}}_{\max }}{L}^{i}dt.$$

$${\mathrm{St}}_{\mathrm{frag}}=\frac{1}{3{\alpha }_{\mathrm{t}}}{\left(\frac{v_{\mathrm{f}}}{c_{\mathrm{s}}}\right)}^2$$

Full text

1. Introduction

The original ideas about planet formation—valid for terrestrial planets as well as gas- and ice-giants—assume that pebble-sized material, which has grown from tiny dust and ice grains via successive collisions and sedimented towards the midplane of the disk, directly collapses into $\sim 10-100\,\mathrm{km}$ sized planetesimals (Safronov, 1969; Goldreich & Ward, 1973). By definition planetesimals are bodies massive enough to be bound by, and able to accrete further solid material via, gravitational attraction. A critical size or mass for boulders to be called planetesimal can be estimated by requiring that pebbles settle toward the planetesimal with a higher velocity compared with the typical headwind the planetesimal experiences in the gaseous environment (Ormel & Klahr, 2010). The precise mass depends on the location in the disk and the size of the pebbles. The fact that their binding force is gravity makes it plausible gravity could also be involved in their formation in the first place. Nevertheless, Weidenschilling (1995) showed that Safronov (1969) and Goldreich & Ward (1973) neglected the effect of friction between particles and gas, and the associated momentum transfer from the solids to the gas. Sedimentation to the midplane leads to dust-to-gas ratios of the order of unity before the critical density for self-gravity among the pebbles is reached. But starting from dust-to-gas ratios of unity the drag from the particles accelerates the gas disk, which is usually pressure supported and sub-Keplerian, to the Keplerian speed.

The feedback from the dust onto the gas modifies the average azimuthal velocity of the gas which leads to resonant drag instabilities (Squire & Hopkins, 2018a, b), for example the streaming instability (Youdin & Goodman, 2005). The resulting turbulence will then prevent further vertical sedimentation and, ultimately, gravitational collapse of pebble clouds to planetesimals, because the necessary critical density in the pebble cloud cannot be reached. Originally this effect was identified as a Kelvin-Helmholtz instability (Sekiya, 1998; Johansen et al., 2006), yet in modern 3D simulations (Johansen et al., 2009b) one finds streaming instability to dominate the dynamical evolution.

After this important discovery, research was focused on the growth from grains to planetesimals via sticking collisions. Yet this model was also facing several problems like the drift and fragmentation barriers, because either the drift time becomes shorter than the growth time (drift limit) (Klahr & Bodenheimer, 2006; Birnstiel et al., 2012) or the turbulence in the nebula induces collision speeds among the solids that lead to bouncing (Güttler et al., 2010; Zsom et al., 2010) and fragmentation (Blum & Münch, 1993; Blum & Wurm, 2008; Gundlach & Blum, 2014). Aggregates charging can also act as a growth barrier (Okuzumi, 2009). Neither bouncing, nor charging is considered in this paper. In the early stages of circumstellar disks, micrometer-sized dust grows to pebble-size (typically $\sim\mathrm{mm-dm}$) via collisions. Particle growth models have been developed and improved in the past (Zvyagina et al., 1974; Weidenschilling, 1980; Dullemond & Dominik, 2005; Brauer et al., 2008a; Birnstiel et al., 2010; Windmark et al., 2012; Krijt et al., 2016a, b; Stammler et al., 2017). Birnstiel et al. (2012) have shown that the radial dynamics are mostly determined by particles with sizes close to the growth barriers.

Following the path of fluffy particle growth (Okuzumi et al., 2012; Kataoka et al., 2013), it becomes feasible to form planetesimals, at least within $10\,\mathrm{au}$ in the disks. Yet, Blum et al. (2017) have shown that the resulting planetesimals would have too large a tensile strength to explain the properties of comets, which are believed to be remnants from the planetesimal formation episode in the solar nebula.

The idea of Safronov (1969) and Goldreich & Ward (1973) was extended by Youdin & Shu (2002), showing that this gravitational instability can develop despite Kelvin-Helmholtz instability for sufficiently high dust-to-gas column density ratios at least in a laminar gas disk.

The work by Johansen et al. (2006) started a revival of gravity assisted planetesimal formation, because it was the first work that incorporated turbulence generated via the magneto rotational instability [MRI] (Balbus & Hawley, 1991). They found that turbulence not only leads to diffusion and prevents sedimentation, but also causes solid particles to get trapped and concentrated in non-laminar flow features reaching locally a critical Hill density to undergo collapse. These were created as by-products of the magneto rotational instability. They were able to show that the over-densities are in fact sufficient to trigger gravitational collapse.

This idea is different from the gravitational collapse envisioned by Cuzzi et al. (2008), in which grains get stochasticaly concentrated locally between vortices for a short time. The vortices there are part of a Kolmogoroff cascade and only live for a short time.

In the paper by Johansen et al. (2007) they found that a zonal flow feature arose in a box simulation of the MRI simulation, in which sufficiently many particles got trapped to trigger a gravitational collapse, without being diffused by the streaming instability. In test simulations in the same paper without MRI and thus no zonal flows, the streaming instability prevented the formation of planetesimals. These simulations of planetesimal formation in zonal flows were expanded to higher resolutions in order to study the size of forming planetesimals in Johansen et al. (2012), and the properties of zonal flows were further studied in Dittrich et al. (2013) and Bai & Stone (2014). Zonal flows appear to be transient features in any accretion disk subject to magnetic instabilities (Johansen et al., 2009a). In the simulations of Dittrich et al. (2013) they found a typical separation between the zonal flows of 5 pressure scale heights and a lifetime of about 100 local orbits. Whereas the above mentioned simulations were ideal magneto-hydrodynamic (MHD) runs, non ideal MHD simulations, especially those including the Hall-term, also showed the phenomenon of creating zonal flows (Bai & Stone, 2014; Béthune et al., 2016).

In situations where magnetic effects are weak because of the typically low ionization of disks around young stars, dust absorbs free electrons and the field-lines decouple from the gas. In this situation potentially hydrodynamical instabilities could drive turbulence and structure formation. For instance a radial temperature gradient leads to vertical shear in the disk, which is unstable for short thermal relaxation times (Urpin & Brandenburg, 1998; Arlt & Urpin, 2004; Nelson et al., 2013; Stoll & Kley, 2014; Barker & Latter, 2015; Stoll & Kley, 2016; Latter & Papaloizou, 2017; Manger & Klahr, 2018; Lyra & Umurhan, 2018; Klahr et al., 2018). This was demonstrated to trigger zonal flows, long lived vortices (Raettig et al., 2013, 2015; Manger & Klahr, 2018), and instability based on radial buoyancy. The same can happen in disks which are radially buoyant, but for longer thermal relaxation times. Then, convective overstability (Klahr & Hubbard, 2014; Lyra, 2014) and its non-linear version, Subcritical Baroclinic Instabilities (Klahr & Bodenheimer, 2003; Petersen et al., 2007a, b; Lesur & Papaloizou, 2010; Lyra & Klahr, 2011), can produce and sustain vortices.

Johansen & Youdin (2007) performed 2D simulations on the Kelvin-Helmholtz instability of sedimenting particles and found that while on average particles were still diffused out of the midplane, the turbulence also led to strong fluctuations.

In order to form planetesimals grains have to circumvent the growth barriers. One possibility is a jump from a clump of condensed pebbles directly to planetesimals as first shown in the above mentioned shearing-box simulations by Johansen et al. (2007) in the gravo-turbulent scenario.

Also in the complete absence of turbulence in the disk, which seems unlikely given the observation of accretion in these systems, sedimentation of pebbles in disks with high enough dust to gas ratios, that is $2-3$ times solar metalicity, can lead to gravitational instability and planetesimal formation, regulated by the streaming instability as originally explored by Youdin & Shu (2002) and since then demonstrated in a series of numerical simulations (Johansen et al., 2009b; Simon et al., 2016).

Independent from the source of turbulence (S.I. or external) there may exist a typical planetesimal diameter regulated by particle diffusion (Klahr & Schreiber, 2015; Schreiber, 2018). Grain evolution models, including radial transport, leading to planetesimal formation, have been performed by various authors. Brauer et al. (2008b) have shown that particles with a radius of $\sim 10^{2}\,\mathrm{m}$ can be produced at the edges of regions from low turbulence (called the “dead zone”) to higher turbulence levels. A particle velocity distribution can lead to “lucky particles” which are able to circumvent the growth barriers via low relative speed collisions (Windmark et al., 2012; Drążkowska et al., 2014b). But this effect may be suppressed by lower breakup speeds (Blum & Münch, 1993; Blum & Wurm, 2008) within the (water) ice line, and in the outer disk by rapid removal due to radial drift as indicated by Estrada et al. (2016). We will discuss this in more detail in section 5.2.2. Drążkowska et al. (2013) have shown that planetesimal formation can occur in a pressure bump at the inner edge of the dead zone via sweep-up growth. Another possibility would be that particles become fluffy, and grow continuously to planetesimal size. This model will be further discussed in section 5.2.1. But particles can also pile up due to back-reaction of particles onto gas and evaporation and re-condensation at the water ice line (Drążkowska et al., 2016; Drążkowska & Alibert, 2017) which we discuss in section 5.2.3.

We will follow a different path where a plethora of various instabilities can create turbulent structures in almost the entire disk either by hydrodynamical (Pfeil & Klahr, 2019) or magnetic instabilities in the non ideal regime (Béthune et al., 2016). For a review see Klahr et al. (2018). Indeed observations find a lot of structure in protoplanetary disks, see e. g. the Disk Substructures at High Angular Resolution Project (DSHARP; e. g., Andrews et al., 2018; Huang et al., 2018; Dullemond et al., 2018). Hence, one can assume that trap structures will appear everywhere in the disk. These can then trap the incoming particles, which are concentrated to particle clouds collapsing to planetesimals.

This paper is structured as follows. In section 2 we present a new model for the planetesimal formation rate. Section 3 summarizes the numerical disk model and explains how we implemented our parameterized planetesimal formation within it. In section 4 we show the first results. We discuss limitations of the model and compare with other planetesimal formation models in section 5. Particular attention should be paid to section 5.2.3, where we explain that streaming instability is not mandatory to form planetesimals. We conclude in section 6.

2. Rational of Pebble Flux limited Planetesimal formation: Trapping Mode.

In order to mimic planetesimal formation in our simulations, we use findings from fluid dynamical simulations which include particles. To keep it simple and feasible we use a probabilistic toy model based on parameterization. The important disk model parameters for this model are the lifetime of pebble traps $\tau_{\mathrm{l}}$ (e. g. vortices or zonal flows) and their radial separation $d$. Strictly speaking, one doesn’t necessarily need a pressure bump with a local maximum. If the pressure profile is flattened a bit without a local maximum, this would slow down the particles since smaller pressure gradients lead to slower drift speeds (Whipple, 1972) and conservation of mass flux leads to denser particle accumulations. Instead, the density is increased which may lead to critical $\Sigma_{\mathrm{d}}/\Sigma_{\mathrm{g}}$ values to trigger streaming instability. In this scenario, the $d$ parameter has to be interpreted as the distance between such flat pressure structures. Or in general, $d$ is the distance between structures leading to streaming instability conditions with sufficient density increase to reach Hill density. In our case we just assume that the flux concentrators are vortices or zonal flows with typical numerically measured radial distance $5h_{\mathrm{g}}$ (Dittrich et al., 2013).

The Stokes number $\mathrm{St}$ describes the aerodynamic behavior of particles surrounded by gas. It is the ratio of the friction time (timescale of coupling to the gas) and the dynamical timescale of the gas (here the inverse of the Keplerian frequency $1/\Omega$). Particles with $\mathrm{St}\ll 1$ are well coupled to the gas whereas $\mathrm{St}\gg 1$ are decoupled from the gas. From the particle side, the parameters of our model are the Stokes number of the smallest ($\mathrm{St}_{\mathrm{min}}$) and largest ($\mathrm{St}_{\mathrm{max}}$) particles that are able to participate in the streaming instability to facilitate gravitational collapse of the particle heaps. And finally an efficiency parameter $\varepsilon$ that defines how much of the actual radial mass flux in pebbles can be trapped and transformed into planetesimals. Then one can convert the column density of drifting pebbles into column density of planetesimals with the following recipe:

[TABLE]

where we defined the mass flux of pebbles for a full circumference as

[TABLE]

Here, $v_{\mathrm{drift}}$ is the radial velocity with which particles drift and $\Sigma_{\mathrm{d}}(r,\mathrm{St})$ is the column density in particles with Stokes number $\mathrm{St}$ within a given bin in particle size. The units of this column density are still $\mathrm{g/cm^{2}}$. The cylindrical radius, i. e., distance to the star in the midplane, is given by $r$. Eq. (2) is written in discretized form, as we treat it in our simulations.

The conversion length over which pebbles are transformed into planetesimals is given by

[TABLE]

If, for example, there is only one pebble species with drift speed $v_{\mathrm{drift}}$ and column density $\Sigma_{\mathrm{peb}}$, we distribute the new planetesimal column density over the trap distance $d$. $\ell/v_{\mathrm{drift}}$ gives the timescale of conversion, thus $v_{\mathrm{drift}}/\ell\cdot\Sigma_{\mathrm{peb}}$ is the rate at which the transformation from $\Sigma_{\mathrm{peb}}$ to $\Sigma_{\mathrm{p}}$ occurs. If this is added up for all particle species, one obtains the column density formation rate of planetesimals as shown in Eq. (1). We barely resolve the expected radial trap structures with a few radial bins at maximum.

To give another explanation, $1/d$ is the radial trap density, meaning that the closer the traps are packed radially, the more planetesimals can be produced. For the flux, only $v_{\mathrm{drift}}$ without the gas velocity is considered to take the relative radial velocity between particles and gas. Hereby we assume that the spatial pressure structures (e. g., real pressure-bump traps like vortices) move radially with the gas of the smooth profile.

We assume theoretically infinitesimal structures as traps even though they are practically extended. Since we assume a non-zero $\alpha_{\mathrm{t}}$ value everywhere anyway, we also assume that turbulent structures can be found everywhere. For our numerical radial grid, the grid cells are smaller than the trap distance parameter $d$. This ensures that we can distribute the mass of newborn planetesimals accordingly (probabilistic ansatz). If the resolution is too low, i. e., $\Delta x\gg d$, a “distribution” of this mass is not physical or one cell traps too much mass. Also, the condition of accumulated particles may not be treated correctly anymore because a bigger cell can also host more mass. At the same time, we are far from resolving the length scale on which planetesimal formation occurs. This length scale is typically on the order of $0.01h_{\mathrm{g}}$ (Schreiber, 2018, see his Eq. (3.40)), depending on particle diffusion on these small scales and on particle Stokes number. Even though we may be able to resolve trap structures of the order of a gas pressure scale height with only a few radial bins, it is not enough to reach a proper resolution of the substructure. The properties of the traps in terms of particle trapping and planetesimal formation is embedded in our efficiency parameter $\varepsilon$.

The material density is assumed to be $\rho_{\mathrm{s}}=1.2\,\mathrm{g/cm^{3}}$ for all solid particles, according to asteroid data (e. g., Carry, 2012). For simplicity we again assume particles to be spherical objects with constant material density $\rho_{\mathrm{s}}$, such that the representative planetesimal mass with a diameter of $100\,\mathrm{km}$ (Klahr & Schreiber, 2015; Schreiber, 2018) is given by

[TABLE]

where ${M_{\oplus}}$ is the mass of the earth. The idea of planetesimals of diameter around $100\,\mathrm{km}$ is also supported by data from our solar system (Morbidelli et al., 2009; Delbo et al., 2017).

Planetesimal formation will only occur if, within one (average) lifetime of a trap $\tau_{\mathrm{l}}$, enough mass can be accumulated to form at least one planetesimal, i. e. if the following condition holds

[TABLE]

Hence, there is a critical flux

[TABLE]

that must be reached to allow planetesimal formation. If condition (5) is fulfilled, we call the flux critical, otherwise we call it sub-critical. We give more detailed information on that criterion in Appendix G.

We further assume that the relative speed between the gravoturbulence triggering structures and the particles is the drift speed $v_{\mathrm{drift}}$. But in reality, e.g., zonal flows can have radial velocities which would change this relative velocity. As long as the radial speed of the pressure bump structure is much smaller than the drift speed, the error is small. In this paper, again for simplicity, we assume that the relative speed is always given by $v_{\mathrm{drift}}$. Additionally, since we don’t know how $\varepsilon$ would change with $\mathrm{St}$ and $r$, we assume that it is a constant. As long as condition (5) is fulfilled, only the value of the conversion length $\ell$ matters. That is, one can change the value of $\varepsilon$ and $d$ by the same factor, leading to the same result for planetesimals.

3. Numerical Model

3.1. The disk and dust model

In the following, we will summarize the disk and dust model of the dust and gas evolution code from Birnstiel et al. (2010) in which we have implemented our planetesimal formation model.

We assume that turbulence is described by an effective kinematic viscosity of (Shakura & Sunyaev, 1973)

[TABLE]

The dimensionless parameter $\alpha_{\mathrm{t}}$ describes the strength of turbulence, since it determines the turbulent velocities. Most likely, the turbulent velocities of the largest eddies is roughly given by $\sqrt{\alpha_{\mathrm{t}}}c_{\mathrm{s}}$ as argued by Cuzzi et al. (2001).

We assume further that the gas is vertically in hydrostatic balance (Weizsäcker, 1948) and that consequently the particle and gas density follow the Gaussian profile (Safronov, 1969; Pringle, 1981)

[TABLE]

where $i$ is either dust/particles ($\mathrm{d}$) or gas ($\mathrm{g}$). To eliminate the $z$-dependence, we use column densities111In Eq. (9) $i$ can also be planetesimals (p).

[TABLE]

The mass density of the gas at the mid-plane ($z=0$) can then be expressed as

[TABLE]

The gas pressure scale height is given by

[TABLE]

and the dust/particle scale height for $\mathrm{St}<1$ by (Dubrulle et al., 1995)

[TABLE]

Here we assumed that turbulence is isotropic such that the dimensionless vertical diffusion coefficient is just given by $\alpha_{\mathrm{t}}$. In the code, another expression for $h_{\mathrm{d}}$ is used which is essentially equivalent to the one in this equation (see Eq. (51) of Birnstiel et al. (2010) for details). Since we use a one-dimensional model with $z$-integrated values, and the coagulation equation scales with the density squared, the Stokes number is calculated at mid-plane values where densities are the highest. Furthermore, we assume we always stay in the Epstein drag regime (Epstein, 1924) which leads to (Birnstiel et al., 2010, 2012)

[TABLE]

Since we use $\mathrm{St}_{\mathrm{max}}=10$ throughout in this paper, one might think that in this case one would be in the Stokes drag regime (see Eq. (126) of Stokes, 1851). The transition between Epstein and Stokes drag occurs if the gas mean free path is around $4/9$ the radius of the particle (Weidenschilling, 1977). Since Stokes numbers of $10$ are only reached in the very outer part of the disk, where the gas density is sufficiently low, the gas mean free path is many orders of magnitude larger than the size of any particle. Dependent on the gas density, dust grains don’t have to grow by much, if at all, to reach these large Stokes numbers. The Stokes drag regime is only reached within a few $\mathrm{au}$. But, again, we will ignore this regime here to avoid overlapping effects with our model.

For growth and fragmentation we solve the Smoluchowski equation (von Smoluchowski, 1916). It considers all binary combinations of particle size bins. It uses a probability for fragmentation and coagulation, depending on the collision speed. This method is identical to the one of Birnstiel et al. (2010), but more details on the Smoluchowski equation can be found in Appendix A.

The relative speed of collisions, $\Delta{v}$, can have different contributions. There are systematic velocities from radial drift (Weidenschilling, 1977; Nakagawa et al., 1986), azimuthal drift (Nakagawa et al., 1986), and vertical settling (Safronov, 1969, p. 26). These scale with the particles Stokes number and vanish for equally sized grains. On the other hand, random speeds such as turbulence (Ormel & Cuzzi, 2007) and Brownian motion lead to non-vanishing relative velocities between equally sized particles. The relative turbulent velocities are proportional to $\sqrt{\alpha_{\mathrm{t}}}$, thus $\alpha_{\mathrm{t}}$ can be interpreted as a measure not only of turbulence strength in terms of typical turbulent gas velocities, but also of particle collision speed.

The radial gas velocity is given by (Lynden-Bell & Pringle, 1974)

[TABLE]

Since the gas motion is determined by a balance of forces among stellar gravity, the gas pressure gradient force, and the centrifugal force, it orbits the star with a different speed than solids, which don’t feel a pressure gradient force. Due to a lack in centrifugal force ($\mathrm{St}<1$), single particles drift radially with a velocity of (Adachi et al., 1976; Weidenschilling, 1977)

[TABLE]

with $P=c_{\mathrm{s}}^{2}{\rho_{\mathrm{g}}}$, i. e., assuming an isothermal case with adiabatic index $1$. Since radial gas motion drags dust along, the total radial dust velocity is given by (Takeuchi & Lin, 2002)

[TABLE]

We model the fragmentation probability as a smooth transition from $1$, for collision speeds above the breakup speed $v_{\mathrm{f}}$, to [math] for collision speeds below,

[TABLE]

where the radius dependent fragmentation velocity is modeled by the smooth function

[TABLE]

The coagulation probability is given via $p_{\mathrm{c}}+p_{\mathrm{f}}=1$. Experiments with silicate dust grains measured velocities of $1\,\mathrm{m/s}$ for the onset of fragmentation (Blum & Münch, 1993; Blum & Wurm, 2008), and theoretical studies (Leinhardt & Stewart, 2009) found similar values. It was found numerically (Wada et al., 2009) that the fragmentation velocity for $10\,\mathrm{\umu m}$ sized icy dust aggregates can reach $50\,\mathrm{m/s}$. Laboratory studies found threshold speeds about $10\,\mathrm{m/s}$ for icy grains (Gundlach & Blum, 2014). But this speed depends on the monomer size. Experiments show that $v_{\mathrm{f}}$ decreases with grain size for silicate dust grains (Beitz et al., 2011), which could be partly attributed to the increasing influence of inhomogeneities with growing grain size. For ice, measured erosion threshold velocities are around $15\,\mathrm{m/s}$ (Gundlach & Blum, 2014). However, we choose $v_{\mathrm{f}}=10\,\mathrm{m/s}$ outside the water ice line. As already mentioned by Birnstiel et al. (2010), the decrease of the threshold collision speed within the ice line leads to a traffic jam effect since inward drifting particles are forced to fragment to smaller sizes, which in turn drift slower.

We follow Youdin & Lithwick (2007) and model the radial particle diffusivity as

[TABLE]

and we assume that the radial gas diffusivity $D_{\mathrm{g}}$ equals the turbulent viscosity of Eq. (7).

Our initial but also fixed gas column density profile follows the self-similar profile as found by Lynden-Bell & Pringle (1974),

[TABLE]

The normalization constant is given by

[TABLE]

At the so-called characteristic radius ${r_{\mathrm{c}}}$, the transition between the power-law and the exponential law occurs. The initial profile of the total dust amount is assumed to follow

[TABLE]

The initial dust-to-gas ratio is set to $Z(t_{0})=0.01$ for all simulations, as found for the interstellar medium (e. g., Savage & Jenkins, 1972). We also fix $\gamma=1$ and ${r_{\mathrm{c}}}=35\,\mathrm{AU}$ in this paper. Furthermore, we keep the static gas column density profile and do not let the disk (viscously) evolve. This has the advantage that we can see the effects of our model parameters for planetesimal formation more clearly without overlapping effects from gas mass transport. Even though we don’t let $\Sigma_{\mathrm{g}}$ evolve, we allow a radial gas velocity according to Eq. (14). This leads to outward transport of small grains in the outer disk. For our model parameters, this occurs roughly from $20\,\mathrm{au}$ to $200\,\mathrm{au}$ and below a particle Stokes number of ($\mathrm{St}\ll 1$)

[TABLE]

as derived in Appendix B. For the initial size distribution we assume the dust to follow the distribution in the interstellar medium found by Mathis et al. (1977). I. e., the number density $n$ is described by the power-law

[TABLE]

with dust radii ranging from $a=0.1\,\mathrm{\umu m}$ to $1\,\mathrm{\umu m}$. For a review on dust evolution in circumstellar disks see, e. g., Birnstiel et al. (2016).

We fix the temperature to a profile similar to that of a radiative equilibrium disk (Chiang & Goldreich, 1997, they derived $T\propto r^{-3/7}$)

[TABLE]

The ice line is defined as the position where $T=170\,\mathrm{K}$ in this paper. Within the ice line we reduce the mass of particles to $1/3$ of their value for planetesimal formation, i. e. we assume an ice mass fraction of $2/3$. For this, we define the rocky-to-(ice and rocky) mass parameter

[TABLE]

Lodders (2003) summarizes her findings in Table 11 where the total rock/(total rock+H2O ice) ratio gives $0.46$. In section 3.3 of Min et al. (2011) they find values around $0.54$. Hayashi (1981) uses a value of $0.24$ for this ratio. We choose a value in between by assuming $1/3$ but in the results it is only important to notice the kink in the final planetesimal distribution at the water ice line. For the particle dynamics and growth processes this reduction in mass is not considered. Within the water ice line our simulation thus has too much mass in particles, but the maximum particle size is only weakly affected by this since, in the fragmentation limit, particles grow to the limit again faster than they drift inward. To be able to analyze the impact of our planetesimal model parameters, we do not include accretion heating (Pringle, 1981), opacity effects, or radiative evolution.

3.2. Parameterized Planetesimal Formation

Do we need to restrict the particle species contributing to planetesimal formation? This is the question we are facing in this section. Since particle dynamics depend on the Stokes number, one could think that there exists an interval $[\mathrm{St}_{\mathrm{min}},\mathrm{St}_{\mathrm{max}}]$ between which particles have the right properties for trapping and building clumps of Hill density. The pebble flux regulates the rate at which material is delivered, so this does not need to be restricted by a Stokes number interval. The question is whether there exist particles which cannot participate in trapping and planetesimal formation. We know that both trapping, either in zonal flows (Dittrich et al., 2013) or vortices (Fu et al., 2014; Raettig et al., 2015), and collapse to planetesimals (Schreiber, 2018) work for particles with $\mathrm{St}=10^{-2}$. From Carrera et al. (2015) we learned that $\mathrm{St}_{\mathrm{max}}\approx 10$. Yang et al. (2017) have shown that streaming instability is possible for $\mathrm{St}=10^{-3}$ particles if $Z\geq 4\cdot 10^{-2}$. We know that particles with $\mathrm{St}=10^{-2}$ still contribute to trapping and gravitational collapse to planetesimals. But so far, we cannot entirely exclude that smaller particles can take part in this process as well—especially for mixtures of different particle species with some mass distribution depending on their Stokes number. However, we use $\mathrm{St}_{\mathrm{min}}=10^{-2}$ and $\mathrm{St}_{\mathrm{max}}=10$ in this paper.

The velocity field of a vortex or zonal flow is modified around pressure bumps, and thus gas motion here differs from the gas motion in a smooth pressure structure. We assume that this motion can be neglected in comparison to the drift velocity.

We further assume that these structures have an average lifetime $\tau_{\mathrm{l}}$ proportional to the orbit time. It has been shown that both zonal flows (Dittrich et al., 2013) and vortices (Manger & Klahr, 2018) can endure up to hundreds of orbits. Hence, we assume $\tau_{\mathrm{l}}=100\,t_{\mathrm{orb}}$ in this paper. Based on work by Dittrich et al. (2013) we assume an average radial separation of $d=5h_{\mathrm{g}}$.

Our analytical description of planetesimal formation arises in the code as a sink term in the advection-diffusion equation for the column density of a single particle species of mass $m_{i}$, given by

[TABLE]

The Heaviside functions $\theta(\cdot)$ represent our conditions that the pebble flux must be critical and that only particles with $\mathrm{St}_{\mathrm{min}}\leq\mathrm{St}\leq\mathrm{St}_{\mathrm{max}}$ are allowed to build planetesimals. We use the flux conserving implicit donor-cell scheme from Appendix A.1 of Birnstiel et al. (2010), where we set

[TABLE]

in their Eq. (A.1) in order to have an implicit loss term due to planetesimal formation. The number of new planetesimals created within the time step $\operatorname{d\!}{t}$ is then given, in terms of column density, by

[TABLE]

The pebble flux $\dot{M}_{\mathrm{peb}}$ is estimated from the last time step, but the resulting error from this is very small.

4. Results

We start our analysis with the pebble flux because it is directly linked to our planetesimal formation. Particle dynamics is mostly determined by the growth barriers (Birnstiel et al., 2012), i. e. the fragmentation barrier

[TABLE]

and the drift limit (here in the Epstein regime only)

[TABLE]

As can be seen in Fig. 1, for stronger turbulence levels (here $\alpha_{\mathrm{t}}=10^{-2}$) the particle flux stays high for longer times because particles are smaller and thus drift slower, see Fig. 2. Here we defined the dust column density distribution per logarithmic bin of grain radius $a$ as

[TABLE]

with $n_{a}$ being the number density per grain size bin. The total dust column density can then be written as

[TABLE]

In the early stages, the flux is smaller for lower $\alpha_{\mathrm{t}}$ for the same reason. The pebble growth front, which is the radius at which $\mathrm{St}_{\mathrm{min}}$ particles form the first time, moves from the inside to the outside, while particles drift inward from the outer disk because their Stokes number is high. When these two fronts clash, the density suddenly rises and a particle wave propagates inward which is why one can see a small local peak in the curves of the lower left panel. In the lower right this is not visible since fragmentation is damping this effect.

At late times ($\sim 10^{6}\,\mathrm{yr}$) the pebble to total mass flux ratio drops significantly in the $\alpha_{\mathrm{t}}=10^{-3}$ case. This is due to the drift barrier, which drops to smaller sized particles since $\Sigma_{\mathrm{d}}$ is decreasing, see Fig. 3. The gray zone in the upper panels mark the sub-critical flux values, i. e., the fluxes that are too low to allow formation of planetesimals within a trap lifetime. One could change $\tau_{\mathrm{l}}$ or $m_{\mathrm{p}}$ by a factor of $10$ without changing the results significantly because the flux is orders of magnitude larger then the critical value.

Since the mass reservoir is in the outer disk where gas densities are very low, the pebble flux in the outer region around $\sim 200-400\,\mathrm{au}$ already reaches critical values, and leads to planetesimal formation, in the early phase ($\lesssim 10^{4}\,\mathrm{yr}$). When the slope of $\Sigma_{\mathrm{p}}$ is steeper than $r^{-2}$ (see Fig. 4), the mass contribution to the total planetesimal population in the disk will not rise with $r$. I. e., the contribution becomes less the further out in the disk the considered regions are. For a fixed logarithmic binning, this can be seen by writing the planetesimal mass between disk radius $r_{1}$ and $r_{2}$ as

[TABLE]

Hence, the probably unphysical outer bump in $\Sigma_{\mathrm{p}}$ doesn’t contribute significantly if the slope is steeper than $r^{-2}$. We think this outer formation region would not exist in the early evolution of the disk, since the formation time of structures such as vortices means they need a few orbits to form, which would take much longer in the outer than in the inner disk.

If turbulence is strong (here $\alpha_{\mathrm{t}}=10^{-2}$), fragmentation will be the growth limiting process. In that case not only is the maximum grain size smaller, but there are also more small grains in general. Hence, in comparison to simulations with weaker turbulence (here $\alpha_{\mathrm{t}}=10^{-3}$), particles in the outer disk stay there for longer and can form planetesimals. Fig. 4 and 5 show the resulting planetesimal column density $\Sigma_{\mathrm{p}}$ and formation rate $\dot{\Sigma}_{\mathrm{p}}$, respectively, both for different snapshots showing the entire disk and for local evolutions. For the latter, we chose $3\,\mathrm{au}$ ($\sim$ position of the main Asteroid Belt), $40\,\mathrm{au}$ ($\sim$ position of the Kuiper Belt), and $100\,\mathrm{au}$ to show the behavior in the outer disk. Planetesimal formation basically starts shortly after the first pebbles have formed. We summarize the effects as follows:

• •

$\alpha_{\mathrm{t}}=10^{-2}$: higher $\varepsilon$ generally lead to higher planetesimal formation rates but these also decline faster. In the end this nevertheless ends in a more abundant planetesimal population for higher $\varepsilon$ in the entire disk.

• •

$\alpha_{\mathrm{t}}=10^{-3}$: in the outer disk $\dot{\Sigma}_{\mathrm{p}}$ is also higher for larger $\varepsilon$, such as in the higher $\alpha_{\mathrm{t}}$ case. Radial drift removes the largest particles of the top-heavy size distribution. The fewer particles were converted into planetesimals in the outer disk, the more material is available at later times in the inner regions (conservation of mass). This can lead to more planetesimals in the inner disk for smaller $\varepsilon$. This effect cannot be seen in the results of the $\alpha_{\mathrm{t}}=10^{-2}$ simulations since fragmentation events slow down the radial motion and the feeding zones of planetesimals are much narrower.

If the mass flux of pebbles is roughly $r$-independent, one can estimate the slope of the outcoming planetesimal profile as we show in section 4.2.

In Fig. 6 we show the column densities of gas, dust and planetesimals for different snapshots and locations in the disk. This way one can not only see that for $\alpha_{\mathrm{t}}=10^{-3}$ plantesimals form at $2-3\,\mathrm{au}$ after a few thousand years, but also after a few ten to hundred thousand years at $20-30\,\mathrm{au}$. Fig. 6 also shows the ratios of dust to gas, planetesimals to dust and planetesimals to initial dust. For $\alpha_{\mathrm{t}}=10^{-2}$, the zone where planetesimals are formed moves faster outward due to faster growth from dust to pebbles. We summarize the different ratios as follows:

dust-to-gas ratio $\Sigma_{\mathrm{d}}/\Sigma_{\mathrm{g}}$

This ratio is constantly decreasing since we kept the gas density constant. In the early stages of evolution ($\lesssim 10^{6}\,\mathrm{yr}$) $\Sigma_{\mathrm{d}}/\Sigma_{\mathrm{g}}$ decreases faster for smaller $\alpha_{\mathrm{t}}$ because particles drift faster.

planetesimals-to-dust ratio $\Sigma_{\mathrm{p}}/\Sigma_{\mathrm{d}}$

• •

$\alpha_{\mathrm{t}}=10^{-3}$: from $10^{4}\,\mathrm{yr}$ on, values above unity are reached within $2\,\mathrm{au}$. At $10^{6}\,\mathrm{yr}$ this ratio is $>1$ almost everywhere in the disk.

• •

$\alpha_{\mathrm{t}}=10^{-2}$: only at late times ($\gtrsim 10^{6}\,\mathrm{yr}$) values above unity are reached.

planetesimals-to-initial dust ratio $\Sigma_{\mathrm{p}}/\Sigma_{\mathrm{d}}(t_{0})$

• •

$\alpha_{\mathrm{t}}=10^{-3}$: at $10^{5}-10^{6}\,\mathrm{yr}$ within $10\,\mathrm{au}$ values above unity are reached, indicating a larger feeding zone of planetesimals, i. e. material significantly contributing to planetesimals can come from far outside.

• •

$\alpha_{\mathrm{t}}=10^{-2}$: the ratio never becomes greater than $1$ which indicates that fragmentation allows more time to form planetesimals out of the available material, i. e. the drift timescale of single particles including destruction events is slower than the timescale of planetesimal formation.

4.1. Necessity for lower turbulence ($\alpha_{\mathrm{t}}$) to form Earth

If our assumptions are correct, no planetesimals are formed within $2\,\mathrm{au}$. The reasons for that are on the one hand the value of $\mathrm{St}_{\mathrm{min}}$ and on the other hand the temperature profile. The latter dictates the radial shape of the fragmentation barrier. For a more realistic temperature model, accretion heating would lead to higher temperatures in the inner disk and a steeper slope. As a result of that the water ice line would be shifted further out. Hotter temperatures lead to a larger gas scale height and thus also to larger $d$. At the same time particles would be smaller due to the lower fragmentation barrier. Both effects would cause a planetesimal free zone, or at least a much smaller planetesimal population in the inner disk ($\lesssim 3\,\mathrm{au}$). I. e., these problems we see here would be even worse, which makes it very likely that our conclusion still holds if more physics is taken into account. As we show in Appendix E, the total mass of planetesimals in the entire disk is also problematic if $\alpha_{\mathrm{t}}$ is high.

Yang et al. (2018) measured in shearing box simulations with non-ideal MHD and a dead zone values of the order of $\alpha_{\mathrm{t}}\sim 10^{-4}$ to $10^{-3}$. In layered accretion theory values of $\alpha_{\mathrm{t}}\sim 10^{-6}$ to $10^{-4}$ seem to be preferred (Johansen et al., 2015). Hence, our finding is congruent with their findings.

4.2. Planetesimal column density steeper than that of initial gas and dust

Fig. 6 shows that for lower turbulence (here $\alpha_{\mathrm{t}}=10^{-3}$) the slope of the planetesimal column density can be steeper than that of initial gas and dust. This result stems mostly from the $r$-dependency of $d\propto h_{\mathrm{g}}$. On top of this, transport phenomena can change it as well. For $\varepsilon=0.1$, this can even be shown analytically since the pebble flux does not depend on $r$ over a larger region of the inner disk (Fig. 1). We find that the final slope of the planetesimal column density is steeper than that of the gas if the effective efficiency $\varepsilon/d$ is sufficiently small (here $0.02/h_{\mathrm{g}}$, $\varepsilon=0.1$). If this effective efficiency is large (e. g., $0.16/h_{\mathrm{g}}$, $\varepsilon=0.8$), the feeding zone of planetesimals is much smaller, i. e., planetesimal material stems from a local region, leading to a slope much closer to the initial dust density, which is the same as for the gas in our simulations.

If the pebble flux is constant for all $r$, i. e., for size distributions dominated by particles with Stokes numbers $\mathrm{St}_{\mathrm{min}}\leq\mathrm{St}\leq\mathrm{St}_{\mathrm{max}}$ and sufficiently small $\varepsilon/d$, then we can conclude that

[TABLE]

where $v_{1}$ is the pebble drift speed. For a constant pebble flux, $\dot{M}_{\mathrm{peb}}$, the pebble column density is given by

[TABLE]

leading to

[TABLE]

For the temperature power-law index we used throughout this paper, $q=0.4$, this gives a planetesimal column density profile of $\Sigma_{\mathrm{p}}\propto r^{-2.3}$, independent of the gas column density profile. In this case, the slope of the planetesimal column density $\Sigma_{\mathrm{p}}$ does not depend on $\varepsilon$, but the total value of $\Sigma_{\mathrm{p}}$ scales with it linearly. For regions limited by fragmentation, the condition that the flux must be dominated by pebbles ($\mathrm{St}_{\mathrm{min}}\leq\mathrm{St}\leq\mathrm{St}_{\mathrm{max}}$) is almost always violated since fragmentation grinds particles to sizes with Stokes numbers smaller than $\mathrm{St}_{\mathrm{min}}$. Due to fragmentation events small particles are constantly replenished. But the result on the right side of proportionality (37) does not, in principal, depend on whether growth is limited by drift or fragmentation—as long as the conditions for a radius-independent pebble flux are fulfilled. Furthermore, as can be seen in Fig. 1, the pebble flux may only vary with $r$ within roughly one order of magnitude, which is still good enough to estimate the slope of $\Sigma_{\mathrm{p}}$.

5. Discussion

5.1. Limitations of the current model

In this paper, the gas density is kept constant, i. e., we neither allow the disk to evolve viscously (Lüst, 1952; Pringle, 1981) nor do we include sink terms such as photoevaporation (e. g. Owen et al., 2012). We also don’t include a disk build up phase as tested by Birnstiel et al. (2010) or Drążkowska & Dullemond (2018). In order to isolate the influence of the new parameters of the model for analysis, we assume a constant temperature profile according to Eq. (25) and ignore viscous heating (Pringle, 1981). We keep $v_{\mathrm{drift}}$ as for a smooth gas profile, but regions with shallow pressure profiles or even pressure bumps would slow down the drift speed. This effect was not considered here since we have a lack of information about how exactly $v_{\mathrm{drift}}$ would be affected by zonal flows, vortices, or other structures leading to a significant slow-down of particles and subsequent gravoturbulent planetesimal formation. Also, planetesimals stay in the annulus where they are born. Their mass and radius are not allowed to change, i. e., we ignore collisions between them and also growth due to pebble accretion (Ormel & Klahr, 2010), as well as radial motion. Furthermore, we can only give a first impression of how planetesimal formation is described and behaves in this model, since we have not explored the full parameter space ($M_{\mathrm{disk}}$, ${r_{\mathrm{c}}}$, $Z(t_{0})$,…).

The turbulence parameter $\alpha_{\mathrm{t}}$ mostly controls the relative velocities between grains (Brauer et al., 2008a). Though this parameter from Eq. (7) can be given by disk winds which cause the transport of angular momentum (see e. g. Papaloizou & Lin, 1995; Béthune et al., 2017), it can also be entangled in, e. g., vortex formation. Hence, a trap formation time, together with a $\alpha_{\mathrm{t}}(r,t)$ model, would make the results more realistic—especially in the outer disk region where it may take a long time to form vortices. Due to conservation of momentum, there is not only the influence of gas on dust but also a back-reaction of particles on gas (Tanaka et al., 2005; Nakagawa et al., 1986). This back-reaction becomes important only for $Z\gtrsim 0.1$, which is not reached in the simulations presented in this paper. But without planetesimal formation, this value can be reached in the vicinity of the ice line where particles are forced to fragment and, thus, to slow down. The sink term of pebbles due to planetesimal formation confines this effect. Furthermore, we do not consider the bouncing or charging barrier and we do not trace volatiles.

As pointed out by Drążkowska et al. (2014a), Smoluchowski solvers need a high resolution in mass to mimic the upper end of the size distribution well enough. Especially in the high $\alpha_{\mathrm{t}}=0.01$ case this upper end can be important for planetesimals formation.

5.2. Comparison to other models

In this section we will compare with the most prominent models for planetesimal formation.

5.2.1 Continuous fluffy growth

Okuzumi et al. (2012) and Kataoka et al. (2013) describe continuous growth to planetesimals. They assume that the threshold velocity for fragmentation is never reached, thus avoiding disruptive events. Without replenishment of small particles these shouldn’t be seen in observations. This "fluffy path" of planetesimal formation is highly sensitive to the disk parameters (e. g. $\alpha_{\mathrm{t}}$) and works roughly up to $10\,\mathrm{au}$ for the conditions presented in these papers. Another problem would be that the size distribution of initial planetesimals cannot explain the kink feature we see in the size distribution of the Asteroid belt (Bottke Jr et al., 2005), the Kuiper belt (Fuentes & Holman, 2008; Fraser & Kavelaars, 2008), Jupiter trojans (Jewitt et al., 2000), and of Neptune trojans (Sheppard & Trujillo, 2010, see their Fig. 4 for an overview) today. Even though runaway growth may start at the size of the kink feature (Kobayashi et al., 2016) which depends on the turbulence level $\alpha_{\mathrm{t}}$ and the mass which is initially available (see their Fig. 6). But this feature is covered by the theory of planetesimal formation via gravitational instability (Morbidelli et al., 2009; Klahr & Schreiber, 2015; Schreiber, 2018), almost independent of disk radius and very robust against changes in $\alpha_{\mathrm{t}}$. However, one cannot dismiss the idea of fluffy growth entirely, especially in the early stages of particle growth. Even some planetesimal could be formed this way without a conflict with today’s data.

However, as shown by Krijt et al. (2015), erosion, where particles lose mass due to collisions with smaller projectiles, may prevent growth through the radial drift barrier. They also show that only for high erosion velocity values of $60\,\mathrm{m/s}$ porous particles can overcome this barrier within $10\,\mathrm{au}$ for a minimum mass solar nebula (Weidenschilling, 1977; Hayashi, 1981).

5.2.2 Lucky particle growth

By including a Maxwell-Boltzmann like velocity distribution, lucky particle effects can occur as described by Windmark et al. (2012) and tested by Drążkowska et al. (2014b) and Estrada et al. (2016). I. e., growth barriers are smoothed out or can be overcome by low velocity collisions leading to growth where the mean collisional speed would lead to bouncing or fragmentation. Although these low velocity collisions have a low probability, by experiencing multiple lucky collisions in a row, the growth barriers can be overcome. In the case of the bouncing barrier, this effect is even stronger because higher collision velocities lead only to bouncing, i. e. no disruption but possible further compactification, but low velocity collisions let the grains grow. Therefore, the bouncing barrier may not be a strict solid barrier as long as the growth timescale (considering bouncing) is sufficiently small. By overcoming the bouncing and fragmentation barrier, these particles present the seeds for a sweep-up scenario which can cause a growth process toward planetesimal size.

But this path is quite inefficient in forming planetesimals and it did not take drift into account. The key here is to have most collision partners leading to mass transfer or sticking. But a lower fragmentation threshold speed within the ice line would basically make it impossible to grow beyond the barriers (or very unlikely). And in the outer disk, drift is removing the lucky particles faster than they can grow beyond the bouncing barrier (Drążkowska et al., 2014b; Estrada et al., 2016).

5.2.3 Planetesimal formation directly regulated via streaming instability

Drążkowska et al. (2014a) and Drążkowska et al. (2016) proposed a parameterized formation model of planetesimals by asking when and where the pebble density ($\mathrm{St}>0.01$) in the disk midplane reaches a critical value of ${\rho_{\mathrm{d}}}/{\rho_{\mathrm{g}}}=1$, assuming a balance of sedimentation and vertical diffusion of pebbles. At this critical dust-to-gas ratio it can be argued that the streaming instability will be triggered (Youdin & Goodman, 2005). But as was shown by Johansen et al. (2009b), Carrera et al. (2015), and Yang et al. (2017) triggering the streaming instability at ${\rho_{\mathrm{d}}}(\mathrm{St}>10^{-2})/{\rho_{\mathrm{g}}}=1$ is not sufficient to create strong enough over-densities to lead to gravitational collapse. The additional criterion to have the total local column density dust-to-gas ratio ($\Sigma_{\mathrm{d}}/\Sigma_{\mathrm{g}}$) increased by a factor of 2–4 seems to be implicitly tested by having sufficient sedimentation even at $\alpha_{\mathrm{t}}=10^{-4}$ as vertical diffusion. But, as argued by Drążkowska et al. (2016), and also shown by Carrera et al. (2017, see their Fig. 2), reaching the critical midplane dust-to-gas ratio is typically the stronger condition. Similar to our ansatz, Drążkowska et al. (2016) convert pebbles into planetesimals with a certain efficiency per orbit, where the column density formation rate $\dot{\Sigma}_{\mathrm{p}}\propto\Sigma_{\mathrm{d}}(\mathrm{St}>10^{-2})$ in their model. I. e., they use a local approach whereas we focus on the pebble delivery via the local pebble flux.

In the work of the follow-up paper by Drążkowska & Alibert (2017), they follow a trap that will actually show up in 1D radial disk models, e.g. pebble pile beyond the (water) ice line due to traffic jam and recondensation of water vapor. This means planetesimal formation will mainly happen outside of and close to the ice line in their model. As a result, planetesimals will be confined to a relatively narrow region between 2 and 5$\,\mathrm{au}$. Besides that they likely will be too water rich to explain the composition of terrestrial planets. Not only are planetesimals only formed where the water ice line is moving (relatively small annulus compared to the disk size), but the time for planetesimal formation is also very limited ($\sim 2-3\cdot 10^{5}\,\mathrm{yr}$).

Since in our model planetesimals are also formed in the outer disk, less material will reach the inner part. Hence, we may not reach their condition ${\rho_{\mathrm{d}}}(\mathrm{St}>10^{-2})/{\rho_{\mathrm{g}}}\geq 1$, which we illustrate in Fig. 7. Even though the back-reaction of particles on gas and the effect of collective drift may not be important here, tracing water vapor could make a difference. Nevertheless, we don’t reach the condition used by Drążkowska et al. (2016) since the midplane dust-to-gas ratio is more than one order of magnitude below the critical value of unity over almost the entire disk. In principle, a combined formation model of both pebble flux regulated formation in pressure bumps and the local conversion could coexist. Maybe photoevaporation can reduce the gas content at late times ($\gtrsim 10^{6}\,\mathrm{yr}$) so much that the condition can be reached in some parts of the disk.

We would like to emphasize that streaming instability per se is not mandatory for gravitational assisted planetesimal formation, and in fact it prevents planetesimal formation for low enough overall dust-to-gas ratios in the disk (Johansen et al., 2007). As in star formation (Klessen et al., 2005), turbulence plays a dual role. On global scales it promotes planetesimal formation by leading to clumps and enhancing the local dust-to-gas ratio, while at the same time it can prevent collapse on smaller scales by diffusing material away. I. e., turbulence stabilizes dust against collapse on small scales (Schreiber, 2018). Also the resulting planetesimal size is a question of particle diffusion, which does not necessarily have to stem from streaming instability. We summarize the work that has been done on fundamental research concerning planetesimal formation with and without turbulence in Table 2. If disk models assume some non-zero value for $\alpha_{\mathrm{t}}$ as in Drążkowska et al. (2016), it means that large scale turbulence is active. Whether there is a critical $Z$, ${\rho_{\mathrm{d}}}/{\rho_{\mathrm{g}}}$ or streaming instability criterion for these turbulent cases (magneto-rotational instability (MRI) etc.) in order to form planetesimals has not been studied so far. The simulations by Johansen et al. (2007) and Johansen et al. (2011) were in fact for MRI turbulent cases, and found planetesimal formation in zonal flows with streaming instability playing a minor role in enhancing perturbations on smaller scales. Planetesimals formed with and without particle feedback, i. e., streaming instability. These simulations unfortunately had to use quite large particles ($\mathrm{St}=0.25$), which is due to the problem of covering the large scales of MRI turbulence on scales of $\sim h_{\mathrm{g}}$, whereas the streaming instability unstable wavelength for smaller Stokes number particles of $\mathrm{St}=0.1-0.01$ is more like $0.01h_{\mathrm{g}}$. This makes 3D simulations of MRI and streaming instability leading to planetesimal formation for small Stokes numbers too challenging at the moment.

Higher $Z$ damps streaming instability (Schreiber, 2018), which means lower diffusion. And when there is less diffusion, planetesimals form (i. e., for weaker streaming instability), as seen in Bai & Stone (2010, see their table 2 and their section 7.2.4.). In our Table 2 we collect evidence that streaming instability is not needed for planetesimal formation, but it is a controlling effect—it can even prevent the formation. Hence, we suggest following Johansen et al. (2006) and calling the process "gravoturbulent planetesimal formation". Remember, in a disk without a local radial pressure gradient, there is no streaming instability and the Goldreich & Ward (1973) mechanism will work perfectly as recently demonstrated by Abod et al. (2018).

Auffinger & Laibe (2017) find streaming instability to be active in local pressure maxima, despite the pressure gradient vanishing locally. Whether this triggering of streaming instability has an impact on the formation of planetesimals has not been studied in detail. Yet, Johansen et al. (2007) and Johansen et al. (2011) find indeed the formation of planetesimals in these bumps, leaving the role of streaming instability undetermined so far. Furthermore, there can also be streaming instability in the $r$-$\varphi$ (radial-azimuthal) plane (Schreiber & Klahr, 2018).

6. Conclusion and Outlook

We presented a new model for the planetesimal formation rate that is directly linked to the local pebble flux—and not local dust-to-gas ratio or local density—relying on parameters such as lifetime of pressure-bump structures, contributing particle Stokes numbers ($\mathrm{St}_{\mathrm{min}}$ & $\mathrm{St}_{\mathrm{max}}$), the radial particle trap density ($1/d$), and an efficiency parameter $\varepsilon$ describing the conversion from pebbles to planetesimals. We have implemented the presented model into a sophisticated dust and gas evolution code, where we switched off gas evolution. This model was investigated for two different values of the turbulence parameter, $\alpha_{\mathrm{t}}=10^{-3}$ (moderate turbulence) and $10^{-2}$ (strong turbulence), as well as two different efficiency values, $\varepsilon=0.1$ and $0.8$. The results can be summarized as follows:

• •

This model is the first one creating planetesimals everywhere in the disk from a few $\mathrm{au}$ to a few hundred $\mathrm{au}$.

• •

The resulting planetesimal distribution and the total mass of initial planetesimals highly depends on the level of turbulence ($\alpha_{\mathrm{t}}$ value), i. e. whether particle growth is limited by fragmentation or drift.

• •

If the disk is fragmentation limited, higher $\varepsilon/d$ leads to higher planetesimal formation rates, but these also decline faster. Nevertheless, higher efficiencies yield more planetesimals.

• •

If the disk is mostly limited by drift, material drifts faster and particles which were not converted into planetesimals in the outer disk can build planetesimals in the inner part. That is, if the level of turbulence is not too high (here moderate turbulence for $\alpha_{\mathrm{t}}=10^{-3}$), the slope of the final planetesimal profile can be steeper than that of initial dust and gas. The overall planetesimal profile $\Sigma_{\mathrm{p}}$ looks also more like a power-law in the inner disk part, with an exponential roll off in the outer part for lower turbulence.

• •

A few planetesimals should have formed in the terrestrial region. Thus, strong turbulence ($\alpha_{\mathrm{t}}=10^{-2}$) seems to be unlikely for the Solar Nebula in order to form the Solar System, especially Earth— at least as long as the minimum Stokes number contributing to planetesimal formation is as high as $\mathrm{St}_{\mathrm{min}}=10^{-2}$ and the efficiency parameter $\varepsilon$ is not very high ($\gtrsim 0.8$ for $d=5h_{\mathrm{g}}$).

• •

The late stage planetesimal-to-dust ratios indicate a larger feeding zone of planetesimals, i. e. material forming planetesimals can originate from far away if $\varepsilon/d$ is small. Here, it is difficult to tell how small exactly, since the feeding zone will stretch with smaller values of that ratio. In this paper, $0.1/(5h_{\mathrm{g}})$ is small compared to $0.8/(5h_{\mathrm{g}})$. But one has to investigate in the future how exactly this feeding zone scales with $\varepsilon/d$.

Future work will analyze the influence of disk evolution and other disk parameters such as the initial disk mass and disk size. Implementing of a proper temperature model with accretion heating, photoevaporation, planetesimal collisions, and pebble accretion will be the focus of future work. For an upcoming study, one could compare the outcome of our simulations for the initial planetesimal population with constraints in special regions of the solar nebula to exclude certain pairings of disk parameters, or even extreme cases for specific parameters such as $\alpha_{\mathrm{t}}\geq 10^{-2}$. The efficiency parameter $\varepsilon$ of the combination of trapping and planetesimal formation can depend on both the particle Stokes number and the disk radius. Fluid dynamics codes have to show how this efficiency depends on particle and disk properties. In this paper, we assumed that particle traps are present from the beginning of the simulation. In reality these traps will take some time to form. Thus, one could introduce another parameter which describes the trap formation time. During this time, particles can drift inward, reducing the mass reservoir and the pebble flux available for transformation into planetesimals once particle traps are considered to be active.

C.L. is thankful to Andreas Schreiber, Joanna Drążkowska, Cornelis Dullemond, Akimasa Kataoka, Joachim A. Maruhn, Andrew Youdin, Konstantin Gerbig, Katherine Kretke, Hans Baehr, and Vincent Carpenter for helpful discussions. We also thank the anonymous referee for a useful report. This work was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) as part of the Schwerpunktprogramm (SPP, Priority Program) SPP 1833 “Building a Habitable Earth” and in part at KITP Santa Barbara by the National Science Foundation under Grant No. NSF PHY17-48958. Part of this work was performed at the Aspen Center for Physics, which is supported by National Science Foundation grant PHY-1607761. This research was supported by the Munich Institute for Astro- and Particle Physics (MIAPP) of the DFG cluster of excellence “Origin and Structure of the Universe”. T.B. acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program under grand agreement No. 714769 and funding from the DFG Ref no. FOR 2634/1.

Appendix A The coagulation approach

The first mathematical expression (in discrete form) for the coagulation process was obtained by von Smoluchowski (1916). It was first written in integral form by Schumann (1940), which is appropriate for a continuous mass spectrum. Consider a dust grain distribution $n_{\mathrm{m}}(m,r,z)$ as the number of particles per spatial volume and per particle mass interval $[m,m+\operatorname{d\!}m]$. In protoplanetary disks, $n_{\mathrm{m}}$ is a function of the solid particle mass $m$ as well as the distance $r$ and height $z$ (cylindrical coordinates). Since we assume a cylindrically symmetric system, the quantities do not depend on the azimuthal angle $\varphi$. Thus, the total dust mass density is

[TABLE]

A further generalized coagulation equation including fragmentation would be (see e. g., Barrow, 1981; Birnstiel et al., 2010)

[TABLE]

In this expression

[TABLE]

is the coagulation kernel,

[TABLE]

the fragmentation kernel, and

[TABLE]

the geometrical cross section , where $a_{i}$ is the particle radius of the colliding particle with mass $m_{i}$. $p_{\mathrm{c}}$ and $p_{\mathrm{f}}$ denote the probability that a collision between two particles of mass $m_{1}$ and $m_{2}$ with a relative velocity $\Delta{v}$ at collision leads to coagulation or fragmentation ($p_{\mathrm{c}}+p_{\mathrm{f}}=1$), respectively. Colliding particles with $m_{1}$ and $m_{2}$ give a gain term for particles with mass $m$, determined by the distribution of fragments $\mathcal{D}(m,m_{1},m_{2})$, where we assume that the distribution is a power-law according to

[TABLE]

We follow Brauer et al. (2008a) and use $\xi=1.83$, but this power-law index may depend on the collision speed (Husmann, 2017, his section 4.2). The Dirac delta-distribution is denoted by $\delta\left(\cdot\right)$.

The first term in Eq. (39b) corresponds to an increasing number of particles with mass $m$ due to grain growth, i. e., particles with masses $m_{1}$ and $m-m_{1}$ coagulate. The second term represents the loss of particles with mass $m$ ($m_{1}$ coagulates with $m$). The third term stands for the fragmentation due to a collision between particles with mass $m_{1}$ and $m_{2}$, and includes the fact that subsequently a distribution of some of their mass will occur via fragmentation to smaller sizes. The fourth term describes the fragmentation of particles with mass $m_{1}$ and $m$ and thus represents the loss of particles with mass $m$. The factors $1/2$ eliminate double counting of the collisions increasing the number of particles of mass $m$. We use 27 bins per size decade. The radial grid has 300 points spanning from $0.2$ to $800\,\mathrm{au}$. The grid points are logarithmically equally spaced, both in mass and disk radius. In order to optimize computational time, we only numerically consider mass cells up to twice the mass of the current maximum mass in the system at each disk radius, similar to section 2.2 of Lee (2000).

In general, all these quantities can also depend on other material properties such as composition, porosity and charge. To keep it simple we consider just the mass of the colliding particles. As Eq. (39a) is proportional to the product of two densities and particle densities become larger toward the mid-plane due to vertical settling, we follow Birnstiel et al. (2010) and approximate the kernels by the mid-plane values. However, the $z$-dependence is eliminated by vertical integration, assuming a Gaussian $z$-profile for the particle density. In the code, the discrete Smoluchowski equation is solved numerically. See Appendix A of Birnstiel et al. (2010) for details.

Appendix B Radial Particle velocities dominated by gas flow

It is worthwhile to estimate the Stokes number below, which the first velocity term in Eq. (16) dominates. Therefor we use

[TABLE]

assuming a falling temperature power-law with exponent $q$. Further, we assume a gas column density profile following Eq. (20) with $\gamma=1$, ${r_{\mathrm{c}}}=35\,\mathrm{au}$. This leads to the logarithmic pressure gradient

[TABLE]

which is obtained by using $P=c_{\mathrm{s}}^{2}{\rho_{\mathrm{g}}}$ for mid-plane values, see Eq. (10). There are regions in the $r$-$\mathrm{St}$ space where the total particle velocity is dominated by the gas velocity, i. e., by the left term in Eq. (16). We will now determine the condition on the Stokes number for this term to dominate the radial velocity. In the most cases the dust-to-gas ratio is small on large scales, thus, we assume ${\rho_{\mathrm{d}}}/{\rho_{\mathrm{g}}}\ll 1$. This way we obtain

[TABLE]

which can be written in the form

[TABLE]

By inserting the turbulent viscosity from Eq. (7) as well as Eqns. (41) and (42) we get, for a constant $\alpha_{\mathrm{t}}$ (after rearranging)

[TABLE]

which finally yields Eq. (23). If ${\rho_{\mathrm{d}}}/{\rho_{\mathrm{g}}}\ll 1$ is not fulfilled, this value has to be multiplied by $[1+({\rho_{\mathrm{d}}}/{\rho_{\mathrm{g}}})^{2}]$ for mid-plane values. Since $\Sigma_{\mathrm{g}}$ decreases exponentially for $r/{r_{\mathrm{c}}}>1$, whereas $\nu\sqrt{r}$ increases with $r$, there exists a radius from which gas is flowing outward. This outflow of gas can drag dust with sufficiently small Stokes numbers along. The transition radius where the gas velocity $u_{\mathrm{r}}$ switches from an inflow ($u_{\mathrm{r}}<0$) to an outflow ($u_{\mathrm{r}}>0$) can be computed the same way:

[TABLE]

leads with Eq. (14) to

[TABLE]

giving

[TABLE]

Since in our simulations the pressure gradient is always negative, drift is also always pointing inward. In this case, the region of outflowing solid matter in the $r$-$\mathrm{St}$ space is given by Eq. (23), giving the Stokes number of the velocity sign flip, and Eq. (47), giving the disk radius beyond which gas, and thus also these particles, are flowing outward.

Appendix C Averaged Size of planetesimal forming material

What is the average size of the material which is building planetesimals? To answer this question we use flux averaging such that

[TABLE]

where $X$ is either Stokes number, $\mathrm{St}$, or particle radius, $a$, and the sum goes over all particles with Stokes number in the range between $\mathrm{St}_{\mathrm{min}}$ and $\mathrm{St}_{\mathrm{max}}$. The Heaviside function reflects the condition (5). We use a flux average, since in our model the contribution to planetesimals scales with it. Fig. 8 shows that the flux averaged particle size building planetesimals is roughly a constant over time, as long as the dust-to-gas ratio is high enough to keep the drift limit close to the fragmentation barrier. The large peak beyond $300\,\mathrm{au}$ stems from the exponential drop in gas density. This yields planetesimals which are built from $\sim\mathrm{\umu m}$-sized particles. The $r$-dependence of $\overline{a}$ is mostly determined by the gas column density $\Sigma_{\mathrm{g}}$, as can be seen from the almost $r$-independent $\overline{\mathrm{St}}$ values at $10^{6}\,\mathrm{yr}$, which transforms $\overline{\mathrm{St}}$ via Eq. (13) into $\overline{a}\propto\Sigma_{\mathrm{g}}$.

Appendix D Necessity of $\mathrm{St}_{\mathrm{min}}$ and $\mathrm{St}_{\mathrm{max}}$

If the flux is not dominated by particles in the range $\mathrm{St}_{\mathrm{min}}\leq\mathrm{St}\leq\mathrm{St}_{\mathrm{max}}$, these parameters become very important. Though, $\mathrm{St}_{\mathrm{max}}$ is only important for the outer disk part where the gas density becomes very low, leading to high Stokes numbers $>1$ even for the smallest particles in the simulation ($0.1\,\mathrm{\umu m}$). Fluid dynamical simulations have to show which values of the Stokes number represent the boundaries at which both particle trapping and planetesimal formation can occur in a gravo-turbulent scenario. One prerequisite would be the onset of streaming instability which can further regulate planetesimal formation after enough mass has been accumulated. In Fig. 1 we plot both the particle flux fulfilling $10^{-2}=\mathrm{St}_{\mathrm{min}}\leq\mathrm{St}\leq\mathrm{St}_{\mathrm{max}}=10$ and the flux of all particles. In the extreme case where $\mathrm{St}_{\mathrm{min}}=0$ and $\mathrm{St}_{\mathrm{max}}=\infty$, the difference between the total flux (dashed lines) and the pebble flux shown (solid lines, $\mathrm{St}_{\mathrm{min}}=10^{-2}$, $\mathrm{St}_{\mathrm{max}}=10$) would also yield a difference in planetesimals. Since the particles’ contribution to planetesimals is weighted by the contribution of their flux to the overall pebble flux, smaller particles also contribute less. Smaller values of $\mathrm{St}_{\mathrm{min}}$, say $10^{-3}$, would lead to earlier planetesimal formation, which also persists longer. The edges of the planetesimal zone would also stretch out.

Appendix E Total dust and planetesimal mass in the disk

Fig. 9 shows the integrated mass of planetesimals and smaller particles. Reaching the total solid mass of the MMSN does not necessarily imply that the result is compatible with planetesimal formation in the Solar Nebula. It is more like a necessary condition and, therefore, knowledge about the planetesimal distribution as a function of radius is mandatory. For $\alpha_{\mathrm{t}}=10^{-2}$ values over the MMSN mass can be reached only for high $\varepsilon$ (here 0.8), which is unlikely since the efficiency in the planetesimal formation process itself is already lower (Schreiber, 2018, see $\kappa$ in his Table 7.3).

High $\alpha_{\mathrm{t}}$ leads to wider regions which are fragmentation limited. Hence, more mass is in small, slow drifting particles. The total mass strongly depends on $\varepsilon$. Whereas for mostly drift limited disks, the mass reservoir in the outer disk is transported to the inner part much faster. Over this long distance, particles which were not transformed into planetesimals in the outer part may turn into one in the inner part. This can lead to similar disk masses of the planetesimal population but with a different $r$-dependent distribution.

Appendix F Origin of particles beyond growth barriers

Physical and numerical diffusion mix particles radially, which can lead to particles larger than the growth barriers. They can remain there for a long time if the particle density is low enough to keep collision rates small, see figures 10 and 11. Converting particle radius into Stokes number, this effect looks more extreme in the outer disk (figures 2 and 3) because the gas density decreases with disk radius.

Appendix G More detailed derivation of the planetesimal formation condition

In order to be able to form planetesimals, there exists a critical cloud diameter $l_{\mathrm{c}}$ that has to be reached in order for the cloud to be able to contract while an underlying turbulent particle diffusion is trying to dissolve it (Schreiber, 2018, his section 3.3). This critical length scale can be derived by equating the diffusion and collapse timescale leading to

[TABLE]

Here, $\delta$ is the dimensionless diffusion parameter (diffusion constant over $c_{\mathrm{s}}h_{\mathrm{g}}$) which acts on the scale of $l_{\mathrm{c}}$. Particles have to concentrate and reach Hill density $\rho_{\mathrm{Hill}}$ in a volume $l_{\mathrm{c}}^{3}$. I. e., under self-gravity particles will settle to an effective scale height $\ell_{z}=l_{\mathrm{c}}$ in the $z$-direction and concentrate azimuthally out of a full $2\pi r$ ring to a length of $\ell_{\varphi}=l_{\mathrm{c}}$. In the radial direction the particles that have been in the trap structure of radial extent $d_{r}$ will be concentrated within the length $\ell_{r}=l_{\mathrm{c}}$, and the particles that would drift the radial distance $\left|{v_{\mathrm{peb}}(t_{0},\mathrm{St})}\right|\tau$ during the time $\tau$ are also concentrated within $\ell_{r}=l_{\mathrm{c}}$. We add the factor $4\pi/3$ because these particles are more likely to concentrate into a spherical shape than a cubic one. The condition for planetesimal formation then reads

[TABLE]

The second term will mostly dominate the first term. The expression

[TABLE]

gives the resulting planetesimal mass (Klahr & Schreiber, 2015; Schreiber, 2018) and, by neglecting the first term, we obtain the same condition as in Eq. (5).

The timescale over which the concentration will occur is assumed to be the lifetime of these structures, $\tau=\tau_{\mathrm{l}}$. Since we need a local criterion to make the model work, we have to assume that the particle flux does not change significantly over $\tau$, or if it does, that it will still lead to the same result of either reaching or not reaching Hill density.