Evidence for Universality in the Initial Planetesimal Mass Function

The streaming instability (Youdin & Goodman 2005) leads to clustering of aerodynamically coupled solids across a broad range of protoplanetary disk conditions (Johansen & Youdin 2007; Bai & Stone 2010a). Because the physical origin of the instability is closely tied to that of radial drift (Whipple 1972; Weidenschilling 1977)—the most widespread barrier to growth beyond small macroscopic sizes—there is a compelling circumstantial argument that streaming plays a major role in planetesimal formation. The simplest scenario is that coagulation and radial drift lead to local conditions that trigger the streaming instability, which then forms dense particle clumps that collapse gravitationally to form planetesimals (Johansen et al. 2007). More involved variants, in which persistent or transient disk structures (ice lines, zonal flows, vortices, dead zone edges, etc.) are prerequisites for streaming-initiated collapse, are also possible (Johansen et al. 2014; Armitage 2015).

The parameters that determine the operation of the streaming instability include the particle size (measured via the dimensionless stopping time τ), the ratio of the solid to gas surface density Z, which we colloquially refer to as "metallicity," and the degree of pressure support in the gas. These parameters vary with radius in the disk, and a successful theory of planetesimal formation must therefore work across a range of starting conditions. This requirement is readily satisfied by the streaming instability at the linear level, where a broad array of gas/particle mixtures are linearly unstable (Youdin & Goodman 2005). At the nonlinear level, simulations for $\tau \sim 1$ show that the metallicity Z needs to exceed the nominal dust-to-gas ratio of 0.01 before the instability produces the strong clumping that precedes collapse (Johansen et al. 2009), but provided high Z can be reached, particles with ${10}^{-3}\lesssim \tau \lesssim 5$ are viable progenitors (Carrera et al. 2015; Yang et al. 2016).⁶ The allowable range of τ may be more strongly restricted at the low end by intrinsic turbulence in the gas, though recent theoretical results (Simon et al. 2015) and observations (Flaherty et al. 2015) suggest that at least some disks may be less turbulent than was previously thought.

In this Letter, we investigate whether the outcome of streaming-initiated collapse is universal, in the sense of forming an initial mass function of planetesimals whose shape is independent of the aerodynamic properties of the particles that participate in the instability. High-resolution simulations have shown that the prediction for the initial mass function is a power law, with a cutoff at high masses that is set by the local mass of solids that participates in the instability (Johansen et al. 2012; Simon et al. 2016; Schäfer et al. 2017). The existing simulations, however, have focused on a range of τ that is much smaller than that realized in actual disks. Here, we present results from simulations that span a range of τ between 0.006 and 2. We analyze the nonlinear particle structures and the mass function of collapsed objects produced in the simulations, and we show that any deviations from universality across this range of parameters must be small.

Our results are based on supplementing the highest-resolution simulation (512³ gas zones, $1.536\times {10}^{8}$ particles) from Simon et al. (2016) with two further calculations initialized with differing values of the stopping time and metallicity. We work in a locally Cartesian "shearing box" with radial, azimuthal, and vertical coordinates $(x,y,z)$ and model the fluid as an isothermal gas with pressure $P=\rho {c}_{s}^{2}$. In the rotating frame, the hydrodynamics of the gas is described by

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}+{\rm{\nabla }}\cdot (\rho {\boldsymbol{u}})=0,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}\displaystyle \frac{\partial \rho {\boldsymbol{u}}}{\partial t}+{\rm{\nabla }}\cdot (\rho {\boldsymbol{uu}}+P{\boldsymbol{I}}) & = & 3\rho {{\rm{\Omega }}}^{2}x\hat{{\boldsymbol{x}}}-\rho {{\rm{\Omega }}}^{2}z\hat{{\boldsymbol{z}}}\\ & & -2{\boldsymbol{\Omega }}\times \rho {\boldsymbol{u}}+{\rho }_{p}\displaystyle \frac{{\boldsymbol{v}}-{\boldsymbol{u}}}{{t}_{\mathrm{stop}}}.\end{eqnarray} \tag{ 2 }$$

Here, ${\boldsymbol{u}}$ and ${\boldsymbol{v}}$ are velocities of the gas and particle fluids, respectively; ${\boldsymbol{I}}$ is the identity matrix; ${\rho }_{p}$ is the particle mass density; ${t}_{\mathrm{stop}}$ is the dimensional stopping time; and Ω the angular velocity that is assumed to be Keplerian.

The solids are represented as discrete super-particles (Youdin & Johansen 2007) i = 1, ..., N. In the shearing box frame, they are subject to the fictitious forces arising from the rotating coordinate system, the vertical component of stellar gravity, a force representing the dynamical effect of the radial pressure gradient, and self-gravity,

$$\begin{eqnarray}\displaystyle \frac{d{{\boldsymbol{v}}}_{i}}{{dt}} & = & 2{{\boldsymbol{v}}}_{i}\times {\boldsymbol{\Omega }}+3{{\rm{\Omega }}}^{2}x\hat{{\boldsymbol{x}}}-{{\rm{\Omega }}}^{2}z\hat{{\boldsymbol{z}}}\\ & & -\displaystyle \frac{{{\boldsymbol{v}}}_{i}-{\boldsymbol{u}}}{{t}_{\mathrm{stop}}}-{{\boldsymbol{F}}}_{{\rm{p}}}+{{\boldsymbol{F}}}_{{\rm{g}}}.\end{eqnarray} \tag{ 3 }$$

A shearing box representation necessarily has no radial pressure gradient across the domain. To induce radial drift, we instead impose a constant inward force on the particles ${{\boldsymbol{F}}}_{{\rm{p}}}=-2\eta {v}_{K}{\rm{\Omega }}\hat{{\boldsymbol{x}}}$, where $\eta {v}_{K}$ is the deviation from Keplerian orbital velocity due to radial pressure gradients in the physical system. The term ${{\boldsymbol{F}}}_{{\rm{g}}}$ is the particle self-gravity term, derived from a solution to Poisson's equation:

$$\begin{eqnarray}&&{{\rm{\nabla }}}^{2}{{\rm{\Phi }}}_{{\rm{p}}}=4\pi G{\rho }_{{\rm{p}}}\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{{\boldsymbol{F}}}_{{\rm{g}}}=-{\rm{\nabla }}{{\rm{\Phi }}}_{{\rm{p}}}.\end{eqnarray} \tag{ 5 }$$

Where necessary, interpolation is used to map fluid quantities defined on a fixed grid to the locations of individual particles, and vice versa (Bai & Stone 2010a; Simon et al. 2016).

The coupled particle–gas system is solved using the athena code (Stone et al. 2008), in practice in a slightly different form that subtracts off the local orbital advection velocity. Established methods are used to implement the shearing boundary conditions (Stone & Gardiner 2010), particle integration (Bai & Stone 2010a), and particle self-gravity (Koyama & Ostriker 2009; Simon et al. 2016), which is calculated using a Particle–Mesh scheme.

The streaming instability is characterized by the dimensionless stopping time of the participating particles,

$$\begin{eqnarray}&&\tau \equiv {t}_{\mathrm{stop}}{\rm{\Omega }},\end{eqnarray} \tag{ 6 }$$

the local metallicity, defined via the ratio of particle to gas surface densities,

$$\begin{eqnarray}&&Z\equiv \displaystyle \frac{{{\rm{\Sigma }}}_{p}}{{{\rm{\Sigma }}}_{g}},\end{eqnarray} \tag{ 7 }$$

a pressure gradient parameter,

$$\begin{eqnarray}&&{\rm{\Pi }}\equiv \eta {v}_{K}/{c}_{s},\end{eqnarray} \tag{ 8 }$$

and a parameter describing the the relative strength of self-gravity versus tidal shear,

$$\begin{eqnarray}&&\tilde{G}\equiv \displaystyle \frac{4\pi G{\rho }_{0}}{{{\rm{\Omega }}}^{2}}.\end{eqnarray} \tag{ 9 }$$

Here, ${\rho }_{0}$ is the midplane gas density. We fix ${\rm{\Pi }}=0.05$, and $\tilde{G}=0.05$, equivalent to a gaseous Toomre $Q\approx 32$. Our three runs sample a range of stopping times and metallicities.

• 1.  
  $\tau =0.3$, Z = 0.02. This is the highest-resolution run from Simon et al. (2016).
• 2.  
  $\tau =2$, Z = 0.1. This is a distinct parameter set that remains relatively easy to simulate. Radial drift means that it is hard to attain the $\tau \gt 1$ regime in real disks (Birnstiel et al. 2012), but it could be physically relevant if solids grow in particle traps in the outer disk (Pinilla et al. 2012).
• 3.  
  $\tau =0.006$, Z = 0.1. These parameters approach those expected for millimeter- to centimeter-sized solids that drift and pile up in the dense region interior to the snow line (Youdin & Chiang 2004).

The results of Carrera et al. (2015) and Yang et al. (2016) imply that all three runs ought to result in strong clustering that is unstable to gravitational collapse, and this expectation is confirmed.

All other numerical details follow those described in Simon et al. (2016). The simulations use a cubical box of size ${(0.2H)}^{3}$, where $H={c}_{s}/{\rm{\Omega }}$, 512³ grid zones for the gas, and $N\,=1.536\times {10}^{8}$ particles.

In common with most prior work, the simulations are run in two steps. Initially, we evolve the aerodynamically coupled system in the absence of self-gravity. When the system has attained a saturated state (defined as when the maximum particle mass density is statistically constant in time, with moderate stochastic fluctuations about this constant value) self-gravity is turned on and collapse to planetesimals proceeds. For the $\tau =0.006,0.03,$ and 2 runs, self-gravity is switched on at $346.3\ {{\rm{\Omega }}}^{-1}$, $110\ {{\rm{\Omega }}}^{-1}$, and $37\ {{\rm{\Omega }}}^{-1}$ (at which point the rms scale heights for the particles are $0.01H,0.006H$ and $0.002H$), respectively. This two-step procedure is followed to reduce the computational expense of the runs. It is justified by tests (admittedly at lower resolution) that find no evidence that the outcome depends on the timing of self-gravity onset (Simon et al. 2016).

In addition, we have chosen to stop each simulation when ∼100 planetesimals have formed. This number is a compromise between statistical precision and computational cost. The chosen stopping times are $468.2\ {{\rm{\Omega }}}^{-1}$, $117.6\ {{\rm{\Omega }}}^{-1}$, and $45.8\ {{\rm{\Omega }}}^{-1}$ for the $\tau =0.006,0.03,$ and 2 runs, respectively. Upon examination of the particle mass surface density, it is possible that, at least for the $\tau =0.006$ and 0.3 runs, further planetesimal formation will occur. However, based on the evolution of the power-law index and minimum and maximum planetesimal masses (as discussed further below), we believe that further planetesimal formation will not appreciably change the mass functions.

Figure 1 shows projections of the surface density of solid material for the three simulation runs in the orbital plane after dense clumps have formed. Visually, they are quite distinct. The prominence of axisymmetric bands in the solid surface density decreases with increasing τ, and notably less material collapses promptly into bound clumps for the $\tau =0.006$ run. In the smaller τ runs, the planetesimals form primarily in two azimuthally extended bands. In contrast, the planetesimals in the $\tau =2$ run fill the box much more uniformly.

Zoom In Zoom Out Reset image size

From these snapshots, we use a clump-finding routine to identify and measure the bound masses ${M}_{{\rm{p}},i}$ of collapsed objects (we use an identical algorithm to Simon et al. 2016). For visual purposes, we compute an unsmoothed estimate of the differential mass function at masses corresponding to each planetesimal:

$$\begin{eqnarray}&&{\left.\displaystyle \frac{{dN}}{{{dM}}_{{\rm{p}}}}\right|}_{i}=\displaystyle \frac{2}{{M}_{{\rm{p}},i+1}-{M}_{{\rm{p}},i-1}}.\end{eqnarray} \tag{ 10 }$$

The resulting mass functions are plotted in the left panel of Figure 2. The shape of the mass functions is generally consistent with a single power law for all three runs, though, as we discuss in more detail below, there is likely a cutoff at high mass. The single power-law fit is clearest in the intermediate- and high-τ runs that form planetesimals with a broader range of masses. Assuming that the data are indeed drawn from a power-law distribution, ${dN}/{{dM}}_{p}\propto {M}_{{\rm{p}}}^{-p}$, we proceed to estimate the index p and error σ using a maximum likelihood estimator (Clauset et al. 2009). Given n planetesimals with masses ${M}_{{\rm{p}},i}\geqslant {M}_{{\rm{p}},\min }$, we have

$$\begin{eqnarray}&&p=1+n{\left[\displaystyle \sum _{i=1}^{n}\mathrm{ln}\left(\displaystyle \frac{{M}_{{\rm{p}},i}}{{M}_{{\rm{p}},\min }}\right)\right]}^{-1},\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&\sigma =\displaystyle \frac{p-1}{\sqrt{n}}.\end{eqnarray} \tag{ 12 }$$

The derived slopes and their associated errors are given in Table 1. Within 1–1.5σ all three runs are consistent with p = 1.6, as found in our previous work (Simon et al. 2016) and in Johansen et al. (2015). It is also worth noting that despite the consistent value of p, the total fraction of mass converted to planetesimals varies strongly with stopping time; $\sim 10 \%$, 40%, and 70% of the total mass in solids is converted to planetesimals for the $\tau =0.006,0.3,$ and 2 simulations, respectively. For the smaller millimeter- to centimeter-sized solids present in the inner regions of disks, our results suggest a relatively low planetesimal formation efficiency in these regions.

Zoom In Zoom Out Reset image size

Table 1.  Summary of Simulation Results

τ	Z	pa	mb
0.006	0.1	1.73 ± 0.11	1.50 ± 0.02
0.3	0.02	1.61 ± 0.08	1.65 ± 0.10
2.0	0.1	1.54 ± 0.04	1.62 ± 0.11

Notes.

^aSlope of the differential mass function. ^bSlope of the power spectrum of solids within the range $100\lt {kH}/2\pi \lt 1250$ prior to turning on self-gravity.

Download table as:  ASCIITypeset image

We also calculate the cumulative mass function, as shown in the right panel of Figure 2. A power-law index of −0.6 clearly agrees with the mass function slope for low-mass objects. This suggests that the differential mass functions are well fit by a single power law because there is a significantly larger number of planetesimals at small masses, thus weighting the fit toward the small-mass end. Therefore, while in reality the mass distribution will have a cutoff at some high-mass value (which depends on τ), the differential distribution can be well fit by a single power law that is representative of masses away from the cutoff.

Finally, we have also examined the evolution of p and the maximum and minimum planetesimal masses (in all three simulations) for times after which a significant number of planetesimals have formed; the results are shown in Figure 3. Note that due to limitations with the clump finder algorithm (as described more fully in Simon et al. 2016), we smoothed the curves (with a boxcar average) to remove noise associated with this algorithm. As the figure shows, p and the minimum planetesimal mass are relatively constant in time,⁷ whereas the maximum planetesimal mass grows by a factor of order unity in each simulation, likely due to accretion of smaller particles and/or mergers with smaller planetesimals. Despite this growth, the mass function does not evolve appreciably once planetesimals have formed.

Zoom In Zoom Out Reset image size

Our ability to directly measure any potential dependence of the mass function on particle properties is limited by the relatively small number of collapsed clumps, but across the range of τ considered, we can bound deviations at the level of ${\rm{\Delta }}p\lesssim 0.1\mbox{--}0.2$. We can therefore exclude, with moderately high confidence, the possibility that the streaming instability might result in a steep mass function ($p\gt 2$) with most of the mass in the smallest planetesimals.

The similarity in the mass functions from the different runs is somewhat surprising given the visual differences in the large-scale particle structures that are collapsing (Figure 1). To identify possible differences on smaller spatial scales, we compute the power spectra of the pre-collapse particle density fields. The power spectrum maps uniquely to the mass function in the case where the density is a Gaussian random field (Press & Schechter 1974). In the more complex case relevant here, we use the power spectra only to test whether the nonlinear structures produced by the streaming instability prior to the onset of self-gravity are indifferent to the particle size.

Figure 4 shows the three-dimensional power spectra of particle mass density computed from a time slice just before self-gravity is switched on. The thickness of the particle layer at this stage varies substantially with τ (higher τ allows for greater settling). To minimize artifacts in the power spectra created by the different thicknesses, we compute ${\tilde{\rho }}_{p}(k)$ from the interpolated density field $\rho ({\boldsymbol{x}})$ in a midplane slice whose thickness is chosen to be smaller than the scale height of the thinnest particle layer (for $\tau =2$). The three-dimensional power spectra are then averaged over shells of constant $| {\boldsymbol{k}}|$.

Zoom In Zoom Out Reset image size

Up to normalization differences the power spectra for all three runs display a high level of similarity. From calculating the largest Hill radius in each run, we expect scales of ${kH}/2\pi \gtrsim 100$ to seed the collapse.⁸ Fitting a power law, ${\tilde{\rho }}_{p}(k)\propto {k}^{-m}$, to the data at ${kH}/2\pi \gt 100$, we find that $m\approx 1.5$–1.7 with reasonably large errors for the larger two τ values,⁹ with the precise values and their associated errors shown in Table 1. All three runs yield statistically consistent slopes, suggesting that the near-identical mass functions formed from those runs result from commonality in the small-scale particle structures formed in the non-self-gravitating nonlinear evolution of the streaming instability. It should be noted, however, that from a visual inspection (Figure 4) the shape and slope of the power spectrum for $\tau =0.006$ is significantly different (and better fit with a single power law) than for the higher τ cases, even though their fitted power-law slopes are statistically equal. It is possible that any such differences in the power spectra translate into differences in the mass function at a level that is smaller than our current measurement uncertainties.

We have presented results from a small number of high-resolution simulations of the streaming instability (Youdin & Goodman 2005) that model the gravitational collapse of overdense structures toward planetesimals. The simulations were initialized with values of the dimensionless stopping time and local metallicity that promote strong clustering and prompt collapse (Carrera et al. 2015; Yang et al. 2016). For $0.006\leqslant \tau \leqslant 2$ we find that the power-law part of the resulting mass function (${dN}/{{dM}}_{p}\propto {M}_{p}^{-p}$) has an approximately universal slope, $p\,\simeq 1.6$, consistent with that measured previously for particle sizes in the middle of this range (Johansen et al. 2012; Simon et al. 2016).

The similar planetesimal mass distributions in the different runs is somewhat surprising, given differences in the larger-scale geometry of the particle clustering. However, this similarity is consistent with the approximately equal slopes in the power spectra of particle clustering on the smaller scales relevant to collapse. It is possible that the observed differences in the shape of the power spectra would translate into small deviations from universality, but within the uncertainties outlined here, our results strongly support a top-heavy mass function for planetesimals. This finding is consistent with previous streaming instability simulations, but is now shown over a broader range of τ–Z parameter space.

Our ability to directly determine the predicted initial mass function is limited by small number statistics. The statistics could be improved by co-adding samples derived from independent runs, by increasing the spatial resolution, or by increasing the domain size. In Simon et al. (2016), we showed that the planetesimal mass distribution is essentially independent of the time at which self-gravity is turned on in the simulation. However, this independence still remains to be tested at higher resolutions and across the broader ranges of parameters considered here.

Where gravitational collapse is encountered in other astrophysical settings, notably in the hydrodynamic formation of stars (Bastian et al. 2010) and in the collisionless collapse of dark matter halos (Navarro et al. 1997), it is known to exhibit universal features. Planetesimals may in principle form from gravitational collapse via other routes, for example, from the direct collapse of dense particle layers (Goldreich & Ward 1973; Youdin & Shu 2002; Shi & Chiang 2013), secular gravitational instability (Youdin 2011; Takahashi & Inutsuka 2014), or when vortices accumulate solids (Barge & Sommeria 1995; Raettig et al. 2015). It is of interest to explore whether these processes lead to similar or identical top-heavy power-law mass functions to those found here, and hence whether constraints on planetesimal formation from the asteroid (Morbidelli et al. 2009) and Kuiper Belt (Nesvorný et al. 2010) test specifically the streaming instability, or rather a broader class of gravitational collapse scenarios. On the other hand, if the mass function is indeed intimately coupled to the nonlinear state of the streaming instability, turbulence driven by other means (e.g., the magnetorotational instability; Balbus & Hawley 1998) and imposed onto the streaming instability may fundamentally alter the mass function of planetesimals.

J.B.S. and P.J.A. acknowledge support from NASA through grants NNX13AI58G and NNX16AB42G, and from the NSF through grant AST 1313021. A.N.Y. acknowledges support from the NSF through grant AST 1616929. We thank Doug Lin and Katherine Kretke for useful discussions regarding this work. We also thank the anonymous referee whose comments greatly improved the quality of this Letter. The computations were performed on Stampede and Maverick at the Texas Advanced Computing Center through XSEDE grant TG-AST120062. We thank the Kavli Institute for Theoretical Physics, supported in part by the National Science Foundation under grant No. NSF PHY-1125915, for hospitality during the completion of the Letter.