Demographics of Planetesimals Formed by the Streaming Instability

TL;DR

This study uses high-resolution simulations of the streaming instability with self-gravity to analyze the initial mass distribution of planetesimals, revealing non-universality and features like a low-mass turnover consistent with gravitational collapse.

Contribution

It introduces a new high-resolution simulation approach with a novel clump-finding tool and applies rigorous model selection to characterize planetesimal mass distributions.

Findings

Mass distribution is non-universal and varies with model parameters.

Fits suggest a steeper distribution above 100 km, similar to Kuiper Belt observations.

Evidence of a low-mass turnover in the distribution, indicating gravitational collapse effects.

Abstract

The streaming instability (SI) is a mechanism to aerodynamically concentrate solids in protoplanetary disks and facilitate the formation of planetesimals. Recent numerical modeling efforts have demonstrated the increasing complexity of the initial mass distribution of planetesimals. To better constrain this distribution, we conduct SI simulations including the self-gravity with hitherto the highest resolution. To subsequently identify all of the self-bound clumps, we develop a new clump-finding tool, PLanetesimal ANalyzer (\texttt{PLAN}). We then apply a maximum likelihood estimator to fit a suite of parameterized models with different levels of complexity to the simulated mass distribution. To determine which models are best-fitting and statistically robust, we apply three model selection criteria with different complexity penalties. We find that the initial mass distribution of clumps is not universal regarding both the functional forms and parameter values. Our model selection criteria prefer models different from those previously considered in the literature. Fits to multi-segment power law models break to a steeper distribution above masses close to 100 km collapsed planetesimals, similar to observed Kuiper Belt size distributions. We find evidence for a turnover in the low mass end of the planetesimal mass distribution in our high resolution run. Such a turnover is expected for gravitational collapse, but had not previously been reported.

Tables

Table 1. Table 1: Simulation Parameters

Run	Domain Size	Number of Cells	$N_{\mathrm{par}}$$*$$*$The number of particles. For reference, $2^{27}={512}^3=134,217,728$.	${\tau }_s$	$Z$	$t_0$${\dagger}$${\dagger}$The time when the particle self-graivty is switched on.	$N_{\mathrm{tot}}$${\ddagger}$${\ddagger}$The number of clumps identified by PLAN at the snapshot where we perform analyses and fitting in Section 4.
	$\left(L_X\times L_Y\times L_Z\right)H^3$	$N_X\times N_Y\times N_Z$				$\left({\Omega }_0^{-1}\right)$
I	$0.1\times 0.1\times 0.2$	$512\times 512\times 1024$	$134,217,728$	$2.0$	$0.1$	$36.0$	$284$
II	$0.2\times 0.2\times 0.2$	$512\times 512\times 512$	$153,600,000$	$0.3$	$0.02$	$110.0$	$174$

Table 2. Table 2: Interpretation Guidelines for BIC

${\Delta }_{\mathrm{BIC}}$	Evidence against the preferred model
$0-2$	Not worth more than a bare mention
$2-6$	Positive
$6-10$	Strong
$>10$	Very Strong

Table 3. Table 3: Interpretation Guidelines for AIC

${\Delta }_{\mathrm{AIC}}$	Level of empirical support for a model
$0-2$	Substantial
$2-4$	Strong
$4-7$	Considerably less
$>10$	Essentially none

Table 4. Table 4: Model Fitting Results for Run I ( t sg = 7.5 / Ω subscript 𝑡 sg 7.5 Ω t_{\rm sg}=7.5/\Omega )

Models	Best-fit Parameters	Mass Scales [$M_G$]	$-\ln ℒ$	${\Delta }_{\mathrm{BMS}}$	${\Delta }_{\mathrm{BIC}}$	${\Delta }_{\mathrm{AIC}}$
Simply Tapered Power Law K=2, $𝜽=(\alpha ,x_{\exp })$	$\alpha ={0.208}_{-0.014}^{+0.011}$ $x_{\exp }={8.905}_{-0.464}^{+0.323}$	$M_{\exp }={0.1385}_{-0.0515}^{+0.0529}$	660.443	63.7	53.1	60.4
Variably Tapered Power Law K=3, $𝜽=(\alpha ,\beta ,x_{\exp })$	$\alpha ={0.036}_{-0.041}^{+0.041}$ $\beta ={0.298}_{-0.040}^{+0.061}$ $x_{\exp }={4.734}_{-1.128}^{+1.022}$	$M_{\exp }={0.0021}_{-0.0014}^{+0.0038}$	633.695	10.6	5.3	8.9
Broken Cumulative Power Law K=3, $𝜽=({\alpha }_1,{\alpha }_2,x_{\mathrm{br}})$	${\alpha }_1={0.162}_{-0.019}^{+0.009}$ ${\alpha }_2={0.589}_{-0.071}^{+0.042}$ $x_{\mathrm{br}}={4.550}_{-0.607}^{+0.003}$	$M_{\mathrm{br}}={0.0018}_{-0.0008}^{+0.0000}$	631.683	10.1	1.2	4.9
Truncated Power Law K=2, $𝜽=(\alpha ,x_{\mathrm{tr}})$	$\alpha ={0.140}_{-0.029}^{+0.012}$ $x_{\mathrm{tr}}={10.880}_{-0.000}^{+0.000}$	$M_{\mathrm{tr}}={0.9981}_{-0.0000}^{+0.0000}$	651.846	47.5	35.9	43.2
Broken Power Law K=3, $𝜽=({\alpha }_1,{\alpha }_2,x_{\mathrm{br}})$	${\alpha }_{\mathrm{𝟏}}=-{0.079}_{-0.049}^{+0.043}$ ${\alpha }_2={0.628}_{-0.055}^{+0.078}$ $x_{\mathrm{br}}={5.620}_{-0.285}^{+0.329}$	$M_{\mathrm{br}}={0.0052}_{-0.0013}^{+0.0020}$	631.946	7.2	1.8	5.4
Truncated Broken Power Law K=4, $𝜽=({\alpha }_1,{\alpha }_2,x_{\mathrm{br}},x_{\mathrm{tr}})$	${\alpha }_{\mathrm{𝟏}}=-{0.102}_{-0.043}^{+0.058}$ ${\alpha }_2={0.468}_{-0.070}^{+0.082}$ $x_{\mathrm{br}}={5.066}_{-0.144}^{+0.519}$ $x_{\mathrm{tr}}={10.880}_{-0.000}^{+0.000}$	$M_{\mathrm{br}}={0.0030}_{-0.0004}^{+0.0020}$ $M_{\mathrm{tr}}={0.9981}_{-0.0000}^{+0.0000}$	628.241	0.0	0.0	0.0
Three-segment Power Law K=5, $𝜽=({\alpha }_1,{\alpha }_2,{\alpha }_3,x_{\mathrm{br1}},x_{\mathrm{br2}})$	${\alpha }_{\mathrm{𝟏}}=-{0.102}_{-0.043}^{+0.057}$ ${\alpha }_2={0.468}_{-0.071}^{+0.081}$ ${\alpha }_3=5.78{\mathrm{e5}}_{-3.9\mathrm{e4}}^{+6.07\mathrm{e7}}\text{<span><sup>*</sup><span><span><sup>*</sup>The best-fit third slope, <xmath><xmapp><xmtok><xmtok>α</xmtok><xmtok>3</xmtok></xmtok></xmapp></xmath>, of the Three-segment Power Law model is extremely large because this model essentially reverts to the Truncated Broken Power Law model (see also Section <a><span>3.2</span></a> and Fig. <a><span>3</span></a>).</span></span></span>}$ $x_{\mathrm{br1}}={5.066}_{-0.142}^{+0.511}$ $x_{\mathrm{br2}}={10.880}_{-0.632}^{+0.000}$	$M_{\mathrm{br1}}={0.0030}_{-0.0004}^{+0.0020}$ $M_{\mathrm{br2}}={0.9981}_{-0.4674}^{+0.0000}$	628.241	0.0	5.6	2.0

Table 5. Table 5: Model Fitting Results for Run II ( t sg = 7.6 / Ω subscript 𝑡 sg 7.6 Ω t_{\rm sg}=7.6/\Omega )

Models	Best-fit Parameters	Mass Scales [$M_G$]	$-\ln ℒ$	${\Delta }_{\mathrm{BMS}}$	${\Delta }_{\mathrm{BIC}}$	${\Delta }_{\mathrm{AIC}}$
Simply Tapered Power Law K=2, $𝜽=(\alpha ,x_{\exp })$	$\alpha ={0.388}_{-0.039}^{+0.030}$ $x_{\exp }={6.397}_{-0.721}^{+0.473}$	$M_{\exp }={11.2531}_{-5.7787}^{+6.8113}$	311.520	16.1	10.2	13.3
Variably Tapered Power Law K=3, $𝜽=(\alpha ,\beta ,x_{\exp })$	$\alpha ={0.304}_{-0.060}^{+0.063}$ $\beta ={0.527}_{-0.088}^{+0.233}$ $x_{\exp }={5.178}_{-0.777}^{+0.743}$	$M_{\exp }={3.3255}_{-1.7965}^{+3.6627}$	309.274	11.6	10.8	10.8
Broken Cumulative Power Law K=3, $𝜽=({\alpha }_1,{\alpha }_2,x_{\mathrm{br}})$	${\alpha }_1={0.348}_{-0.033}^{+0.035}$ ${\alpha }_2={0.866}_{-0.080}^{+0.143}$ $x_{\mathrm{br}}={3.226}_{-0.065}^{+0.035}$	$M_{\mathrm{br}}={0.4719}_{-0.0298}^{+0.0169}$	303.864	1.4	0.0	0.0
Truncated Power Law K=2, $𝜽=(\alpha ,x_{\mathrm{tr}})$	$\alpha ={0.360}_{-0.051}^{+0.029}$ $x_{\mathrm{tr}}={7.891}_{-0.000}^{+0.000}$	$M_{\mathrm{tr}}={50.126}_{-0.000}^{+0.000}$	310.769	14.9	8.7	11.8
Broken Power Law K=3, $𝜽=({\alpha }_1,{\alpha }_2,x_{\mathrm{br}})$	${\alpha }_1={0.240}_{-0.068}^{+0.064}$ ${\alpha }_2={0.895}_{-0.103}^{+0.265}$ $x_{\mathrm{br}}={4.431}_{-0.180}^{+0.322}$	$M_{\mathrm{br}}={1.5754}_{-0.2594}^{+0.5992}$	309.058	11.0	10.4	10.4
Truncated Broken Power Law K=4, $𝜽=({\alpha }_1,{\alpha }_2,x_{\mathrm{br}},x_{\mathrm{tr}})$	${\alpha }_1={0.249}_{-0.077}^{+0.065}$ ${\alpha }_2={0.729}_{-0.200}^{+0.152}$ $x_{\mathrm{br}}={4.330}_{-0.691}^{+0.261}$ $x_{\mathrm{tr}}={7.891}_{-0.000}^{+0.000}$	$M_{\mathrm{br}}={1.4241}_{-0.7103}^{+0.4254}$ $M_{\mathrm{tr}}={50.1262}_{-0.0000}^{+0.0000}$	307.718	8.4	12.9	9.7
Three-segment Power Law K=5, $𝜽=({\alpha }_1,{\alpha }_2,{\alpha }_3,x_{\mathrm{br1}},x_{\mathrm{br2}})$	${\alpha }_1={0.771}_{-0.159}^{+0.291}$ ${\alpha }_2=-{0.362}_{-0.303}^{+0.181}$ ${\alpha }_3={0.833}_{-0.075}^{+0.164}$ $x_{\mathrm{br1}}={1.827}_{-0.314}^{+0.270}$ $x_{\mathrm{br2}}={3.501}_{-0.079}^{+0.198}$	$M_{\mathrm{br1}}={0.1166}_{-0.0314}^{+0.0362}$ $M_{\mathrm{br2}}={0.6216}_{-0.0474}^{+0.1359}$	303.711	0.0	10.0	3.7

Table 6. Table 6: Mass Distribution Models

Name	Mass Distribution Function	Mass Scale	PDF in Likelihood Estimator
K: # of parameters	CDF[$P_{>}\left(M\right)$] or PDF[$\xi \left(M\right)$]	$x\equiv \ln \left(M/M_{\min }\right)$	$p(x;𝜽)$
Simply Tapered Power Law (Abod et al., 2019) K=2, $𝜽=(\alpha ,x_{\exp })$	$P_{>}\left(M\right)={\left(\frac{M}{M_{\min }}\right)}^{-\alpha }\exp \left[-\frac{M-M_{\min }}{M_{\exp }}\right]$	$x_{\exp }\equiv \ln \left(\frac{M_{\exp }}{M_{\min }}\right)$	$\frac{\left(\alpha +e^{-x_{\exp }}{e}^x\right)}{\exp \left[\alpha x+e^{-x_{\exp }}\left(e^x-1\right)\right]}$
Variably Tapered Power Law (Schäfer et al., 2017) K=3, $𝜽=(\alpha ,\beta ,x_{\exp })$	$P_{>}\left(M\right)={\left(\frac{M}{M_{\min }}\right)}^{-\alpha }\exp \left[-\frac{M^{\beta }-M_{\min }^{\beta }}{M_{\exp }^{\beta }}\right]$	$x_{\exp }\equiv \ln \left(\frac{M_{\exp }}{M_{\min }}\right)$	$\frac{\left(\alpha +\beta e^{\beta \left(x-x_{\exp }\right)}\right)}{\exp \left[\alpha x+e^{-\beta x_{\exp }}\left(e^{\beta x}-1\right)\right]}$
Broken Cumulative Power Law K=3, $𝜽=({\alpha }_1,{\alpha }_2,x_{\mathrm{br}})$		$x_{\mathrm{br}}\equiv \ln \left(\frac{M_{\mathrm{br}}}{M_{\min }}\right)$
Truncated Power Law (Schäfer et al., 2017) K=2, $𝜽=(\alpha ,x_{\mathrm{tr}})$		$x_{\mathrm{tr}}\equiv \ln \left(\frac{M_{\mathrm{tr}}}{M_{\min }}\right)$
Broken Power Law K=3, $𝜽=({\alpha }_1,{\alpha }_2,x_{\mathrm{br}})$		$x_{\mathrm{br}}\equiv \ln \left(\frac{M_{\mathrm{br}}}{M_{\min }}\right)$
Truncated Broken Power Law K=4, $𝜽=({\alpha }_1,{\alpha }_2,x_{\mathrm{br}},x_{\mathrm{tr}})$
Three-segment Power Law K=5, $𝜽=({\alpha }_1,{\alpha }_2,{\alpha }_3,x_{\mathrm{br1}},x_{\mathrm{br2}})$

Equations

$$\frac{\partial {\rho }_{\mathrm{g}}}{\partial t}+\nabla \cdot ({\rho }_{\mathrm{g}}\mathbf{u})$$

$$\displaystyle\begin{split}\frac{\partial{(\rho_{\rm g}\bm{u})}}{\partial{t}}+\nabla\cdot(\rho_{\rm g}\bm{u}\bm{u}+P\bm{I})&=\\ \rho_{\rm g}\biggl{[}2\bm{u}\times\bm{\Omega}_{0}+3{\Omega}_{0}^{2}\bm{x}&-{\Omega}_{0}^{2}\bm{z}\biggr{]}+\rho_{\rm p}\frac{\bar{\bm{v}}-\bm{u}}{t_{\mathrm{stop}}},\end{split}$$

$$\begin{array}{rl}\frac{d{\mathbf{v}}_i}{dt}=2{\mathbf{v}}_i\times {\mathbf{\Omega }}_0+3{\Omega }_0^2{\mathbf{x}}_i& -{\Omega }_0^2{\mathbf{z}}_i\\ -\frac{{\mathbf{v}}_i-\mathbf{u}}{t_{\mathrm{stop},\mathrm{i}}}& -\nabla {\Phi }_{\mathrm{sg}}-2\eta v_{\mathrm{K}}{\Omega }_0\hat{x},\end{array}$$

$${\nabla }^2{\Phi }_{\mathrm{sg}}=4\pi G{\rho }_{\mathrm{p}}$$

$${\text{uptau}}_{\mathrm{s}}={\Omega }_0{t}_{\mathrm{stop}},$$

$$Z=\frac{{\Sigma }_{\mathrm{p}}}{{\Sigma }_{\mathrm{g}}},$$

$$\Pi \equiv \frac{\eta v_{\mathrm{K}}}{c_{\mathrm{s}}}\equiv -\frac{1}{2}\frac{c_{\mathrm{s}}}{v_{\mathrm{K}}}\frac{\partial \ln P}{\partial \ln R},$$

$$\tilde{G}\equiv \frac{4\pi G{\rho }_0}{{\Omega }_0^2}=\frac{4}{\sqrt{2\pi }Q},$$

$${\rho }_{\mathrm{g}}$$

$${\rho }_{\mathrm{p}}$$

$$t_{\mathrm{sg}}=t-t_0,$$

$$M_{\mathrm{G}}=\pi {\left(\frac{{\lambda }_{\mathrm{G}}}{2}\right)}^2{\Sigma }_{\mathrm{p}}=4{\pi }^5\frac{G^2{\Sigma }_{\mathrm{p}}^3}{{\Omega }_0^4}=\frac{\sqrt{2}}{2}{\pi }^{\frac{9}{2}}{Z}^3{\tilde{G}}^2({\rho }_0{H}^3),$$

$$p(x;\mathbf{\theta })=\frac{1}{N_{\mathrm{tot}}}\frac{dN(x;\mathbf{\theta })}{dx}=C(\mathbf{\theta })g(x;\mathbf{\theta }),$$

$$P_{>}(x;\mathbf{\theta })=\underset{x}{\overset{+\infty }{\int }}p(x^{\prime };\mathbf{\theta })d{x}^{\prime }=C(\mathbf{\theta })\underset{x}{\overset{+\infty }{\int }}g(x^{\prime };\mathbf{\theta })d{x}^{\prime },$$

$$C(\mathbf{\theta })=\frac{1}{\underset{0}{\overset{+\infty }{\int }}g(x;\mathbf{\theta })dx},$$

$$\mathcal{L}(\mathbf{x}|\mathbf{\theta })\equiv \prod _{i=1}^{N_{\mathrm{tot}}}p(x_i;\mathbf{\theta })=\prod _{i=1}^{N_{\mathrm{tot}}}C(\mathbf{\theta })g(x_i;\mathbf{\theta }).$$

$$\ln \mathcal{L}(\mathbf{x}|\mathbf{\theta })=\sum _{i=1}^{N_{\mathrm{tot}}}\ln p(x_i;\mathbf{\theta })=N_{\mathrm{tot}}\ln C(\mathbf{\theta })+\sum _{i=1}^{N_{\mathrm{tot}}}\ln g(x_i;\mathbf{\theta }).$$

$$\{\begin{array}{rl}\frac{\partial \ln \mathcal{L}(\mathbf{x}|\mathbf{\theta })}{\partial {\theta }_j}& =0,\\ \frac{{\partial }^2\ln \mathcal{L}(\mathbf{x}|\mathbf{\theta })}{\partial {{\theta }_j}^2}& <0,\end{array}$$

$$\Delta {\mathbf{\theta }}_{\mathrm{bs}}^{+}$$

$$\Delta {\mathbf{\theta }}_{\mathrm{bs}}^{-}$$

$$P_{>}(M;\mathbf{\theta })=c_1{M}^{-\alpha }\exp \left(-\frac{M}{M_{\exp }}\right),$$

$$P_{>}(M;\mathbf{\theta })=c_2{M}^{-\alpha }\exp \left[-{\left(\frac{M}{M_{\exp }}\right)}^{\beta }\right],$$

$$P_{>}(M;\mathbf{\theta })=\{\begin{array}{rlr}& c_{31}{M}^{-{\alpha }_1}& M⩽M_{\mathrm{br}}\\ & c_{32}{M}^{-{\alpha }_2}& M>M_{\mathrm{br}}\end{array},$$

$$\xi (M;\mathbf{\theta })=\{\begin{array}{rlr}& c_4{M}^{-\alpha -1}& M⩽M_{\mathrm{tr}}\\ & 0& M>M_{\mathrm{tr}}\end{array},$$

$$\xi (M;\mathbf{\theta })=\{\begin{array}{rlr}& c_{51}{M}^{-{\alpha }_1-1}& M⩽M_{\mathrm{br}}\\ & c_{52}{M}^{-{\alpha }_2-1}& M>M_{\mathrm{br}}\end{array}.$$

$$\xi (M;\mathbf{\theta })=\{\begin{array}{rlr}& c_{61}{M}^{-{\alpha }_1-1}& M⩽M_{\mathrm{br}}\\ & c_{62}{M}^{-{\alpha }_2-1}& \text{otherwise}\\ & 0& M>M_{\mathrm{tr}}\end{array}.$$

$$\xi (M;\mathbf{\theta })=\{\begin{array}{rlr}& c_{71}{M}^{-{\alpha }_1-1}& M⩽M_{br1}\\ & c_{72}{M}^{-{\alpha }_2-1}& \text{otherwise}\\ & c_{73}{M}^{-{\alpha }_3-1}& M>M_{br2}\end{array}.$$

$$\text{BIC}=K\ln (N_{\mathrm{tot}})-2\ln \mathcal{L},$$

$$\text{AIC}=2K-2\ln \mathcal{L}.$$

$${\Delta }_{\mathrm{BIC}}$$

$${\Delta }_{\mathrm{AIC}}$$

$$\text{BMS}=-2\times \text{median}\left(\ln \mathcal{L}({\mathbf{x}}_k|{\mathbf{\theta }}_{\mathrm{MLE}})\right).$$

$${\Delta }_{\mathrm{BMS}}=\text{BMS}-{\text{BMS}}_{\min },$$

Full text

protoplanetary disks—hydrodynamics—instabilities—planets and satellites: formation—methods: data analysis—methods: statistical

††journal: ApJ††software: ATHENA (Stone et al., 2008; Bai & Stone, 2010c), Matplotlib (Hunter, 2007), Numpy & Scipy (Jones et al., 2001), emcee v3.0rc2 (Foreman-Mackey et al., 2013), PLAN (Li, 2019a), Pyridoxine (Li, 2019b).

1 Introduction

An indispensable step in planet formation is to build planetesimals—super-kilometer objects bound by self-gravity—in protoplanetary disks (Chiang & Youdin, 2010; Johansen et al., 2014). One of the compelling pathways to planetesimal formation is the efficient concentration of solids by the streaming instability (SI) followed by their gravitational collapse (Youdin & Goodman, 2005; Johansen et al., 2007).

The SI is an aerodynamic instability arising from the relative drift and the mutual drag forces between gas and solids (Youdin & Goodman, 2005). The SI is one example of a broader class of drag instabilities in protoplanetary disks (Goodman & Pindor, 2000; Lin & Youdin, 2017; Squire & Hopkins, 2018).

Strong particle clumping can be induced by the SI under the right conditions, that is, when the midplane dust-to-gas volume density ratio exceeds unity, which further depends on the particle size, dust-to-gas surface density ratio (often referred as metallicity), and gas pressure gradient (Johansen & Youdin, 2007; Bai & Stone, 2010a, b; Carrera et al., 2015; Yang et al., 2017). In a typical smooth disk with cm-sized pebbles, slightly super solar metallicity is required to trigger the strong SI (Johansen et al., 2009). Higher metallicity may be reached by gas removal (e.g., via photoevoporation) or dust pile-up (e.g., at pressure bumps, snow lines), where smaller solids can also lead to strong particle clumping (Carrera et al., 2017; Dr\każkowska et al., 2016; Dr\każkowska & Alibert, 2017; Schoonenberg & Ormel, 2017). Previous studies have also shown that high-resolution SI simulations with self-gravity can produce a broad initial top-heavy mass distribution of planetesimals (Johansen et al., 2015; Simon et al., 2016, 2017).

Recently, Yang et al. (2018) suggest that the midplane dust-to-gas density ratio exceeding order unity may not be necessary for the strong SI when the radial diffusion of particles driven by turbulence near the midplane is weak. Lin (2019) find that particle back-reaction on to the gas can lead to self-sustained dust sedimentation against the vertical shear instability turbulence at a modest super solar metallicity, which may trigger the SI. Furthermore, Krapp et al. (2019) investigate the linear growth of the SI with multiple dust species for the first time(Bai & Stone, 2010a; Schaffer et al., 2018).

Observations of asteroids and Trans-Neptunian objects support models in which planetesimals were born big (Morbidelli et al., 2009), with evidence of a drop-off in planetesimal numbers below $\sim$1–50 kilometers, depending on the population (Delbo’ et al., 2017; Singer et al., 2019). Nesvorný et al. (2019) find that SI simulations correctly predict the primarily prograde mutual inclinations of the abundant binaries in the Cold Classical Kuiper Belt (Grundy et al., 2019). Further studies on the demographics of planetesimals formed via the SI, and by other mechanisms, offer the promise of more detailed observational comparisons and tests.

Quantifying the mass distribution of planetesimals formed by the SI is of significant interest. Due to the high computational cost of SI simulations, a parameterized mass function can be used as the input for global studies of disk evolution and planet formation (Dr\każkowska & Dullemond, 2014). Furthermore, the shape of the mass distribution offers insights to the physical processes of particle clumping and gravitational collapse.

Previous work has fit the mass distribution to a simple power law (Simon et al., 2016, 2017) or to a power law with an exponential cutoff or truncation (Schäfer et al., 2017; Abod et al., 2019). These work suggested that the initial planetesimal mass function might be near-universal. However, it is not trivial to determine the best parameterization of the initial planetesimal mass function in SI simulations. Moreover, it is not clear if a single functional form can describe planetesimal formation with different physical conditions, i.e. simulation parameters.

Motivated by these issues, our goal in this paper is to better understand and constrain the broad initial mass distribution of planetesimals with robust statistical analyses. This work will fit many different parameterizations to simulated planetesimal mass distributions. To determine which models best describe the data, model selection techniques weigh the goodness of fit against a complexity penalty, intended to avoid the overfitting of data features that might be spurious. Since there is no universal agreement on complexity penalties either, we apply different model selection techniques, including a bootstrap method that we developed independently.

The paper is organized as follows. In Section 2, we begin with an overview of the numerical models and our simulations. Section 2.3 then introduce our newly-developed clump-finding tool, PLAN. Section 3 lays out all the statistical models and our fitting procedure as well as the model selection criteria. In Section 4, we show the fitting results and the model selection results. Section 5 discusses the implications of our statistical understanding, with a summary and conclusions in the end.

2 Method

To simulate the formation of planetesimals, we use the ATHENA code with a similar setup to Simon et al. (2017). In Section 2.1, we briefly introduce the numerical methods employed in ATHENA for modeling the coupled dynamics of gas and particles, including the self-gravity of solids, in a protoplanetary disk (see Bai & Stone, 2010c; Simon et al., 2016, for more details). Section 2.2 then summarizes the numerical setup and parameters used in our simulations. Section 2.3 explains how PLAN identifies and characterizes all the self-bound clumps in the output particle data.

2.1 Planetesimal Formation Simulations

We use ATHENA to simulate a small three-dimensional vertically-stratified patch of the protoplanetary disk with the local shearing box approximation (Stone et al., 2008; Stone & Gardiner, 2010; Hawley et al., 1995). This approximation—which is justified by the small length scales of the SI compared to the radial position in the disk—maps the global disk geometry $(R,\phi,z^{\prime})$ onto a local Cartesian coordinate system $(x,y,z)$ (Goldreich & Lynden-Bell, 1965). The local box is centered at a fiducial disk radius ($R_{0}$) in the midplane, where $(x,y,z)\equiv(R-R_{0},R_{0}\phi,z^{\prime})$, where the Keplerian frequency and velocity are $\Omega_{0}$ and $v_{\rm K}=\Omega_{0}R_{0}$, respectively.

In this non-inertial computational domain, ATHENA solves the equations of gas dynamics and the equation of motion for each particle (indexed by $i$)

[TABLE]

where $\rho_{\rm g}$, $\bm{u}$ and $P$ are density, velocity and pressure of gas, $\bm{I}$ is the identity matrix, $\bm{\Omega}_{0}=\Omega_{0}\hat{z}$, $\rho_{\rm p}$ and $\bar{\bm{v}}$ are the average density and velocity of the particles in a hydrodynamic grid cell, $t_{\rm stop}$ is the dimensional stopping time, $\bm{v}_{i}$ is the velocity of the $i$-th particle, $\Phi_{\rm sg}$ is the potential field of the self-gravity of solids, and $\eta$ denotes the relative difference between the gas orbital velocity and the Keplerian velocity due to the radial pressure gradient in the disk.

Our model calculates the Coriolis forces, radial and vertical tidal gravity, and the particle feedback exerted on the gas, as in the right hand side of Eq. 2. The equation of state for the gas is assumed to be isothermal, $P=c_{\rm s}^{2}\rho_{\rm g}$, where the constant $c_{\rm s}$ is the isothermal sound speed. We neglect the self-gravity of the gas because the gas density fluctuations are relatively negligible.

For solids, ATHENA adopts the super-particle treatment, where each particle in our simulations statistically represents a large number of pebbles in terms of mass. The acceleration of each particle is governed by Eq. 3 with the Coriolis and tidal forces (similar to those in Eq. 2), and also the gas drag as well as the force due to particle self-gravity. The gravitational potential field, $\Phi_{\rm sg}$, is obtained by solving Poisson’s equation

[TABLE]

with the Fast Fourier Transform (FFT) method of Simon et al. (2016), where $G$ is the gravitational constant. The accuracy of particle self-gravity from such a method depends on the grid resolution. The last source term in Eq. 3, $-2\eta v_{\rm K}\Omega_{0}\hat{x}$, is a constant radial force in ATHENA to implement the effective global radial pressure gradient under the restrictions of the local model.

The radial and azimuthal boundary conditions (BCs) for our model are the standard shearing-periodic BCs (Stone & Gardiner, 2010). In the vertical direction, we use a modified outflow BCs that extrapolates the gas density into the ghost zones exponentially and prohibits any gas inflow (Simon et al., 2011; Li et al., 2018). These vertical BCs maintain hydrostatic equilibrium and reduce artificial gas motions at vertical boundaries, which is beneficial for simulations in short boxes. Furthermore, the total gas mass is renormalized to compensate the gas outflow at each time step to ensure the mass conservation.

The physical behavior of our simulations are dominated by four key dimensionless parameters. The SI is characterized by the first three of them: the dimensionless particle stopping time

[TABLE]

which represents the ratio of a particle’s aerodynamic ($t_{\rm stop}$) and orbital ($\Omega_{0}^{-1}$) timescales, increases with a particle’s size, and decreases with the local gas density; the surface density ratio between the solids ($\Sigma_{\rm p}$) and the gas ($\Sigma_{\rm g}$)

[TABLE]

which is sometimes called the total solid-to-gas mass ratio; the global radial pressure gradient parameter

[TABLE]

which accounts for the strength of the headwind on the particles. The fourth key parameter controls the relative strength of the particle self-gravity compared with the tidal shear

[TABLE]

where $\rho_{0}$ is the midplane gas density, $Q$ is the Toomre’s Q (Toomre, 1964).

The code units of our simulations are set to the natural units of the shearing box. The density unit and the time unit are $\rho_{0}$ and $\Omega_{0}^{-1}$, respectively. While $\rho_{0}$ and $\Omega_{0}$ are code units, their allowed physical values and thus the choice of disk model are constrained by the $\tilde{G}$ parameter. The length unit is $H=c_{\rm s}/\Omega_{0}$, the vertical scale height of the gas.

2.2 Numerical Setup

All of our simulations are initiated with Gaussian vertical density profiles for both the gas and particles

[TABLE]

where $H_{\rm p}$ is the particle scale height and is set to $0.02H$ in the beginning. The choice of such an initial particle scale height matches our previous work in Simon et al. (2017), which let particles naturally sediment to a pseudo-equilibrium scale height of the order of $0.01H$ (Yang & Johansen, 2014; Li et al., 2018). Both gas and particles are then initialized with the Nakagawa–Sekiya–Hayashi (NSH) equilibrium drift velocities (Nakagawa et al., 1986).

In this work, we fix $\Pi=0.05$ and $\tilde{G}=0.05$ ($Q\simeq 32$), which are typical values in protoplanetary disks. Table 1 lists other physical and numerical parameters for all two of our simulations. Run I has $(\uptau_{\rm s},Z)=(2.0,0.1)$ and a higher resolution in a smaller domain (one grid cell $\Delta x$ is $H/5120$ wide, the highest resolution to date). Run II has $(\uptau_{\rm s},Z)=(0.3,0.02)$ and a lower resolution in a larger domain. The data of Run II are directly taken from Simon et al. (2017), and our Run I is a higher resolution version of another simulation in Simon et al. (2017). The relation between particle size and $\uptau_{\rm s}$ depends on uncertain properties of particles and gas disk. The range of $\uptau_{\rm s}$ adopted here may correspond to solids with any size from millimeter to decimeter. These physical parameters are chosen known to produce strong particle clumping that triggers gravitational collapse. The resulting bound clumps are the subject of our statistical analyses below.

Following previous studies, we start our simulations first without particle self-gravity. Only after the SI has fully developed and saturated, do we switch on the self-gravity. This approach significantly reduces the computational expense and has little influence on the final properties of planetesimals (Simon et al., 2016; Abod et al., 2019). For convenience, we define a self-gravity time

[TABLE]

where $t$ is the simulation time and $t_{0}$ denotes the time when the self-gravity is turned on (see the next to last column in Table 1). Also, we present all planetesimal masses in units of the dimensional mass for a self-gravitating particle disk

[TABLE]

where $\lambda_{\rm G}$ is the critical unstable wavelength from the standard Toomre dispersion relation.

To translate $M_{\rm G}$ into a physical mass unit and planetesimal size, additional assumptions about the disk model are required. For instance, if we assume a fiducial disk radius of $R_{0}=10$ au in the modified minimum mass solar nebula (MMSN) model of Chiang & Youdin (2010, adapted from , hereafter CY10 model), where $\Sigma_{\rm g}\propto R^{-3/2}$ and the temperature $T\propto R^{-3/7}$, then our $\tilde{G}=0.05$ parameter implies the gas mass in the CY10 model is about half the original MMSN values. For these parameters, the CY10 model has $\Pi=0.068$, slightly higher than $\Pi=0.05$ in our simulations. Nevertheless, a smaller $\Pi$ value might arise from a weak pressure bump, a common substructure in protoplanetary disks (Pinilla & Youdin, 2017; Dullemond et al., 2018). Moreover, $\Pi$ values do not significantly affect the planetesimal mass distribution, as studied in Abod et al. (2019). The mass unit for Run I then equates to $M_{\rm G}=1.82\times 10^{23}$ g $=0.19\ M_{\rm Ceres}$. With $1$ g cm*-3* as the mean density, the physical radius of a $1M_{\rm G}$ is $\simeq 350$ km in this model. For Run II, the same gas disk model and location gives $M_{\rm G}=1.45\times 10^{21}$ g $=0.0015\ M_{\rm Ceres}$, and a physical radius of $\simeq 70$ km.

Fig. 1 shows a snapshot of the particle surface density from Run I at $t_{\rm sg}=4/\Omega_{0}$, where solids already collapse into self-bound clumps. In the following section, we describe how PLAN finds these clumps in detail.

2.3 Clump-Finding with PLanetesimal ANalyzer

To identify and further characterize the properties of planetesimals produced in our simulations, we develop a new clump-finding tool, PLanetesimal ANalyzer (PLAN, Li (2019a)). It is designed to work with the 3D particle output of ATHENA and find self-bound clumps robustly and efficiently. PLAN is scalable to analyze billions of particles and many snapshots simultaneously because of its massively parallelized scheme written in C**+****+** with OpenMP/MPI.

We now briefly present the workflow of PLAN. The approach is based on the dark matter halo finder HOP developed by Eisenstein & Hut (1998), which is able to quickly group physically related particles. PLAN first builds a memory-efficient linear Barnes-Hut tree representing all the particles in the Morton order (Barnes & Hut, 1986). Each particle is then assigned a density computed from the nearest $N_{\rm den}$ particles ($N_{\rm den}=64$ by default). For particles with densities higher than a threshold, $\delta_{\rm outer}=8\rho_{0}/\tilde{G}$, PLAN chains them up towards their densest neighbors recursively until a density peak is reached. The value of $\delta_{\rm outer}$ is physically motivated to be slightly smaller than the Roche density ($9\rho_{0}/\tilde{G}$) such that PLAN can quickly find most relevant particles. All the particle chains leading to the same density peak are combined into a group.

PLAN then merges those groups by examining their boundaries to construct a list of bound clumps. Based on the total kinematic and gravitational energies, deeply intersected groups are merged if bound. However, two particle groups with a saddle point less dense than $\delta_{\rm saddle}=2.5\delta_{\rm outer}$ remain separated (Eisenstein & Hut, 1998). Next, PLAN goes through each group—or raw clump—to unbind any contamination (i.e., passing-by and not bound) particles and gather possibly unidentified member particles within its Hill sphere. After discarding those clumps (if any) with Hill radii ($R_{\rm Hill}$) smaller than one hydrodynamic grid cell ($\Delta x$) or density peaks less than $\delta_{\rm peak}=3\delta_{\rm outer}$, PLAN outputs the final list of clumps with their physical properties derived from particles.

Most clumps in our high-resolution simulations are highly-concentrated, where particles often collapse into regions much smaller than the Hill radius (see Fig. 1). For small clumps, these regions are comparable to one cell size. The particle-mesh method of calculating self-gravity does not resolve scales below $\Delta x$, which is a primary motivation for our high-resolution simulation and the reason why PLAN compares $R_{\rm Hill}$ to $\Delta x$. While this work was in progress, PLAN was already used in the analyses of Abod et al. (2019), and Nesvorný et al. (2019).

Our former clump-finding tool analyzed the surface density of solids in the ($x,y$) cells of the hydrodynamic grid. This technique identifies prominent clumps and calculates the clump masses as the encapsulated column mass within their projected Hill radii. Such a treatment is limited by the grid resolution and has difficulty detecting small planetesimals, especially when a massive clump is nearby. Consequently, this method tends to overestimate the clump mass by including the mass of surrounding small planetesimals as well as other solids that are vertically far away. PLAN overcomes those issues by diagnosing the particle data instead, as described above.

Fig. 2 shows that the PLAN results agree with our previous results at large masses, though the previous analyses slightly overestimated masses, as explained above. The significant advantage of the PLAN analysis is that we identify gravitationally bound clumps at much lower masses than before, which improves our ability to statistically characterize the resulting mass distributions.

3 Statistical Modeling of the Mass Distribution

This section details our statistical methodology for analyzing the planetesimal mass distribution. We present a maximum likelihood estimator (MLE) for estimating the parameters of a given model, and the uncertainty on those parameters in Section 3.1. Section 3.2 lists the models that we fit, which vary in complexity from 2 to 5 parameters. Finally, Section 3.3 describes our model selection criteria, and how they apply a penalty to more complex models.

3.1 Maximum Likelihood Parameter Estimation

We assume that the masses, $M$, of planetesimals (strictly, protoplanetesimal clumps in the simulations) are drawn from a probability density function (PDF), $\xi(M)$, parameterized by a vector $\bm{\theta}$, where $\xi(M;\bm{\theta})dM$ represents the probability that a given clump forms in the mass interval $M$ to $M+dM$.

In practice, it is easier to work with a logarithmic mass coordinate, $x=\ln{(M/M_{\rm min})}$, referenced to the minimum mass of the distribution, $M_{\rm min}$. With this transformation, the functional form of the PDF is different and written as

[TABLE]

where the first equality simply relates the PDF to the total number of bodies $N_{\rm tot}$ and the number in a given logarithmic mass interval, $dN$. In the second equality, the normalization factor $C(\bm{\theta})$ is introduced for later convenience (see also Youdin 2011).

Accordingly, the cumulative distribution function (CDF) is

[TABLE]

which denotes the expected fraction of clumps with masses larger than $M(=e^{x}M_{\rm min})$ in the distribution111Note the minus sign when relating the PDF to the CDF: $\displaystyle{p(x;\bm{\theta})=-\frac{d{}}{d{x}}P_{>}(x;\bm{\theta})}$.. The normalization of the CDF, $P_{>}(x;\bm{\theta})|_{x=0}=1$, gives

[TABLE]

which requires that $g(x;\bm{\theta})$ does not diverge as $x\to+\infty$.

To fit a model to the data, i.e. the simulated mass distribution of planetesimals, we consider the likelihood of the data given the model

[TABLE]

It is usually easier to consider the log-likelihood

[TABLE]

The maximum likelihood estimator (MLE) of $\bm{\theta}$ estimates the best-fit parameters, $\bm{\theta}_{\rm MLE}$, by maximizing the log-likelihood (within certain physical bounds if necessary). For some simple log-likelihood functions, $\bm{\theta}_{\rm MLE}$ can be solved as the root(s) of

[TABLE]

where $\theta_{j}$ means the $j$-th parameter in $\bm{\theta}$. However, constraints on the allowed values of parameters sometimes confound traditional root-finding methods.

In this work, we apply numerical techniques to maximize the likelihood of the trial PDF (see Section 3.2 for our choices). In practice, we first use the python package emcee to explore the parameter space with a Markov-chain Monte Carlo (MCMC) approach (Foreman-Mackey et al., 2013) to obtain an initial guess of parameters, $\bm{\theta}_{\rm MCMC}$. We then use the minimize method222A complicated PDF may lead to a non-convex or non-smooth log-likelihood function, which is known to be difficult to minimize. We always test different algorithms (e.g., “Powell”, “Newton-CG”, “L-BFGS-B”, etc.) provided by minimize and run a set of optimizations with initial guesses selected in a mesh grid centered on $\bm{\theta}_{\rm MCMC}$. We then take the solution leading to the lowest $-\ln\mathcal{L}$. in the scipy.optimize package (Jones et al., 2001) to find the most likely $\bm{\theta}$ that minimizes $-\ln\mathcal{L}(\bm{x}|\bm{\theta})$.

To quantify the uncertainties of the best-fit parameters, $\bm{\theta}_{\rm MLE}$, we adopt the nonparametric bootstrap method (Efron & Tibshirani, 1994; Burnham & Anderson, 2002). By repeatedly taking a random sample of size $N_{\rm tot}$ with replacement from the actual mass data, we first generate $N_{\rm bs}$ independent bootstrap samples. They serve as a proxy for a set of $N_{\rm bs}$ independent real samples from the same mass distribution, because taking extra data (from additional simulations) is too costly. The MLE is then employed to fit the model PDF to each bootstrap sample to obtain the best-fit parameters, $\bm{\theta}_{\mathrm{bs},k}$ ($k=1,\cdots,N_{\rm bs}$). Parameter uncertainties expected from real samples are estimated by calculating the distance between $\bm{\theta}_{\rm MLE}$ and the 84th and 16th percentiles of the distribution of $\bm{\theta}_{\mathrm{bs},k}$, i.e. $\bm{\theta}_{\mathrm{bs},k}^{84\%}$ and $\bm{\theta}_{\mathrm{bs},k}^{16\%}$, as

[TABLE]

Efron & Tibshirani (1994) have shown that this bootstrap method works reasonably well if $N_{\rm bs}$ is large (e.g., $>$1000). In this work, we fix $N_{\rm bs}=10000$. Appendix A gives a model fitting example in detail.

Our maximum likelihood estimator is similar to what used in Simon et al. (2016, 2017); Abod et al. (2019) but we use bootstrapping to estimate parameter uncertainties. A different method of parameter estimation is used by Johansen et al. (2015); Schäfer et al. (2017), etc., who applied curve fitting routines to the CDF. While maximizing the likelihood functions seems more fundamental to us, we make no claim that our method is actually superior. We conduct tests to recover the slope of a single power law distribution from randomly generated data using these two methods. This test yields non-identical results, which confirms the methods are not equivalent, but offers no evidence to clearly favor either method. A more rigorous investigation of this statistical issue could be warranted.

3.2 Statistical Models

We now describe the seven statistical models that are used to fit the mass distribution of planetesimals. To focus on the shapes of these models, this section only gives their basic functional forms. All the normalization coefficients and the full functional expressions are put in the Appendix B. For convenience, we define $K$ as the number of parameters in a model.

The first three models below are presented as CDFs. We simply convert them to their corresponding PDFs to apply our maximum likelihood estimator. However, it is natural to expect, if the planetesimal masses arise probabilistically, that a continuous PDF would be a more physical description. The reason to consider the CDFs is that some appeared previously in the literature (the first two models) and one of our runs shows visual evidence for a kink in the CDF (the third model). Because this kink gives a discontinuity in the PDF, it is arguably unphysical, but in this work we only examine statistical robustness, as no comprehensive physical theory for the distribution of masses exists.

Simply Tapered Power Law The concave downward profile of the CDFs of clump masses (see Fig. 2) suggests a power law distribution with exponentially tapering (Abod et al., 2019)

[TABLE]

where $M_{\rm exp}$ is the characteristic mass scale and $c_{1}$ is the renormalization coefficient (see Appendix B, same for coefficients below). This model has two free parameters ($K=2$) and $\bm{\theta}=(\alpha,M_{\rm exp})$, with $M_{\rm min}\leqslant M_{\rm exp}\leqslant M_{\rm max}$ but no constraints on $\alpha$. 2. 2.

Variably Tapered Power Law In addition to the first model, this model frees the tapering power by adding one more free parameter $\beta$ inside the exponential function (Schäfer et al., 2017),

[TABLE]

where $\bm{\theta}=(\alpha,\beta,M_{\rm exp})$ and $K=3$. This model requires that at least one of $\alpha$ and $\beta$ is positive, and again $M_{\rm min}\leqslant M_{\rm exp}\leqslant M_{\rm max}$. 3. 3.

Broken Cumulative Power Law A broken cumulative power law distribution connects two power law segments with different slopes in the cumulative distribution. It also manifests different behaviors at different mass ranges

[TABLE]

where $M_{\rm br}$ denotes the characteristic mass scale at which the slope breaks. This model has three free parameters ($K=3$) and $\bm{\theta}=(\alpha_{1},\alpha_{2},M_{\rm br})$. There is no constraints on $\alpha_{1}$, but $\alpha_{2}$ needs to be positive, and $M_{\rm min}\leqslant M_{\rm br}\leqslant M_{\rm max}$. 4. 4.

Truncated Power Law Schäfer et al. (2017) also tested a truncated power law model

[TABLE]

where an upper bound $M_{\rm tr}$ truncates the PDF (and CDF). In this model, $\bm{\theta}=(\alpha,M_{\rm tr})$, $K=2$, $\alpha>0$, and $M_{\rm tr}\geqslant M_{\rm max}$. Here the power law exponent in PDF becomes $-\alpha-1$ because the exponent in the corresponding CDF is $-\alpha$. Furthermore, it is easy to show that the PDF monotonically increases with increasing $M_{\rm tr}\leq M_{\rm max}$, and hence $-\ln\mathcal{L}(\bm{x}|\bm{\theta})$ minimizes if and only if $M_{\rm tr}=M_{\rm max}$. 5. 5.

Broken Power Law Another compelling possibility is the broken power law distribution333Not to be confused with the broken cumulative power law.. The corresponding PDF consists of two different power law segments, leading to a smooth transition in the CDF near the breaking point

[TABLE]

This model has three free parameters ($K=3$), $\bm{\theta}=(\alpha_{1},\alpha_{2},M_{\rm br})$, where $\alpha_{1}$ has no limits, $\alpha_{2}>0$, and $M_{\rm min}\leqslant M_{\rm br}\leqslant M_{\rm max}$. When $\alpha_{2}\to+\infty$, this model reverts to the Truncated Power Law model. 6. 6.

Truncated Broken Power Law This model complicates the last model by adding a truncation to the PDF

[TABLE]

This model has four free parameters ($K=4$) and $\bm{\theta}=(\alpha_{1},\alpha_{2},M_{\rm br},M_{\rm tr})$, where $\alpha_{1}$ and $\alpha_{2}$ has no limits, $M_{\rm min}\leqslant M_{\rm br}\leqslant M_{\rm max}\leqslant M_{\rm tr}$. Similar to the Truncated Power Law model, the PDF monotonically decreases with $M_{\rm tr}$ and $-\ln\mathcal{L}(\bm{x}|\bm{\theta})$ minimizes when $M_{\rm tr}=M_{\rm max}$. 7. 7.

Three-segment Power Law We take a step further to consider another broken power law distribution but with three segments in the PDF,

[TABLE]

This model has five free parameters ($K=5$) and $\bm{\theta}=(\alpha_{1},\alpha_{2},\alpha_{3},M_{\rm br1},M_{\rm br2})$. Both $\alpha_{1}$ and $\alpha_{2}$ have no boundaries, but $\alpha_{3}>0$ and $M_{\rm min}\leqslant M_{\rm br1}\leqslant M_{\rm br2}\leqslant M_{\rm max}$. When $\alpha_{3}\to+\infty$, this model reverts to the Truncated Broken Power Law model.

We choose these seven statistical models as candidates since they have been previously used to fit the planetesimal mass function or are commonly applied to fit top-heavy mass distributions. Other models are also certainly possible, but are beyond the scope of this paper. Note that all the models above are transformed to a PDF function of $x$ (see Table 6) to be used in our MLE.

3.3 Model Selection Criteria

Out next goal is to select the statistical models that best represent simulation data. Models with more parameters (larger $K$) have the flexibility to provide closer fits to the data, i.e. higher likelihood values. Often, a well-chosen function with fewer parameters can provide a better fit than a different function with more parameters. The much larger concern is the opposite case, where a more complex model does not better represent reality, but merely overfits statistical fluctuations in the data.

For the problem of planetesimal formation by the streaming instability, this statistical concern is relevant. The high mass tail of the planetesimal distribution is very significant, but with low numbers of high mass clumps in any simulation, the risk of statistical fluctuations impacting model fitting is potentially high.

To address these issues, we first review two of the most commonly-used model selection criteria and then introduce a selection criterion that we develop independently, motivated by the nonparametric bootstrap method.

3.3.1 Information Criteria

The most commonly used model selection criteria are (i) the Bayesian Information Criterion (BIC) (Kass & Raftery, 1995)

[TABLE]

and (ii) the Akaike Information Criterion (AIC) (Akaike, 1974)

[TABLE]

Both the BIC and the AIC involve the calculations of the log-likelihood $-2\ln{\mathcal{L}}$, which are affected by arbitrary constants and the sample size. Thus, the individual BIC/AIC values are not significant and the relative differences between models

[TABLE]

are more important, where BICmin/AICmin is the minimum of the BIC/AIC values of all the model candidates. In this way, the preferred model naturally has $\Delta_{\rm BIC}=0$ and other models have positive $\Delta_{\rm BIC}$’s (similar for AIC). To interpret $\Delta_{\rm BIC}$ and $\Delta_{\rm AIC}$ quantitatively in model selection, we follow the conventional categorical guidelines in Tables 2 and 3.

Formally, the value of $\Delta_{\rm BIC}$ represents the complexity-corrected likelihood ratio in a natural logarithmic scale, or the evidence provided by the data in favor of the preferred statistical model over another model (Kass & Raftery, 1995). The value of $\Delta_{\rm AIC}$ measures the Kullback–Leibler distance, or the information lost when a less preferred model is used to approximate the true distribution (Burnham & Anderson, 2002). For further discussions on the differences between the BIC and the AIC, we refer the reader to Burnham & Anderson (2002) and Burnham & Anderson (2004).

These two criteria put different weights on the penalty on the number of parameters, $K$, which becomes quite significant for large $N_{\rm tot}$, and which can lead to different results in model selection. It is difficult (for us) to determine which information criterion is more appropriate, or indeed if either is reliable. More complex and computationally methods exist to assess the complexity penalty based not simply on the number of free parameters and/or data points, but the actual geometry of the model space (e.g., Fisher Information Approximation, Ly et al., 2017). However these methods were beyond the scope of this work. Instead we describe an alternate model selection method below which we compare to the conventional AIC/BIC methods.

3.3.2 Bootstrap Model Selection

Motivated by concerns about the applicability of standard model selection techniques (BIC and AIC, discussed above), we consider an alternative method where the complexity penalty is not given as a fixed, simple function of the number of parameters but instead is generated automatically by bootstrapping.

Inspired by the nonparametric bootstrap method for uncertainty estimation (described in Section 3.1), we again consider all the bootstrap samples a good proxy for mass distributions from $N_{\rm bs}$ independent simulations, which in reality are too computationally costly expensive to be conducted. Through such a proxy, the median likelihood of all the bootstrap samples given the best-fit parameters can be used as a model selection criterion

[TABLE]

where BMS stands for Bootstrap Model Selection, $\bm{x}_{k}$ is the $k$-th bootstrap sample, and the factor of $2$ is chosen for similarity to AIC/BIC. This empirical criterion represents to what extent the best-fit parameters can explain/describe other samples drawn from the same mass distribution. Also, it naturally penalizes more complex models that tend to overfit data because they yield poorer fits to those bootstrap samples that deviate farther from the original data. In the following work, we therefore also consider

[TABLE]

as one of our model selection metrics and follow the similar conventional categorical guidelines. The comparison between BMS and the commonly-used BIC/AIC are discussed in the following sections.

4 The Initial Planetesimal Mass Distribution

In this section, we analyze and compare our two high resolution simulations (Run I and II, see Table 1) that have different physical parameters $(\uptau_{\rm s},Z)$. We fit seven statistical models (described in Section 3.2 and Appendix B) to the simulated mass distribution of planetesimals identified by PLAN (described in Section 2.3). In this section, we first present the fitting results and then the preferred models by model selection criteria. Section 4.3 describes an interesting turnover in the fitted mass distribution of Run I. In the end, Section 4.4 compares our results to recent studies.

4.1 Maximum Likelihood Fitting Results

Tables 4 and 5 list the best-fit likelihood and parameters for all the statistical models for Run I and II, respectively. The mass distribution used in fitting is extracted at $t_{\rm sg}=7.5/\Omega_{0}$ for Run I and at $t_{\rm sg}=7.6/\Omega_{0}$444Run I is the same simulation snapshot presented in Fig. 1, panel 2 in Simon et al. (2017), but re-analyzed here by PLAN. for Run II. By this time hundreds of planetesimals have formed, and the time evolution of the mass distribution has slowed. As shown in the two tables, for each model, the best-fit parameters for the two simulations are statistically different and are outside the uncertainty ranges of each other, indicating two different mass distributions are produced.

Fig. 3 visualizes all the resulting model CDFs for both of our simulations. They produce a broad mass distribution spanning more than three orders of magnitude in mass. However, Run I and II cover different mass regimes due to the different choices of physical parameters and hence the disk conditions, which may also require differently shaped distribution functions to describe. We emphasize that a physical understanding of these differences is elusive, hence our focus on statistics.

For Run I, we find the best-fit first-segment power law slopes ($-\alpha_{1}$) for the last three models in Table 4 become positive, which is of great interest for understanding which planetesimal masses dominate in number counts, and will be discussed in Section 4.3. In addition, the third-segment power law slope ($-\alpha_{3}$) for the Three-segment Power Law model is extremely steep. As mentioned in Section 3.2, when $\alpha_{3}$ approaches $+\infty$, this model reverts to the Truncated Broken Power Law, with the other four parameters being identical between these two models.

In this paper, we do not further consider the mass distributions in other snapshots and the possible variations with time, which belongs to future work. Simon et al. (2017) used a single power law model to fit the mass distribution and found that the power law index remains relatively constant in time after an initially rapid collapse. Schäfer et al. (2017) used the Variably Tapered Power Law model and also found that the power law index, characteristic mass scale, and tapering exponent all remain approximately constant in time for five orbital periods. Based on these results, we do not expect the mass distribution in our simulations evolve rapidly, i. e. on dynamical timescales, after the snapshots.

4.2 Model Selection

Our model selection analyses are based on the information criteria described in Section 4.2, i.e. $\Delta_{\rm BMS}$, $\Delta_{\rm BIC}$, and $\Delta_{\rm AIC}$ values, presented in Tables 4 and 5. We find that more complex models ($K\geqslant 3$) are generally much preferred than the other two simpler models ($K=2$), regardless of the model selection criteria. In other words, the Simply Tapered Power Law and the Truncated Power Law models are consistently less favored.

In the case of Run I, all the model selection criteria reach a consensus on choosing the $K=4$ Truncated Broken Power Law as the preferred model. However, there is disagreement on the strength of preference, as we now explain. In this specific case, the $K=5$ Three-segment Power Law model reverts to the preferred $K=4$ model, because the high mass power law is very steep (effectively a truncation) with the other four parameters identical. Since the BMS does not count parameters, it does not distinguish between these identically shaped distributions. However the AIC and BIC both apply a penalty for the 5th parameter, which is much more significant for the BIC. While all model selection metrics prefer the $K=4$ Truncated Broken Power Law, the BIC method (only) finds that the evidence against a pair of $K=3$ models – the Broken Cumulative Power Law model ($\Delta_{\rm BMS}=1.2$) and the Broken Power Law model ($\Delta_{\rm BMS}=1.8$) – is not significant.

In the Run II case, both the BIC and the AIC designate the Broken Cumulative Power Law model as the preferred model, which only has a moderate complexity level among all the model candidates. However, the BMS prefers the Three-segment Power Law models slightly more but also substantially supports the Broken Cumulative Power Law model. The overall preference for a broken CDF model may not be surprising given that the mass distribution shows a kink in the CDF (see Fig. 3), but is less physical intuitive since the planetesimal masses are expected to arise probabilistically and broken PDFs have been observed in the size-frequency distributions of asteroids and Kuiper Belt Objects (KBOs) (e.g. Morbidelli et al., 2009; Fraser et al., 2010; Shankman et al., 2013; Singer et al., 2019). More work is needed to further understand whether our case is special or not.

Fig. 3 allows a visual inspection of our model fitting and selection results. Not surprisingly, more complex models generally better fit the mass distributions. The effect of the logarithmic number axis is worth emphasizing. A small deviation from the (logarithmically plotted) CDF at the low mass end is more statistically significant than a larger deviation at the high mass end, where numbers, especially cumulative numbers, are much higher. The advantage of CDF plots is that no binning choices are required. While no data features are lost to binning, a disadvantage is the difficulty in interpreting slopes at the low mass end (where the CDF is near unity).

Fig. 4 provides an alternate view of binned planetesimal numbers (masses) which are compared to the PDFs of the best fit models. The PDF of the Broken Cumulative Power Law model – fit to Run II – has a discontinuous value at the CDF bread. As noted earlier, a physical explanation for such break does not exist. Overall, the selected models are excellent fits to the binned data, especially when accounting for the error bars, which represent the Poisson noise on the expected number (and mass) of bodies per bin. For Run II, the two preferred models seem to represent the mass distribution equally well.

4.3 A Turnover in the Mass Distribution

In SI simulations to date, the low mass end of the PDF is described by a power law with $\alpha>0$, where $dN/d\ln(M)\propto M^{-\alpha}$ (most of these works, described below, use $p\equiv\alpha+1$). Such a slope cannot extend to arbitrarily low mass, because the total mass of planetesimals would diverge. Sufficiently high resolution simulations should solve this problem, by revealing a low mass tail with $\alpha<0$. Such a measurement would determine the mass of planetesimals formed by SI that initially dominate by number. We present here the first evidence of such a turnover.

In Fig. 4, the mass frequency distribution of Run I presents a turnover below $\sim 0.003M_{\rm G}$ (roughly $100$ km-sized objects for the disk model in Section 2.2). The number of clumps drops with decreasing mass, except for the lowest mass bin (more on this below).

Our preferred statistical model confirms this visual evidence. The $K=4$ Truncated Broken Power Law model has a positive low mass slope of $-\alpha_{1}=0.102$ (as does the in practice identical $K=5$ Three-segment Power Law model). Our bootstrap error estimates (see Table 4) confirm the significance of the positive slope. Moreover, the simpler $K=3$ Broken Power Law model (while not the most preferred model) also has a positive low mass slope ($-\alpha_{1}=0.079$), which is flatter, but agrees within statistical uncertainties.

The evidence for a mass turnover (i.e. a positive slope at low masses) in Run I is fairly compelling. However, Run I also has an increase in the number of bodies in the lowest mass bin in Fig. 4. It is unclear if this increase has a physical origin, though we suspect that insufficient resolution on the smallest scales is an issue. We note that the lower resolution Run II shows an uptick in the planetesimal numbers in the two lowest mass bins.

Future studies with higher resolution are required to better resolve the low mass planetesimals and better characterize the inevitable turnover in the gravitational collapse mass function at low masses.

By contrast, for Run II resolution is not sufficient for evidence of a low mass turnover. However, one of the preferred statistical model (the $K=5$ Three-segment Power Law model) reveals a positive slope at intermediate masses ($-\alpha_{2}=0.362$). Whether this slope extends to the low-mass end and whether the uptick of the two lowest mass bins is numerical again require higher resolution studies. Moreover, the binned mass distribution of Run II shows a much flatter slope at the high mass regime than that of Run I, which is also shown by the CDFs in Fig. 3. This comparison further demonstrates that the different physical conditions of our two simulations produce different mass distributions.

4.4 Comparison with Previous Studies

In this section, we compare our fitting results to two recent works on planetesimal mass distribution that considered models beyond a single power law.

Schäfer et al. (2017) ran a suite of simulations with different resolutions and box sizes, fixing the physical parameters, $(\uptau_{\rm s},Z,\Pi,\tilde{G})=(0.314,0.02,0.05,0.318)$, i.e. similar to our Run II with weaker self-gravity. They considered a two-parameter Truncated Power Law model and a three-parameter Variably Tapered Power Law model and found the latter provides a better fit. Our analysis did not favor either of these models, especially for the similar run II. For the tapering exponent (of the Variably Tapered Power Law), they fit $\beta\simeq 0.3$ – $0.4$, similar to our results ($0.298$ and $0.527$ for Run I and II, respectively). Schäfer et al. (2017) explained that with limited resolution ($\Delta x=H/640$ at best), their simulations did not produce enough planetesimals to constrain the shape of the CDF in the power law part, and thus the values of $\alpha$ and $M_{\rm exp}$. Their work used a different code, PENCIL, and also used a sink-particle algorithm to handle bound clumps. These differences are a useful check on numerical robustness.

Abod et al. (2019) used high-resolution ($\Delta x=H/2560$) simulations with $(\uptau_{\rm s},Z,\tilde{G})=(0.05,0.1,0.02)$ – i.e. smaller solids in a lower mass disk than our runs – to study how the initial mass distribution of planetesimals depends on the pressure gradient, $\Pi$, with values from [math] to $0.1$. They found that the planetesimal mass function depends at most weakly on $\Pi$. Abod et al. (2019) used a two-parameter Simply Tapered Power Law model to fit the simulation data. Our analysis did not prefer this model, though it does have an advantage of simplicity. They fit the power law exponent $\alpha\approx 0.3$ and the characteristic mass scale of $\sim 0.3M_{\rm G}$ when $\Pi=0.05$. Our results give similar power law slopes ($\alpha=0.208$ for Run I and $0.388$ for Run II) and, for Run I a similar characteristic mass scale ($M_{\rm exp}=0.1385M_{\rm G}$). Our Run II fit, $M_{\rm exp}(=11.2531M_{\rm G})$, differs by a factor of $\sim 81$, for reasons that are not yet clear.

Since there are always more than one difference in the physical conditions and only limited model candidates are considered, these comparisons are somewhat inconclusive. Though costly, more parameter studies are needed to understand how the initial planetesimal mass function varies with each of the four physical parameters $(\uptau_{\rm s},Z,\Pi,\tilde{G})$ and eventually how these parameter dependencies couple.

5 Discussions And Conclusions

We investigate the mass distribution of planetesimals formed in high-resolution SI simulations. This mass distribution is of great astrophysical interest since it provides insights for the observations of small bodies in the Solar System (e.g., Cold Classical Kuiper Belt Objects, Nesvorný et al., 2019, etc.) as well as for the modeling of subsequent protoplanet formation (e.g. Liu et al., 2019).

In this work, we conduct SI simulations including particle self-gravity with the highest resolution to date, which produce broad mass distributions of planetesimals. While such distributions are top-heavy in mass for all numerical resolution choices, higher resolution simulations probe the lower-mass tail that dominates planetesimal numbers. We also develop and apply a new clump-finding tool, PLAN (described in Section 2.3), to accurately identify self-bound clumps in our simulations and extract their mass distributions. PLAN was used in previous work (Abod et al., 2019; Nesvorný et al., 2019), but the details of the method are presented here.

We fit the mass distribution to statistical models with different parameterizations (described in Section 3.2. To determine and select the preferred model from simulation data is a difficult statistical art, especially when different model candidates have different numbers of parameters. Thus, this work considers a variety of model selection criteria: the commonly-used BIC and AIC, as well as a bootstrap model selection method that we call BMS.

Based on our analyses, we find that

• •

Run I is best described by a four-parameter model, the Truncated Broken Power Law.

• •

For Run II (with smaller solids and a lower solid abundance than Run I) the preferred model depends on the model selection criterion used. The AIC and BIC prefer a three-parameter Broken Cumulative Power Law. The BMS selects a five-parameter model, the Three-Segment Power Law.

The interpretations and conclusions are drawn as follows:

The initial mass distribution of planetesimals formed by the streaming instability is shown to be numerically robust for the high mass regimes, and is most probably also robust at lower masses. Simulations with different numerical resolution (Run I and an equivalent run with lower resolution) show a similar mass distribution at the high mass end. Higher resolution gives a correction at intermediate masses and an extension to lower masses. 2. 2.

For different physical conditions, the initial mass distribution is not universal. While all cases produce top-heavy mass distributions with similar overall shapes, simulations with different physical parameters produce statistically different mass distributions. Fitting the same model to different runs often yields different best-fit parameters, e.g. power law slopes and characteristic mass scales. Moreover, the preferred models for different runs have different functional forms. More work and more high-resolution simulations are needed to better understand the initial mass distribution. 3. 3.

Our preferred models were not previously considered in the literature. We analyze the models that were used in previous studies, and find alternate models which rank higher by all model selection criteria. We make no claim to have found the optimal model, which we may not have considered and which may change as simulation data improves. 4. 4.

We find evidence for a turnover in the mass frequency distribution at the low mass end. This evidence is most prominent for Run I, where the PDF of the logarithmic masses transitions to a positive slope below $M\sim 0.003M_{\rm G}$ at roughly 2-$\sigma$ significance in the estimated slope. There is also statistical evidence for a turnover in the case of Run II, but only at intermediate masses. To better characterize the turnover of initial planetesimal mass distributions, higher resolution simulations are required. 5. 5.

The most complex model is not always selected as the preferred model. This result emphasizes the importance of applying complexity penalties for model selection. 6. 6.

Different model selection criteria disagree on both the absolute and relative rankings of different models. It is often difficult to rigorously justify a single model selection criteria for a given (astrophysical) application. Absent this justification, we generally recommend that multiple selection criteria be considered to increase confidence in model selection analyses.

Nesvorný et al. (2019) recently find that the clumps formed by the SI possess excess angular momenta and are likely to form binaries or multiple systems. In that case, the mass distributions from our simulations may describe the mass function of binaries/systems, not individual objects. This finding introduces corrections to the overall mass distribution and also some uncertainties to our modeling, which are beyond the scope of this work. However, those corrections and uncertainties would be minor if all clumps eventually form equal-size binaries as proposed in Nesvorný et al. (2010).

Acknowledgements

We thank Kaitlin Kratter, Paola Pinilla, Philip Pinto, and Peter Behroozi for useful discussions. RL acknowledges support from NASA headquarters under the NASA Earth and Space Science Fellowship Program grant NNX16AP53H. ANY acknowledges partial support from NASA Astrophysics Theory Grant NNX17AK59G and from the NSF through grant AST-1616929.

Appendix A A Model Fitting Example

In this appendix section, we take the fitting of the Variably Tapered Power Law model to the data from Run I as an example and describe it in detail. First, we assume a uniform prior distribution of parameters and use emcee to explore the posterior likelihood sampling (see Fig. 5). In this example, our chains consist of 32 walkers and 6000 steps with 1000 burn-in steps, while the estimated autocorrelation time returned by emcee is only about 75 steps and is well blow the number of burn-in steps. From this MCMC result, the best-fit parameters and the log-likelihood are

[TABLE]

We now feed a set of initial guesses generated on a mesh grid centered on the $\bm{\theta}_{\rm MCMC}$ to the minimize method provided by the scipy.optimize package. Various minimization algorithms are then employed for further minimization, including “Nelder-Mead”, “Powell”, “CG”, “BFGS”, “Newton-CG”, “L-BFGS-B” and “TNC”. The latter five algorithms compute the gradient vector, $\partial\ln\mathcal{L}/\partial\bm{\theta}$, to converge more quickly to the solution. We have implemented the methods to compute the gradient vector and Hessian matrix for all seven statistical models and make them available on GitHub. In this specific example, all algorithms but “L-BFGS-B” converge at the final best-fit parameters with an acceptable tiny gradient vector

[TABLE]

where Eq. 17 is satisfied with a good numerical precision.

In the next step, we use the nonparametric bootstrap method to estimate the uncertainties of $\bm{\theta}_{\rm MLE}$. Fig. 6 shows all the bootstrap samples in the left panel. We again apply our MLE to obtain the best-fit parameters, $\bm{\theta}_{\mathrm{bs},k}$, for the $k$-th bootstrap sample ($k=1,\cdots,N_{\rm bs}$). All these best-fits are plotted in the right panel of Fig. 6. The uncertainties are then calculated based on Eq. 18

[TABLE]

Furthermore, the uncertainty of the characteristic mass scale, $M_{\rm exp}$, can be derived as

[TABLE]

Appendix B Model Coefficients and Full Functional Forms

In this section, we list all the renormalization coefficients for the statistical models in Section 3.2 and show their full functional forms in Table 6.

Simply Tapered Power Law

[TABLE] 2. 2.

Variably Tapered Power Law

[TABLE] 3. 3.

Broken Cumulative Power Law

[TABLE] 4. 4.

Truncated Power Law

[TABLE] 5. 5.

Broken Power Law

[TABLE] 6. 6.

Truncated Broken Power Law

[TABLE] 7. 7.

Three-segment Power Law

[TABLE]