# Planetesimal Initial Mass Functions following Diffusion Regulated
Gravitational Collapse

Konstantin Gerbig Department of Astronomy, Yale University, 52 Hillhouse Ave,
New Haven, CT 06511, USA  Rixin Li (李日新)  Center for Astrophysics and
Planetary Science, Department of Astronomy, Cornell University, Ithaca, NY
14853, USA

(Received January 27, 2023; Revised March 19, 2023; Accepted April 2, 2023)

###### Abstract

The initial mass function (IMF) of planetesimals is of key importance for
understanding the initial stages of planet formation, yet theoretical
predictions so far have been insufficient in explaining the variety of IMFs
found in simulations. Here, we connect diffusion-tidal-shear limited
planetesimal formation within the framework of a Toomre-like instability in
the particle mid-plane of a protoplanetary disk to an analytic prediction for
the planetesimal IMF. The shape of the IMF is set by the stability parameter
$Q_{\mathrm{p}}$, which in turn depends on the particle Stokes number, the
Toomre $Q$ value of the gas, the local dust concentration and the local
diffusivity. We compare our prediction to high-resolution numerical
simulations of the streaming instability and planetesimal formation via
gravitational collapse. We find that our IMF prediction agrees with numerical
results, and is consistent with both the ‘planetesimals are born big’ paradigm
and the power law description commonly found in simulations.

††journal: The Astrophysical Journal††software: Athena (Stone et al., 2008),
PLAN (Li, 2019), NumPy (Harris et al., 2020), Matplotlib (Hunter, 2007),
CMasher (van der Velden, 2020)

##  1 Introduction

Planetesimals are the initial building blocks of planets and their Initial
Mass Function (IMF) is consequently of great interest. The planetesimal
formation process in protoplanetary disks, which ultimately dictates the
planetesimal IMF, is connected to an ensemble of challenges, one of the most
prominent of which is the so-called meter barrier (see the reviews by Chiang &
Youdin, 2010; Liu & Ji, 2020; Drazkowska et al., 2022). Particles exceeding
sizes of order meters are expected to be limited in their capacity to grow via
continued coagulation due to both rapid radial drift inwards and increased
relative velocities that result in preferentially destructive collisions
(Birnstiel et al., 2012). One well-received solution to this growth barrier is
to rapidly form planetesimals via gravitational collapse of over-dense
particle clouds. Initial ideas focused on the gravitational instability of the
entire particle mid-plane (Safronov, 1969; Goldreich & Ward, 1973). However,
Kelvin-Helmholtz stirring prevents razor-thin particle settling and thus
renders such a global gravitational instability challenging (Weidenschilling,
1980; Sekiya, 1998; Youdin & Shu, 2002; Johansen et al., 2006a). On the other
hand, local patches have been shown to collapse if local particle
concentrations are sufficiently high (Johansen et al., 2006b). In particular,
instabilities energized by the relative dust-gas streaming velocity (Youdin &
Goodman, 2005; Squire & Hopkins, 2018) can concentrate particles to densities
sufficient for gravitational collapse to trigger and planetesimals to form
(e.g., Johansen et al., 2015; Gerbig et al., 2020). In the past, the IMF of
planetesimals has been obtained by performing numerical simulations of this
setup and then counting produced planetesimals, a process that resulted in
power-laws of various kinds (Simon et al., 2016; Schäfer et al., 2017; Li et
al., 2019).

Recently, the formation of planetesimals has been connected to the diffusion
of particles as well as their stability to stellar tidal forces (Klahr &
Schreiber, 2016; Gerbig et al., 2020; Klahr & Schreiber, 2020, 2021). This
framework envisions particles to be subject to a diffusive flux away from the
concentration maximum, much how pressurized gas clouds resist collapse in star
formation. In addition, any particle cloud on the verge of gravitational
collapse must be stable to stellar tidal gravity. These effects limit
planetesimal formation on small scales and large scales respectively, thus
together prescribe a characteristic scale on which planetesimal formation is
expected to occur. This in turn, can be translated to a characteristic
planetesimal mass, which Klahr & Schreiber (2020) hypothesized to be the
center of a Gaussian-shaped IMF, in agreement with IMFs of primordial
asteroids (Delbo et al., 2019), yet seemingly implying a mismatch to
numerically obtained IMFs.

In this paper, we connect these two paradigms by deriving planetesimal IMFs
within the framework of diffusion-tidal-shear limited planetesimal formation.
Thereby, we argue that the probability density function of a given scale to
collapse scales with the scale’s Toomre-like growth rate. We also directly
test our prediction by conducting numerical simulations using proven setups in
ATHENA (Stone et al., 2008), that produce the streaming instability (Youdin &
Goodman, 2005) and planetesimal formation. In the process, we for the first
time conduct diffusion measurements in large-scale stratified streaming
instability simulations, as well as develop a method for obtaining local
particle concentrations that are appropriate for characterizing the onset of
planetesimal formation.

The paper is structured as follows. In Sect. 2 we review the Toomre-like
instability for planetesimal formation. In Sect. 3, we derive planetesimal
IMFs, which we compare to numerical simulations in Sect. 4. We discuss our
results, namely applicability, caveats and implications, in Sect. 5.

##  2 Toomre-like instability

###  2.1 Stability parameter

The stability of self-gravitating particles subject to a diffusive flux
induced by coupling to turbulent gas velocities has been studied in numerous
occasions (Youdin, 2011; Gerbig et al., 2020; Klahr & Schreiber, 2020, 2021),
and can be assessed using a Toomre-like value $Q_{\mathrm{p}}$, such as

$\displaystyle
Q_{\mathrm{p}}\equiv\sqrt{\frac{\delta_{r}}{\mathrm{St}}}\frac{1}{Z}\frac{c_{\mathrm{s}}\Omega}{\pi
G\Sigma_{\mathrm{g}}}=\sqrt{\frac{\delta_{r}}{\mathrm{St}}}\frac{Q}{Z}$ (1)

where $Z=\Sigma_{\mathrm{p}}/\Sigma_{\mathrm{g}}$ is the (local) dust
concentration, $Q=c_{\mathrm{s}}\Omega/(\pi G\Sigma_{\mathrm{g}})$ is the
standard Toomre value (Toomre, 1964), $\mathrm{St}=t_{\mathrm{s}}\Omega$ is
the dimensionless stopping time, and
$\delta_{r}=D_{\mathrm{p,}r}/(c_{\mathrm{s}}H)$ is the dimensionless (radial)
diffusion coefficient for particles. $\Omega$ is the orbital frequency,
$c_{\mathrm{s}}$ the sound-speed of the gas, $H=c_{\mathrm{s}}/\Omega$ the
disk pressure scale-height, $\Sigma_{\mathrm{p}}$ and $\Sigma_{\mathrm{g}}$
particle and gas surface densities respectively, and $G$ is the gravitational
constant. In the Epstein regime, the stopping time relates to particle size
$a$, and dust material density $\rho_{\bullet}$ via
$t_{\mathrm{s}}=\rho_{\bullet}a/(\sqrt{2\pi}\Sigma_{\mathrm{g}})$. We consider
dust and pebbles where $t_{\mathrm{s}}$ is such that $\mathrm{St}$ remains
below unity (e.g., Birnstiel et al., 2012). This corresponds to well or at
least marginally coupled particles.

If $Q_{\mathrm{p}}<1$, the system is unstable and expected to collapse and
form planetesimals. We review the corresponding instability analysis in the
following section.

###  2.2 Dynamical equations and dispersion relation

Figure 1: Dispersion relation $\omega_{0}^{2}$ in Eq. (5) (left panel) and
growth rates $\gamma$ in Eq. (11) (right panel) for different values of
$Q_{\mathrm{p}}$. Parameters are $\delta=10^{-5}$ and $\mathrm{St}=0.4$. The
gray shaded region in the left panel indicates stability against axisymmetric
perturbations. As $Q_{\mathrm{p}}$ decreases, more modes become unstable, and
the fastest growing mode (dashed lines in right panel) shifts to smaller
scales. For reference, we plot the scale $2\pi\sqrt{\delta/\mathrm{St}}H$
(dotted line right panel).

Following (Klahr & Schreiber, 2021), we start with the a set of dynamical
equations for dust particles in the shearing sheet approximation and under the
assumption of a razor-thin disk. We adopt the coordinates $r,\phi,z$ for
radial, azimuthal and vertical directions respectively, and consider a patch
at distance $R$ from a solar-mass star. The linearized set of equation is
given by

$\displaystyle\frac{1}{\Sigma_{\mathrm{p}}}\frac{\partial\Sigma_{\mathrm{p}}^{\prime}}{\partial
t}+\frac{\partial v_{r}^{\prime}}{\partial r}=0,$ (2)
$\displaystyle\frac{\partial v_{r}^{\prime}}{\partial t}-2\Omega
v_{\phi}^{\prime}=-\frac{1}{\Sigma_{\mathrm{p}}}\frac{D_{\mathrm{p,r}}}{t_{\mathrm{s}}}\frac{\partial\Sigma_{\mathrm{p}}^{\prime}}{\partial
r}-\frac{\partial\Phi^{\prime}}{\partial
r}-\frac{v_{r}^{\prime}}{t_{\mathrm{s}}},$ (3) $\displaystyle\frac{\partial
v_{\phi}^{\prime}}{\partial t}+\frac{\Omega v_{r}^{\prime}}{2}=0,$ (4)

where the prime denotes perturbed quantities. They describe mass continuity,
and conservation of radial and azimuthal momenta. Notably, we ignore an
explicit azimuthal drag term as in e.g., Youdin (2011), an assumption that is
justified at dust-to-gas ratios of order unity. Instead, coupling to the gas
is assumed to be wholly described by gas pressure counteracting radial
contraction and the diffusive flux in the radial momentum equation, which can
be understood as an effective pressure flux induced by turbulent gas motions.
As a consequence, there is also no mass diffusion term in the continuity
equation.

We continue by introducing axisymmetric WKB waves scaling with
$\Sigma_{\mathrm{p}}^{\prime}\propto\exp(-i(kr-\omega t))$.
$\Phi^{\prime}=-2\pi G\Sigma^{\prime}/|k|$ is the potential for a perturbed
disk assuming $\rho_{\mathrm{p}}(k,z)=(k\Sigma_{\mathrm{p}}/2)\exp(-|k|z)$.
The resulting dispersion relation is given by (Eq. (B22) in Klahr & Schreiber,
2021),

$\displaystyle\omega^{2}_{0}=\frac{\delta_{r}}{\mathrm{St}}c_{\mathrm{s}}^{2}k^{2}-2\pi\Sigma_{\mathrm{p}}G|k|+\Omega^{2},$
(5)

and can be expressed in terms of the stability parameter $Q_{\mathrm{p}}$

$\displaystyle\frac{\omega^{2}_{0}}{\Omega^{2}}=\frac{\delta_{r}}{\mathrm{St}}(kH)^{2}-\frac{2}{Q_{\mathrm{p}}}\sqrt{\frac{\delta_{r}}{\mathrm{St}}}|k|H+1.$
(6)

Here, $\omega_{0}^{2}$ is defined as
$\omega_{0}^{2}=\omega(\omega-i/t_{\mathrm{s}})$, and represents the complex
frequency without the drag term in Eq. (3). Note that the herein used solution
to the Poisson equation does not take into account a softening term caused by
the particle layer thickness (see Eq. (12) in Youdin, 2011).

###  2.3 Fastest growing mode and growth rates.

Figure 2: Probability density functions (left panel), cumulative probability
density functions (center panel) and initial mass functions (right panel) for
the growth rates in Fig. 2 where $\mathrm{St}=0.4$. Diffusion is isotropic and
set to $\delta=10^{-6}$ (solid lines) and $\delta=10^{-5}$ (dashed lines).
Different colors correspond to different values of $Q_{\mathrm{p}}$, all of
which chosen to be unstable. We use Eq. (16) for converting unstable scale to
seed mass, which assumes collapse at Hill density, and requires knowledge of
the aspect ratio which we set to $h=0.05$. We assume the system to produce
$N=1000$ planetesimals total to calculate the normalized IMF
$N(<m_{\mathrm{p}})$.

The fastest growing mode is found by solving $\partial\omega_{0}^{2}/\partial
k=0$, which yields

$\displaystyle
k_{\mathrm{fgm}}=\sqrt{\frac{\mathrm{St}}{\delta_{r}}}\frac{1}{Q_{\mathrm{p}}H}.$
(7)

The complex frequency of the fastest growing mode is simply

$\displaystyle\frac{\omega^{2}_{0,\mathrm{fgm}}}{\Omega^{2}}=1-\frac{1}{Q_{\mathrm{p}}^{2}},$
(8)

which highlights that only for $Q_{\mathrm{p}}<1$ exist $k$ for which
$\omega_{0}^{2}<0$ and the instability can grow, thus justifying the
definition of $Q_{\mathrm{p}}$ as a stability parameter. We can also calculate
largest and smallest unstable scales $\lambda=2\pi/k$ by solving
$\omega^{2}_{0}=0$, resulting in

$\displaystyle\frac{\lambda_{\mathrm{min/max}}}{\lambda_{\mathrm{fgm}}}=\frac{1}{Q_{\mathrm{p}}}\left(\frac{1}{Q_{\mathrm{p}}}\pm\sqrt{\frac{1}{Q_{\mathrm{p}}^{2}}-1}\right),$
(9)

with $\lambda_{\mathrm{fgm}}=2\pi/k_{\mathrm{fgm}}$. The range of unstable
scales is
$\lambda_{\mathrm{max}}-\lambda_{\mathrm{min}}=2\lambda_{\mathrm{fgm}}/Q_{\mathrm{p}}\sqrt{1/Q_{\mathrm{p}}^{2}-1}$,
which increases for decreasing $Q_{\mathrm{p}}$.

Note, that the fastest growing mode relates to the critical cloud radius
$l_{\mathrm{c}}$ in Klahr & Schreiber (2020) and Gerbig et al. (2020) via

$\displaystyle\lambda_{\mathrm{fgm}}=2\pi\sqrt{\frac{\delta_{r}}{\mathrm{St}}}Q_{\mathrm{p}}H=6\pi
Q_{\mathrm{p}}l_{\mathrm{c}}.$ (10)

The growth rate $\gamma(k)$ can be found solving $\gamma=i\omega$ or
equivalently $\gamma(\gamma+1/t_{\mathrm{s}})=-\omega^{2}_{0}$ (Klahr &
Schreiber, 2021):

$\displaystyle\frac{\gamma(k)}{\Omega}=-\frac{1}{2\mathrm{St}}+\sqrt{\frac{1}{4\mathrm{St}^{2}}-\frac{\omega^{2}_{0}(k)}{\Omega^{2}}}.$
(11)

Hence, the fastest growing mode $k_{\mathrm{fgm}}$ grows with

$\displaystyle\frac{\gamma(k=k_{\mathrm{fgm}})}{\Omega}=-\frac{1}{2\mathrm{St}}+\sqrt{\frac{1}{4\mathrm{St}^{2}}+\frac{1}{Q_{\mathrm{p}}^{2}}-1}$
(12)

Dispersion relation and growth rates are shown for a set of unstable
$Q_{\mathrm{p}}$ in Fig. 1. In the right panel, we also plot
$\lambda_{\mathrm{fgm}},l_{\mathrm{c}}$ as well as the scale
$2\pi\sqrt{\delta_{r}/\mathrm{St}}H$ which is of order the radial extent of
small-scale particle structures (Gerbig et al., 2020).

##  3 Planetesimal initial mass functions

###  3.1 Planetesimal masses at Hill density

In order for a region in the particle disk to be stable against stellar tidal
gravity, its mass must be be contained within its Hill-radius. In other words,
the region’s density must be at least Hill density

$\displaystyle\rho_{\mathrm{H}}=\frac{9\Omega^{2}}{4\pi G}.$ (13)

Assuming a given scale $k=2\pi/\lambda$ is unstable under the previously
discussed dispersion relation, at Hill density $\rho_{\mathrm{H}}$, the mass
available to the produced planetesimal111The mass $m_{\mathrm{p}}$ is only
equal to the initial planetesimal mass if the entire particle cloud collapses
to material density of the planetesimal, and thus should be understood as an
approximate mass scale. As such, in Klahr & Schreiber (2020, 2021), this mass
scale is called ‘equivalent mass’, and associated with some conversion
efficiency that quantifies the fraction that ultimately ends up in the formed
planetesimal. can be estimated with

$\displaystyle
m_{\mathrm{p}}=\frac{\pi}{4}\left(\frac{\lambda}{2\pi}\right)^{2}\Sigma_{\mathrm{p}}(\rho_{\mathrm{p}}=\rho_{\mathrm{H}})$
(14)

We assumed that the seed mass has access to a region of size $\lambda/2\pi$.
The surface density if the mid-plane has Hill-density can be estimated with
$\Sigma_{\mathrm{p}}=\sqrt{2\pi}H_{\mathrm{p}}\rho_{\mathrm{p}}$, where the
particle scale height relates to the vertical diffusion coefficient via
(Youdin & Lithwick, 2007)

$\displaystyle H_{\mathrm{p}}=\sqrt{\frac{\delta_{z}}{\mathrm{St}}}H.$ (15)

Due to anisotropic diffusion where $\delta_{r}\neq\delta_{z}$, contraction
does not necessarily proceed spherically symmetric (compare to Klahr &
Schreiber, 2021).

Assuming a solar mass star, the mass associated with an unstable scale
$\lambda$ is of order,

$\displaystyle
m_{\mathrm{p}}=\frac{9}{64}\sqrt{\frac{2}{\pi^{3}}}h^{3}\sqrt{\frac{\delta_{z}}{\mathrm{St}}}\left(\frac{\lambda}{H}\right)^{2}M_{\odot},$
(16)

and equivalently, the unstable scale $\lambda$ that is expected to produce a
planetesimal of mass of order $m_{\mathrm{p}}$ is

$\displaystyle\frac{\lambda}{H}=\left(\frac{64}{9}\sqrt{\frac{\pi^{3}}{2}}\frac{1}{h^{3}}\sqrt{\frac{\mathrm{St}}{\delta_{z}}}\frac{m_{\mathrm{p}}}{M_{\odot}}\right)^{1/2}$
(17)

Here $h=H/R$ is the disk aspect ratio, typically of order $0.03<h<0.1$. The
mass associated with the fastest growing mode is given by

$\displaystyle
m_{\mathrm{p,fgm}}=\frac{9}{8}\sqrt{\frac{\pi}{2}}\sqrt{\frac{\delta_{z}}{\mathrm{St}}}\frac{\delta_{r}}{\mathrm{St}}h^{3}Q_{\mathrm{p}}^{2}M_{\odot}$
(18)

which scales as

$\displaystyle\begin{split}m_{\mathrm{p,fgm}}&\approx
0.37M_{\mathrm{Ceres}}\cdot\left(\frac{\delta_{z}}{10^{-5}}\right)^{\frac{1}{2}}\left(\frac{\delta_{r}}{10^{-5}}\right)\\\
&\cdot\left(\frac{0.1}{\mathrm{St}}\right)^{\frac{3}{2}}\left(\frac{h}{0.05}\right)^{3}\left(\frac{Q_{\mathrm{p}}}{1}\right)^{2}\end{split}$
(19)

Under isotropic diffusion, this equals the characteristic mass in Klahr &
Schreiber (2020) if $Q_{\mathrm{p}}\approx 0.28$. We can also compare this
mass scaling to Liu et al. (2020) who comprise numerical results Johansen et
al. (2015); Simon et al. (2016); Schäfer et al. (2017); Abod et al. (2019); Li
et al. (2019) into a characteristic mass which, for their fiducial parameters
and a solar mass star likewise yields $\sim 0.3M_{\mathrm{Ceres}}$ and is thus
consistent with Eq. (19).

Note, that the Toomre instability a priori does not require Hill-density to
operate. Indeed, as pointed out by Klahr & Schreiber (2021), if the vertical
density structure is not set by vertical diffusion, but instead stringently
follows $\rho_{\mathrm{p}}(k,z)=(k\Sigma_{\mathrm{p}}/2)\exp(-|k|z)$, then the
fastest growing linear instability would be achieved at a mid-plane density of
$\rho_{\mathrm{p}}(k=k_{\mathrm{fgm}},z=0)=2\rho_{\mathrm{H}}/9$. By
presupposing Hill density, our approach detaches from this assumption of mode-
dependent stratification, and in the process excludes Toomre unstable clouds
that fail to withstand tidal gravity and are thus of little physical
importance for the IMF.

###  3.2 IMF from Toomre growth rate

Given this context of Toomre-like instability and Hill density, we proceed by
providing predictions for the IMF. Our ansatz is to take the mode-depended
instability growth rates as a probability density function $p(\lambda)$ for
unstable scales, and then convert unstable modes to planetesimal masses via
the recipes discussed in Sect. 3.1. This yields a probability density function
$p(m_{\mathrm{p}})$ such that the probability of a seed mass within
$[m_{\mathrm{p}},m_{\mathrm{p}}+\mathrm{d}m_{\mathrm{p}}]$ is
$p(m_{\mathrm{p}})\mathrm{d}m_{\mathrm{p}}$. This first order approach agrees
with intuition in that fastest growing modes should most preferentially
collapse, slowly growing modes only sometimes, stable modes never.

We write the probability density function $p(\lambda)$ for scale $\lambda$ to
collapse as

$\displaystyle
p(\lambda)\propto\mathrm{max}[\hat{\gamma}^{-1}\gamma(\lambda),0]$ (20)

$\hat{\gamma}^{-1}$ is a normalization constant given by
$\hat{\gamma}=\int_{\lambda_{\mathrm{min}}}^{\lambda_{\mathrm{max}}}\gamma(\lambda)\mathrm{d}\lambda$.
Since $\lambda\propto m_{\mathrm{p}}^{1/2}$ via Eq. (17), we can map
$p(\lambda)\mathrm{d}\lambda$ onto the probability density function for seed
planetesimal masses $p(m_{\mathrm{p}})\mathrm{d}m_{\mathrm{p}}$.

The probability density function $p(m_{\mathrm{p}})$, the cumulative
probability function $P(m_{\mathrm{p}}^{\prime}\leq m_{\mathrm{p}})$, and the
resulting IMF are shown in Fig. 2, for $\mathrm{St}=0.4$ and different values
of $\delta$ and $Q_{\mathrm{p}}$. The smaller $Q_{\mathrm{p}}$, the flatter
the IMF, as more scales become unstable. The steepest IMF is achieved for
$Q_{\mathrm{p}}\rightarrow 1$. Indeed, for $Q_{\mathrm{p}}=1$, the probability
density function becomes a Dirac delta function, i.e.,
$p(m_{\mathrm{p}})=\delta(m_{\mathrm{p,H}}-m_{\mathrm{p}})$. Increasing
diffusivity shifts the IMF to larger masses, provided $Q_{\mathrm{p}}$ remains
constant, which would require a corresponding decrease in $Q/Z$.

###  3.3 Comparison to statistical approaches

Our ansatz is distinct from past means of deriving planetesimal IMFs. Cuzzi et
al. (2008, 2010); Hartlep & Cuzzi (2020) approach the problem statistically,
and consider turbulent clustering of particles. In particular, the assumed
scale invariance of the turbulent spectrum implies the statistically
appearance of regions of highly enhanced particle density. The argument is
that only sufficiently dense clumps can withstand ram pressure disruption and
thus contract to planetesimals in a process called primary accretion, thus
limiting the formation of planetesimals on the low-mass end. Our mechanism
likewise prohibits the formation of arbitrarily small planetesimals, yet the
physical intuition differs in that (1) the system is, in fact, gravitationally
unstable under the Toomre-like instability discussed in Sect. 2 resulting in
(2) ram pressure being negligible compared to diffusion (also see Klahr &
Schreiber (2020)), and (3), that the smallest planetesimals are those
resulting from the smallest scale that can be gravitationally unstable under
the dispersion relation in Eq. (5).

Another statistical ansatz was taken by Hopkins & Christiansen (2013), where
turbulent density fluctuations can render local regions of the (gas) disk
unstable to gravitational collapse. The initial mass function of collapsing
clumps therein depends on the critical density for self-gravitating clumps and
as well as the properties of the ambient turbulence. In our ansatz, the gas
remains stable throughout. Moreover, we take all clumps to collapse exactly at
Hill density in Eq. (13). Different masses are the result of collapse of
differently sized regions, i.e. on different scales $\lambda$.

##  4 Numerical Tests

We perform high-resolution numerical tests using ATHENA (Stone et al., 2008;
Bai & Stone, 2010a) to directly test our predictions for the IMF. The
procedure is similar to that in e.g. Li et al. (2019); Gerbig et al. (2020),
in that we let a small patch around a protoplanetary disk mid-plane evolve
into some turbulent state. We then turn on self-gravity with different values
for $Q_{\mathrm{p}}$ by varying the self-gravity parameter $\hat{G}$ in code
units. The measurement of $Q_{\mathrm{p}}$, as well as the calculation of the
predicted IMF requires measurement of radial and vertical diffusion, prior to
turning on-self gravity. Since the streaming instability will concentrate
particles into filaments, we must also determine the local particle
concentration.

###  4.1 Numerical Setup

Figure 3: Map of the particle surface density at after 60 $\Omega^{-1}$, which
is when self-gravity is turned on. Top panel is vertically integrated, wheares
bottom panel is azimuthally integrated. The white dashed lines mark the border
of the physical simulation domain and the ghost cells. Since, one of the two
filaments at 60 $\Omega^{-1}$ is located right at the radial simulation
boundary, we chose to include the ghost cells in this figure, which results in
that filament being depicted twice.

Figure 4: Maps of the columns that, at 60 $\Omega^{-1}$, contain at least one
cell at or above Hill-density for $Q=16$, $Q=8$, and $Q=4$ (top panels, from
left to right), and corresponding column density histograms (bottom panel).
The bottom panel also indicates the chosen values for $Z$ in dashed lines,
i.e. $Z=0.137$, $Z=0.100$, and $Z=0.0075$ for $Q=16$, $Q=8$, and $Q=4$
respectively. The full surface density map of this snapshot is shown in Fig.
3.

We employ ATHENA (Stone et al., 2008) to solve the hydrodynamic equations on
an Eulerian grid including Lagrangian super-particles (Bai & Stone, 2010a).
Our numerical setup is similar to simulations of streaming instability
regulated planetesimal formation such as Johansen et al. (2007); Simon et al.
(2016); Li et al. (2019); Gerbig et al. (2020), in that we use the local
shearing box approximation (Goldreich & Lynden-Bell, 1965) with coordinates
$(x,y,z)$. We consider a non-magnetized gas with an isothermal equation of
state. Gas is initialized in hydrostatic equilibrium. Particles are forced
into a super-Keplerian rotation by an external pressure gradient, which we
parameterize using

$\displaystyle\Pi=\frac{\eta v_{\mathrm{K}}}{c_{\mathrm{s}}},$ (21)

where $\eta=-(1/2)h^{2}\mathrm{d}\ln\rho/(\mathrm{d}\ln r)$. This is
mathematically equivalent to real disks where particles orbit Keplerian, and
gas experiences sub-Keplerian forcing (see Bai & Stone 2010a for details on
the pressure gradient implementation in ATHENA). Note, that $\Pi$ relates to
the pressure gradient parameter used in Schreiber & Klahr (2018); Gerbig et
al. (2020) simply via $\Pi=\beta/2$. Particles are initialized with a narrower
Gaussian, however at an scale height of $\eta r/2$, which is the
characteristic scale of Kelvin-Helmholtz instability in protoplanetary disks
(Gerbig et al., 2020). As such, particles are not expected to undergo
significant settling or lofting due to vertical stellar gravity or Kelvin-
Helmholtz stirring respectively.

The pressure gradient and the resulting relative velocity between particle and
gas flow energize the linear streaming instability (Youdin & Goodman, 2005;
Squire & Hopkins, 2018), which then saturates non-linearly and in the process
concentrates particles into high density regions (Johansen & Youdin, 2007).
The degree of which the streaming instability operates, and as a consequence
the strength of particle diffusion and concentration prior to gravitational
collapse and planetesimal formation, is largely set by two quantities. First,
the Stokes number $\mathrm{St}$, and second, the ratio of metallicity and
pressure gradient $Z_{0}/\Pi$ (Sekiya & Onishi, 2018). The former quantifies
the coupling of particles to the gas flow, in particular the gas turbulence.
The latter traces the ratio of dust abundance relative to gas and dust layer
scale-height, and thus maps onto the mid-plane dust-to-gas ratio. Note, that
the metallicity $Z_{0}$ is the global (as in simulation domain averaged)
particle concentration and thus not necessarily equal to the local enhancement
$Z$ we introduced in Eq. (1). We elaborate on this important difference in
Sect. 4.2.

In this work, we choose $\mathrm{St}=0.4$, $Z_{0}=0.02$ and $\Pi=0.05$. This
setup is specifically designed to provide favorable conditions for the
streaming instability and, provided $Q_{\mathrm{p}}<1$, planetesimal
formation, in order to produce a large number of planetesimals which allows
for a more statistically robust determination of the numerical IMF. Our
simulations use a computational domain of $0.2H\times 0.2H\times 0.15H$, with
a resolution of $2560/H$ (i.e., $\Delta x\approx 3.9\cdot 10^{-4}H$) and
$N_{\rm par}=2^{26}\approx 6.71\cdot 10^{7}$ particles. The vertical extent is
slightly reduced from a cube to mitigate computational costs and is still tall
enough because the particle layer remains thin all the time. Moreover, in the
vertical direction, we adopt outflow boundary conditions that are known to
reduce boundary artifacts, especially in shorter boxes (Li et al., 2018). In
the radial and azimuthal directions, the standard shearing-periodic boundary
conditions are imposed.

Following the precedent set by many past works of streaming instability, we
express the results of our scale-free simulations in the dimensionless unit
system of dynamical timescale $\Omega^{-1}$, $H$, and $\rho_{0}$. The
simulation is run for $t=60\Omega^{-1}$ without self-gravity. This allows for
sufficient amount of time in the non-linear phase of the streaming instability
to measure diffusion.

Figure. 3 shows the vertically (top panel) and azimuthally (bottom panel)
integrated particle densities at $60\Omega^{-1}$, i.e. just before self-
gravity is turned on. The streaming instability in a stratified disk collects
particles into two azimuthally elongated filaments, which are enhanced in
particle density relative to the prescribed average of $Z_{0}=0.02$.

Following the snapshot at 60 $\Omega^{-1}$ depicted in Fig. 3 we turn on self-
gravity for the particles, the strength of which is parameterized by the self-
gravity parameter $\tilde{G}\equiv 4\pi G\rho_{0}/\Omega^{2}$, which relates
to $Q$ via $Q=4/(\sqrt{2\pi}\tilde{G})$. Self-gravity is required for the
concept of Hill-density to be meaningful. Specifically, the Hill-density
depends on the self-gravity parameter via

$\displaystyle\frac{\rho_{\mathrm{H}}}{\rho_{0}}=\frac{9}{\tilde{G}}=9\sqrt{\frac{\pi}{8}}Q$
(22)

We conduct three self-gravity runs with $Q=16,Q=8$ and $Q=4$ which correspond
to Hill densities of $\rho_{\mathrm{H}}/\rho_{0}\sim
90,\rho_{\mathrm{H}}/\rho_{0}\sim 45$ and $\rho_{\mathrm{H}}/\rho_{0}\sim 23$
respectively. Note that we can also now associate a Hill-radius
$r_{\mathrm{H}}$ with a given mass $m$, i.e.,

$\displaystyle
r_{\mathrm{H}}=R\left(\frac{m}{M_{\odot}}\right)^{\frac{1}{3}}=\left(\frac{\tilde{G}m}{4\pi\rho_{0}}\right)^{\frac{1}{3}}=\left(\frac{9}{4\pi}\frac{m}{\rho_{\mathrm{H}}}\right)^{\frac{1}{3}}.$
(23)

###  4.2 Local particle concentration

The stability parameter $Q_{\mathrm{p}}$ depends on the local particle
concentration $Z=\Sigma_{\mathrm{p}}/\Sigma_{\mathrm{g}}$. This is importantly
not equal to the initial, global particle concentration $Z_{0}$, which in our
case is set to $Z_{0}=0.02$. The reason for this lies within predominanetly
radial concentration of particle surface density within two filaments as
evident in Fig. 3, where the typical particle column density exceeds $Z=0.02$
by between one and two orders of magnitude. This challenge in applying the
diffusion-limited collapse criterion to proven shearing box simulation
simulations of the streaming instability was also recognized in Gerbig et al.
(2020), where the radial extent of streaming instability filaments was related
to a radial enhancement. While this argument was sufficient to demonstrate the
general applicability of a diffusion-limited collapse criterion to
simulations, it lacks the precision to reliably predict the $Q_{\mathrm{p}}$
as it by construction cannot account for additional azimuthal enhancement
within filaments.

Because of this, in this work, we choose a different pathway and utilize the
requirement of Hill-density for planetesimal formation to measure the average
particle concentration of the cells that are expected to participate in
planetesimal formation. More precisely, in order for a given vertical column
to be counted towards the metallicity, it must contain at least one cell at or
above Hill density.

Figure. 4 shows maps of the columns that satisfy this requirement for all
three $Q$-values. Note, we are here depicting the same snapshot at
$60\Omega^{-1}$ as in Fig. 3, i.e. before self-gravity has impacted the
simulation. We are merely evaluating which columns satisfy our requirement for
a given Hill-density. A comparison to Fig. 3 confirms that this scheme
constrains the metallicity measurement to the over-dense filaments only. In
the bottom panel of Fig. 4, we show histograms of the surface density of both
the entire simulation domain in grey, as well only those columns that contain
at least one cell above Hill density. The dashed lines indicate the average
particle concentration for the three runs given, and correspond to $Z=0.137$,
$Z=0.100$, and $Z=0.075$ for $Q=16$, $Q=8$, and $Q=4$ respectively.

###  4.3 Radial diffusion

Figure 5: Evolution of an initial particle distribution at different times
$\Delta t$ for $t_{0}=18\Omega^{-1}$ (top panel) and $t_{0}=50.25\Omega^{-1}$
(bottom panel). Due to the two over-dense filaments the particle distribution
develops bi-modality — an effect not seen in unstratified streaming
instability simulations (e.g., Johansen & Youdin, 2007) Figure 6: Radial
diffusion depends on particle density. Left panel: Evolution of the variance
of the particle distribution for different densities. Colors correspond to
densities. Right panel: diffusion vs particle density. Colors correspond to
density and map onto the left panel. Different line styles indicate different
averaging times when calculating the gradient in Eq. (24). We overlay a
histogram of the density distribution of the particles, where vertical solid
lines indicate Hill density for chosen $Q$-values. Vertical dashed-dotted
lines correspond to the weighted average density
$\tilde{\rho}_{\mathrm{p}}(Q)$, which pinpoint the diffusivities to
$\delta_{r}=2.4\cdot 10^{-6}$, $\delta_{r}=2.0\cdot 10^{-6}$, and
$\delta_{r}=2.6\cdot 10^{-6}$ for $Q=16$, $Q=8$, and $Q=4$ respectively

We measure radial particle diffusion via (Youdin & Lithwick, 2007; Johansen &
Youdin, 2007; Schreiber & Klahr, 2018; Baehr et al., 2022)

$\displaystyle\frac{\delta_{r}}{c_{\mathrm{s}}H}=D_{\mathrm{p,}r}=\frac{1}{2}\frac{\partial\sigma_{r}^{2}}{\partial
t},$ (24)

The idea is, that as the turbulent state evolves particles, the underlying
diffusion will widen the variance
$\sigma_{r}^{2}=\langle|x_{i}(t)-x_{i}(t_{0})|^{2}\rangle_{i}$ of the radial
particle distribution, and as such allow for a calculation of a value for
particle diffusion. Note, that since we are investigating the streaming
instability in stratified disks, this method cannot be used to determine
vertical diffusivity which we will measure in the subsequent section.

Figure. 5 shows the evolution of the initial particle distribution from
$18\Omega^{-1}$ to $60\Omega^{-1}$, when we turn on self-gravity. Our particle
distribution depicts much stronger non-Gaussianity than Johansen & Youdin
(2007) Fig. 17, since our simulation is stratified. This leads the streaming
instability to produce as azimuthally-extended filaments, where particles
drift slower due to the enhanced local dust-to-gas ratio. This can lead to
multiple peaks in the particle distribution — in our specific case there are
two peaks as there are two dominant filaments.

If we take all particles into account when evaluating Eq. (24), we will
determine a global value for the radial diffusivity. However, the strength of
diffusion is expected to vary locally and correlate with particle density (see
e.g., Schreiber & Klahr, 2018). As such, we sort all particles by local
density and measure the evolution of the distribution up to $60\Omega^{-1}$
for each density bin, and then calculate the corresponding diffusion
coefficients via Eq. (24). Fig. 6 shows that denser regions are less diffusive
than less dense regions. More specifically, up to a dust-to-gas ratio of order
$10^{1}$, the diffusion is approximately constant at $\delta_{r}=1.3\cdot
10^{-4}$. This is a value consistent with that found in Li & Youdin (2021).

For more dense regions, the diffusion decreases significantly down to
$\delta\sim 10^{-6}$ for $\rho_{\mathrm{p}}/\rho_{0}>4\cdot 10^{2}$. This
turbophoretic behavior of the particles (Caporaloni et al., 1975; Belan et
al., 2014) makes it challenging to pinpoint an exact value of $\delta$ that is
appropriate for the determination of $Q_{\mathrm{p}}$ and the subsequent IMF.
Additionally, calculation of $\delta$ requires temporal averaging. In Fig. 6
we show multiple averaging times in order to highlight the trend of longer
averaging times yielding larger diffusivity.

In this work, we choose diffusivities by evaluating
$\delta(\rho_{\mathrm{p}})$ at the particle density weighted average

$\displaystyle\tilde{\rho}_{\mathrm{p}}(Q)=\frac{\sum_{i,\rho_{\mathrm{p},i}>\rho_{\mathrm{H}}(Q)}N_{i}\rho_{\mathrm{p},i}}{\sum_{i,\rho_{\mathrm{p},i}>\rho_{\mathrm{H}}(Q)}}.$
(25)

Note that we only count those particles in Hill-stable regions. The resulting
densities are indicated in the right panel of Fig. 6. We evaluate the radial
diffusion coefficient $\delta_{r}$ (see Eq. 24) with a second-order one-side
derivative at $t=60/\Omega$ (by supplying edge_order=2 to numpy.gradient).
This results in diffusivities of $\delta_{r}=2.4\cdot 10^{-6}$,
$\delta_{r}=2.0\cdot 10^{-6}$, and $\delta_{r}=2.6\cdot 10^{-6}$ for $Q=16$,
$Q=8$, and $Q=4$ respectively.

We acknowledge that due to the steepness of
$\delta_{x}(\rho_{\mathrm{H}}/\rho_{0})$ at high dust-to-gas ratios, as well
as the dependence on averaging time $\Delta t_{\mathrm{ave}}$, there remains a
relatively large ambiguity about the most appropriate values for $\delta_{r}$.
Indeed, we find any diffusivity within $1\cdot
10^{-6}\lesssim\delta_{x}\lesssim 1\cdot 10^{-5}$ justifiable, and note that
this full range is consistent with diffusivities obtained by Schreiber & Klahr
(2018) at similar dust-to-gas ratios (although with smaller particles).

###  4.4 Vertical diffusion

Figure 7: Map of particle scale height $H_{\mathrm{p}}$(left panel) and
vertical diffusion coefficient $\delta_{z}$(center panel) at 60 $\Omega^{-1}$.
The right panel shows the mean vertical diffusion $\overline{\delta}_{z}$ vs
the local particle concentration (compare to top panel of Fig. 3). The shaded
region characterizes $\overline{\delta}_{z}\pm\sigma_{\delta}$, where
$\sigma_{\delta}$ is the standard deviation of the vertical diffusivities at a
given $\Sigma_{\mathrm{p}}/\Sigma$. The vertical lines indicate the the
determined values for $Z$ (see Sect. 4.2, Fig. 4), which set the chosen values
for vertical diffusivities to $\delta_{z}=2.3\cdot 10^{-6}$,
$\delta_{z}=3.1\cdot 10^{-6}$, and $\delta_{z}=3.3\cdot 10^{-6}$ for $Q=16$,
$Q=8$, and $Q=4$ respectively.

While vertical diffusion is not required to calculate the stability parameter
$Q_{\mathrm{p}}$ which arises from a purely 2D consideration, it is required
to calculate planetesimal masses by setting the surface density necessary to
achieve tidal-stability. We first measure the particle scale-height in each
vertical column of grid cells by calculating the standard deviation of the
vertical particle positions, and then convert to vertical diffusivity
$\delta_{z}$ using Eq. (15). Similar procedures were also done by e.g., Bai &
Stone (2010b); Yang et al. (2018); Li & Youdin (2021) when analyzing the
vertical diffusivity within a stratified particle layer subject to the
streaming instability. The resulting maps are seen in the left and center
panel of Fig. 7. The right panel of Fig. 7 depicts the mean vertical
diffusivity plotted vs local particle concentration, where the gray shaded
region indicates one standard deviation from the mean. While there is
significant spread of diffusivity for a given particle concentration, the
vertical diffusivity behaves similar to the radial diffusivity in that it
decreases in more concentrated regions. While in the plateau region, the
diffusion anisotropy is of order $\delta_{r}/\delta_{z}\sim 10$ which is
consistent with e.g., Fig. 6 in Li & Youdin (2021), in denser regions the
diffusion is of order isotropic. Indeed, we evaluate
$\delta_{z}(\Sigma_{\mathrm{p}}/\Sigma)$ at the particle concentrations
identified in Sect. 4.2, which yields $\delta_{z}=2.3\cdot 10^{-6}$,
$\delta_{z}=3.1\cdot 10^{-6}$, and $\delta_{z}=3.3\cdot 10^{-6}$ for $Q=16$,
$Q=8$, and $Q=4$ respectively.

This relatively low, and isotropic diffusion is consistent with high particle
concentrations damping diffusion from gas, regardless of direction (Ida et
al., 2021).

###  4.5 Planetesimals and mass scalings

Table 1: Simulation parameters for the three self-gravity runs. For all simulations, $\mathrm{St}=0.4,Z_{0}=0.02$, $\Pi=0.05$, and self-gravity is turned on after $60\Omega^{-1}$. Calculated masses assume $h=0.05$ and collapse at Hill-density. Our simulations have a radial and azimuthal domain size of $0.2H$, and the smallest resolved scale is $\Delta x=3.9\times 10^{-4}H$. The first two rows are prescribed values. Rows three, four and five correspond to quantities measured at $60\Omega^{-1}$ as described in Sects. 4.2, 4.3 and 4.4. Subsequent rows are properties calculated within the diffusion-limited collapse framework outlined in Sects. 2 and 3.  $Q$ | 16 | 8 | 4
---|---|---|---
$\rho_{\mathrm{H}}/\rho_{0}$ | 90.2 | 45.1 | 22.6
$Z$ | 0.137 | 0.100 | 0.075
$\delta_{r}$ | $2.42\cdot 10^{-6}$ | $2.04\cdot 10^{-6}$ | $2.60\cdot 10^{-6}$
$\delta_{z}$ | $2.35\cdot 10^{-6}$ | $3.07\cdot 10^{-6}$ | $3.28\cdot 10^{-6}$
$Q_{\mathrm{p}}$ | 0.287 | 0.180 | 0.136
$\lambda_{\mathrm{min}}[H]$ | $2.26\cdot 10^{-3}$ | $1.29\cdot 10^{-3}$ | $1.09\cdot 10^{-3}$
$\lambda_{\mathrm{fgm}}[H]$ | $4.42\cdot 10^{-3}$ | $2.56\cdot 10^{-3}$ | $2.18\cdot 10^{-3}$
$\lambda_{\mathrm{max}}[H]$ | $1.05\cdot 10^{-1}$ | $1.56\cdot 10^{-1}$ | $2.34\cdot 10^{-1}$
$m_{\mathrm{p,min}}[M_{\mathrm{Ceres}}]$ | $1.16\cdot 10^{-4}$ | $4.32\cdot 10^{-5}$ | $3.22\cdot 10^{-5}$
$m_{\mathrm{p,fgm}}[M_{\mathrm{Ceres}}]$ | $4.44\cdot 10^{-4}$ | $1.70\cdot 10^{-4}$ | $1.27\cdot 10^{-4}$
$m_{\mathrm{p,max}}[M_{\mathrm{Ceres}}]$ | 0.253 | 0.630 | 1.47
Figure 8: Snapshots of the three self-gravity simulations at $62\Omega^{-1}$,
which is $2\Omega^{-1}$ after self-gravity is turned on. All three simulations
formed gravitationally bound clumps indiciated by red circles.

Table 1 summarizes the simulation parameters, the measured quantities
$\delta_{r}$, $\delta_{z}$ and local particle concentration $Z$, as well as
the resulting predicted properties, namely $Q_{\mathrm{p}}$ as well as
characteristic size and mass scales. All three simulations have
$Q_{\mathrm{p}}<1$ and are thus expected to be gravitationally unstable. This
is unsurprising as the simulation parameters were chosen in line with past
setups proven to produce planetesimals (e.g., Li et al., 2019). Our smallest
resolved scale is $\Delta x=3.9\cdot 10^{-4}H$, which is about an order of
magnitude smaller than $\lambda_{\mathrm{min}}$ for our simulations. We thus
expect to capture the small scales of diffusion regulated collapse in these
simulations. On the other hand, the domain size is $0.2H$, which is of order
$\lambda_{\mathrm{max}}$. Moreover, the largest scale is further limited by
the maximum extent of over-dense regions, which radially and vertically is of
order $\eta r=\Pi H$ (Gerbig et al., 2020), which for $\Pi=0.05$ is just
$5\cdot 10^{-2}H$, and thus almost an order of magnitude less than
$\lambda_{\mathrm{max}}$. Hence, we do not expect to the simulation to
collapse on the largest unstable scales and produce planetesimals on the very
high mass end.

Fig. 8 shows snapshots of all three runs at $62\Omega^{-1}$ which is
$2.0\Omega^{-1}$ after self-gravity is turned on. In addition to the practical
benefit of reducing computational costs, we chose this short time for self-
gravity to act, in order to capture a relatively pristine mass distribution
that is devoid of mergers. We used the clump finder algorithm PLAN (Li, 2019;
Li et al., 2019) to identify gravitationally bound clumps , which are produced
in all three runs. Indeed, all runs collapse rapidly enough for diffusion
measurements to be unfeasible post 60 $\Omega^{-1}$. Thus, we cannot confirm
whether or not diffusion increases in this gravo-turbulent state as it does in
Klahr & Schreiber (2021).

Note, that the requirement of Hill-density imposes a physical mass unit onto
the analytic prediciton, and only implicitly depends on distance to the star
(see Klahr & Schreiber, 2020, for a discussion on the radial dependence of the
diffusion-limited collapse criterion). On the other hand, the numerical
results, i.e. the densities of bound clumps in Fig. 8 are fully scale-free in
gas surface density $\rho_{\mathrm{0}}$. However, the two can be connected
using the self-gravity parameter. More specifically, the mass unit of the
simulation $M_{0}$ can be expressed as

$\displaystyle M_{0}$
$\displaystyle=\rho_{0}H^{3}=h^{3}\tilde{G}\frac{M_{\odot}}{4\pi}=\frac{1}{4\sqrt{2\pi^{3}}}\frac{h^{3}}{Q}M_{\odot}$
(26) $\displaystyle=5.21\cdot
10^{2}\left(\frac{Q}{16}\right)^{-1}\left(\frac{h}{0.05}\right)^{3}M_{\mathrm{Ceres}},$
(27)

which introduces the same cubic dependence on aspect ratio as the analytic
predictions in e.g., Eq. (16). The aspect ratio relates to the pressure
gradient parameter via

$\displaystyle h=\left(-2\eta\frac{\mathrm{d}\ln\rho}{\mathrm{d}\ln
r}\right)^{1/2}=-2\Pi\frac{\mathrm{d}\ln\rho}{\mathrm{d}\ln r}.$ (28)

Throughout this work, we choose $h=0.05$, which together with our simulation
parameter of $\Pi=0.05$ implicitly assumes $\mathrm{d}\ln\rho/\mathrm{d}\ln
r=-1/2$. These choices are relatively generic for disk models (compare to
e.g., Dullemond et al., 2007; Gerbig et al., 2019; Gerbig & Laughlin, 2022).

The smallest mass above which a planetesimal is well resolved is the mass
associated with a Hill-radius that equals the cell-size, i.e., $m_{r_{\rm
H}=\Delta x}$. This limitation is caused by the accuracy of the self-gravity
solver set by the finite grid scale (Simon et al., 2016). PLAN automatically
discards clumps, if any, less massive than $m_{r_{\rm H}=\Delta x}$. Using Eq.
(23) and the mass scaling in Eq. (26), this mass limit is independent of $Q$,
i.e.

$\displaystyle\begin{split}m_{r_{\rm H}=\Delta
x}&=\frac{9}{4\pi}\rho_{\mathrm{H}}\Delta
x^{3}=\frac{9}{4\pi}\left(\frac{\rho_{\mathrm{H}}}{\rho_{0}}\right)\left(\frac{\Delta
x}{H}\right)^{3}M_{0}\\\ &=\frac{81}{64}h^{3}\left(\frac{\Delta
x}{H}\right)^{3}M_{\odot}=1.97\cdot 10^{-5}M_{\mathrm{Ceres}}\end{split}$ (29)

for $\Delta x=3.9\cdot 10^{-4}H$. As shown in Table 1, $m_{r_{\rm H}=\Delta
x}$ is smaller than $m_{\mathrm{p,min}}$, and $\lambda_{\rm min}$ is larger
than $2\Delta x$ for all $Q$, suggesting we sufficiently resolve the smallest
expected planetesimals.

###  4.6 IMFs and K-S test

Figure 9: Initial mass functions $N(>m_{\mathrm{p}})$ for the three self-
gravity runs. The numerically obtained IMFs are shown in solid lines and are
based on the snapshots seen in Fig. 8. The IMF resulting from the growth rate-
based PDF are shown in dashed lines. The dark and light shaded regions
indicate one and three standard deviations from the distribution of IMFs that
originates from only drawing $N\sim 100$ planetesimals from the distribution.
The vertical dotted lines indicate the masses associated with the fastest
growing mode $m_{\mathrm{p,fgm}}$. Lastly, the $p$-value from the K-S-Test is
shown in the top right corner of each panel. Figure 10: Differential initial
mass functions for the three self-gravity runs (compare to Figs. 8 and 9). We
show the number of clumps formed in each logarithmic mass bin as in similar
figures in Li et al. (2019) in comparison to the analytical PDF we calculate
(solid lines). The mass associated with the fastest growing mode (vertical
dotted line) approximately corresponds to the turnover in the differential
initial mass function.

Figure 9 shows the IMFs obtained from the snapshot depicted in Fig. 8 in solid
lines. We compare these numerically obtained IMFs, to the predicted IMF from
the Toomre-like instability paradigm, evaluated for the properties shown in
Tab. 1. Hereby, we draw $N$ planetesimals from the theoretical PDF, where $N$
is the number of bound clumps found in the respective simulation. We conduct
this random draw 1000 times in order to obtain an average theoretical IMF
(dashed curves), together with one and three $\sigma$ intervals (dark and
light shaded regions respectively). In addition, we draw a dotted vertical
line, corresponding to the mass associated with the fastest growing mode of
the IMF prediction.

In order to assess the quality of the analytical predictions, we perform two
sample Kolmogorov-Smirnov (K-S) tests, which quantify the likelihood $p$ of
two samples (in our case the numerical and analytically obtained planetesimal
masses) originating from the same underlying distribution. The null
hypothesis, that is the two samples indeed coming from the same distribution,
can be rejected if $p<0.05$. For the growth-rate based PDF, the K-S test
yields $p$-values of $0.501,0.395$, and $0.0629$ for $Q=16,Q=8$, and $Q=4$
respectively. The numerical IMFs are therefore consistent with the analytical
IMFs obtained from the Toomre-like analysis.

In Fig. 10, we compare the differential mass distribution for our three runs
to the analytical prediction for the underlying PDF. The turnover of the
differential mass distributions approximately corresponds to the masses
associated with the fastest growing mode, i.e. the peaks of the PDFs. For
$Q=4$, the distribution is bimodal with another, smaller peak toward the
higher mass end, presumably due to the stronger self-gravity where large clump
mergers may already be happening when the planetesimals were counted. This
deviation is reflected in the relative low $p$-value for $Q=4$ (see Fig. 9).

##  5 Discussion

We developed a framework of analytically obtaining planetesimal IMFs, based on
the Toomre-like instability in the particle layer in protoplanetary disks. In
Sect. 4 we presented a simulation of the streaming instability, measured
diffusivity and particle concentration for three self-gravity runs, and then
compared resulting $Q_{\mathrm{p}}$-values and corresponding analytical IMFs
to the numerically obtained IMFs. We find that the IMFs are consistent with
each other, as seen Fig. 9. Specifically, our prediction can explain a number
of key properties of the planetesimal IMF, which we will outline in the
following.

Our IMF prediction explains the somewhat counter-intuitive property of the
least massive clumps being present in the simulation with the most mass, i.e.
$Q=4$. Since this simulation has the smallest $Q_{\mathrm{p}}$ with
$Q_{\mathrm{p}}=0.14$, the fastest growing mode shifts to smaller scales (see
Fig. 1) and thus the most likely planetesimal becomes less massive.

Next, we observe that the IMF becomes flatter if $Q$ and therefore
$Q_{\mathrm{p}}$ decreases, which is due to the increase of the range of
unstable scales. Notably, all three numerical IMFs fall short of the
analytical prediction on the high mass end. As alluded to in Sect. 4.5, this
is expected as these planetesimals would have to form from scales that exceed
the maximum scale of particle over-densities $\eta r$ (see Gerbig et al.,
2020). Indeed, the idea that the most massive planetesimal is not set by the
largest Toomre-scale $\lambda$ but by $\eta r$, is consistent with our
simulations where the most massive clump has $\sim 0.1M_{\mathrm{Ceres}}$ in
all three simulations, and also with Fig. 10 in Simon et al. (2016), where the
most massive planetesimal shortly after self-gravity is turned on is
relatively independent of $\tilde{G}$ (or equivalently $Q$).

Lastly, the fact that $\lambda_{\rm fgm}$ is much closer to $\lambda_{\rm
min}$ when $Q_{\rm p}$ is small, provides an explanation for why only high-
resolution simulations can sufficiently capture the turnover in the mass
frequency distribution at the low-mass end of planetesimal IMFs (Li et al.,
2019).

###  5.1 Caveats

In this section, we outline a number of caveats that are to be kept in mind
when relating the herein presented theory and analysis to real disks.

First, unstable scales $\lambda$ were converted to planetesimal masses by
assuming the contraction of a circular region of radius $\lambda/(2\pi)$ at
Hill-density. While such a conversion agrees with timescale arguments
presented in e.g. Gerbig et al. (2020); Klahr & Schreiber (2020), it certainly
is an order of magnitude calculation. In particular, as pointed out in Polak &
Klahr (2022), a sphere at Hill-density is not fully tidally stable. The Hill-
sphere is derived by considering the tidal forces on either end of the Roche
lobe in the restricted three-body potential. This has the Hill-sphere extend
beyond the tidally-stable Roche lobe, and implies Hill-density is insufficient
for tidal stability by an order of unity factor. Moreover, it is not the case
that all unstable scales are at Hill-density exactly. In fact, as evidenced by
Fig. 6, many particles are located within regions which exceed Hill-density by
up to an order of magnitude.

Additionally, our model assumed planetesimal formation to be perfectly
efficient. Fig. 8 however demonstrates that there are still particles in the
filaments that are not bound to any clump. Moreover, previous works showed
that higher resolution simulations may extend the mass distribution to lower
masses (Li et al., 2019). Another possibility that is not considered in our
model is a single clump forming binary planetesimals upon contraction to
material density (Nesvorný et al., 2019, 2021).

We note, that while multiplicity and inefficient contraction would shift the
IMF towards smaller masses, collapse beyond Hill-density and double counting
would shift the IMF towards larger masses. We expect that this allows the
numerically obtained IMFs to remain consistent with our analytical
predictions.

Other simplifying assumptions include gas density being considered constant.
For stable, yet small $Q$ values, the gas develops a mid-plane cusp with a
non-neglible vertical density gradient (see e.g., Armitage, 2015, for a
review). Our simulations are conducted in quiescent disks without an external
source of turbulence in the gas phase. Gas turbulence can increase particle
diffusivity which is manifested in turbulent thickening of the particle mid-
plane (Gole et al., 2020). This increases planetesimal masses (see Eq. (16)),
but at the same time, imposes higher demands on the particle concentration $Z$
for achieving an unstable $Q_{\mathrm{p}}<1$ in the first place. We also
assumed a mono-disperse dust population with a constant Stokes number. This
choice is likely reasonable to first order as in real disks, most solid mass
is contained in the largest grains (e.g., Birnstiel et al., 2011; McNally et
al., 2021), and the largest grains dominate dust-gas interactions in the mid-
plane and participate most vigorously in planetesimal formation (Yang & Zhu,
2021). Still, a multi-species dust fluid can alter the turbulent behavior of
the system itself (see e.g., Krapp et al., 2019; Schaffer et al., 2021, who
pointed out the damping effect multi-species dust can have on the streaming
instability). Therefore, the assumption of single-species dust is to be kept
in mind when applying our results to real disks.

Lastly, we acknowledge that our numerical setup does not perfectly replicate
the assumptions the Toomre-like instability is based on. More specifically,
the Toomre analysis, assumes a thin, two-dimensional disk of constant density
which is then linearly perturbed in form of axisymmetric waves. This picture
evidently differs from the turbulent state, i.e. non-linear streaming
instability in vertically-stratified disks, our simulations are in when self-
gravity is turned on. We attempted to mitigate this discrepancy by evaluating
diffusivities and concentration locally, rather than globally, which, as
discussed in Sect. 4.6, yielded consistent IMFs.

###  5.2 The shape of the initial mass function in real disks

The diffusion-tidal-shear collapse paradigm suggests that the shape of
planetesimal IMF depends on the gravitational, Toomre-like instability and on
the stability parameter $Q_{\mathrm{p}}$ only. While independent information
on diffusivities $\delta_{r}$ and $\delta_{z}$, as well as Stokes number
$\mathrm{St}$, and disk aspect ratio $h$ are required to match the IMF to the
appropriate mass range, they do not affect the steepness of the IMF (for a
fixed $Q_{\mathrm{p}}$). As such, as a disk region collects more massive
particles, it experiences an increase in average grain size, cools or becomes
less turbulent, its $Q_{\mathrm{p}}$-value decreases, thus flattening the
planetesimal IMF. Within this context, the role of the streaming instability
(or other processes affecting particle dynamics) for the planetesimal IMF in
this work, is to set the initial condition at the onset of gravitational
collapse. The key question within our framework therefore is to determine what
value for $Q_{\mathrm{p}}$ takes on in real disks.

Our simulations, much like many in previous works (e.g., Simon et al., 2016;
Schäfer et al., 2017; Li et al., 2019; Gerbig et al., 2020) show collapse as
soon as self-gravity is turned on. The turbulent state before self-gravity is
turned on is therefore not physical, and thus the $Q_{\mathrm{p}}$-values and
IMFs we calculate not necessarily reflective of real disks. Indeed, as
discussed in Gerbig et al. (2020), if a system softly evolves towards
instability — for example by collecting more particle mass — its
$Q_{\mathrm{p}}$-value likewise softly evolve from $Q_{\mathrm{p}}>1$ to
$Q_{\mathrm{p}}<1$. If $Q_{\mathrm{p}}\sim 1$ is sufficient for the system to
fully collapse and form planetesimals as assumed in Klahr & Schreiber (2020,
2021), then the resulting IMF is expected to be rather steep (Polak & Klahr,
2022). On the other hand, the system may be unable to fully contract at
$Q_{\mathrm{p}}\sim 1$ if it is below Hill-density. To a similar effect, as
the contraction time, given by (Gerbig et al., 2020),

$\displaystyle t_{\mathrm{c}}=\frac{\Omega}{4\pi
G\rho_{\mathrm{p}}\mathrm{St}}=\sqrt{\frac{\pi}{8}}\sqrt{\frac{\delta_{z}}{\delta_{r}}}\frac{Q_{\mathrm{p}}}{\mathrm{St}}\Omega^{-1},$
(30)

only falls below the dynamical timescale $\Omega^{-1}$ once $Q_{\mathrm{p}}$
is substantially less than unity, it is plausible that the system can collect
relevant amounts of mass during contraction. Both situations would lead to
smaller values for $Q_{\mathrm{p}}$ and thus flatter IMFs like the ones in our
simulations in Fig. 9.

One observational constraint may be provided by the Asteroids size
distribution, which, at least in part, constitutes pristine planetesimals from
the Solar System’s formation era (see Klahr et al., 2022, for a review).
Taking, the collisional evolution into account, it is possible to reconstruct
model of the primordial size distribution of Asteroids (see e.g., Delbo et
al., 2019). Within the here presented framework, the “knee” at $\sim 100$ km
(e.g., Gladman et al., 2001) typically found in the resulting mass functions
can be interpreted as an indication of a marginally unstable origin system
with $Q_{\mathrm{p}}$ close to unity — that is assuming the asteroids in
question indeed formed via gravitational collapse of locally over-dense
regions, which very well may not have been the case. Cold Classical Kuiper
Belt Objects (KBOs) are believed to be even more primordial with little to no
collisional evolution (Morbidelli & Nesvorný, 2020). Kavelaars et al. (2021)
found an exponentially tapered power law as the best fit for the Cold
Classical KBO mass distribution, and connected the lack of large planetesimals
($>$ 400 km in size) to streaming instability regulated planetesimal
formation. Our work is consistent with this interpretation, as
$Q_{\mathrm{p}}$ imposes a limit on the maximum planetesimal mass.

One often invoked pathway of forming planetesimals and circumventing the meter
barrier, that is worth discussing in the context of our work, is the existence
of long-lived pressure bumps associated with disk sub-structures. While such
structures are robustly confirmed observationally (e.g., Andrews et al.,
2018), their role in producing the first generation of planetesimals is
elusive, as already existing planets seem to most promisingly explain the sub-
structures in the first place (Teague et al., 2021). Either way, a pressure
bump is characterized by very small pressure gradients, which diminishes the
relative velocity between dust and gas. If $\Pi=0$, then particles can
trivially settle razor-thin and go gravitationally unstable (as seen in e.g.,
Abod et al., 2019). In this state, both vertical and radial diffusivities are
expected to trend towards zero, implying $Q_{\mathrm{p}}\rightarrow 0$ as
well, which our work associates with the planetesimal IMF shifting towards
very small masses. Unless hierarchical mergers between small planetesimal
immediately after formation is very efficient, such a bottom-heavy mass
function seems unlikely. We therefore suspect that either the majority of
planetesimal formation does not occur in $\Pi=0$ regions; or, that other
mechanisms, possibly gravito-turbulence (see e.g., Riols et al., 2017, for the
gas disk), provide a lower bound on diffusivities, thus also setting a lower
limit on $Q_{\mathrm{p}}$.

To conclude, our analysis that connects the PDF of formed planetesimals to the
growth rates of the Toomre-like instability of a particle layer subject to a
diffusive flux, unifies both the flat, power-law shaped, IMFs previously
obtained numerically (e.g., Simon et al., 2016; Li et al., 2019), and the
‘Asteroids are born big’ (Morbidelli et al., 2009) paradigm that arises when
investigating marginally unstable systems (Klahr & Schreiber, 2020). The
analytically obtained IMFs are consistent with our numerical setups. Further
work is required to test the predictions for different numerically setups, and
in particular, to assess how the value of $Q_{\mathrm{p}}$ at the on-set of
planetesimal formation, and thus the steepness of the resulting IMF, depend on
disk properties and radii. Such an analysis informs initial conditions for
planet formation models (e.g., Emsenhuber et al., 2021), and conversely, if
the formation of an observed exoplanetary population presupposes a specific
planetesimal IMF(e.g., Batygin & Morbidelli, 2023, for rocky Super-Earths),
provides constraints on the disk state during the era of planetesimal
formation.

##  6 Acknowledgements

We thank the referee for insightful comments that improved the quality of the
manuscript. This work and KG and RL benefited from the 2022 Exoplanet Summer
Program in the Other Worlds Laboratory (OWL) at the University of California,
Santa Cruz, a program funded by the Heising-Simons Foundation. The authors
thank Greg Laughlin, Hubert Klahr, Andrew Youdin, Ruth Murray-Clay and Malena
Rice for insightful comments and discussions.

## References

  * Abod et al. (2019) Abod, C. P., Simon, J. B., Li, R., et al. 2019, ApJ, 883, 192, doi: 10.3847/1538-4357/ab40a3
  * Andrews et al. (2018) Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJ, 869, L41, doi: 10.3847/2041-8213/aaf741
  * Armitage (2015) Armitage, P. J. 2015, arXiv e-prints, arXiv:1509.06382.  https://arxiv.org/abs/1509.06382
  * Baehr et al. (2022) Baehr, H., Zhu, Z., & Yang, C.-C. 2022, ApJ, 933, 100, doi: 10.3847/1538-4357/ac7228
  * Bai & Stone (2010a) Bai, X.-N., & Stone, J. M. 2010a, ApJS, 190, 297, doi: 10.1088/0067-0049/190/2/297
  * Bai & Stone (2010b) —. 2010b, ApJ, 722, 1437, doi: 10.1088/0004-637X/722/2/1437
  * Batygin & Morbidelli (2023) Batygin, K., & Morbidelli, A. 2023, arXiv e-prints, arXiv:2301.04680.  https://arxiv.org/abs/2301.04680
  * Belan et al. (2014) Belan, S., Fouxon, I., & Falkovich, G. 2014, Phys. Rev. Lett., 112, 234502, doi: 10.1103/PhysRevLett.112.234502
  * Birnstiel et al. (2012) Birnstiel, T., Klahr, H., & Ercolano, B. 2012, A&A, 539, A148, doi: 10.1051/0004-6361/201118136
  * Birnstiel et al. (2011) Birnstiel, T., Ormel, C. W., & Dullemond, C. P. 2011, A&A, 525, A11, doi: 10.1051/0004-6361/201015228
  * Caporaloni et al. (1975) Caporaloni, M., Tampieri, F., Trombetti, F., & Vittori, O. 1975, Journal of Atmospheric Sciences, 32, 565, doi: 10.1175/1520-0469(1975)032<0565:TOPINA>2.0.CO;2
  * Chiang & Youdin (2010) Chiang, E., & Youdin, A. N. 2010, Annual Review of Earth and Planetary Sciences, 38, 493, doi: 10.1146/annurev-earth-040809-152513
  * Cuzzi et al. (2010) Cuzzi, J. N., Hogan, R. C., & Bottke, W. F. 2010, Icarus, 208, 518, doi: 10.1016/j.icarus.2010.03.005
  * Cuzzi et al. (2008) Cuzzi, J. N., Hogan, R. C., & Shariff, K. 2008, ApJ, 687, 1432, doi: 10.1086/591239
  * Delbo et al. (2019) Delbo, M., Avdellidou, C., & Morbidelli, A. 2019, A&A, 624, A69, doi: 10.1051/0004-6361/201834745
  * Drazkowska et al. (2022) Drazkowska, J., Bitsch, B., Lambrechts, M., et al. 2022, arXiv e-prints, arXiv:2203.09759, doi: 10.48550/arXiv.2203.09759
  * Dullemond et al. (2007) Dullemond, C. P., Hollenbach, D., Kamp, I., & D’Alessio, P. 2007, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil, 555.  https://arxiv.org/abs/astro-ph/0602619
  * Emsenhuber et al. (2021) Emsenhuber, A., Mordasini, C., Burn, R., et al. 2021, A&A, 656, A69, doi: 10.1051/0004-6361/202038553
  * Gerbig & Laughlin (2022) Gerbig, K., & Laughlin, G. 2022, ApJ, 930, 68, doi: 10.3847/1538-4357/ac6500
  * Gerbig et al. (2019) Gerbig, K., Lenz, C. T., & Klahr, H. 2019, A&A, 629, A116, doi: 10.1051/0004-6361/201935278
  * Gerbig et al. (2020) Gerbig, K., Murray-Clay, R. A., Klahr, H., & Baehr, H. 2020, ApJ, 895, 91, doi: 10.3847/1538-4357/ab8d37
  * Gladman et al. (2001) Gladman, B., Kavelaars, J. J., Petit, J.-M., et al. 2001, AJ, 122, 1051, doi: 10.1086/322080
  * Goldreich & Lynden-Bell (1965) Goldreich, P., & Lynden-Bell, D. 1965, MNRAS, 130, 125, doi: 10.1093/mnras/130.2.125
  * Goldreich & Ward (1973) Goldreich, P., & Ward, W. R. 1973, ApJ, 183, 1051, doi: 10.1086/152291
  * Gole et al. (2020) Gole, D. A., Simon, J. B., Li, R., Youdin, A. N., & Armitage, P. J. 2020, ApJ, 904, 132, doi: 10.3847/1538-4357/abc334
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Hartlep & Cuzzi (2020) Hartlep, T., & Cuzzi, J. N. 2020, ApJ, 892, 120, doi: 10.3847/1538-4357/ab76c3
  * Hopkins & Christiansen (2013) Hopkins, P. F., & Christiansen, J. L. 2013, ApJ, 776, 48, doi: 10.1088/0004-637X/776/1/48
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Ida et al. (2021) Ida, S., Guillot, T., Hyodo, R., Okuzumi, S., & Youdin, A. N. 2021, A&A, 646, A13, doi: 10.1051/0004-6361/202039705
  * Johansen et al. (2006a) Johansen, A., Henning, T., & Klahr, H. 2006a, ApJ, 643, 1219, doi: 10.1086/502968
  * Johansen et al. (2006b) Johansen, A., Klahr, H., & Henning, T. 2006b, ApJ, 636, 1121, doi: 10.1086/498078
  * Johansen et al. (2015) Johansen, A., Mac Low, M.-M., Lacerda, P., & Bizzarro, M. 2015, Science Advances, 1, 1500109, doi: 10.1126/sciadv.1500109
  * Johansen et al. (2007) Johansen, A., Oishi, J. S., Mac Low, M.-M., et al. 2007, Nature, 448, 1022, doi: 10.1038/nature06086
  * Johansen & Youdin (2007) Johansen, A., & Youdin, A. 2007, ApJ, 662, 627, doi: 10.1086/516730
  * Kavelaars et al. (2021) Kavelaars, J. J., Petit, J.-M., Gladman, B., et al. 2021, ApJ, 920, L28, doi: 10.3847/2041-8213/ac2c72
  * Klahr et al. (2022) Klahr, H., Delbo, M., & Gerbig, K. 2022, in Vesta and Ceres. Insights from the Dawn Mission for the Origin of the Solar System (Cambridge University Press), 199, doi: 10.1017/9781108856324.017
  * Klahr & Schreiber (2016) Klahr, H., & Schreiber, A. 2016, in Asteroids: New Observations, New Models, ed. S. R. Chesley, A. Morbidelli, R. Jedicke, & D. Farnocchia, Vol. 318, 1–8, doi: 10.1017/S1743921315010406
  * Klahr & Schreiber (2020) Klahr, H., & Schreiber, A. 2020, ApJ, 901, 54, doi: 10.3847/1538-4357/abac58
  * Klahr & Schreiber (2021) —. 2021, ApJ, 911, 9, doi: 10.3847/1538-4357/abca9b
  * Krapp et al. (2019) Krapp, L., Benítez-Llambay, P., Gressel, O., & Pessah, M. E. 2019, ApJ, 878, L30, doi: 10.3847/2041-8213/ab2596
  * Li (2019) Li, R. 2019, PLAN: A Clump-finder for Planetesimal Formation Simulations, Astrophysics Source Code Library, record ascl:1911.001.  http://ascl.net/1911.001
  * Li & Youdin (2021) Li, R., & Youdin, A. N. 2021, ApJ, 919, 107, doi: 10.3847/1538-4357/ac0e9f
  * Li et al. (2018) Li, R., Youdin, A. N., & Simon, J. B. 2018, ApJ, 862, 14, doi: 10.3847/1538-4357/aaca99
  * Li et al. (2019) —. 2019, ApJ, 885, 69, doi: 10.3847/1538-4357/ab480d
  * Liu & Ji (2020) Liu, B., & Ji, J. 2020, Research in Astronomy and Astrophysics, 20, 164, doi: 10.1088/1674-4527/20/10/164
  * Liu et al. (2020) Liu, B., Lambrechts, M., Johansen, A., Pascucci, I., & Henning, T. 2020, A&A, 638, A88, doi: 10.1051/0004-6361/202037720
  * McNally et al. (2021) McNally, C. P., Lovascio, F., & Paardekooper, S.-J. 2021, MNRAS, 502, 1469, doi: 10.1093/mnras/stab112
  * Morbidelli et al. (2009) Morbidelli, A., Bottke, W. F., Nesvorný, D., & Levison, H. F. 2009, Icarus, 204, 558, doi: 10.1016/j.icarus.2009.07.011
  * Morbidelli & Nesvorný (2020) Morbidelli, A., & Nesvorný, D. 2020, in The Trans-Neptunian Solar System, ed. D. Prialnik, M. A. Barucci, & L. Young, 25–59, doi: 10.1016/B978-0-12-816490-7.00002-3
  * Nesvorný et al. (2021) Nesvorný, D., Li, R., Simon, J. B., et al. 2021, , 2, 27, doi: 10.3847/PSJ/abd858
  * Nesvorný et al. (2019) Nesvorný, D., Li, R., Youdin, A. N., Simon, J. B., & Grundy, W. M. 2019, Nature Astronomy, 3, 808, doi: 10.1038/s41550-019-0806-z
  * Polak & Klahr (2022) Polak, B., & Klahr, H. 2022, arXiv e-prints, arXiv:2211.13318.  https://arxiv.org/abs/2211.13318
  * Riols et al. (2017) Riols, A., Latter, H., & Paardekooper, S. J. 2017, MNRAS, 471, 317, doi: 10.1093/mnras/stx1548
  * Safronov (1969) Safronov, V. S. 1969, Evolution of the protoplanetary cloud and formation of the earth and planets (Nauka Press)
  * Schäfer et al. (2017) Schäfer, U., Yang, C.-C., & Johansen, A. 2017, A&A, 597, A69, doi: 10.1051/0004-6361/201629561
  * Schaffer et al. (2021) Schaffer, N., Johansen, A., & Lambrechts, M. 2021, A&A, 653, A14, doi: 10.1051/0004-6361/202140690
  * Schreiber & Klahr (2018) Schreiber, A., & Klahr, H. 2018, ApJ, 861, 47, doi: 10.3847/1538-4357/aac3d4
  * Sekiya (1998) Sekiya, M. 1998, Icarus, 133, 298 , doi: https://doi.org/10.1006/icar.1998.5933
  * Sekiya & Onishi (2018) Sekiya, M., & Onishi, I. K. 2018, ApJ, 860, 140, doi: 10.3847/1538-4357/aac4a7
  * Simon et al. (2016) Simon, J. B., Armitage, P. J., Li, R., & Youdin, A. N. 2016, ApJ, 822, 55, doi: 10.3847/0004-637X/822/1/55
  * Squire & Hopkins (2018) Squire, J., & Hopkins, P. F. 2018, MNRAS, 477, 5011, doi: 10.1093/mnras/sty854
  * Stone et al. (2008) Stone, J. M., Gardiner, T. A., Teuben, P., Hawley, J. F., & Simon, J. B. 2008, ApJS, 178, 137, doi: 10.1086/588755
  * Teague et al. (2021) Teague, R., Bae, J., Aikawa, Y., et al. 2021, ApJS, 257, 18, doi: 10.3847/1538-4365/ac1438
  * Toomre (1964) Toomre, A. 1964, ApJ, 139, 1217, doi: 10.1086/147861
  * van der Velden (2020) van der Velden, E. 2020, The Journal of Open Source Software, 5, 2004, doi: 10.21105/joss.02004
  * Weidenschilling (1980) Weidenschilling, S. J. 1980, Icarus, 44, 172, doi: 10.1016/0019-1035(80)90064-0
  * Yang et al. (2018) Yang, C.-C., Mac Low, M.-M., & Johansen, A. 2018, ApJ, 868, 27, doi: 10.3847/1538-4357/aae7d4
  * Yang & Zhu (2021) Yang, C.-C., & Zhu, Z. 2021, MNRAS, 508, 5538, doi: 10.1093/mnras/stab2959
  * Youdin (2011) Youdin, A. N. 2011, ApJ, 731, 99, doi: 10.1088/0004-637X/731/2/99
  * Youdin & Goodman (2005) Youdin, A. N., & Goodman, J. 2005, The Astrophysical Journal, 620, 459, doi: 10.1086/426895
  * Youdin & Lithwick (2007) Youdin, A. N., & Lithwick, Y. 2007, Icarus, 192, 588, doi: 10.1016/j.icarus.2007.07.012
  * Youdin & Shu (2002) Youdin, A. N., & Shu, F. H. 2002, ApJ, 580, 494, doi: 10.1086/343109