The Mass and Size Distribution of Planetesimals Formed by the Streaming Instability. II. The Effect of the Radial Gas Pressure Gradient

The Mass and Size Distribution of Planetesimals Formed by the Streaming Instability. II. The Effect of the Radial Gas Pressure Gradient

Charles P. Abod11affiliation: Department of Astrophysical and Planetary Sciences, University of Colorado, Boulder, CO 80309-0391 22affiliation: JILA, University of Colorado and NIST, 440 UCB, Boulder, CO 80309-0440 , Jacob B. Simon22affiliation: JILA, University of Colorado and NIST, 440 UCB, Boulder, CO 80309-0440 33affiliation: Department of Space Studies, Southwest Research Institute, Boulder, CO 80302 , Rixin Li44affiliation: Department of Astronomy and Steward Observatory, University of Arizona, 933 North Cherry Avenue, Tucson, AZ 85721 , Philip J. Armitage55affiliation: Center for Computational Astrophysics, Flatiron Institute, NY 10010 66affiliation: Department of Physics and Astronomy, Stony Brook University, Stony Brook, NY 11794 22affiliation: JILA, University of Colorado and NIST, 440 UCB, Boulder, CO 80309-0440 ,
Andrew N. Youdin44affiliation: Department of Astronomy and Steward Observatory, University of Arizona, 933 North Cherry Avenue, Tucson, AZ 85721 , Katherine A. Kretke33affiliation: Department of Space Studies, Southwest Research Institute, Boulder, CO 80302

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 self-gravity to study how the initial mass function of planetesimals depends on the radial pressure gradient. Fitting our results to a power-law, , 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 power-law 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 streaming-derived 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 top-heavy mass function but with a significantly different shape. We discuss the consistency of the theoretically predicted mass function with observations of small-body 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 disks

1. 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 m-mm scale particles (Blum, 2018), no viable path is known that leads to planetesimal formation via two-body 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 mm-cm 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 self-gravity 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 (Moro-Martín et al., 2009; Raymond et al., 2018).

Simple considerations suggest that the streaming-derived 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 solid-to-gas 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 self-gravity (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 over-densities 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 high-resolution 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 streaming-initiated 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 self-gravity 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 aaPressure gradient parameter, defined by Equation 7 bbTime (in units of ) at which particle self-gravity is initiated ccBest fit power law slope to single power law model ddBest fit power law slope to truncated power law model eeBest fit truncation mass in truncated power law model ffThe mass at which , where is the total planetesimal mass greater than normalized to the total mass in planetesimals. Mass Fraction
in Planetesimals
P0-SG0 0 0 1.75 0.03 1.52 0.05 0.78 0.21 0.57 0.81
P0.0375-SG0 0.0375 0 1.64 0.02 1.28 0.03 0.26 0.02 0.31 0.66
P0.05-SG0 0.05 0 1.60 0.03 1.29 0.03 0.35 0.05 0.42 0.55
P0.0625-SG0 0.0625 0 1.59 0.03 1.31 0.04 0.32 0.04 0.35 0.29
P0.075-SG0 0.075 0 1.58 0.05 1.31 0.06 0.56 0.17 0.61 0.22
P0.0875-SG0 0.0875 0 1.58 0.07 1.33 0.09 0.82 0.34 1.0 0.12
P0.1-SG0-Lz0.4H 0.1 0 1.49 0.07 1.21 0.07 0.44 0.12 0.52 0.073
P0.0375-SG72 0.0375 72 1.56 0.02 1.19 0.03 0.16 0.02 0.22 0.53
P0.05-SG80 0.05 80 1.64 0.03 1.30 0.05 0.24 0.04 0.28 0.28
P0.075-SG66 0.075 66 1.51 0.06 1.21 0.06 0.26 0.05 0.39 0.062
Table 1Streaming Instability Simulations

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 self-gravity 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 self-gravity is implemented using a particle-mesh 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 power-law, rather than the simple power-law 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 co-rotating 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 co-rotates 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:


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 left-hand side of the momentum equation (i.e., the motion of the gas in the absence of source terms) is solved using a second-order accurate Godunov flux-conservative method, with the dimensionally unsplit corner transport upwind method of Colella (1990) and the third-order 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 non-inertial terms), including orbital advection (the background Keplerian velocity is subtracted and integrated analytically; Masset 2000; Johnson et al. 2008) and Crank-Nicholson 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 super-particle approach (i.e., each super-particle is a statistical representation of a number of smaller particles). The equation of motion for super-particle (hereafter, simply “particle” for simplicity) is given by


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 self-gravity. 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 semi-implicit 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 self-gravity is activated, there is an additional force in equation 3, , which is found by solving Poisson’s equation for particle self-gravity. 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 non-periodic 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 shearing-periodic 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; shearing-periodic 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

Figure 1.— Self-gravity parameter (solid, black curve) and dimensionless stopping time (dashed, red curve) as a function of radius (in AU) for a protoplanetary disk with surface density , disk mass 0.05, and temperature profile . is calculated for the Epstein drag regime (this may not apply at small radii, but we are only interested in relatively large distances from the star) assuming a constant particle size of 2.5 mm. The higher (red) dotted line corresponds to the value used in all of our simulations, and the lower (black) dotted line corresponds to used in our simulations. The vertical dotted line marks 10 AU, and corresponds to the radius at which and . At this radius, which is reasonable for planet formation to occur, our simulation parameters are physically realistic for typical protoplanetary disk conditions.

All of our simulations are initialized as follows. The gas is in a hydrostatic (Gaussian) vertical profile,


where is the mid-plane 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,


which characterizes the aerodynamic interaction between a single particle species and the gas, the particle concentration,


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 sub-Keplerian gas in real disks


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 self-gravity is included in the simulations: the relative importance of self-gravity and tidal shear, which we quantify as


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 two-dimensional 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.111Of 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 Kelvin-Helmholtz 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 over-densities 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.222There 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 thus-far, we have also carried out three additional simulations in which self-gravity remained off until the streaming instability was fully non-linear. We then chose a point during the saturated state of the instability to switch on self-gravity, 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 self-gravity is switched on. As mentioned above, for each simulation in which the self-gravity is not on from the beginning of the run, a simulation without self-gravity 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 clump-finding algorithm, PLanetesimal ANalyzer (Plan). Briefly, Plan is based on the halo finder HOP (Eisenstein & Hut, 1998) and is designed to find self-bound 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 user-defined smooth kernel using a memory-efficient Barnes-Hut 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 sub-structures within clumps because (i) bound planetesimals in our simulations are highly-concentrated without subhalo-rich 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 clump-finding 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.

Figure 2.— Surface density of solids in the (orbital plane) for four of our SG0 simulations. The color map is the logarithm of the surface density perturbations relative to the average surface density. Time increases from left to right, with the rightmost panel representing the fiducial time at which we analyze planetesimal properties. The top row corresponds to , second row , third row , and bottom row . As increases, the pre-collapse structure of solids moves towards larger scale, more axisymmetric structures, and the number of planetesimals decreases. In the rightmost column, white circles depict Hill spheres for a subset of planetesimals.

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,


as in Simon et al. (2016) and Simon et al. (2017). The differential mass function can be written as


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


where is the minimum value of the planetesimal mass in our data, and is the total number of planetesimals. The error in is


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,


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,


where is the minimum mass planetesimal. The normalization factor is simply set (following equation 2.15 of Meerschaert et al. 2012) by


In fitting the cumulative distributions below, we solve equations 1415 using the Newton-Rhapson method, after which we determine from equation 16. To calculate the standard errors on these parameters, we use the simple bootstrap method of sampling with replacement.

2.3.3 Mass Scaling and Fiducial Times

In this work, we normalize all planetesimal masses to the dimensional mass for a self-gravitating particle disk. From balancing tidal and self-gravitational forces, one gets, from the standard Toomre dispersion relation, a critical unstable wavelength of


For length-scales 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 ;


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 large-scale 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 P0-SG0, P0.0375-SG0, P0.05-SG0, P0.0625-SG0, P0.075-SG0, P0.0875-SG0, and P0.1-SG0-Lz0.4H, respectively. For the late gravity runs, these times are , , and for P0.0375-SG72, P0.05-SG80, and PG0.075-SG66, respectively.

Figure 3.— Fraction of total solid mass within the domain converted to planetesimals versus for all simulations with particle self-gravity always on. The values plotted were taken at the fiducial times. The planetesimal formation efficiency strongly depends on , and for , only 7.3% of the solid mass is converted to planetesimals.

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 self-gravity 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 Kolmogorov-Smirnov 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.0-SG0, P0.0375-SG0, P0.075-SG0, and P0.1-SG0-Lz0.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 pre-collapse structure on . First, for , the pre-collapse structure is largely non-axisymmetric, with small filaments that appear “web-like” being azimuthally stretched (presumably by shear). As we describe further in Section 3.3 below, this non-axisymmetric structure is likely related to the effect of particle self-gravity. For larger , the pre-collapse 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.

Figure 4.— Differential mass function for three of our simulations with (blue crosses), (black asterisks), and (red x’s). Over-plotted on each set of points is a line representing the best-fit power law index. Note that the normalization of this line was chosen to improve visibility and is not part of the fit. The masses are given in units of the gravitational mass (Equation 18). As increases, fewer planetesimals are produced, but the power law index is approximately constant.

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.333As 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 cut-off 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 best-fit slope, from Equation 11. While the number of planetesimals differs drastically, as already discussed, the power law index remains approximately the same.

Figure 5.— Best fit power law index as a function of for all of our simulations. The top panel shows at the fiducial times in the SG0 simulations. The bottom panel combines these points with before the fiducial time (blue triangles), before the fiducial time (red asterisks), and simulations with self-gravity switched on at late times (green squares). Error bars denote uncertainty from equation 12. Apart from , all data points are consistent with to within .

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.

Figure 6.— Cumulative mass function for our simulations with (yellow triangles), 0.0375 (blue crosses), 0.05 (green squares), 0.0625 (purple diamonds), 0.075 (black asterisks), 0.0875 (brown circles), and 0.1 (red x’s). The masses are given in units of the gravitational mass (Equation 18). Over plotted on each distribution is a best fit line of the same color (using equations 1416) assuming an exponentially truncated power law. The exponentially truncated power law model provides a good fit to the data for , though there remain moderate discrepancies at the highest masses.

3.1.3 Initial Mass Function: Truncated Power-law

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 1416 to fit the mass function assuming an exponentially truncated power law in the cumulative distribution (i.e., equation 13). The best fit curves are over-plotted 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.

Figure 7.— Best fit power law index , calculated from fitting to an exponentially truncated power law, as a function of for all of the SG0 runs (black circles) and the late gravity runs (green squares). Error bars denote uncertainty. Apart from , nearly all data points are consistent with to within . Run P0.0375-SG72 displays a slightly lower value, which is still consistent with but to within . While is consistently less than its counterpart fitted to a single power law () (see Fig. 5), is remarkably independent of for .

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 cut-off 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,


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 well-defined linear prediction for the unstable modes of the physical system that we are simulating (which is not in initial equilibrium, and which has self-gravity on from the start). However, to the extent that the system behaves similarly to the model linear problem (unstratified, and without self-gravity), 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 length-scale of the streaming instability ,


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 self-gravitating simulations.

Figure 8.— Best fit characteristic planetesimal mass as a function of for all SG0 simulations. Error bars denote uncertainty. Apart from , all data points are statistically consistent with either a constant (dotted line) or a linear trend with (; dashed line). The data clearly rules out a dependence (dot-dashed line), as is expected if purely linear modes determined .

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 self-gravity 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 non-linear streaming-induced 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 self-gravity, which is only turned on at a later time. This choice allows the structure of two-phase 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 self-gravity.

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 .

Figure 9.— Cumulative mass functions at the fiducial snapshots for simulations with (left) and (right). In the left (right) panel, the black crosses (asterisks) correspond to the SG0 simulations, whereas the green triangles (diamonds) correspond to the late gravity runs. In both cases, the shape of the distributions are quite similar (particularly the slope at small masses), which suggests that the mass function is largely independent of the time at which self-gravity is activated. However, there is evidence that the mass scale (; equivalently where the mass function becomes exponentially steepened) is larger in the SG0 runs compared to the late-gravity runs, though given uncertainties associated with continued accretion and planetesimal formation, as discussd in the main text, such evidence remains tentative.

Our fiducial runs for this work all initiate self-gravity at . We supplement these with three additional simulations in which we first let the streaming instability develop into its fully non-linear state without any particle self-gravity enabled. We then activated self-gravity at a time during the non-linear state for each simulation. In particular, for runs P0.0375-SG72, P0.05-SG80, and P0.075-SG66, we switched on self-gravity 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 self-gravity is initiated, and we expect the system to be fully nonlinear when self-gravity 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.075-SG0 compared to P0.075-SG66. 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 self-gravity start time.

Our conclusion is that the interplay between the initial conditions and the time at which self-gravity 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 self-gravity is activated from the beginning of the run. Physically, it may suggest that the clumping of solids due to the non-self-gravitating streaming instability alone, and that caused by the combined action of streaming and secular self-gravity, are similar. More in-depth work, however, will be needed to determine precisely how the non-linear turbulence driven by the streaming instability in the presence of self-gravity 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 gas-particle 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.0375-SG P0.05-SG0 P0.0625-SG0 P0.075-SG0 P0.875-SG0 P0.1-SG0-Lz0.4H
P0-SG0 0.007
P0.0375-SG0 x 0.29 0.03 0.26 0.52
P0.05-SG0 x x 0.12 0.59 0.84 0.002
P0.0625-SG0 x x x 0.65 0.81 0.05
P0.075-SG0 x x x x 0.69 0.05
P0.0875-SG0 x x x x x 0.13
Table 2Kolmogorov-Smirnov probabilities when comparing mass distributions between runs.

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 pre-collapse structure seen in this simulation is quite similar to that in P0.0375-SG0 (as well as in P0.05-SG0 and P0.0625-SG0; not shown in the figure).

On the other hand, for the late-gravity runs, P0.0375-SG72 and P0.05-SG80, the particle surface density at early times (before self-gravity 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 “web-like” shape demonstrated by the SG0 runs. Evidently, self-gravity plays a role in determining the pre-collapse 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 sub-dominant role in P0.0375-SG0; the value of is slightly larger than the other runs and the late gravity simulation. However, whether P0.0375-SG0 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 444Figure 5 shows that at before the fiducial time, for P0-SG0. 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 well-connected to the small-scale 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 top-heavy 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 well-resolved. We have briefly explored how the results change if we lower the low mass cut-off 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 cut-off, 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 model-independent 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. Model-Independent 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 model-independent test. Namely, we determine the probability that any two mass functions are derived from the same distribution via a Kolmogorov-Smirnov (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 P0-SG0 to the other simulations are mostly (with the exception of 0.007 for P0.0875-SG0 compared to P0-SG0). 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 Trans-Neptunian 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 power-laws in mass , diameter , and absolute magnitude (taking the band as a relevant band for TNOs),


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 ,


Fitting the streaming instability results to a single power-law 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 power-law 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 late-time 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 power-law fit to streaming simulations but somewhat steeper than the low mass limit of the truncated power-law derived in this work.

4.2. Kuiper Belt size distribution

Figure 10.— Comparison of theoretical estimates for the power-law index or of the initial planetesimal mass function with observational determinations for populations of Kuiper Belt Objects. The blue boxes show estimates for a dynamically cold population from Fraser et al. (2014), assuming that the inferred break in the slope at corresponds to . The red boxes (with arbitrarily assigned errors of ) show estimates for a population of scattered TNOs from Lawler et al. (2018), assuming the inferred break in the slope at corresponds to . The grey bands shows the best-fitting single power-law fit () and best-fitting low mass limit of the truncated power-law fit () derived from this work. The green curves show the running of for the truncated power-law assuming a break radius of , 200, or 400 km. The dashed line shows the classical prediction for a strengthless collisional cascade (Dohnanyi, 1969).

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 faint-end slope (), and a bright-end 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 faint-end slope of that is similar to that found for the cold population, and a bright-end 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 best-fitting single power-law derived from our (and other groups) simulation results. The agreement with the low mass limiting slope of the truncated power-law 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 high-mass 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 solid-to-gas 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 Athena-based work (Simon, 2016; Simon et al., 2017) with one important addition: the use of a new clump-finding routine (Plan) to robustly identify bound objects in three-dimensions. We also use a maximum likelihood estimator to fit an exponentially truncated power-law model to the simulation data, consistently extending our earlier power-law 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 power-law, yields an index to within .

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

  • The exponentially truncated power-law 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 self-gravity is operative from the start of the simulations, or instead turned on at a (mostly arbitrary) later time after non-linear 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 over-densities 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 top-heavy) than the streaming unstable cases, irrespective of whether the data is fit by a single power-law or an exponentially truncated power-law 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 non-self-gravitating streaming instability produces over-densities that would collapse rapidly in the presence of self-gravity (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 faint-end KBOs that is somewhat too flat ( versus the that is inferred observationally)555It 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 sharp666The 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 top-heavy 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 sub-km scale planetesimals, a top-heavy 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 streaming-based models for forming planetary systems is the requirement for at least modestly super-Solar solid-to-gas 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 pile-up (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 over-densities 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 large-scale 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.

Author contributions. CPA ran the simulations and analyzed the results. JBS initiated and directed the project, and contributed to the analysis. RL developed the PLAN code. All authors contributed to the interpretation and write-up of the results. We thank Daniel Wik for his insight into statistics and interpretation of our results. CPA, JBS and PJA acknowledge support from NASA under awards NNX13AI58G, NNX16AB42G and 80NSSC18K0640. JBS acknowledges support from NASA under Emerging Worlds through grant 80NSSC18K0597. RL acknowledges support from NASA under grant NNX16AP53H, and ANY acknowledges support from NSF AAS grant 1616929. The numerical simulations and analyses were performed on Stampede 2 through XSEDE grant TG-AST120062.


  • 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
  • Moro-Martín et al. (2009) Moro-Martí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
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description