The Mass and Size Distribution of Planetesimals Formed by the
Streaming Instability.
I. The role of SelfGravity
Abstract
We study the formation of planetesimals in protoplanetary disks from the gravitational collapse of solid overdensities generated via the streaming instability. To carry out these studies, we implement and test a particlemesh selfgravity module for the Athena code that enables the simulation of aerodynamically coupled systems of gas and collisionless selfgravitating solid particles. Upon employment of our algorithm to planetesimal formation simulations, we find that (when a direct comparison is possible) the Athena simulations yield predicted planetesimal properties that agree well with those found in prior work using different numerical techniques. In particular, the gravitational collapse of streaminginitiated clumps leads to an initial planetesimal mass function that is wellrepresented by a powerlaw, , with , which equates to a differential size distribution , with . We find no significant trends with resolution from a convergence study of up to grid zones and particles. Likewise, the powerlaw slope appears indifferent to changes in the relative strength of selfgravity and tidal shear, and to the time when (for reasons of numerical economy) selfgravity is turned on, though the strength of these claims is limited by small number statistics. For a typically assumed radial distribution of minimum mass solar nebula solids (assumed here to have dimensionless stopping time ), our results support the hypothesis that bodies on the scale of large asteroids or Kuiper Belt Objects could have formed as the highmass tail of a primordial planetesimal population.
Subject headings:
planets and satellites: formation — hydrodynamics — instabilities — turbulence — protoplanetary disks1. Introduction
The discovery of thousands of exoplanet systems from Kepler and other missions has confirmed the ubiquity and diversity of planetary systems in our galaxy. The dominant physical processes that lead to the observed configurations of exoplanet systems, however, remain unclear. A central question, for both high and lowmass planets, is whether what we see reflects the in situ growth of planets from a population of planetesimals, or is instead determined largely by migration at a later stage. Answering this question robustly requires determining, from protoplanetary disk initial conditions, where and when planetesimals form.
In the “bottomup” model for planet formation, planets are built through several key stages, beginning with the coagulation of small particles into larger particles and particle aggregates, the formation of planetesimals from these solids, and the growth of planetesimals from accretion into larger bodies that become terrestrial planets and the cores of gas giants. In the first step, the collisions of particles through Brownian and turbulent motions lead to their growth into larger bodies. This process is efficient at producing solids of sizes mm (e.g., Brauer et al., 2008; Birnstiel et al., 2010; Zsom et al., 2010). Recent work has shown that it may be possible for solids to grow to cmm sizes as they drift through radially varying turbulence (Dra̧żkowska et al., 2013), though interior to the snow line, such growth is difficult and solids likely remain stunted at mm sizes (Dra̧żkowska & Dullemond, 2014).
The propensity of these solids to grow in size only goes so far, however, and forming larger bodies of km size scales or larger (i.e., planetesimal scales) faces two theoretical difficulties. First, if particles grow to sufficiently large sizes ( meter sizes at 1 AU), they experience a significant headwind from the subKeplerian rotating gas, which causes them to lose angular momentum and rapidly spiral into the star (Weidenschilling, 1977). Several mechanisms have been invoked to counteract this drift. For example, pressure bumps due to zonal flows (Johansen et al., 2009a; Simon & Armitage, 2014), ice lines (Kretke & Lin 2007; but see Yang & Menou 2010 and Bitsch et al. 2014), and abrupt transitions in ionization fraction (Dzyurkevich et al., 2010; Dra̧żkowska et al., 2013) may stop or slow this radial drift and may even be required given observations showing the presence of these particles at large disk radii (e.g., Andrews et al., 2012).^{1}^{1}1Though, as shown by Birnstiel et al. (2012), the presence of small grains at large radii may be the result of very long timescales over which these grains are swept up by larger solids. However, even if these particle traps are efficient at slowing radial drift, particles still cannot easily grow beyond mmcm sizes. A combination of laboratory experiments (primarily on silicates) and modeling shows that particles within this size range do not easily stick together, and instead fragment or bounce (Brauer et al., 2008; Blum & Wurm, 2008; Birnstiel et al., 2010; Zsom et al., 2010). Whether this is also the case for icy particles is less clear (e.g., Blum & Wurm, 2008; Okuzumi et al., 2012; Kataoka et al., 2013; Wada et al., 2013; Krijt et al., 2015; Musiolik et al., 2016).
A promising route toward surpassing these difficulties involves instabilities of aerodynamically coupled systems of gas and solid particles. As the gas removes angular momentum from the particles via a headwind, the particles experience inward radial drift. However, the backreaction of the particle momentum on the gas causes particle overdensities to experience less gas drag as these overdensities act to boost the gas and reduce their own headwind. Therefore, regions of higher particle density will experience slower radial drift, leading to a pile up of solids at various radii. This mechanism is the essence of what is known as the streaming instability (Youdin & Goodman, 2005; Youdin & Johansen, 2007; Johansen & Youdin, 2007; Bai & Stone, 2010c, a).
The streaming instability has been studied analytically in the linear regime (Youdin & Goodman, 2005; Youdin & Johansen, 2007) and in the nonlinear regime via numerical simulations (e.g., Johansen et al., 2007; Johansen & Youdin, 2007; Bai & Stone, 2010a). The initial growth of the instability does not involve selfgravity, and in this limit, two and threedimensional calculations have been used to quantify the strength of clumping as a function of the physical (Johansen et al., 2009b; Bai & Stone, 2010c; Carrera et al., 2015), and numerical parameters, e.g., resolution (Bai & Stone, 2010b). In particular, Johansen et al. (2009b) and Bai & Stone (2010c) carried out numerical calculations of the streaming instability without selfgravity and found that the instability generally leads to efficient particle clumping when the dimensionless stopping time (i.e., the stopping time multiplied by the orbital frequency) is and the heightintegrated solidtogas ratio is supersolar. Furthermore, the instability occurs for a wide range in background gas pressure gradients, but for a given metallicity, it favors smaller gradients (Bai & Stone, 2010c).
Threedimensional calculations that include the mutual gravitational forces between particles have been used to study how planetesimals form (Johansen et al., 2007, 2009b, 2011, 2012, 2015). While many characteristics of planetesimal formation have yet to be fully explored, these studies have generally shown that planetesimals with masses consistent with dwarf planets and planetesimals in the main asteroid belt and Kuiper belt can form after the streaming instability generates sufficiently strong clumps that selfgravity can take over (e.g., Johansen et al., 2007, 2011, 2012, 2015).
These studies have established the streaming instability as the leading candidate mechanism for the efficient formation of planetesimals. Many critical questions, however, remain open. These include whether the streaming instability — in isolation, or in concert with largescale structures such as zonal flows (Johansen et al., 2009a) — can form planetesimals in the inner disk, where the dimensionless stopping time of mmsized particles . Furthermore, the dependence of the initial planetesimal mass and size distributions on physical properties, such as , the metallicity, and gasphase turbulence is an open issue. Indeed, characterizing these distributions is of substantial importance to explaining properties of both the main asteroid belt and the Kuiper belt planetesimal populations. Todate, simulations of the streaming instability in the presence of particle selfgravity (Johansen et al., 2015) have revealed a shallower planetesimal size distribution than what is inferred from main belt and Kuiper belt observations (Jedicke et al., 2002; Morbidelli et al., 2009; Fraser et al., 2010, 2014). However, an extensive systematic survey of the dependence of this distribution on physical parameters has yet to be carried out.
Purely numerical issues are also of interest. The relative strengths and weakness of different numerical schemes are wellcharacterized in a variety of situations, including low and highMach number turbulence, but much less is known in the more complex situation where we have multiple phases, aerodynamic coupling and selfgravity. All previous simulations of planetesimal formation that have included selfgravity have used highorder finite difference methods (implemented in the Pencil code). Here, we instead combine particle selfgravity with a higher order Godunov hydrodynamic scheme (implemented in the Athena code; Stone et al., 2008) to study planetesimal formation and again address the dependence of planetesimal mass and size distributions on numerical effects.
In this paper, the first of a series, we describe our implementation of particle selfgravity within the Athena code, and present results from a baseline set of 3D simulations. Our primary focus is on independently verifying (and in some cases extending) prior numerical simulations of streaminginitiated planetesimal formation, and hence we start with parameter choices that match those in the existing literature. We specifically address the convergence of the initial planetesimal mass function with numerical resolution, and the effect of varying the strength of selfgravity as compared to tidal shear. We also study whether the results depend on when in the simulation particle selfgravity is turned on, given that it is common (and computationally expedient) to do so at a relatively late time when the nonselfgravitating streaming instability is already strongly nonlinear. There has been some study of the impact of this approximation (Johansen et al., 2011), but a detailed investigation of the effects of pregravity conditions on the outcome of planetesimal formation has not yet been performed.
The outline of the paper is as follows. Our methodology is described in detail in Section 2, which includes a description of the numerical algorithm and the implementation of particle selfgravity, two test problems to check this implementation, and a description of our streaming instability calculations. In Section 3, we present our results from each set of calculations, and we discuss these results in Section 4. We wrap up with a summary and conclusions in Section 5.
Run  Domain Size  Resolution  Comments  

()  
SI64G0.05  300,000  0.3  0.02  0.05  400  –  
SI128G0.05  2,400,000  0.3  0.02  0.05  170  Fiducial Run  
SI256G0.05  19,200,000  0.3  0.02  0.05  150  –  
SI512G0.05  153,600,000  0.3  0.02  0.05  110  –  
SI128G0.02  2,400,000  0.3  0.02  0.02  170  –  
SI128G0.02_tm20  2,400,000  0.3  0.02  0.02  150  Restarted earlier  
SI128G0.02_tm10  2,400,000  0.3  0.02  0.02  160  Restarted earlier  
SI128G0.02_tp10  2,400,000  0.3  0.02  0.02  180  Restarted later  
SI128G0.02_tp20  2,400,000  0.3  0.02  0.02  190  Restarted later  
SI128G0.05  2,400,000  0.3  0.02  0.05  170  Fiducial Run  
SI128G0.05_tm20  2,400,000  0.3  0.02  0.05  150  Restarted earlier  
SI128G0.05_tm10  2,400,000  0.3  0.02  0.05  160  Restarted earlier  
SI128G0.05_tp10  2,400,000  0.3  0.02  0.05  180  Restarted later  
SI128G0.05_tp20  2,400,000  0.3  0.02  0.05  190  Restarted later  
SI128G0.1  2,400,000  0.3  0.02  0.1  170  –  
SI128G0.1_tm20  2,400,000  0.3  0.02  0.02  150  Restarted earlier  
SI128G0.1_tm10  2,400,000  0.3  0.02  0.1  160  Restarted earlier  
SI128G0.1_tp10  2,400,000  0.3  0.02  0.1  180  Restarted later  
SI128G0.1_tp20  2,400,000  0.3  0.02  0.1  190  Restarted later  
SI128G0.05no_clump  2,400,000  0.3  0.02  0.05  0  –  
SI128G0.05low_clump  2,400,000  0.3  0.02  0.05  40  –  
SI128G0.05med_clump  2,400,000  0.3  0.02  0.05  170  Same as Fiducial Run  
SI128G0.05high_clump  2,400,000  0.3  0.02  0.05  240  – 
2. Method
In this section, we explain our methodology. We first describe the algorithmic details of Athena in the shearing box approximation with the inclusion of aerodynamic coupling between the gas and the particles and the mutual gravitational attraction between particles solved via a Fast Fourier Transform (FFT) method. We then present tests of our particle selfgravity module, followed by a description of our set up, parameters, and diagnostics for the local, streaming instability simulations that we carry out in this paper.
2.1. Numerical Algorithm
Our simulations use Athena, a secondorder accurate Godunov fluxconservative code for solving the equations of hydrodynamics and magnetohydrodynamics. The simulations we have carried out here neglect magnetic fields, and we use the hydrodynamics module of Athena. We use the Athena configuration that includes the dimensionally unsplit corner transport upwind method of Colella (1990) coupled with the thirdorder in space piecewise parabolic method of Colella & Woodward (1984). We use the HLLC Riemann solver to calculate the numerical fluxes (Toro, 1999). A detailed description of the base Athena algorithm and the results of various test problems are given in Gardiner & Stone (2005), Gardiner & Stone (2008), and Stone et al. (2008).
The simulations employ a local shearing box approximation. The shearing box models a corotating disk patch whose size is small compared to the radial distance from the central object, . This allows the construction of a local Cartesian frame that is defined in terms of the disk’s cylindrical coordinates via , , and . The local patch corotates with an angular velocity corresponding to the orbital frequency at , the center of the box; see Hawley et al. (1995).
There are two sets of equations to solve. Comprising the first set are the continuity and momentum equation for the gas dynamics:
(1) 
(2)  
where is the mass density, is the momentum density, is the gas pressure, is the identity matrix, is the (dimensional) stopping time of the particles, and is the shear parameter, defined as lnln. We use , appropriate for a Keplerian disk. For simplicity and numerical convenience, we assume an isothermal equation of state , where is the isothermal sound speed. 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. The feedback term consists of the local particle mass density , the difference between the particle velocity and gas velocity, and the stopping time of particles 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.
The second set of equations describes the particle evolution. The equation of motion for particle is given by
(3)  
where the prime denotes a frame in which the background shear velocity has been subtracted. This is part of the orbital advection scheme (Stone & Gardiner, 2010; Masset, 2000; Johnson et al., 2008), which has been implemented for the gas dynamics as well. The term accounts for the inward radial drift of particles resulting from a gas headwind, where is the fraction of the Keplerian velocity by which the orbital velocity of particles is reduced (see Section 2.3). In real disks, this headwind results from a radial pressure gradient that causes the gas to orbit at subKeplerian speeds while the particles continue to orbit at Keplerian velocities. However, such a radial pressure gradient is inconsistent with the radial (shearing) periodic boundary conditions in the shearing box model. Following Bai & Stone (2010b), we circumvent this issue by imposing 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); in particular, we use the 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 ().
The force due to the particle selfgravity is denoted by , which is found by first solving Poisson’s equation for particle selfgravity,
(4) 
where is the gravitational potential of the particle selfgravity. The force is then calculated via,
(5) 
In solving Equation (4), we follow the methodology outlined in Section 2.3 and the Appendix of Koyama & Ostriker (2009), which consists of a discrete 3D FFT adapted to handle the shearingperiodic boundaries in and the vacuum boundaries in . Solving the Poisson equation with FFT requires periodicity in all three dimensions. In the azimuthal dimension, this is trivially satisfied. For the radial boundaries, the particle mass density (which is distributed to the gas grid points using the same interpolation method as used to calculate the momentum exchange terms; TSC) is mapped to the nearest time in which the shearingperiodic boundaries are purely periodic. From Hawley et al. (1995), the times at which the radial boundaries are purely periodic are with . For each value of , the particle density is reconstructed along via a conservative remap; the density is calculated by differencing numerical “fluxes” along . These fluxes are calculated via third order reconstruction coupled with the extremum preserving algorithm of Sekora & Colella (2009). Note that this reconstruction is the same method employed to remap the fluid quantities in the radial ghost zones as part of the shearingperiodic boundary conditions (Stone & Gardiner, 2010).
In the vertical direction, the solution to Equation (4) is found using the Green’s function of the Poisson equation for a horizontal sheet of sinusoidal source mass, as described in detail in the Appendix of Koyama & Ostriker (2009). The authors of that work developed this method to calculate the potential of selfgravitating gas, but we have trivially extended this same algorithm to use the particle mass density instead of the gas density.
With these methods in hand, the potential is calculated as follows. The particle mass density is reconstructed along via the same thirdorder conservative remap as is used for the boundary conditions, multiplied by appropriate functions of , as described above (see Equation (A11) of Koyama & Ostriker, 2009), and then Fourier transformed via a 3D FFT. The density in Fourier space is then multiplied by the appropriate coefficients and transformed back to real space via another 3D FFT (see Equation (A8) of Koyama & Ostriker, 2009). In calculating the 3D FFT’s, a domain of size is used; twice the vertical domain size is required for the use of Green’s function in solving the potential with vacuum boundaries. Finally, the gravitational potential is mapped back to the original time using the thirdorder conservative remap that was used on the particle mass density.
The resulting is a cellcentered quantity, and we calculate the forces at the cell center using a central finite difference method over three grid cells in every direction. Thus, the force (Equation 5) in the direction in cell is calculated by
(6) 
and similarly for the and directions. The forces then have to be interpolated back to the location of the particles, and we again use the TSC method. We add the selfgravity force to the particles simultaneously with the drag force.
The boundary conditions are the shearingperiodic boundaries in the radial direction, 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 grid zones (Simon et al., 2011, Li et al. 2016 in prep). This latter boundary condition is not standard in shearing box setups for the streaming instability as the small domain size makes it difficult to prevent substantial outflow in the vertical direction. However, in a forthcoming publication by the authors (Li et al. 2016 in prep), we have experimented with the vertical boundary conditions and find that coupled with a routine to renormalize the total mass in the domain to make it constant in time, the outflow boundaries work reasonably well.
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. Particle SelfGravity Tests
In this section, we carry out two tests of our numerical algorithm for particle selfgravity. The first test is the collapse of a uniform density sphere of particles. The second test consists of the selfgravitating, shearing wave of particles described in the Supplementary Material of Johansen et al. (2007).
2.2.1 Spherical Collapse
The collapse of a uniform density sphere under its own gravity has a simple analytic solution with which to compare the numerical solution. The equation of motion for a test particle at radius starting at the outer edge of the sphere is
(7) 
where is the total mass within the sphere and remains constant in time. Parameterizing the time dependence of via , we assume that , where is the initial radius of the uniform density sphere, and find that
(8) 
The analytical solutions to the radius and mass density are given by
(9) 
(10) 
where is given by Equation (8).
To test our selfgravity solver against this analytic solution, we initialize a sphere of particles with radius on a cubic domain of resolved by zones. Every zone within the radius of the sphere has one particle per grid cell initially. We turn off gas drag on the particles in order to only test the particle selfgravity module. We use TSC interpolation, and the semiimplicit particle integrator. The boundary conditions are open in all three dimensions, and consequently we use the Green’s function approach of Koyama & Ostriker (2009) to solving for the potential with vacuum boundaries; here applied to all three dimensions instead of only the vertical as described above. We arbitrarily set the value of G to .
The solution to this test problem is shown in Fig. 1, which depicts the evolution of and versus time in units of the free fall time where is the initial particle mass density. As the figure shows, the numerical solution matches the analytic solution quite well. The differences between the two curves that show up at late times is a result of the FFT solver becoming less accurate as the particles are concentrated into fewer and fewer grid cells. At such particle densities, the discretization error of the FFT solution starts to dominate, as we discuss more below. Even with these errors, the numerical solution does not deviate too much from the analytic solution. This limitation should be kept in mind when considering the formation of planetesimals. In fact, the natural softening length of our FFT solver is where is the length of a grid cell.^{2}^{2}2The grid cell is uniform along all three dimensions; .
To further test our gravity module, we have also implemented a direct summation method for calculating the mutual gravitational forces between particles. This method is exact, but scales as , where is the number of particles. For large numbers of particles the direct summation method becomes prohibitively expensive. However, for small numbers of particles we can use this method to compare the forces calculated by the FFT method with the forces from the direct summation, which will be exact to within machine precision. The comparison also tests the accuracy of our interpolation scheme since the direct sum method doesn’t require any interpolation of the gravitational forces to the particle locations. We calculate the error in the force in the th direction as
(11) 
where is the force due to selfgravity on particle calculated via the FFT method and is the same but calculated with the direct summation method. We run the particle sphere collapse problem at a resolution of (and again, one particle per grid cell within the sphere initially) and find that the errors are typically very small early on (on the order of or less). As time progresses, the errors grow; as the particles get more concentrated within a smaller number of grid cells, the truncation level error from the particle selfgravity solver will become more and more dominant. At , the average errors are approximately % or less, and near , the average errors are on order of a few percent.
2.2.2 SelfGravitating, Shearing Wave
While the spherical collapse problem tests the core of our Poisson solver, we require a problem to test the implementation of this solver in the shearing box setup. To this end, we employ the linearized selfgravitating, shearing particle wave with gas drag described in detail in Section 1.3.1 of the Supplementary Material of Johansen et al. (2007). Following that work, we linearize the continuity equation, equation of motion, and Poisson’s equation for a particle “fluid” of density , velocity , and selfgravity potential where prime represents the linear perturbation on the background. The perturbation is assumed to be of the form . Note that we assume uniformity along the vertical direction. The evolution of the density and velocity perturbations are then governed by the following equations,
(12) 
(13) 
(14) 
where the radial wavenumber is timedependent,
(15) 
We set , , , , and . We solve Equations (12)(14) using a simple finite difference algorithm to obtain our linear solution to the selfgravitating, shearing wave problem in the linear limit. For the numerical test of our algorithm, we initialize a grid of size at resolution and with one particle per cell initially; the small size/resolution in the vertical direction makes the problem effectively twodimensional while still testing the full 3D FFT algorithm. In this setup, the boundary conditions are shearingperiodic in and purely periodic in both and . Note that we turn off particle feedback to the gas for this test problem.
The time evolution of the amplitude of the perturbations are shown in Fig. 2. The numerical solution agrees very well with the linear solution up to amplitudes of . At this point, the numerical solution approaches the nonlinear regime, and agreement between the linear and numerical solutions is not expected. This test problem, which includes several essential ingredients relevant to the streaming instability calculations in this paper, (i.e., shear, selfgravity, and drag), demonstrates the validity of our selfgravity implementation in the shearing box setup.
2.3. Simulation setup
All of our streaming instability simulations use a shearing box domain of size = . Here, is the vertical scale height of the gas, which is initialized with a hydrostatic (Gaussian) vertical profile,
(16) 
where is the midplane gas density. We choose code units so that the standard gas parameters are unity; . Given the limited vertical extent of our domain, the gas density does not vary significantly.
The particles are distributed uniformly in and and with a Gaussian profile in . The scale height of this profile is . While we add random noise to the particle locations to seed the streaming instability, it is worth mentioning that the particles are not initially in an equilibrium state. There is radial motion induced by the radial drift term as described above, and the lack of a pressure gradient to counteract vertical gravity means that the particles will settle towards the disk midplane.
Four dimensionless quantities characterize the evolution of the streaming instability and the ability of dense clumps to gravitationally collapse. The streaming instability depends upon the dimensionless stopping time,
(17) 
which characterizes the aerodynamic interaction between a single particle species and the gas, the metallicity,
(18) 
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
(19) 
This radial pressure gradient produces a headwind on the particles, which has velocity a fraction of the Keplerian velocity. For all of our runs, , metallicity is moderately superSolar, , and .
With the inclusion of selfgravity, an additional parameter is needed to describe the relative importance of selfgravity and tidal shear. We define a parameter ,
(20) 
which describes the strength of selfgravity in the simulation. Physically, varying is equivalent to changing the gas density (and thus through assuming the same , the particle mass density) or the strength of tidal stretching (i.e., changing ) or both. Thus, in a real disk, will vary with radius, and in a minimum mass solar nebula model (MMSN; Hayashi, 1981), this parameter increases very gradually with radius.
In what follows, we convert code units to physical units using the mass unit and assuming a radius of 3 AU in a MMSN. For , this equates to . This physical unit conversion depends on , as we discuss further below.
Our fiducial simulation is a relatively low resolution run with = gas zones and 2,400,000 particles. For this simulation, we turn on selfgravity at a time and set .
2.4. Parameter variation
We then explore parameter space by varying one parameter in each simulation subset (demarcated in Table 1 by separate boxes). In the first subset, we vary the numerical resolution and the number of particles. We explore a lower resolution with 300,000 particles, and two higher resolutions: with 19,200,000 particles and with 153,600,000 particles. In the next subset, we vary the value of from 0.02 to 0.1; this was chosen to match the test carried out by Johansen et al. (2012) in which the strength of gravity was varied.
Finally, we determine the effect of the state at which selfgravity is initiated by examining two diagnostics in a simulation identical to the fiducial one, but with no selfgravity. The first is the maximum value of the volumeaveraged particle mass density, , and the second is a weighted version of the volumeaveraged particle mass density, defined as . In general, these two quantities behave similarly, but while the maximum density has been used more frequently in the literature to track the degree of clumping via the streaming instability, it can be largely affected by a small number of grid cells (or even one cell), and the latter quantity is more representative of the degree of clumping over the entire domain. Using both of these diagnostics in concert, we chose restart times of (the “low clumping” case), (medium clumping; this is the same as the fiducial simulation), and (high clumping). We also run one case with selfgravity turned on from the beginning state .
All of these runs are labelled by the numerical resolution employed, the strength of selfgravity, and any other modifying information. For example, run SI128G0.05low_clump has 128 zones per dimension, , and is initiated from a “low clumping” state, as described above. All simulations are shown in Table 1 with the parameters outlined and any additional, relevant comments.
2.5. Diagnostics
There are several diagnostics employed in Section 3 that we now describe. First is the maximum particle mass density, as described above, as a proxy for the degree of clumping during the streaming instability. Another often employed quantity is the particle mass surface density, which is simply the integral of over the vertical extent of the domain.
Most of our diagnostics require that we locate gravitationally bound mass clumps (after selfgravity has been included) and calculate their mass. Following the basic arguments in Johansen et al. (2011), we calculate at a given time and use a peakfinding algorithm to locate local maxima in the surface density above a certain threshold. The mass of any given clump is determined by calculating a circular region surrounding the clump and iteratively increasing the radius of this region until the mass enclosed within it equals the mass of a clump that is bound by selfgravity compared to tidal effects. Mathematically, we seek a radius such that the tidal term, , and the gravitational acceleration of a test particle at the edge of the clump’s Hill sphere, are in balance;
(21) 
when the mass enclosed within our test circle equals , we set as the Hill radius and as the clump mass.
There are two limitations to our algorithm. First, it has difficulty determining the masses of clumps when two or more clumps are very close together in the plane. This is alleviated somewhat by setting a minimum distance, equivalent to the Hill radius of a 0.1 mass (which, as we will see, is consistent with masses at the higher end of the resulting mass distributions), below which the algorithm counts two clumps as being one. Furthermore, we subtract contributions of neighboring clumps that are beyond this minimum distance but still within each clump’s Hill radius. This limitation with overlapping clumps is particularly troublesome at early times when gravitationally bound clumps have yet to fully form, yet there are still peaks in . The error associated with overlapping clumps commonly manifests itself as very sharp peaks or troughs in the time evolution of the clump masses. While this limitation causes issues on short timescales when clumps temporarily interact, it does not affect the clump mass on long timescales. As a first order approach to remove this effect, we apply a simple boxcar smoothing technique to the time evolution of clump masses in all that follows.
The second limitation arises from applying a twodimensional approach to a fully threedimensional problem. However, even though clumps exist in three dimensions, the particle layer is vertically thin enough to make the twodimensional clump finding approach a good approximation. Indeed, the algorithm appears to do a sufficient job at finding clumps that remain gravitationally bound as the simulation progresses and does not find too many false positives. Furthermore, as shown in Section 3, the diagnostics that employ this algorithm return results roughly consistent with Johansen et al. (2012), lending support to the idea that our analysis works reasonably well.
From this clump finding algorithm, we calculate both the time history of the mass of the formed clumps and the mass distribution function for these clumps. Many of the simulations form only small numbers of clumps, so we seek a simple functional form to fit to the mass distribution. Consistent with prior work, a powerlaw differential mass distribution,
(22) 
works well. Here, is a constant that for the purposes of this paper is arbitrary. We also consider the cumulative mass distribution,
(23) 
where is another arbitrary constant. We can easily translate Equation (22) to a size distribution to find where .
Under the assumption that the data is drawn from a powerlaw distribution, the power law index of the differential distribution can be determined directly from the set of measured masses using a maximum likelihood estimator (MLE; Clauset et al., 2009). From Equation (3.1) in Clauset et al. (2009), we estimate the value of the power law index by
(24) 
where is the minimum value of the planetesimal mass in our data, and is the total number of planetesimals. The error is given by their Equation (3.2),
(25) 
Using the MLE for avoids any binning step, which can introduce bias into the estimate. To visually represent the differential distribution, however, we make a local estimate of by taking to be half the distance in the mass coordinate between a given planetesimal and its nearest neighbors. A least squares linear fit in loglog space to vs gives an alternate measure of , that is very roughly consistent but generally slightly smaller.
Irrespective of the statistical methods used, two cautions are in order. First, the underlying distribution of planetesimal masses formed in the simulations cannot in reality be a simple powerlaw, because there is zero probability of forming a body with a mass greater than the total mass of solids in the domain. If the actual distribution is a truncated powerlaw, then our naive fit will overestimate the value of . Second, our lower resolution runs form rather small number of planetesimals, and any method for estimating in this situation is noisy, possibly biased, and liable to return a nongaussian distribution of slope estimates. Simple experiments suggest that only the higher resolution runs, for which there are on order 100 planetesimals, can be expected to return robust estimates of the distribution.
3. Results
In this section, we present the results from our various parameter studies. We first describe the properties of the fiducial simulation SI128G0.05 in detail, followed by quantifying the effects of changing the strength of gravity, the numerical resolution, and the degree of clumping when selfgravity is initiated. Each study is described in its own subsection and the simulations corresponding to each study are presented in Table 1.
3.1. Fiducial Run
In this section, we describe some properties of our fiducial run, SI128G0.05. Selfgravity is switched on at , after which the mutual gravitational attraction between particles causes the particle density to increase; in units of the initial midplane gas density, the maximum particle density rapidly increases to , and then slowly increases to afterwards.
A time progression of the particle surface density is shown by a series of snapshots in Fig. 3. The streaming instability produces several azimuthally extended structures at early times, which eventually form into one large clumping structure. After selfgravity turns on, some of the high density regions collapse and become gravitationally bound structures with a variety of masses. Near the end of the simulation, there are on order 10 of these bound structures, to which we refer from now on as planetesimals.
The evolution of the mass of these planetesimals is shown in Fig. 4; both the total mass and maximum mass are shown. This evolution can be compared to Fig. 13 of Johansen et al. (2012). Although, most of the simulations in Johansen et al. (2012) contain collisional microphysics, they find that the mass evolution of planetesimals only depends very weakly on this, if at all. We find that the total and maximum masses in units of the mass of Ceres are approximately the same by in our fiducial simulation and the simulation of Johansen et al. (2012). Our simulations show a faster growth rate for the masses of planetesimals. Given the difficulties associated with the clump finding analysis early on, some differences are not surprising. However, that we find approximately the same values for the maximum and total masses is encouraging.
Finally, we analyze the mass and size distribution of these planetesimals. We find that there are 13 clumps, and roughly, the number of these clumps increase towards the low mass end.
Choosing the point at which to calculate the mass distribution is a bit subtle. Our general approach is to visually inspect the particle mass density (e.g., Fig. 3) and choose a time shortly after planetesimals have formed (so as to not sample later times when these planetesimals have substantially grown in mass due to accretion of smaller particles and/or mergers) and have become separate objects (i.e., no significant overlap with the over density of mass from which they formed). We have calculated the mass distribution at three different times; , and . While the large scatter in the differential mass makes a comparison between the different snapshots and a precise quantification of the power law slope difficult, we find approximate “byeye” agreement in the values of and . To reduce the noise inherent in the differential mass distributions, we also checked the cumulative distributions at these three times and found excellent agreement. Thus, we are confident that there is not a substantial variation in the clump properties over short timescales.
The significant scatter present in the mass distribution brings into question any robust quantification of the value of . The best fit value of from this simulation is ; clearly there is a large error on . In order to improve the statistics, we have run four additional simulations with the same parameters as this fiducial run, but restarted from the “parent” nonselfgravitating simulation at different times. Specifically, we initiated selfgravity at before, before, after, and after the fiducial run. These runs are included in Table 1; the run appended with “tm20” (“tp20”) means selfgravity is initiated at the fiducial case restart time minus (plus) . We chose these relatively large time displacements to reduce the likelihood that we would be sampling planetesimals that occur from quite nearly the same initial conditions; in such a case, the formation of planetesimals would not be independent. The result is shown in Fig. 5. We fit this data with a power law and now find ; the error has been reduced by roughly a factor of three.
There is a question of whether or not the properties of the distribution of planetesimal masses depend on the time at which selfgravity is activated. We carry out such an investigation below, and as described further in that section, we find that these properties do not appear to have any strong correlation with the initial state, justifying this approach.
3.2. Effect of Resolution
Figure 6 shows a snapshot of the particle mass surface density for each resolution. As with the fiducial simulation, the time chosen in each case is sufficiently long after selfgravity has been turned on for bound planetesimals to form, but sufficiently early so as to avoid the effect of merger events and additional mass accretion onto the formed planetesimals. The times corresponding to these snapshots are also chosen in the calculation of the mass distribution function below. As expected, the size of the largest planetesimals decrease as the grid scale decreases; the minimum radius follows the grid zone size. Furthermore, as resolution is increased, the number of planetesimals increases and the smallest bound mass decreases. From resolution to , the number of formed planetesimals at this time are 2, 13, 30, and 53, respectively.
In order to show the development of these small scale structures in the highest resolution run, we have plotted the various stages of the streaming instability before and after selfgravity was turned on in Fig. 7. As both Fig. 6 and Fig. 7 show, there is significant structure present on small scales, though the large scale structure remains consistent with the lower resolution runs. In particular, there are large scale axisymmetric enhancements in the particle mass density. Once selfgravity has been activated, the smaller structure available at the higher resolution allows smaller mass planetesimals to condense out of the high densities induced by the streaming instability. Despite the increase in the number of planetesimals at small scales, larger mass planetesimals still form.
This resolution effect is also demonstrated via the differential mass distribution, shown in Fig. 8. We do not include a power law fit to either the SI64G0.05 planetesimals since there are only two objects formed at that resolution or to the SI128G0.05 planetesimals because we only include the single standard run here (in order to more clearly see the effect of the resolution on the number of planetesimals). Furthermore, the time chosen for the highest resolution run was taken to be the end of the simulation. It is not entirely clear if planetesimals have stopped forming by this point, and the very high computational expense of this simulation makes integrating further not feasible at this time.
With the exception of the lowest resolution, the highest mass planetesimals are for all resolutions. The higher masses (by only a factor of a few) in the run could be related to more rapid growth of the planetesimals that do form since their physical cross section will be larger. However, since there are only two planetesimals, very small number statistics make this notion difficult to test.
The power law slope is for SI128G0.05, for SI256G0.05, and for SI512G0.05. There is substantial scatter around the best fit power law functions, as the figure shows, though the number of formed planetesimals increases with higher resolution, improving the statistics somewhat. In the right panel of Fig. 8, we plot the power law index with the errors given by Equation (25). There does not appear to be any significant trend in the value of with resolution, and –.
Our best fit value of for the highest resolution simulation is the same as that found for the simulation of Johansen et al. (2015); they found a value of by fitting the cumulative distribution to a power law with an exponential tail. Furthermore, comparing the left panel of Fig. 8 to Fig. 3 of Johansen et al. (2015) reveals similar scatter about the best fit line for each resolution.
3.3. Strength of Gravity
We vary to match a similar exploration of this parameter in Johansen et al. (2012); specifically, . Note that with the exception of one of the simulations with , the Johansen et al. (2012) simulations included collision microphysics. However, as discussed above, these authors concluded that the collisions did not have a strong impact on the formation of planetesimals, making a comparison between our work and theirs viable.
Figure 9 shows the mass evolution of the planetesimals for the three values of . Both the total and maximum planetesimal mass increase with increasing gravity. This figure can be compared to Figure 13 of Johansen et al. (2012). In general, we find good agreement with their results; the mass values differ by a factor of less than two. Furthermore, there is no systematic direction in which the masses differ; i.e., some of the mass values in Johansen et al. (2012) are larger than what we find for the same and others are smaller. Given the differences likely present in the clump finding algorithms, we believe that the agreement is pretty strong.
In converting the planetesimal mass to physical mass units (to normalize by ), the gravity parameter is folded into the calculation. Thus, a more meaningful comparison between these simulations is to renormalize the planetesimal mass to the total mass in the grid, which is independent of . This is shown in Fig. 10. While the total and maximum planetesimal masses increase with as was the case in Fig. 9, the differences between the curves are significantly reduced, and the maximum mass of SI128G0.02 and SI128G0.05 appear to be roughly the same.
As with the previous simulations, we calculate the differential mass distribution. The result is shown in Fig. 11. As with the fiducial calculation, we have run an additional four simulations for each value of corresponding to selfgravity activated at and both before and after the time that it is activated in the reference simulations. As before, this approach significantly improves the statistical scatter in the mass distribution plot.
We examine the best fit power law slope as a function of in the right panel of Fig. 11. The power law slope is for , for , and for ; there is some variance in with , but these values fall within the range of –.
Despite the improvement introduced by combining simulation data, there is still significant scatter in the values of the differential mass function. To remove some of this noise, we calculate the cumulative mass distribution for each value of as shown in Fig. 12. As with the differential mass distribution, the slope of the distribution remains roughly constant with , but the distribution shifts towards larger mass.
3.4. Effect of Initial Clumping
In this section, we consider the effect of changing the time at which particle selfgravity is activated. We first consider the maximum particle mass density as a function of time for this series of simulations, which is shown in Fig. 13. There is a clear difference in the initial growth of depending on the initial degree of clumping from which the simulation was activated. At late times, the maximum particle density reaches approximately the same level for all four simulations.
These basic results are corroborated by examination of the total and maximum planetesimal mass evolution, as shown in Fig. 14. We should reiterate that the initial stages of growth should be treated with some caution here due to the possible presence of falsepositives and overlapping clumps that cause errors in our clumpfinding algorithm. However, that the initial evolution is generally consistent with the evolution of the maximum particle density in Fig. 13 is encouraging.
Given the different growth rates for the planetesimals in each of these runs, one has to be careful in choosing the correct times to analyze properties of the planetesimal distribution. Our general method has been the same as described above; visually examine the collapse of planetesimals and then average the mass distribution over a period corresponding to roughly when individual clumps can be identified but before there is significant merging between these clumps.
As shown in Fig. 14, there is a steep growth early on in SI128SGweak_clump, followed by a relatively flat evolution, which is then followed by growth again. This second growth period, which happens roughly between and , is due to a second phase of planetesimal formation. From an examination of the particle surface density evolution, this second phase of planetesimal formation appears to follow the formation of another largely axisymmetric enhancement in the particle density. Apparently, in this particular simulation, the conditions allowed the streaming instability to act on remaining small solids in the disk after the initial formation of planetesimals. The second density enhancement induced by the streaming instability then went gravitationally unstable and added to the total number of planetesimals.
In choosing our analysis time for SI128G0.05weak_clump, we thus chose a time after the second period of planetesimal growth. At the chosen times, SI128G0.05no_clump produces 11 clumps, SI128G0.05weak_clump produces 10, SI128G0.05med_clump produces 13 and SI128G0.05strong_clump produces 11.
In Fig. 15, we plot the differential mass distribution for the four different start times. There is significant scatter in the differential mass function, but the power law indices appear roughly consistent between the different simulations. Specifically, for SI128G0.05no_clump, for SI128G0.05weak_clump, for SI128G0.05med_clump, and for SI128G0.05strong_clump. These values are slightly larger than that found for the previously discussed simulations, which is a result of smaller number statistics biasing the fitted slope value, as we found for the fiducial run in Section 3.1.
As we did in Section 3.3, we calculate the cumulative mass distribution to reduce the noise inherent in the differential distribution. This is shown in Fig. 16.
Roughly speaking, there does not appear to be significant differences between the mass distributions in each of these simulations. They all occupy roughly the same space in the differential distribution plot. Based on all of these results combined, it seems that the properties of the planetesimals do not strongly depend on the initial state from which they collapse.
4. Discussion
The initial size distribution of planetesimals has a number of astrophysical implications. Planetesimals that are too small may be subject to turbulent stirring and collisional destruction (Ormel & Okuzumi, 2013), while planetesimals that are too large will accrete less efficiently on to giant planet core as a consequence of less efficient gravitational focusing. While important, independent uncertainties in the modeling of these processes make it hard to use them for quantitative constraints on planetesimal size. The best prospect for a direct comparison between theory and observations thus remains studies of the size distribution of small bodies in the Solar System. From our simulations we predict that generally lies in the range 1.4–1.8, which equates to a slope of the differential size distribution –. For reasonable choices of parameters, we find that the largest bodies in the region of the asteroid belt have masses , corresponding to a size of several hundred km.
The two largest populations of small bodies in the Solar System with which to compare our results are the main asteroid belt and the Kuiper belt. The asteroid belt has experienced substantial dynamical depletion (which does not alter the sizes of surviving bodies) and collisional evolution in the time since it formed, and detailed modeling is needed to assess whether the current size distribution preserves information about the primordial population. Using such modeling, Morbidelli et al. (2009) have argued that reproducing the current size distribution of asteroids requires primordial planetesimals with diameters , and that the slope of the current size distribution above roughly this scale is similar to the initial slope. To the extent that our work agrees with prior simulations in predicting the prompt formation of some very massive planetesimals, the Morbidelli et al. (2009) analysis supports the scenario of streaminginitiated planetesimal formation. However, the power law index for the size distribution of asteroids larger than 120 km in diameter (which is relevant to our calculations) is (Jedicke et al., 2002; Bottke et al., 2005). This would correspond to , which is significantly steeper than the powerlaw portion of the predicted mass function. We note, however, that if the most massive planetesimals are identified with the most massive (surviving) asteroids, then the observed slope would correspond to the region near the cutoff in the primordial mass function, and would be expected to be steeper than the powerlaw fit seen at lower masses.
Similar arguments apply to the Kuiper Belt. Fraser et al. (2010) calculated the observed luminosity function of the Kuiper belt and found a value of for the cold classical Kuiper belt objects and for the excited, hot Kuiper belt population. Later corrections to this work (Fraser et al., 2014) that removed the model dependence associated with uncertainties in the distance to KBOs found even steeper values of (cold population) and (hot population). Moreover, the largest Kuiper belt objects are bigger than those in the main asteroid belt. This is consistent with a streaming instability model, since a robust inference from our simulations with varying is that for larger gravity parameters, more massive planetesimals tend to form. In a minimum mass solar nebula model (Hayashi, 1981), slowly increases with radius (at the location of the asteroid belt, , and at 40 AU, ), suggesting that larger planetesimals should be formed at further distances from the Sun. However, as with the asteroid belt, the observed slopes are not what we predict for the primordial population.
The above comparisons are suggestive but should be regarded as preliminary. In addition to the uncertainties in deciphering what the current size distributions in the asteroid and Kuiper belts tell us about the primordial population, there are several sources of possible uncertainty in the prediction of the size distribution itself. First, at a purely numerical level, the collapse of planetesimals in our simulations is halted at the grid scale, which is unphysically large. The crosssection for subsequent accretion or mergers is therefore boosted, which could impact the measured size distribution. We note, however, that the recent simulations of Johansen et al. (2015) — in which collapsed planetesimals were replaced with sink particles — yielded roughly consistent values of and as compared to our simulation. Thus, it would seem that there are no systematic errors in the determination of with high resolution simulations.
The largest scales in the box may also play an important role in properties of collapsed planetesimals. Indeed, it has been recently shown through nonselfgravitating simulations that the domain size influences the temporal and spatial properties of particle clumping during the nonlinear state of the streaming instability (Yang & Johansen, 2014). It is entirely conceivable that such effects will influence the outcome of the gravitational collapse phase, and we are currently pursuing such investigations (Li et al., in prep).
Numerical effects aside, however, it is possible that the size distribution of planetesimals formed from gravitational collapse is not universal, and may be different in the physical environment relevant to the asteroid or Kuiper belts. The streaming instability can produce strong clumping of solids across a wide region in the parameter space of metallicity, radial gas pressure gradient and particle size (Bai & Stone, 2010c; Carrera et al., 2015). These parameters are expected to vary substantially with disk radius (for example, the dimensionless stopping time for particles in the asteroid belt region is likely to be much smaller than for the Kuiper Belt region), and this may result in different size distributions after gravitational collapse.
It is also conceivable that the size distribution is a function of how unstable the disk is to the streaming instability and gravitational collapse. We have simulated a system in which a significant fraction of the total mass of solids forms planetesimals on a time scale of just . Meteoritic evidence from the asteroid belt, conversely, suggests a broad spread in the formation times of primitive material (Villeneuve et al., 2009). This can be interpreted as implying that planetesimal precursors persist in the disk over time scales more akin to . It is not known whether a marginally unstable system would form the same planetesimal mass function as the one that we have simulated.
Furthermore, it is worth noting that the size distribution at the end of the gravitational collapse can be subsequently modified by the longer term accretion of planetesimals and smaller solids. The effects of this accretion depend on the evolving size distribution of solids and the evolution of the gas, including its turbulent state (e.g., Johansen et al., 2015). However, the size distribution at the end of the collapse phase, as studied here, is important as input for these longer term accretion studies.
5. Conclusions
We have developed and tested a module for the mutual gravitational interaction of particles in the Athena code and applied it to preliminary studies of the streaming instability under the influence of particle gravity. As this is the first study of planetesimal formation with the Athena code, we have carried out a basic parameter sweep in order to provide a baseline of calculations and results from which to spawn further investigations and with which to compare the results already present in the literature.
In this paper, we have varied the numerical resolution, the relative strength of gravity and tidal effects, and the degree of clumping induced by the streaming instability before the activation of particle selfgravity. While we will follow up this study in future papers to explore more parameters, we can draw some preliminary conclusions from this work.

The streaming instability leads to enhanced particle clumping, after which the mutual gravitational attraction between solid particles leads to the formation of a number of bound planetesimals.

For the choice of metallicity and stopping time used here, the masses of these planetesimals range from to . The typical radii of these planetesimals are 50 km to a few hundred km.

Where a direct comparison is possible, we find excellent agreement between planetesimal properties in our Athena simulations and those carried out with the Pencil code.

As resolution is increased, more planetesimals are produced at lower masses and smaller radii, while the high end mass of the distribution remains approximately the same. There is no significant trend of the power law index with resolution. The values of fall within the range –.

The power law slope of the size distribution for the highest resolution simulation is , which is significantly shallower than that measured for both the main asteroid belt () and the classical Kuiper belt (cold population; hot population; ).

Varying the relative strength of gravity through the parameter changes the mass range of planetesimals produced. The value of does not appear to have any consistent trend with .

The properties of planetesimals appear largely independent of the initial degree of clumping before gravity takes over, justifying the method employed here and in various papers by Johansen and coauthors of activating particle selfgravity after the streaming instability has already produced clumps.
As this is an initial step into planetesimal studies with Athena, we have only explored a few parameters. A more complete study of planetesimal properties will require varying other parameters, both physical (e.g., metallicity, particle size, radial pressure gradient) and numerical (e.g., boundary conditions and domain size). We will address these issues in future publications.
A more serious uncertainty lies in the relatively small mass and size range produced by the moderate resolution simulations, which generally spans only a decade in mass. While we were able to alleviate the resulting issue of small number statistics associated with this limited mass range by combining datasets, the small mass end of the distribution is ultimately determined by the finite grid scale. Going to higher resolution and a larger number of particles is the obvious solution to this problem, and indeed the highest resolution simulation spanned nearly two mass decades. However, such a simulation is very computationally expensive and this expense makes a large exploration of parameter space infeasible.
In a similar vein, at any given resolution, we do not resolve compact planetesimals; i.e., the planetesimals that do form are large compared to the size of planetesimals that would form if we resolved particle selfgravity below the grid scale. These unphysical sizes can cause enhanced accretion of smaller solids onto formed planetesimals and make the power law distributions shallower. The future implementation of sink particles or subgridcell gravity will help to alleviate these issues.
Despite these uncertainties, our first principles calculations provide additional evidence that the streaming instability can form planetesimals with properties consistent with Solar System constraints. The largest remaining source of uncertainty may be in the initial metallicity and particle size distribution (i.e., prior to the onset of the streaming instabilty). If the appropriate initial conditions can be pinned down via astronomical observations, more comprehensive numerical simulations should be able to provide a detailed understanding of the formation and evolution of the planetesimal population.
References
 Andrews et al. (2012) Andrews, S. M., Wilner, D. J., Hughes, A. M., et al. 2012, The Astrophysical Journal, 744, 162
 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
 Birnstiel et al. (2010) Birnstiel, T., Dullemond, C. P., & Brauer, F. 2010, Astronomy and Astrophysics
 Birnstiel et al. (2012) Birnstiel, T., Andrews, S. M., & Ercolano, B. 2012, A&A, 544, A79
 Bitsch et al. (2014) Bitsch, B., Morbidelli, A., Lega, E., Kretke, K., & Crida, A. 2014, A&A, 570, A75
 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
 Brauer et al. (2008) Brauer, F., Dullemond, C. P., & Henning, T. 2008, Astronomy and Astrophysics, 480, 859
 Carrera et al. (2015) Carrera, D., Johansen, A., & Davies, M. B. 2015, arXiv.org
 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
 Dra̧żkowska et al. (2013) Dra̧żkowska, J., Windmark, F., & Dullemond, C. P. 2013, A&A, 556, A37
 Dra̧żkowska & Dullemond (2014) Dra̧żkowska, J., & Dullemond, C. P. 2014, A&A, 572, A78
 Dzyurkevich et al. (2010) Dzyurkevich, N., Flock, M., Turner, N. J., Klahr, H., & Henning, T. 2010, A&A, 515, 70
 Fraser et al. (2014) Fraser, W. C., Brown, M. E., Morbidelli, A., Parker, A., & Batygin, K. 2014, The Astrophysical Journal, 782, 100
 Fraser et al. (2010) Fraser, W. C., Brown, M. E., & Schwamb, M. E. 2010, Icarus, 210, 944
 Gardiner & Stone (2005) Gardiner, T. A., & Stone, J. M. 2005, JCP, 205, 509
 Gardiner & Stone (2008) —. 2008, JCP, 227, 4123
 Hawley et al. (1995) Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1995, ApJ, 440, 742
 Hayashi (1981) Hayashi, C. 1981, Progress of Theoretical Physics Supplement, 70, 35
 Jedicke et al. (2002) Jedicke, R., Larsen, J., & Spahr, T. 2002, Asteroids III, 71
 Johansen et al. (2011) Johansen, A., Klahr, H., & Henning, T. 2011, Astronomy and Astrophysics, 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 & Youdin (2007) Johansen, A., & Youdin, A. 2007, The Astrophysical Journal, 662, 627
 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
 Kataoka et al. (2013) Kataoka, A., Tanaka, H., Okuzumi, S., & Wada, K. 2013, A&A, 557, L4
 Koyama & Ostriker (2009) Koyama, H., & Ostriker, E. C. 2009, The Astrophysical Journal, 693, 1316
 Kretke & Lin (2007) Kretke, K. A., & Lin, D. N. C. 2007, ApJ, 664, L55
 Krijt et al. (2015) Krijt, S., Ormel, C. W., Dominik, C., & Tielens, A. G. G. M. 2015, A&A, 574, A83
 Masset (2000) Masset, F. 2000, A&AS, 141, 165
 Morbidelli et al. (2009) Morbidelli, A., Bottke, W. F., Nesvorný, D., & Levison, H. F. 2009, Icarus, 204, 558
 Musiolik et al. (2016) Musiolik, G., Teiser, J., Jankowski, T., & Wurm, G. 2016, ApJ, 818, 16
 Okuzumi et al. (2012) Okuzumi, S., Tanaka, H., Kobayashi, H., & Wada, K. 2012, ApJ, 752, 106
 Ormel & Okuzumi (2013) Ormel, C. W., & Okuzumi, S. 2013, The Astrophysical Journal, 771, 44
 Sekora & Colella (2009) Sekora, M., & Colella, P. 2009, arXiv:0903.4200
 Simon et al. (2011) Simon, J. B., Hawley, J. F., & Beckwith, K. 2011, ApJ, 730, 94
 Simon & Armitage (2014) Simon, J. B., & Armitage, P. J. 2014, ApJ, 784, 15
 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
 Toro (1999) Toro, E. F. 1999, Riemann solvers and numerical models for fluid dynamics
 Villeneuve et al. (2009) Villeneuve, J., Chaussidon, M., & Libourel, G. 2009, Science, 325, 985
 Wada et al. (2013) Wada, K., Tanaka, H., Okuzumi, S., et al. 2013, A&A, 559, A62
 Weidenschilling (1977) Weidenschilling, S. J. 1977, MNRAS, 180, 57
 Yang & Menou (2010) Yang, C.C., & Menou, K. 2010, MNRAS, 402, 2436
 Yang & Johansen (2014) Yang, C.C., & Johansen, A. 2014, ApJ, 792, 86
 Youdin & Johansen (2007) Youdin, A., & Johansen, A. 2007, The Astrophysical Journal, 662, 613
 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