The Mass and Size Distribution of Planetesimals Formed by the Streaming Instability. II. The Effect of the Radial Gas Pressure Gradient
Abstract
The streaming instability is a mechanism for concentrating solid particles in protoplanetary disks that can lead to gravitational collapse and planetesimal formation. The energy source that drives this instability is the radial pressure gradient in the disk. Despite its key role in producing particle clumping and determining critical length scales in the instability’s linear regime, the influence of the pressure gradient on planetesimal properties has not been examined in detail. In this work, we use high resolution simulations of the streaming instability that include particle selfgravity to study how the initial mass function of planetesimals depends on the radial pressure gradient. Fitting our results to a powerlaw, , we find that a single value of describes simulations in which the pressure gradient varies by more than a factor of two. An exponentially truncated powerlaw provides a better representation of the cumulative mass function, and fitting to this function we obtain a single low mass slope . The characteristic truncation mass is found to be of the order of . We exclude the cubic dependence of the characteristic mass with pressure gradient suggested by linear considerations, finding instead at most a weak and possibly not significant dependence. These results strengthen the case for a streamingderived initial mass function that depends at most weakly on the aerodynamic properties of the disk and participating solids. A simulation initialized with zero pressure gradient—which is not subject to the streaming instability—also yields a topheavy mass function but with a significantly different shape. We discuss the consistency of the theoretically predicted mass function with observations of smallbody populations in the Kuiper Belt, and describe implications for models of early stage planet formation.
Subject headings:
planets and satellites: formation — hydrodynamics — instabilities — turbulence — protoplanetary disks1. Introduction
An outstanding issue in planet formation is how small particles of size mm–cm grow to larger 1–100 km scale planetesimals. While growth through sticking is efficient for mmm scale particles (Blum, 2018), no viable path is known that leads to planetesimal formation via twobody collisions across all radii in the disk. In particular, silicate aggregates that may be representative of solids interior to the snow line fragment for mutual collision velocities (Güttler et al., 2010). For reasonable values of the disk turbulence parameter and sound speed this threshold is substantially smaller than the peak collision speeds, , predicted for a turbulent flow (Ormel & Cuzzi, 2007). Icy materials exhibit higher threshold fragmentation velocities (Wada et al., 2009; Gundlach & Blum, 2015), and given enough time collisions could form larger bodies. In the presence of radial drift (Whipple, 1972; Weidenschilling, 1977), however, even icy particle growth is limited to mmcm sizes by the competition between inward drift and growth (Birnstiel et al., 2012).
A solution to these problems emerges when one considers not only the aerodynamic drag on the particles from the gas, but also the feedback of the particles onto the gas. The resulting equilibrium solution for radial drift in a vertically unstratified disk (Nakagawa et al., 1986) is linearly unstable to the streaming instability (Youdin & Goodman, 2005), which leads to particle clumping (Youdin & Johansen, 2007; Johansen et al., 2007; Bai & Stone, 2010a). Under conditions that can plausibly, but not inevitably, be attained within protoplanetary disks, the clumping becomes strong enough (Carrera et al., 2015; Yang et al., 2017) that selfgravity takes over and a fraction of the solid mass becomes bound objects. The most massive of these objects could have the size of large Solar System asteroids or Kuiper Belt objects (Johansen et al., 2007, 2012, 2015; Simon et al., 2016, 2017; Schäfer et al., 2017).
Given the difficulty of directly observing planetesimals in the process of forming, one of the key routes to testing models is via their predictions for the initial mass or size distribution of planetesimals. The initial mass distribution is a direct input to models of the current size distribution of small bodies in the Solar System (Morbidelli et al., 2009), and it has an important, if indirect, impact on the predicted growth rates of giant planet cores (Pollack et al., 1996; Fortier et al., 2013), on the evolution of debris disks (Krivov & Booth, 2018), and on the abundance of interstellar planetesimals (MoroMartín et al., 2009; Raymond et al., 2018).
Simple considerations suggest that the streamingderived mass function might depend on as many as four independent parameters. The model linear problem (Youdin & Goodman, 2005), which considers aerodynamically coupled dust and gas within an unstratified disk with zero intrinsic turbulence, yields growth rates and unstable modes that are functions of the dimensionless stopping time , the solidtogas ratio , and the radial pressure gradient parameter (we defer to §2 the formal definition of these parameters, but note here that , where is the difference between the gas velocity and the local Keplerian speed). Although not considered as part of the model linear problem the strength of intrinsic disk turbulence may also play a role, both by controlling the thickness of settled particle layers (Dubrulle et al., 1995) and by acting as an effective diffusivity on small scales (Youdin & Goodman, 2005).
The relationship between the instability of the unstratified model system (specified by , , and ) and the properties of planetesimals forming in real disks is unclear. Several additional pieces of physics need to be considered. The vertical component of stellar gravity leads to particle settling, and although this can be balanced by turbulent diffusion (leading to an equilibrium state) the linear theory including stratification has not been worked out. The combination of particle selfgravity (which is also not part of the model streaming problem) and aerodynamic drag leads to secular gravitational instability (Ward, 1976; Youdin, 2011; Takahashi & Inutsuka, 2014), with distinct linear properties. Finally, gravitational collapse requires such high particle overdensities that expectations grounded in the linear physics may be a poor guide to planetesimal properties, even absent the extra physical effects discussed above. Indeed, numerical simulations suggest that some of the important linear parameters do not strongly modify the resulting mass function. Based on three highresolution simulations that spanned a range of stopping time , Simon et al. (2017) found no evidence for significant variations with . This led us to conjecture that the streaminginitiated planetesimal mass function might be universal, in the sense that some of its parameters (specifically, the low mass slope) are independent of the aerodynamic properties of the particles that participate in the instability.
Here, we extend our prior work to study how the planetesimal mass function depends on the radial pressure gradient in the disk. Since the linear scale of unstable modes in an unstratified disk without selfgravity is proportional to , there is a natural expectation that masses might scale as (Youdin & Goodman, 2005; Taki et al., 2016).
In addition to testing this expectation, there are also observational motivations. Although is a weak function of orbital radius for commonly adopted disk profiles, stronger variations are expected near the outer edge of the gas disk, which could correspond to the outer edge of the Kuiper Belt (Trujillo et al., 2001). Variations in would also be expected if a significant fraction of planetesimals form from particles trapped within local pressure maxima (Whipple, 1972; Pinilla et al., 2012; Dittrich et al., 2013), caused for example by zonal flows (Johansen et al., 2009a; Simon & Armitage, 2014; Bai & Stone, 2014; Béthune et al., 2017). Such traps may be responsible for the observed rings in a number of protoplanetary systems (e.g., ALMA Partnership et al., 2015; Andrews et al., 2016; Isella et al., 2016). If these rings are indeed sites for planetesimal formation, building a complete model of planetary system formation motivates us to understand the role of the radial pressure gradient in the birth of planetesimals.
This paper is organized as follows. In §2, we describe our numerical algorithm in detail, define our initial conditions and choice of parameters, and describe the diagnostics we employ to characterize planetesimal formation. We present our results in §3. §4 discusses the implications of those results, particularly in light of the observed size distribution of objects in the Kuiper Belt. We summarize and conclude in §5.
Run  ^{a}^{a}Pressure gradient parameter, defined by Equation 7  ^{b}^{b}Time (in units of ) at which particle selfgravity is initiated  ^{c}^{c}Best fit power law slope to single power law model  ^{d}^{d}Best fit power law slope to truncated power law model  ^{e}^{e}Best fit truncation mass in truncated power law model  ^{f}^{f}The mass at which , where is the total planetesimal mass greater than normalized to the total mass in planetesimals.  Mass Fraction 

in Planetesimals  
P0SG0  0  0  1.75 0.03  1.52 0.05  0.78 0.21  0.57  0.81 
P0.0375SG0  0.0375  0  1.64 0.02  1.28 0.03  0.26 0.02  0.31  0.66 
P0.05SG0  0.05  0  1.60 0.03  1.29 0.03  0.35 0.05  0.42  0.55 
P0.0625SG0  0.0625  0  1.59 0.03  1.31 0.04  0.32 0.04  0.35  0.29 
P0.075SG0  0.075  0  1.58 0.05  1.31 0.06  0.56 0.17  0.61  0.22 
P0.0875SG0  0.0875  0  1.58 0.07  1.33 0.09  0.82 0.34  1.0  0.12 
P0.1SG0Lz0.4H  0.1  0  1.49 0.07  1.21 0.07  0.44 0.12  0.52  0.073 
P0.0375SG72  0.0375  72  1.56 0.02  1.19 0.03  0.16 0.02  0.22  0.53 
P0.05SG80  0.05  80  1.64 0.03  1.30 0.05  0.24 0.04  0.28  0.28 
P0.075SG66  0.075  66  1.51 0.06  1.21 0.06  0.26 0.05  0.39  0.062 
2. Methods
We describe our numerical methods and simulation diagnostics in detail below. For readers familiar with the field the setup is essentially identical to prior work (Simon et al., 2016, 2017). We model a small but vertically stratified local volume of the disk, whose properties are specified by two parameters: the gas pressure gradient and a parameter describing the balance between selfgravity and shear. The solid component is specified by a single dimensionless stopping time and a ratio of the particle to gas surface densities. We solve the coupled system using Athena, modeling the gas as a compressible isothermal fluid (ignoring magentic fields and imposing no intrinsic disk turbulence) and the solids as an ensemble of individual particles. Particle selfgravity is implemented using a particlemesh method with shearing radial boundary conditions. There are two key differences from our previous work (Simon et al., 2016, 2017), both of which affect the analysis rather than the simulations themselves. We now identify gravitationally bound clumps of solids directly in three dimensions, rather than working with the projected surface density. We also apply a maximum likelihood estimator to fit our results to an exponentially truncated powerlaw, rather than the simple powerlaw that we restricted ourselves to before.
2.1. Numerical Algorithm
The simulations use the Athena code to solve the equations of coupled particle motion and hydrodynamics within the local, shearing box approximation. In this approximation we model a corotating disk patch whose size is small compared to the radial distance from the central object. This allows for the construction of a local Cartesian frame , which is defined in terms of disk’s cylindrical coordinate system as , , and . This domain corotates around the central object at angular velocity , which corresponds to the angular velocity at the center of the box (i.e., at ). A detailed description of the shearing box can be found in Hawley et al. (1995).
The equations of gas dynamics in this approximation are comprised of the continuity and momentum equations:
(1) 
(2)  
where is the mass density, the momentum density, the gas pressure, and is the identity matrix. The quantity is the shear parameter, defined as lnln, which we take to be , as appropriate for a Keplerian disk. From left to right, the source terms in equation (2) correspond to radial tidal forces (gravity and centrifugal), vertical gravity, the Coriolis force, and the feedback from the particle momentum onto the gas; is the mass density of particles, is their velocity, and is the timescale over which a particle will lose a factor of of its momentum due to gas drag. As we describe below, this feedback term is calculated at the location of every particle and then distributed to the gas grid points. Finally, these gas dynamics equations are supplemented by an equation of state , where is the isothermal sound speed.
The lefthand side of the momentum equation (i.e., the motion of the gas in the absence of source terms) is solved using a secondorder accurate Godunov fluxconservative method, with the dimensionally unsplit corner transport upwind method of Colella (1990) and the thirdorder in space piecewise parabolic method of Colella & Woodward (1984). A detailed description of these components of Athena along with various test problems is given in Gardiner & Stone (2005), Gardiner & Stone (2008), and Stone et al. (2008). There are several algorithmic additions employed to integrate these equations within the shearing box framework (thus handling the noninertial terms), including orbital advection (the background Keplerian velocity is subtracted and integrated analytically; Masset 2000; Johnson et al. 2008) and CrankNicholson differencing to preserve epicyclic energy to machine precision. Details for the implementation of these algorithms can be found in Stone & Gardiner (2010).
Particles are treated in Athena via the superparticle approach (i.e., each superparticle is a statistical representation of a number of smaller particles). The equation of motion for superparticle (hereafter, simply “particle” for simplicity) is given by
(3)  
where the prime denotes a frame in which the background shear velocity has been subtracted, as part of the orbital advection scheme mentioned above. From left to right, the source terms are the radial acceleration of the particles due to the Coriolis effect, the gravitational + centrifugal force, and radial drift; azimuthal motion due to the Coriolis effect; vertical motion due to the central star’s gravity; the gas drag; and the force due to particle selfgravity. The term accounts for the inward radial drift of particles resulting from a gas headwind, which in turn results from the radial pressure gradient. Thus, , which is a primary parameter of interest in this work, can also be described as the fraction of the Keplerian velocity by which the orbital velocity of particles is reduced (see Section 2.2). An actual radial gradient in the gas pressure is inconsistent with the radial (shearing) periodic boundary conditions in the shearing box model. To circumvent this issue, we follow Bai & Stone (2010b) and impose an inward force on the particles, resulting in the term as described above. As a result, both the particles and gas (azimuthal) velocities are shifted to slightly higher values (by ) than what would be present in a real disk, but the essential physics of differential motion between the gas and particles is accurately captured.
Equation (3) is solved using the algorithms described in Bai & Stone (2010b). We integrate the equations of motion with a semiimplicit integration method combined with a triangular shaped cloud (TSC) scheme to map the particle momentum feedback to the grid cell centers and inversely to interpolate the gas velocity to the particle locations (). A more thorough description of the particle integration algorithm, as well as test problems, can be found in Bai & Stone (2010b).
When particle selfgravity is activated, there is an additional force in equation 3, , which is found by solving Poisson’s equation for particle selfgravity. The details of solving Poisson’s equation can be found in Simon et al. (2016). Briefly, we map the mass density of the particles to the grid cell centers using the TSC method and then calculate the gravitational potential by shifting the radial boundaries to be purely periodic (see Simon et al. 2016 and Hawley et al. 1995 for more details), applying a 3D FFT to calculate the gravitational potential in Fourier space, transforming the potential to real space using another 3D FFT, and then reversing the azimuthal shift to the original nonperiodic frame. In our simulations, the vertical boundaries are open, and as such, we apply a Green’s function approach to the Poisson equation in the vertical direction (see Koyama & Ostriker, 2009; Simon et al., 2016). We then finite difference the potential to calculate the gravitational force at each grid cell center and then apply TSC to map these forces to the precise particle locations. Tests of this algorithm are shown in Simon et al. (2016).
The boundary conditions for the gas and particles are shearingperiodic in the radial direction (Hawley et al., 1995), purely periodic in the azimuthal direction, and a modified outflow condition in the vertical direction in which the gas density is extrapolated via an exponential function into the ghost zones (Simon et al., 2011; Li et al., 2018). This modified boundary condition prevents gas mass loss and spurious effects from developing near the vertical boundaries when the gas is in purely hydrostatic equilibrium. Of course, once the system begins to evolve, this hydrostatic equilibrium will not be maintained, and gas mass will be lost through the vertical boundaries (see Li et al., 2018). To ensure mass conservation, at each time step, we renormalize the gas density in every cell so that the total gas mass is conserved. The boundary conditions for the gravitational potential are essentially the same as the hydrodynamic variables; shearingperiodic in and purely periodic in . The vertical boundary conditions are open, and the potential is calculated in the ghost zones via a third order extrapolation.
2.2. Initial Conditions, Parameters and Simulation Setup
All of our simulations are initialized as follows. The gas is in a hydrostatic (Gaussian) vertical profile,
(4) 
where is the midplane gas density. We choose code units so that the standard gas parameters are unity, .
Four dimensionless quantities characterize the evolution of the streaming instability and subsequent gravitationally collapse: the dimensionless stopping time,
(5) 
which characterizes the aerodynamic interaction between a single particle species and the gas, the particle concentration,
(6) 
which is the ratio of the particle mass surface density to the gas surface density , and a radial pressure gradient parameter that accounts for the subKeplerian gas in real disks
(7) 
This radial pressure gradient produces a headwind on the particles, which has velocity a fraction of the Keplerian velocity.
There is an additional parameter that comes into play when selfgravity is included in the simulations: the relative importance of selfgravity and tidal shear, which we quantify as
(8) 
Physically, varying corresponds to a change in gas density (and thus the particle density at constant ) and/or the strength of tidal stretching (i.e., changing ). Another way to envision the importance of this parameter is to consider how it would change with distance from the central star. At larger radial distances, the tidal effects will be reduced compared to gravity, and thus will increase with radius.
In this work, we are primarily interested in the influence of on the properties of planetesimals. Thus, we fix the other parameters to be , , and (which corresponds to a Toomre parameter for the gas ). These particular choices are chosen to closely mimic reasonable disk parameters as much as possible. In particular, Fig. 1 shows and (assuming 2.5mm sized particles, which is on the order of the maximum size particles can reach in a typical protoplanetary disk; Blum & Wurm 2008; Birnstiel et al. 2011; Zsom et al. 2010) for a typical protoplanetary disk with surface density (assuming the disk has mass 0.05) and temperature profile . As the figures shows, near AU, and corresponds to mm sized particles at this radius. Thus, our values for these parameters are roughly consistent with that expected within the planet forming region given typical disk conditions. In principle, a range of would be more appropriate since disk solids will not be confined to one size. However, for the purposes of this work, in which we only vary , it is important to keep the other parameter choices as simple as possible. As shown by Bai & Stone (2010c), must be larger than for . While these simulations were twodimensional and included a range of values, we use them to further inform our parameter choices. Thus, we choose the conservative for our concentration parameter.
The particles are initially distributed uniformly in and with a Gaussian vertical profile. We choose the initial scale height of this vertical profile to be . We apply the Nakagawa–Sekiya–Hayashi (NSH) equilibrium solution in and (Nakagawa et al., 1986) in order to ensure that the system starts in as close to an equilibrium configuration as possible.^{1}^{1}1Of course, in , no such equilibrium exists. To seed the streaming instability, we add random noise to the particle locations in all three dimensions.
One potential complication from these initial conditions, and in particular from using a factor of 10 larger concentration parameter relative to solar, which is , is that once the particles sufficiently settle, the very thin, high density particle layer may overpower the KelvinHelmholtz instability that can in general prevent solids from settling (Cuzzi et al., 1993; Weidenschilling, 1995). If this is the case, the particle layer will undergo gravitational collapse into planetesimals via the secular gravitational instability (hereafter, SGI; Ward 1976; Youdin 2011; Takahashi & Inutsuka 2014), which here refers simply to a regime of gravitational instability that is substantially modified by aerodynamic coupling to the gas.
To circumvent this issue, our parameters, while physically motivated, were also chosen to make the growth rate of the streaming instability comparable to or larger than the settling rate (in the absence of turbulence) of particles, which is . At the given initial particle scale height, the ratio of particle mass density to gas mass density ; extrapolating very roughly from Fig. 1 in Youdin & Johansen (2007), the growth time of the streaming instability is on the scales which we can resolve. This growth time is comparable to, and slightly less than the timescale over which the particles settle. Therefore, turbulence induced via the streaming instability should prevent the direct gravitational collapse of particles from the dense particle layer.
In practice, however, we have found it difficult to guarantee that overdensities that lead to planetesimal formation result only from the streaming instability; in some of our simulations, the particle layer gets very thin before the streaming instability even sets in.^{2}^{2}2There are a number of simplifications in our estimates of the streaming growth rates, including the neglect of vertical gravity and a very approximate extrapolation from the linear analysis of Youdin & Johansen (2007). To test for this potential complication, we also run a simulation with . In such a simulation, the streaming instability vanishes, and any planetesimal formation must be from SGI. Thus, by comparing the evolution and planetesimal properties in the run with our other simulations, we can determine what role SGI is playing in our simulations.
One of the goals of this work is to determine what role, if any, the linear regime of the streaming instability plays in planetesimal formation. As a comparison with the simulations described thusfar, we have also carried out three additional simulations in which selfgravity remained off until the streaming instability was fully nonlinear. We then chose a point during the saturated state of the instability to switch on selfgravity, after which we characterized the formation of planetesimals.
Most of our simulations are run in a domain size matching our previous work (Simon et al., 2016, 2017), = , at a resolution of with the number of particles . However, for , we find that a large, asymmetric vertical outflow develops that carries away particle mass. We have determined that this effect, which only appears for (likely due to more vigorous turbulence in this case), does not significantly change the shape of the initial planetesimal mass functions. However, because of the decreased particle mass, fewer planetesimals are formed. Given these considerations, our simulation with is run in a taller domain with = , at a resolution of , with the the same particle resolution as in the smaller domain runs; . This taller domain does not experience any asymmetric outflow or particle mass loss.
All of our simulations are listed in Table 1. The naming convention for the simulation is P#–SG#, where P correspond to the value, and SG is the time at which selfgravity is switched on. As mentioned above, for each simulation in which the selfgravity is not on from the beginning of the run, a simulation without selfgravity on is run with otherwise identical parameters until the streaming instability reaches a statistically saturated state. We do not include these “base” simulations in the table. Finally, the tall box simulation described above follows the same naming convention, but with “–Lz0.4H” appended.
2.3. Diagnostics
In what follows, we employ several diagnostics to fully analyze our simulations. In this subsection, we first describe our new numerical algorithm for detecting and characterizing planetesimal properties. We then present two statistical methods (based on maximum likelihood estimation) that we use to characterize the mass function of the produced planetesimals. Finally, we wrap up with a brief description of our mass scaling and the fiducial times at which we analyze our simulations.
2.3.1 Plan
In order to quantify the planetesimal properties, we use a new clumpfinding algorithm, PLanetesimal ANalyzer (Plan). Briefly, Plan is based on the halo finder HOP (Eisenstein & Hut, 1998) and is designed to find selfbound clumps by analyzing 3D particle output of Athena. With its hybrid OpenMP/MPI parallelization scheme written in C++, Plan is efficient and scalable to handle billions of particles as well as multiple snapshots simultaneously.
Following the standard HOP method, Plan first computes a density for each particle with a userdefined smooth kernel using a memoryefficient BarnesHut tree (Barnes & Hut, 1986). Particles with densities higher than a customized threshold are selected and chained up towards densest neighbors to create groups. Plan then merges/decomposes those groups by examining boundaries, saddle points, and total kinematic and gravitational energies to construct a list of bound clumps. An unbinding procedure is then performed to remove contamination particles (e.g., those that are passing by clumps and are not bound). Plan then collects possibly unidentified member particles within each clump’s Hill sphere. Currently, Plan does not further identify substructures within clumps because (i) bound planetesimals in our simulations are highlyconcentrated without subhalorich hierarchical structures as seen in cosmological simulations, (ii) most of their masses collapse into regions comparable/smaller than one hydrodynamic grid cell, invalidating further identification. A more detailed description of Plan will be given in a forthcoming publication, Li et al., in prep.
In practice, Plan can calculate the mass of planetesimals much smaller than our previous clumpfinding method (e.g. Simon et al., 2017). However, we leave the inclusion and analysis of these small planetesimals to future work (Li et al., in prep). Here we instead choose a minimum planetesimal mass of 1% of the maximum planetesimal mass in all of our simulations. In this way, we are sampling the same relative range of masses for each , consistent with our previous papers, which by necessity only examined a maximum planetesimal mass range of approximately two decades.
2.3.2 Statistical Methods for Quantifying the Initial Mass Function
After using Plan to output the mass of bound clumps, we use a simple finite difference to determine an unsmoothed estimate of the differential mass function at masses corresponding to each planetesimal,
(9) 
(10) 
where is an arbitrary normalization constant. As done in our previous two papers (Simon et al. 2016 and Simon et al. 2017), we use the maximum likelihood estimator (MLE) of Clauset et al. (2009) to determine via
(11) 
where is the minimum value of the planetesimal mass in our data, and is the total number of planetesimals. The error in is
(12) 
We will also calculate the cumulative distribution function by integrating equation (10). As we find below, the cumulative distributions are not particularly well fit by a single power law, but instead an exponentially truncated power law (similar to that used in Johansen et al. 2015 and Schäfer et al. 2017) of the form,
(13) 
where is another normalization constant, is the same power law index as above, but calculated via fitting to the truncated power law model, and is the characteristic mass for the exponential truncation.
To calculate , and , we use another MLE approach, following the work by Meerschaert et al. (2012). In particular, from their equations 2.13–2.14, and can be found by solving the following set of nonlinear equations,
(14) 
(15) 
where is the minimum mass planetesimal. The normalization factor is simply set (following equation 2.15 of Meerschaert et al. 2012) by
(16) 
2.3.3 Mass Scaling and Fiducial Times
In this work, we normalize all planetesimal masses to the dimensional mass for a selfgravitating particle disk. From balancing tidal and selfgravitational forces, one gets, from the standard Toomre dispersion relation, a critical unstable wavelength of
(17) 
For lengthscales smaller than , gravity will overpower tidal shear. Thus, we can define a “gravitational mass” as patch of the particle disk with surface density and diameter ;
(18) 
were is a dimensional reference mass.
Finally, in what follows we will analyze planetesimal properties at specific fiducial times. These fiducial times were chosen “by eye” at a time in which bound clumps had recently formed but were not too connected to the largescale structure from which they formed. We then verified that the power law slope did not change appreciably at different snaphsots (see below). For reference, the fiducial times for the SG0 simulations are , and for runs P0SG0, P0.0375SG0, P0.05SG0, P0.0625SG0, P0.075SG0, P0.0875SG0, and P0.1SG0Lz0.4H, respectively. For the late gravity runs, these times are , , and for P0.0375SG72, P0.05SG80, and PG0.075SG66, respectively.
3. Results
As we have done in our previous work (Simon et al., 2016, 2017), our goal is to determine how planetesimal properties vary with the physical parameters relevant to the streaming instability. Here, we focus on the pressure gradient parameter ; we will first quantify the effects of a varying on planetesimal formation and properties, followed by an examination, within the context of varying , of whether or not the time at which selfgravity is activated has any influence on the final properties of planetesimals. We then discuss the special case of , for which the streaming instability is absent. We finish the section with an independent verification of our results using a KolmogorovSmirnov test.
3.1. Effect of the Pressure Gradient in the SG0 Runs
Figure 2 shows the particle surface density for four of our simulations: P0.0SG0, P0.0375SG0, P0.075SG0, and P0.1SG0Lz0.4H. In each case, we show the particle mass surface density in a series of snapshots, with time increasing from left to right. The rightmost panel is the surface density at the fiducial times. As the figure shows, there is a dependence of the precollapse structure on . First, for , the precollapse structure is largely nonaxisymmetric, with small filaments that appear “weblike” being azimuthally stretched (presumably by shear). As we describe further in Section 3.3 below, this nonaxisymmetric structure is likely related to the effect of particle selfgravity. For larger , the precollapse structures are mostly axisymmetric, and the width and separation of the structures increases with , a behavior broadly consistent with previous analyses (Li et al., 2018; Sekiya & Onishi, 2018).
Using Plan to detect and analyze planetesimals in all of our simulations, we now quantify the properties of these planetesimals.
3.1.1 Planetesimal Formation Efficiency
From Fig. 2 it is clear that the number of planetesimals formed decreases with increasing . This result can be further quantified by examining the fraction of solid mass converted to planetesimals (the planetesimal formation efficiency), as shown in Fig. 3 and displayed in Table 1. There is a steep dependence of the planetesimal formation efficiency on , and for , 66% of the mass has been converted into planetesimals. At the other extreme, for , only of solid mass is converted to planetesimals.^{3}^{3}3As described below, since we remove mass below 1% of the maximum mass, these are lower limits to the total number of planetesimals. However, since the cumulative distribution begins to flatten in most cases, it is unlikely that the total number of planetesimals is significantly larger than this. Furthermore, we have verified that this mass cutoff does not strongly change the mass fraction converted to planetesimals.
This strong dependence on planetesimal formation efficiency may have profound implications for the initial population of planetesimals in forming planetary systems, since depends on radius, albeit weakly. We should caution, however, that for , there remain significant axisymmetric structures after planetesimals have already formed. It is entirely possible that planetesimal formation will continue with more planetesimals ultimately spawning from these structures, pushing the planetesimal formation efficiency higher for large . Running these simulations further to test this hypothesis is very computationally expensive, and as such, we leave addressing this question to future work.
3.1.2 Initial Mass Function: Single Power Law
A key diagnostic in planetesimal formation studies is the distribution of masses, as it allows for a more direct comparison with Solar System planetesimal populations and can thus constrain planetesimal formation mechansisms. In Fig. 4, we show the differential mass function for three of our simulations. For each run in the figure, we over plot power law functions with the bestfit slope, from Equation 11. While the number of planetesimals differs drastically, as already discussed, the power law index remains approximately the same.
Figure 5 further solidifies this result and displays for the entire range of values explored here at the fiducial times and at 5 and 10 before the fiducial times (which corresponds to the third and second panels in Fig. 2, respectively). The values at the fiducial times are also given in Table 1. With the exception of the simulation, – to within . This is in agreement with our previous results (Simon et al., 2016, 2017); clearly for a wide range of numerical and physical parameters.
Another common diagnostic (Johansen et al., 2015; Simon et al., 2016, 2017; Schäfer et al., 2017), which removes most of the scatter inherent in the differential mass function, is the cumulative mass function. Figure 6 shows this cumulative mass function for all SG0 simulations. Clearly, a single value power law does not fit the cumulative data over the entire range of masses. Indeed, it appears as though itself depends on the planetesimal mass. As we will show momentarily, these data can be better fit by an exponentially truncated power law. However, we should first note that there may be different behavior below the mass cutoff described in Section 2.3.3, as suggested by the “flattening” of the cumulative mass function towards the smallest masses. The issue of how the power law index may change at very small planetesimal masses is the subject of a forthcoming paper by our group (Li et al., in prep). As such, we will not pursue this issue here.
3.1.3 Initial Mass Function: Truncated Powerlaw
As is apparent from Fig. 6, a single slope power law is too simplistic to accurately characterize the mass function over the entire range of masses probed. Here, we use equations 14–16 to fit the mass function assuming an exponentially truncated power law in the cumulative distribution (i.e., equation 13). The best fit curves are overplotted in Fig. 6. There is good agreement between this truncated power law model and the cumulative data for , though we note that there remains moderate discrepancy at the high mass end of the distributions (), which could result from fewer planetesimals contributing to the fit in that region.
The exponentially truncated power law is fit by two parameters. The first, , is the index of the power law component, and is thus the differential power law slope assuming the truncated power law model. Fig. 7 shows (calculated at the fiducial times) versus for all simulations, and these values are also displayed in Table 1; to within for the SG0 simulations with .
This shallower slope, coupled with an exponential cutoff at a characteristic mass (which we discuss below), clearly produces a better fit to the data than a single power law. The addition of the second parameter, which accounts for the steepening of the mass function towards higher masses, allows for a smaller power law index. Furthermore, while the single power law fit returned a value , which was largely determined by small planetesimals, the larger planetesimals inevitably play a role in making the single power law fit steeper; this is no longer an issue as this steepening is fitted as the exponential tail in the truncated model.
3.1.4 Characteristic Planetesimal Masses
The second parameter used in fitting the mass function to a truncated power law is , the mass above which the mass function is substantially steepened. In fitting the data, is the only mass related parameter, and as such, it represents a characteristic mass for the distribution. Furthermore, since the mass function is exponential for , can be treated as a proxy for the maximum planetesimal mass.
Fig. 8 shows this characteristic mass (with associated error bars) as a function of . The masses (for ) lie within the range – to within and are thus statistically consistent with having no dependence on . However, the data is also consistent with, and indeed marginally favors, a relatively shallow increase of with (the dashed line in the Figure shows a linear fit to the data). Similar conclusions follow from an analysis that does not depend on fitting a specific function to the data. Table 1 shows the values of a quantity we label , which is defined as the mass where,
(19) 
Here, is the total planetesimal mass greater than normalized to the total mass in planetesimals. We find that is on the order (ranging from 0.2–1) in all simulations, and, as with , could be consistent with either a constant or relatively weakly increasing function of .
As we have noted, there is no welldefined linear prediction for the unstable modes of the physical system that we are simulating (which is not in initial equilibrium, and which has selfgravity on from the start). However, to the extent that the system behaves similarly to the model linear problem (unstratified, and without selfgravity), we might expect the linear scales to be proportional to (Youdin & Goodman, 2005). This would result in a scaling of the characteristic planetesimal mass , as suggested by Taki et al. (2016). Such a cubic scaling, normalized to the result, is shown in the Figure. It is inconsistent with the available data.
Another way to see this result is to consider the ratio of the gravitational wavelength (Equation 17) to the characteristic lengthscale of the streaming instability ,
(20) 
If for all simulations, then the largest unstable mass (which is on the order of ) will scale with (this is the essence of the Taki et al. 2016 claim). However, the values of range from 0.42 at to 0.16 for ; is clearly not proportional to , and the size scale of gravitationally bound clumps is not actually set by the value of .
In both of these characteristic masses, there is an uncertainty associated with continued planetesimal formation and continued accretion of solids onto and tidal stripping of already formed planetesimals. Longer duration simulations, possibly at higher resolution (to avoid unphysical merger events or accretion) are needed to measure any dependence of on more robustly. That being said, the scatter in the data should, to some extent, include this uncertainty, and even within the given scatter, a dependence is ruled out. We conclude from this that the linear physics of the model streaming instability problem, specified by , and , does not simply translate into a prediction for the final sizes and masses of planetesimals formed in stratified selfgravitating simulations.
3.2. Late Onset Gravity Simulations
Simulations of planetesimal formation via the streaming instability start from simple initial conditions, and do not model processes such as particle growth and global radial drift that would bring the physical system to the point of instability. In particular, two types of simulation have been used in prior studies.

Simulations in which particle selfgravity is turned on from . The streaming instability (in its vertically stratified version) and SGI may both play a role in the formation of particle clumps and subsequent collapse. This choice is largely physically correct. However, by allowing collapse to occur rapidly—potentially before nonlinear streaminginduced structures have fully formed—the results could be affected by the largely arbitrary initial conditions.

Simulations in which the streaming instability is allowed to fully develop in the absence of selfgravity, which is only turned on at a later time. This choice allows the structure of twophase turbulence in disks to fully develop prior to collapse (as likely occurs in reality), at the expense of excluding some of the physical effects of secular selfgravity.
Since neither approach is fully realistic, comparing the two provides insight into whether the limitations of current simulations affect physical quantities of interest such as the mass function. In earlier lower resolution simulations we found that the slope of the mass function (fit as a single power law) was not significantly altered by this numerical choice (Simon et al., 2016). We revisit the issue here at high resolution in the context of models that vary .
Our fiducial runs for this work all initiate selfgravity at . We supplement these with three additional simulations in which we first let the streaming instability develop into its fully nonlinear state without any particle selfgravity enabled. We then activated selfgravity at a time during the nonlinear state for each simulation. In particular, for runs P0.0375SG72, P0.05SG80, and P0.075SG66, we switched on selfgravity at , and (see Table 1), when = 530, 540, and 193, respectively. Given the more than a factor of 10 enhancement of particle mass density over the initial values (), the system should be highly nonlinear at these times. Furthermore, the streaming instability growth timescale with our chosen parameters is very roughly (see Section 2.2). Thus, each of these simulations have evolved for many growth timescales before selfgravity is initiated, and we expect the system to be fully nonlinear when selfgravity is activated.
Figure 9 shows the cumulative mass functions at the fiducial times for two of these late gravity simulations compared with their SG0 counterparts. In both cases, the mass functions have similar slopes, though there is a difference in the mass values themselves, amounting to a factor of difference in P0.075SG0 compared to P0.075SG66. The corresponding parameters of the single power law and truncated power law fits are plotted in Fig. 5 and Fig. 7. As shown in these figures (and by examining the and values in Table 1), the power law slopes are consistent between the SG0 and late gravity runs to within for most of the runs and to within for . Moreover, when comparing the characteristic mass between the SG0 and late gravity simulations, is consistent with to within (see Table 1). We should note however, that in every case, in the late gravity simulations is less than in the corresponding SG0 run, a feature which is seen in Fig. 9 as well. Nevertheless, given the error bars, the discrepancy between the model and data at masses greater than , and the potential role of continued accretion (and possibly tidal interactions) as described above, there is, at best, a weak dependence of on the selfgravity start time.
Our conclusion is that the interplay between the initial conditions and the time at which selfgravity is initiated do not significantly change the properties of the mass function. This is consistent with our our previous result (Simon et al., 2016), derived from substantially lower resolution runs, that is largely independent of initial conditions and whether or not selfgravity is activated from the beginning of the run. Physically, it may suggest that the clumping of solids due to the nonselfgravitating streaming instability alone, and that caused by the combined action of streaming and secular selfgravity, are similar. More indepth work, however, will be needed to determine precisely how the nonlinear turbulence driven by the streaming instability in the presence of selfgravity leads to the planetesimal masses and mass functions that we see.
3.3. The Special Case of Zero Pressure Gradient
As mentioned previously, we have carried out one simulation with . Since, by definition, there is no gasparticle drift for , the streaming instability is necessarily absent. The effects of SGI are however present. A direct comparison between the planetesimals that form in the case and those that form for thus tests to what extent the streaming instability is present for the simulations versus SGI.
P0.0375SG  P0.05SG0  P0.0625SG0  P0.075SG0  P0.875SG0  P0.1SG0Lz0.4H  

P0SG0  0.007  
P0.0375SG0  x  0.29  0.03  0.26  0.52  
P0.05SG0  x  x  0.12  0.59  0.84  0.002 
P0.0625SG0  x  x  x  0.65  0.81  0.05 
P0.075SG0  x  x  x  x  0.69  0.05 
P0.0875SG0  x  x  x  x  x  0.13 
The top row of Fig. 2 shows the particle mass surface density for for a series of time snapshots. There is considerable small scale structure that develops before planetesimals form, with many narrow filaments that appear elongated in azimuth (most likely due to shear). This mechanism generates a large number of planetesimals that are uniformly distributed throughout the computational domain. The precollapse structure seen in this simulation is quite similar to that in P0.0375SG0 (as well as in P0.05SG0 and P0.0625SG0; not shown in the figure).
On the other hand, for the lategravity runs, P0.0375SG72 and P0.05SG80, the particle surface density at early times (before selfgravity is activated) resembles the larger simulations in that there are largely axisymmetric structures (though of significantly smaller radial scale compared to the larger runs) induced that do not have the same “weblike” shape demonstrated by the SG0 runs. Evidently, selfgravity plays a role in determining the precollapse structure of the particles for .
This result suggests that the mass function in the SG0 runs with may be influenced by the same physics that determines the mass function. Indeed it may very well play some subdominant role in P0.0375SG0; the value of is slightly larger than the other runs and the late gravity simulation. However, whether P0.0375SG0 lies in an intermediate physical regime between SGI and streaming instability remains difficult to quantify, particularly since is in agreement with the larger simulations.
For , we find and values more than – lower than the case, in which and ^{4}^{4}4Figure 5 shows that at before the fiducial time, for P0SG0. However, upon examination of the surface density corresponding to this time in Fig. 2, it is not clear that planetesimal formation is complete, as many of the plantesimals are still wellconnected to the smallscale filaments. Given this observation and the agreement between at before the fiducial time with that at the fiducial time, we believe that is too early to accurately quantify planetesimal properties in this particular run.. is also larger than the mass () with which the simulations are statistically consistent. However, apart from , the is also statistically in agreement with the larger simulations to within .
These differences (particularly in the power law slopes) suggest that a disk that is unstable to the streaming instability will produce a mass function that is significantly different from a disk where gravitational instability is triggered by other physical processes. As we suggested in Simon et al. (2017), this is consistent with turbulence playing a key role in determining the mass function. If the streaming instability produces turbulence with a particular power spectrum in turbulent fluctuations that is absent or different from SGI, then one would expect different mass functions to arise.
Despite the steeper mass function in the run, it is noteworthy that . That is, both the streaming instability simulations and the run where streaming is absent produce a topheavy mass function where most of the mass is in the largest planetesimals. Many of the qualitative implications discussed in the context of planetesimals forming via the streaming instability carry over directly to the case with zero pressure gradient.
3.4. Robustness of Derived Fits
For the results presented above we have fit the different models to planetesimal data that spans two orders of magnitude in mass. This choice was made to ensure uniformity across runs, and because we are confident that the least massive planetesimals remain wellresolved. We have briefly explored how the results change if we lower the low mass cutoff to less than the fiducial 1%. Doing so we found modestly lower values, with dropping the most but remaining distinct from the runs. However, we also found that some of the fits to a truncated power law are not nearly as good (as judged by eye) when lowering the cutoff, making our statements about lower more tentative.
We do not have any strong theoretical reason to expect that the true mass function is a truncated exponential—we have adopted this function because it appears to be a reasonable approximation that can be described by only two parameters. More complex functions would yield better fits, and might justify the introduction of extra parameters. For example, exponentially truncated power laws with an additional exponent parameter () have been fit to planetesimal data by Johansen et al. (2015) and Schäfer et al. (2017). Schäfer et al. (2017) find that an exponentially truncated power law with – fit their data reasonably well, though they were unable to strongly constrain the low mass end of the mass function using their data. Further work is thus needed to determine what is the best description of the entire mass function forming from the streaming instability.
Given the above uncertainties the precise extrapolation of the derived mass function to lower masses (determined by the values of and ) remains unclear. However, we consider that the consistency of the derived parameters across different runs (i.e. for the single power law and for the truncated model) to be more robust and physically significant. Both the fitting exercise described above, and the modelindependent tests presented in §3.5, lead to two conclusions. First, for the “streaming unstable” runs with , very similar mass functions are derived for all of the relevant streaming parameters studied to date (Simon et al., 2016, 2017). We have interpreted this previously as pointing to the primary role played by nonlinear turbulence in determining the initial mass function of planetesimals. Second, the run that is not unstable to streaming is significantly different.
3.5. ModelIndependent Verification of Results
As discussed above, the truncated power law fit (which is already significantly better than the single power law fit) to the mass functions shows the largest discrepancy for ; the fit nearly always leads to fewer planetesimal numbers at the highest masses. To determine how such discrepancies might affect our results, we present here a modelindependent test. Namely, we determine the probability that any two mass functions are derived from the same distribution via a KolmogorovSmirnov (KS) test, in which we normalize the masses by the maximum mass for each simulation and then apply the KS test in logarithmic space.
Table 2 displays the probability value from the KS test for every combination of pairs between the SG0 simulations. While a large range of probabilities exist for the runs with , most of the probabilities are . In contrast, the KS probability from comparing P0SG0 to the other simulations are mostly (with the exception of 0.007 for P0.0875SG0 compared to P0SG0). This large discrepancy supports our claim that the mass function for is different than for .
It is also worth noting that the run shows both moderate and low probabilities, ranging from to 0.13, when compared against other runs. This is consistent with our general finding that both and are lower for compared to the other simulations. However, given the large error bars on the fit parameters and the range in probability values from the KS test, we report that there is, at best, tentative evidence for producing a different mass function than the other runs with . Regardless of the robustness of such evidence, the deviation of from the other runs appears to be modest.
4. Compatibility with observations
The most direct constraints on the initial planetesimal mass function come from observations of populations of small bodies—asteroids and TransNeptunian Objects (TNOs)—in the Solar System. Suppose that we fit (allowing for a varying power law index) the differential number distributions of these objects as powerlaws in mass , diameter , and absolute magnitude (taking the band as a relevant band for TNOs),
(21)  
(22)  
(23) 
Then, if we assume constant albedo and constant density, there are simple relations linking the slope of the mass distribution to the slope of the size distribution and to the observationally accessible distribution of absolute magnitudes ,
(24)  
(25) 
Fitting the streaming instability results to a single powerlaw yields , which implies that and that (Johansen et al., 2012; Simon et al., 2016; Schäfer et al., 2017). The low mass end of the truncated powerlaw fit yields , giving and .
Assessing the consistency of these predictions with observations is not an easy task. The smallest bodies, especially in the asteroid belt, are collisionally evolved (Dohnanyi, 1969), while the largest bodies may have grown via accretion of planetesimals or pebbles (Ormel & Klahr, 2010; Johansen et al., 2015). Both of these modifiers depend on the initial mass in small bodies in the belts, which is a large but uncertain multiple of the current mass.
4.1. Asteroid Belt size distribution
The degree to which the current asteroid belt size distribution is consistent with streaming predictions has been studied by Morbidelli et al. (2009) and by Johansen et al. (2015). Models suggest that asteroids with diameters are primordial—in the sense of not having suffered latetime collisional evolution (Bottke et al., 2005). The slope of the size distribution of these bodies is , which is significantly steeper than the streaming prediction. Smaller bodies are typically collisionally evolved. Modeling by Tsirvoulis et al. (2018) backs out a primordial size distribution for that is characterized by a slope , which is somewhat flatter than the single powerlaw fit to streaming simulations but somewhat steeper than the low mass limit of the truncated powerlaw derived in this work.
4.2. Kuiper Belt size distribution
Turning to the Kuiper Belt, Figure 10 shows estimates of derived from selected recent observational determinations of TNO size distributions. The blue boxes show estimates for a dynamically cold populations of TNOs from Fraser et al. (2014). This population is arguably of the greatest interest theoretically, as modeling suggests that its properties are consistent with in situ formation with a low initial mass (Gomes et al., 2018). For their sample of cold TNOs (Fraser et al., 2014) derive a faintend slope (), and a brightend slope that has large uncertainties but which is clearly much steeper (). The inferred break occurs at an absolute magnitude , which corresponds to for an albedo of 0.04. The figure also shows, as the red boxes, estimates for for scattered TNOs derived by Lawler et al. (2018) using data from the Outer Solar System Origins Survey. (Note that the true uncertainty in for this sample is not clear; we have assigned arbitrary errors of in Fig. 10.) They derive a faintend slope of that is similar to that found for the cold population, and a brightend slope of that is shallower than for dynamically cold TNOs. The break between these populations occurs at , corresponding to a diameter of about 120 km.
From Figure 10 it is clear that the size distribution of TNOs in the range is an excellent match to the bestfitting single powerlaw derived from our (and other groups) simulation results. The agreement with the low mass limiting slope of the truncated powerlaw mass function (equation 13) is markedly worse, and the theoretical and observational results are significantly discrepant according to the purely statistical error bars that we quote. It is much harder, of course, to estimate the magnitude of possible systematic errors, arising for example from the fact that we fit only a limited range of planetesimal masses. The largest and brightest objects, conversely, display a much steeper distribution that appears to differ between dynamically distinct populations. The consistency of the steeper slopes at larger sizes with streaming predictions depends sensitively on the precise form of the highmass cutoff. Using the simple exponential form that we have adopted we would infer a characteristic diameter (at the cutoff mass ) of . This, however, leads to very few large KBOs. A distribution more consistent with observations would require either a softer cutoff function (such as that adopted by Schäfer et al., 2017) or a distribution that is the superposition of several mass functions with differing values of .
5. Conclusions
We have studied the mass distribution of planetesimals that form from gravitational collapse in aerodynamically coupled mixtures of solids and gas within protoplanetary disks. Gravitational collapse is part of several scenarios for planetesimal formation (including the classical picture of fragmentation of a thin dust layer; Safronov, 1972; Goldreich & Ward, 1973), but appears easiest to realize as a secondary phase of particle clumping resulting from the streaming instability (Youdin & Goodman, 2005). Accordingly, we have focused on simulation domains and physical parameters (in particular relatively high solidtogas ratios) that capture the regime where the streaming instability generates strong particle clumping (Carrera et al., 2015; Yang et al., 2017). The numerical setup of the simulations closely follows our earlier Athenabased work (Simon, 2016; Simon et al., 2017) with one important addition: the use of a new clumpfinding routine (Plan) to robustly identify bound objects in threedimensions. We also use a maximum likelihood estimator to fit an exponentially truncated powerlaw model to the simulation data, consistently extending our earlier powerlaw only results.
Our results are based on a set of high resolution simulations ( or gas zones, solid particles) that vary the strength of the background gas pressure gradient. We consider a range of between 0.0375 and 0.1, which is representative of the values in a protoplanetary disk under typical assumptions, and separately model the limit, which is of interest because it is not subject to streaming instability. Our primary results are as follows.

For , the simulations produce a mass function of planetesimals that, when fit as a single powerlaw, yields an index – to within .

The slope of the derived mass functions clearly varies with mass, making a single powerlaw a poor representation of the data. An exponentially truncated powerlaw provides a substantially better fit to the cumulative mass distribution. Using a maximum likelihood estimator we fit a function of the form . The derived powerlaw index is remarkably consistent across the runs, at .

The exponentially truncated powerlaw model has a single characteristic mass parameter . Characterizing the masses of the most massive planetesimals formed in the simulations using , we find that is of the order of the gravitational mass , which scales with disk parameters as (equation 18). There is tentative evidence that increases with , and may depend on whether selfgravity is operative from the start of the simulations, or instead turned on at a (mostly arbitrary) later time after nonlinear clustering has already developed.

As noted by Taki et al. (2016) and Youdin & Goodman (2005) the linear scale of the streaming instability is proportional to , defining a mass scale that is proportional to . Our results exclude a scaling as steep as . This suggests that the linear physics of the unstratified model problem does not simply translate into a prediction for the characteristic masses of planetesimals.

When combined with the results of prior simulations that varied the dimensionless stopping time , the particle concentration (Simon et al., 2017), and the strength of gravity relative to shear (Simon et al., 2016), our results support the hypothesis that the initial mass function resulting from the collapse of overdensities is weakly dependent (and possibly independent) on the physical parameters relevant to the streaming instability.

The special case of yields a mass function that is significantly steeper (but still topheavy) than the streaming unstable cases, irrespective of whether the data is fit by a single powerlaw or an exponentially truncated powerlaw in the cumulative distribution. We caution that the numerical challenges of modeling the limit are somewhat different from the cases, and that more work is needed to study in detail. Nonetheless, the observed differences constitute initial evidence that the streaming instability produces a distinct mass function that is not shared by other mechanisms that result in gravitational collapse.
Our work is subject to both numerical and physical limitations. On the numerical side, solving Poisson’s equation on a uniform grid halts collapse on the grid scale at densities that are much smaller than material densities. This severely restricts our ability to say anything about the possible fragmentation of collapsing clumps into binaries (see Nesvorný et al. 2010), or clusters of smaller planetesimals. Furthermore, the dependence (if any) of the mass function on the level of intrinsic disk turbulence remains to be fully explored (though see Johansen et al., 2007, 2011). We also note that, for reasons of numerical expediency, most work to date has focused on regions of parameter space where the nonselfgravitating streaming instability produces overdensities that would collapse rapidly in the presence of selfgravity (as defined by Carrera et al., 2015; Yang et al., 2017). Precisely because of the speed of gravitational collapse, we might instead expect that physical planetesimal formation occurs near the threshold values of the parameters that allow it (Armitage et al., 2016). There is no evidence that a “slow” approach to planetesimal formation—for example as particles increase their via coagulation or as a disk increases locally in a pressure trap—changes the mass function, but neither can such possibilities be excluded based on current results.
Future work will need to quantify the systematic uncertainties in the predicted mass functions that arise from the above limitations. This is needed, most obviously, in order to assess whether the streaming predictions are compatible with observations of small body populations in the Solar System. Adopting our best fitting model (the truncated power law) we predict a slope for faintend KBOs that is somewhat too flat ( versus the that is inferred observationally)^{5}^{5}5It is intriguing that in the simulation, tentatively suggesting that direct collapse played a role in Solar System planetesimal formation. However, more carefully designed numerical experiments must be carried out to solidify such a claim., and a cutoff that is substantially too sharp^{6}^{6}6The more gradual cutoff proposed by Schäfer et al. (2017) would ameliorate this problem substantially. We defer study of the precise nature of the cutoff, however, to future work.. Exactly how varies with mass is clearly sensitive to whether the observed objects formed in a single burst, or in multiple episodes with differing , and hence the significance of the latter discrepancy is dependent on unknown aspects of Solar System planetesimal formation history. However, given the robustness of the slope at small sizes, as we have found here, this slope is a stronger, and possibly ultimately decisive, test.
From these considerations, the most important implication of a streaming origin for planetesimals remains the prediction of a topheavy initial mass function (e.g. Johansen et al., 2007, 2015; Simon et al., 2016, 2017; Schäfer et al., 2017). Compared to the classical alternative of km or subkm scale planetesimals, a topheavy mass function implies a dynamically hotter population that is accreted by growing protoplanets and giant planet cores more slowly. For model builders the apparent scaling of the characteristic or maximum mass with (implying in turn a steep—nominally cubic—scaling with particle surface density) may also be important. The main challenge in developing streamingbased models for forming planetary systems is the requirement for at least modestly superSolar solidtogas ratios (Johansen et al., 2009b; Carrera et al., 2015; Yang et al., 2017). This result has also been found in all work to date, and it suggests two basic scenarios for successful planetesimal formation. One possibility is that a phase of rapid radial drift leads to a particle pileup (Youdin & Chiang, 2004), with planetesimal formation occurring at an early time primarily interior to some critical radius. Alternatively, the conditions needed to trigger gravitational collapse might be obtained only locally, in overdensities associated with zonal flows (Johansen et al., 2009a; Simon et al., 2012), vortices (Barge & Sommeria, 1995), or ice lines (Stevenson & Lunine, 1988). In some variations of this second scenario planetesimal formation would be directly linked to the mechanisms that slow radial drift, and there would be a direct connection between largescale structure in disks and where and when planetesimals form.
The case for the streaming instability playing a central role in planetesimal formation has often hinged, in part, on the perceived shortcomings of other models. Further advances in the fidelity of numerical simulation, together with work to understand how to combine the results with models of particle growth and collisional evolution, offer the exciting prospect of instead comparing streaming predictions directly against Solar System and exoplanet data.
References
 Andrews et al. (2016) Andrews, S. M., Wilner, D. J., Zhu, Z., et al. 2016, ApJ, 820, L40
 Armitage et al. (2016) Armitage, P. J., Eisner, J. A., & Simon, J. B. 2016, ApJ, 828, L2
 Bai & Stone (2010a) Bai, X.N., & Stone, J. M. 2010a, The Astrophysical Journal, 722, 1437
 Bai & Stone (2010b) —. 2010b, The Astrophysical Journal Supplement, 190, 297
 Bai & Stone (2010c) —. 2010c, The Astrophysical Journal Letters, 722, L220
 Bai & Stone (2014) Bai, X.N., & Stone, J. M. 2014, ApJ, 796, 31
 Barge & Sommeria (1995) Barge, P., & Sommeria, J. 1995, A&A, 295, L1
 Barnes & Hut (1986) Barnes, J., & Hut, P. 1986, Nature, 324, 446
 Béthune et al. (2017) Béthune, W., Lesur, G., & Ferreira, J. 2017, A&A, 600, A75
 Birnstiel et al. (2012) Birnstiel, T., Klahr, H., & Ercolano, B. 2012, Astronomy and Astrophysics, 539, A148
 Birnstiel et al. (2011) Birnstiel, T., Ormel, C. W., & Dullemond, C. P. 2011, Astronomy and Astrophysics, 525, 11
 Blum (2018) Blum, J. 2018, Space Sci. Rev., 214, 52
 Blum & Wurm (2008) Blum, J., & Wurm, G. 2008, Annual Review of Astronomy and Astrophysics, 46, 21
 Bottke et al. (2005) Bottke, W. F., Durda, D. D., Nesvorný, D., et al. 2005, Icarus, 175, 111
 Carrera et al. (2015) Carrera, D., Johansen, A., & Davies, M. B. 2015, A&A, 579, A43
 Clauset et al. (2009) Clauset, A., Shalizi, C. R., & Newman, M. E. J. 2009, SIAM Review
 Colella (1990) Colella, P. 1990, JCP, 87, 171
 Colella & Woodward (1984) Colella, P., & Woodward, P. R. 1984, JCP, 54, 174
 Cuzzi et al. (1993) Cuzzi, J. N., Dobrovolskis, A. R., & Champney, J. M. 1993, Icarus, 106, 102
 Dittrich et al. (2013) Dittrich, K., Klahr, H., & Johansen, A. 2013, ApJ, 763, 117
 Dohnanyi (1969) Dohnanyi, J. S. 1969, J. Geophys. Res., 74, 2531
 Dubrulle et al. (1995) Dubrulle, B., Morfill, G., & Sterzik, M. 1995, Icarus, 114, 237
 Eisenstein & Hut (1998) Eisenstein, D. J., & Hut, P. 1998, ApJ, 498, 137
 Fortier et al. (2013) Fortier, A., Alibert, Y., Carron, F., Benz, W., & Dittkrist, K.M. 2013, A&A, 549, A44
 Fraser et al. (2014) Fraser, W. C., Brown, M. E., Morbidelli, A., Parker, A., & Batygin, K. 2014, The Astrophysical Journal, 782, 100
 Gardiner & Stone (2005) Gardiner, T. A., & Stone, J. M. 2005, JCP, 205, 509
 Gardiner & Stone (2008) —. 2008, JCP, 227, 4123
 Goldreich & Ward (1973) Goldreich, P., & Ward, W. R. 1973, ApJ, 183, 1051
 Gomes et al. (2018) Gomes, R., Nesvorný, D., Morbidelli, A., Deienno, R., & Nogueira, E. 2018, Icarus, 306, 319
 Gundlach & Blum (2015) Gundlach, B., & Blum, J. 2015, ApJ, 798, 34
 Güttler et al. (2010) Güttler, C., Blum, J., Zsom, A., Ormel, C. W., & Dullemond, C. P. 2010, A&A, 513, A56
 Hawley et al. (1995) Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1995, ApJ, 440, 742
 Isella et al. (2016) Isella, A., Guidi, G., Testi, L., et al. 2016, Physical Review Letters, 117, 251101
 Johansen et al. (2011) Johansen, A., Klahr, H., & Henning, T. 2011, A&A, 529, A62
 Johansen et al. (2015) Johansen, A., Mac Low, M.M., Lacerda, P., & Bizzarro, M. 2015, Science Advances, 1, 1500109
 Johansen et al. (2007) Johansen, A., Oishi, J. S., Mac Low, M.M., et al. 2007, Nature, 448, 1022
 Johansen et al. (2009a) Johansen, A., Youdin, A., & Klahr, H. 2009a, The Astrophysical Journal, 697, 1269
 Johansen et al. (2009b) Johansen, A., Youdin, A., & Mac Low, M.M. 2009b, The Astrophysical Journal Letters, 704, L75
 Johansen et al. (2012) Johansen, A., Youdin, A. N., & Lithwick, Y. 2012, Astronomy and Astrophysics, 537, A125
 Johnson et al. (2008) Johnson, B. M., Guan, X., & Gammie, C. F. 2008, ApJS, 179, 553
 Koyama & Ostriker (2009) Koyama, H., & Ostriker, E. C. 2009, ApJ, 693, 1316
 Krivov & Booth (2018) Krivov, A. V., & Booth, M. 2018, MNRAS, 479, 3300
 Lawler et al. (2018) Lawler, S. M., Shankman, C., Kavelaars, J. J., et al. 2018, AJ, 155, 197
 Li et al. (2018) Li, R., Youdin, A. N., & Simon, J. B. 2018, ApJ, 862, 14
 Masset (2000) Masset, F. 2000, A&AS, 141, 165
 Meerschaert et al. (2012) Meerschaert, M. M., Roy, P., & Shao, Q. 2012, Communications in Statistics  Theory and Methods, 41, 1839
 Morbidelli et al. (2009) Morbidelli, A., Bottke, W. F., Nesvorný, D., & Levison, H. F. 2009, Icarus, 204, 558
 MoroMartín et al. (2009) MoroMartín, A., Turner, E. L., & Loeb, A. 2009, ApJ, 704, 733
 Nakagawa et al. (1986) Nakagawa, Y., Sekiya, M., & Hayashi, C. 1986, Icarus, 67, 375
 Nesvorný et al. (2010) Nesvorný, D., Youdin, A. N., & Richardson, D. C. 2010, AJ, 140, 785
 Ormel & Cuzzi (2007) Ormel, C. W., & Cuzzi, J. N. 2007, Astronomy and Astrophysics, 466, 413
 Ormel & Klahr (2010) Ormel, C. W., & Klahr, H. H. 2010, A&A, 520, A43
 ALMA Partnership et al. (2015) ALMA Partnership, Brogan, C. L., Pérez, L. M., et al. 2015, ApJ, 808, L3
 Pinilla et al. (2012) Pinilla, P., Birnstiel, T., Ricci, L., et al. 2012, Astronomy and Astrophysics, 538, 114
 Pollack et al. (1996) Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62
 Raymond et al. (2018) Raymond, S. N., Armitage, P. J., Veras, D., Quintana, E. V., & Barclay, T. 2018, MNRAS, 476, 3031
 Safronov (1972) Safronov, V. S. 1972, Evolution of the protoplanetary cloud and formation of the earth and planets.
 Schäfer et al. (2017) Schäfer, U., Yang, C.C., & Johansen, A. 2017, Astronomy and Astrophysics, 597, A69
 Sekiya & Onishi (2018) Sekiya, M., & Onishi, I. K. 2018, ApJ, 860, 140
 Simon (2016) Simon, J. B. 2016, The Astrophysical Journal Letters, 827, L37
 Simon & Armitage (2014) Simon, J. B., & Armitage, P. J. 2014, The Astrophysical Journal, 784, 15
 Simon et al. (2016) Simon, J. B., Armitage, P. J., Li, R., & Youdin, A. N. 2016, The Astrophysical Journal, 822, 55
 Simon et al. (2017) Simon, J. B., Armitage, P. J., Youdin, A. N., & Li, R. 2017, The Astrophysical Journal Letters, 847, L12
 Simon et al. (2012) Simon, J. B., Beckwith, K., & Armitage, P. J. 2012, Monthly Notices of the Royal Astronomical Society, 422, 2685
 Simon et al. (2011) Simon, J. B., Hawley, J. F., & Beckwith, K. 2011, ApJ, 730, 94
 Stevenson & Lunine (1988) Stevenson, D. J., & Lunine, J. I. 1988, Icarus, 75, 146
 Stone & Gardiner (2010) Stone, J. M., & Gardiner, T. A. 2010, ApJS, 189, 142
 Stone et al. (2008) Stone, J. M., Gardiner, T. A., Teuben, P., Hawley, J. F., & Simon, J. B. 2008, The Astrophysical Journal Supplement, 178, 137
 Takahashi & Inutsuka (2014) Takahashi, S. Z., & Inutsuka, S.i. 2014, ApJ, 794, 55
 Taki et al. (2016) Taki, T., Fujimoto, M., & Ida, S. 2016, Astronomy and Astrophysics, 591, A86
 Trujillo et al. (2001) Trujillo, C. A., Jewitt, D. C., & Luu, J. X. 2001, AJ, 122, 457
 Tsirvoulis et al. (2018) Tsirvoulis, G., Morbidelli, A., Delbo, M., & Tsiganis, K. 2018, Icarus, 304, 14
 Wada et al. (2009) Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2009, ApJ, 702, 1490
 Ward (1976) Ward, W. R. 1976, in Frontiers of Astrophysics, ed. E. H. Avrett, 1–40
 Weidenschilling (1977) Weidenschilling, S. J. 1977, MNRAS, 180, 57
 Weidenschilling (1995) Weidenschilling, S. J. 1995, Icarus, 116, 433
 Whipple (1972) Whipple, F. L. 1972, in From Plasma to Planet, ed. A. Elvius, 211
 Yang et al. (2017) Yang, C.C., Johansen, A., & Carrera, D. 2017, A&A, 606, A80
 Youdin & Johansen (2007) Youdin, A., & Johansen, A. 2007, The Astrophysical Journal, 662, 613
 Youdin (2011) Youdin, A. N. 2011, ApJ, 731, 99
 Youdin & Chiang (2004) Youdin, A. N., & Chiang, E. I. 2004, ApJ, 601, 1109
 Youdin & Goodman (2005) Youdin, A. N., & Goodman, J. 2005, The Astrophysical Journal, 620, 459
 Zsom et al. (2010) Zsom, A., Ormel, C. W., Güttler, C., Blum, J., & Dullemond, C. P. 2010, Astronomy and Astrophysics, 513, A57