THE MASS AND SIZE DISTRIBUTION OF PLANETESIMALS FORMED BY THE STREAMING INSTABILITY. I. THE ROLE OF SELF-GRAVITY

Our simulations use Athena, a second-order accurate Godunov flux-conservative 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 third-order 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, 2008), and Stone et al. (2008).

The simulations employ a local shearing box approximation. The shearing box models a co-rotating disk patch whose size is small compared to the radial distance from the central object, R₀. This allows the construction of a local Cartesian frame $(x,y,z)$ that is defined in terms of the disk's cylindrical co-ordinates $(R,\phi ,{z}^{\prime })$ via $x=(R-{R}_{0})$, $y={R}_{0}\phi$, and $z={z}^{\prime }$. The local patch co-rotates with an angular velocity Ω corresponding to the orbital frequency at R₀, 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:

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

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

where ρ is the mass density, $\rho {\boldsymbol{u}}$ is the momentum density, P is the gas pressure, ${\boldsymbol{I}}$ is the identity matrix, ${t}_{{\rm{stop}}}$ is the (dimensional) stopping time of the particles, and q is the shear parameter, defined as q = −dlnΩ/dlnR. We use $q=3/2$, appropriate for a Keplerian disk. For simplicity and numerical convenience, we assume an isothermal equation of state $P=\rho {c}_{{\rm{s}}}^{2}$, where ${c}_{{\rm{s}}}$ 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 ${\rho }_{{\rm{p}}}$, the difference between the particle velocity ${\boldsymbol{v}}$ and gas velocity, and the stopping time ${t}_{{\rm{stop}}}$ 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 i is given by

$$\begin{eqnarray}\displaystyle \frac{d{{\boldsymbol{v}}}_{i}^{\prime }}{{dt}} & = & 2({v}_{{iy}}^{\prime }-\eta {v}_{{\rm{K}}}){\rm{\Omega }}\hat{{\boldsymbol{x}}}-(2-q){v}_{{ix}}^{\prime }{\rm{\Omega }}\hat{{\boldsymbol{y}}}\\ & & -{{\rm{\Omega }}}^{2}z\hat{{\boldsymbol{z}}}-\displaystyle \frac{{{\boldsymbol{v}}}_{i}^{\prime }-{{\boldsymbol{u}}}^{\prime }}{{t}_{{\rm{stop}}}}+{{\boldsymbol{F}}}_{{\rm{g}}},\end{eqnarray} \tag{ 3 }$$

where the prime denotes a frame in which the background shear velocity has been subtracted. This is part of the orbital advection scheme (Masset 2000; Johnson et al. 2008; Stone & Gardiner 2010), which has been implemented for the gas dynamics as well. The $\eta {v}_{{\rm{K}}}$ 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 sub-Keplerian 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 $\eta {v}_{{\rm{K}}}$ term as described above. As a result, both the particles and gas (azimuthal) velocities are shifted to slightly higher values (by $\eta {v}_{{\rm{K}}}$) 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 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 (${{\boldsymbol{u}}}^{\prime }$).

The force due to the particle self-gravity is denoted by ${{\boldsymbol{F}}}_{{\rm{g}}}$, which is found by first solving Poisson's equation for particle self-gravity,

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

where ${{\rm{\Phi }}}_{{\rm{p}}}$ is the gravitational potential of the particle self-gravity. The force is then calculated via,

$$\begin{eqnarray}&&{{\boldsymbol{F}}}_{{\rm{g}}}=-{\boldsymbol{\nabla }}{{\rm{\Phi }}}_{{\rm{p}}}.\end{eqnarray} \tag{ 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 shearing-periodic boundaries in x and the vacuum boundaries in z. 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 shearing-periodic boundaries are purely periodic. From Hawley et al. (1995), the times at which the radial boundaries are purely periodic are ${t}_{n}={{nL}}_{y}/(q{\rm{\Omega }}{L}_{x})$ with $n=0,1,2,3...$. For each value of x, the particle density is reconstructed along y via a conservative remap; the density is calculated by differencing numerical "fluxes" along y. 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 shearing-periodic 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 self-gravitating 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 y via the same third-order conservative remap as is used for the boundary conditions, multiplied by appropriate functions of z, 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 ${L}_{x}\times {L}_{y}\times 2{L}_{z}$ 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 third-order conservative remap that was used on the particle mass density.

The resulting ${{\rm{\Phi }}}_{{\rm{p}}}$ is a cell-centered 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 x-direction in cell $(i,j,k)$ is calculated by

$$\begin{eqnarray}&&{F}_{x,g}^{i,j,k}=-\displaystyle \frac{{{\rm{\Phi }}}_{{\rm{p}}}^{i+1,j,k}-{{\rm{\Phi }}}_{{\rm{p}}}^{i-1,j,k}}{2{\rm{\Delta }}x},\end{eqnarray} \tag{ 6 }$$

and similarly for the y and z 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 self-gravity force to the particles simultaneously with the drag force.

The boundary conditions are the shearing-periodic 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, R. Li et al. 2016, in preparation). 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 (R. Li et al. 2016, in preparation), 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; shearing-periodic in x and purely periodic in y. The vertical boundary conditions are open, and the potential is calculated in the ghost zones via a third order extrapolation.

In this section, we carry out two tests of our numerical algorithm for particle self-gravity. The first test is the collapse of a uniform density sphere of particles. The second test consists of the self-gravitating, 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 r(t) starting at the outer edge of the sphere is

$$\begin{eqnarray}&&\displaystyle \frac{{d}^{2}r}{{{dt}}^{2}}=-\displaystyle \frac{{GM}}{{r}^{2}},\end{eqnarray} \tag{ 7 }$$

where M is the total mass within the sphere and remains constant in time. Parameterizing the time dependence of r(t) via $\alpha (t)$, we assume that $r(t)={r}_{0}{{\rm{cos}}}^{2}(\alpha (t))$, where r₀ is the initial radius of the uniform density sphere, and find that

$$\begin{eqnarray}&&\alpha +\displaystyle \frac{1}{2}{\rm{sin}}(2\alpha )=\sqrt{\displaystyle \frac{2{GM}}{{r}_{0}^{3}}}t.\end{eqnarray} \tag{ 8 }$$

The analytical solutions to the radius and mass density are given by

$$\begin{eqnarray}&&r(t)/{r}_{0}={{\rm{cos}}}^{2}\alpha \end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{\rho }_{{\rm{p}}}(t)/{\rho }_{{\rm{p}},0}=\displaystyle \frac{1}{{{\rm{cos}}}^{6}\alpha },\end{eqnarray} \tag{ 10 }$$

where α is given by Equation (8).

To test our self-gravity solver against this analytic solution, we initialize a sphere of particles with radius ${r}_{0}=0.25$ on a cubic domain of ${L}_{x}\times {L}_{y}\times {L}_{z}=1\times 1\times 1$ resolved by 96³ 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 self-gravity module. We use TSC interpolation, and the semi-implicit 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 8 × 10⁻³.

The solution to this test problem is shown in Figure 1, which depicts the evolution of $r(t)/{r}_{0}$ and ${\rho }_{{\rm{p}}}(t)/{\rho }_{{\rm{p}},0}$ versus time in units of the free fall time ${t}_{{\rm{ff}}}\equiv \sqrt{3\pi /32G{\rho }_{{\rm{p}},0}}$ where ${\rho }_{{\rm{p}},0}$ 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 $\sim {\rm{\Delta }}x$ where ${\rm{\Delta }}x$ is the length of a grid cell.⁷

Zoom In Zoom Out Reset image size

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 ${ \mathcal O }({N}^{2})$, where N 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 does not require any interpolation of the gravitational forces to the particle locations. We calculate the error in the force in the ith direction as

$$\begin{eqnarray}&&{\epsilon }_{i}=\displaystyle \frac{1}{{N}_{{\rm{par}}}}\displaystyle \sum _{n=1}^{{N}_{{\rm{par}}}}\displaystyle \frac{{F}_{n,i,{\rm{FFT}}}-{F}_{n,i,{\rm{direct}}}}{| {F}_{n,i,{\rm{direct}}}| },\end{eqnarray} \tag{ 11 }$$

where ${F}_{n,i,{\rm{FFT}}}$ is the force due to self-gravity on particle n calculated via the FFT method and ${F}_{n,i,{\rm{direct}}}$ is the same but calculated with the direct summation method. We run the particle sphere collapse problem at a resolution of 32³ (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 10⁻⁵ 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 self-gravity solver will become more and more dominant. At $t=0.5{t}_{{\rm{ff}}}$, the average errors are approximately 1% or less, and near $t\approx {t}_{{\rm{ff}}}$, the average errors are on order of a few percent.

2.2.2. Self-gravitating, 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 self-gravitating, 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 ${\rho }_{p}={\rho }_{p,0}+{\rho }_{p}^{\prime }$, velocity ${\boldsymbol{v}}={{\boldsymbol{v}}}_{0}+{{\boldsymbol{v}}}^{\prime }$, and self-gravity potential ${\rm{\Phi }}={{\rm{\Phi }}}_{0}+{{\rm{\Phi }}}^{\prime }$ where prime represents the linear perturbation on the background. The perturbation is assumed to be of the form ${q}^{\prime }(t,x,y)=\delta q(t){\rm{Exp}}[i({k}_{x}(t)x+{k}_{y}y)]$. Note that we assume uniformity along the vertical direction. The evolution of the density and velocity perturbations are then governed by the following equations:

$$\begin{eqnarray}&&\displaystyle \frac{d\delta {\rho }_{p}}{{dt}}=-{\rho }_{p,0}i[{k}_{x}(t)\delta {v}_{x}+{k}_{y}\delta {v}_{y}],\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d\delta {v}_{x}}{{dt}}=2{\rm{\Omega }}\delta {v}_{y}+\displaystyle \frac{4\pi {{iGk}}_{x}(t)\delta {\rho }_{p}}{{k}_{x}^{2}(t)+{k}_{y}^{2}}-\displaystyle \frac{\delta {v}_{x}}{{t}_{{\rm{stop}}}},\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d\delta {v}_{y}}{{dt}}=-\displaystyle \frac{1}{2}{\rm{\Omega }}\delta {v}_{x}+\displaystyle \frac{4\pi {{iGk}}_{y}\delta {\rho }_{p}}{{k}_{x}^{2}(t)+{k}_{y}^{2}}-\displaystyle \frac{\delta {v}_{y}}{{t}_{{\rm{stop}}}},\end{eqnarray} \tag{ 14 }$$

where the radial wavenumber k_x is time-dependent,

$$\begin{eqnarray}&&{k}_{x}(t)={k}_{x}(0)+q{\rm{\Omega }}{{tk}}_{y}.\end{eqnarray} \tag{ 15 }$$

We set $G={t}_{{\rm{stop}}}={\rm{\Omega }}={\rho }_{p,0}=1$, ${{\boldsymbol{v}}}_{0}={{\rm{\Phi }}}_{0}=\delta {\rm{\Phi }}\;=\delta {v}_{x}=\delta {v}_{y}={k}_{x}(0)=0$, k_y = 1, q = 1.5, and $\delta {\rho }_{p}={10}^{-4}$. We solve Equations (12)–(14) using a simple finite difference algorithm to obtain our linear solution to the self-gravitating, shearing wave problem in the linear limit. For the numerical test of our algorithm, we initialize a grid of size ${L}_{x}\times {L}_{y}\times {L}_{z}=2\pi \times 2\pi \times 0.2$ at resolution $128\times 128\times 6$ and with one particle per cell initially; the small size/resolution in the vertical direction makes the problem effectively two-dimensional while still testing the full 3D FFT algorithm. In this setup, the boundary conditions are shearing-periodic in x and purely periodic in both y and z. 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 Figure 2. The numerical solution agrees very well with the linear solution up to amplitudes of $\sim {10}^{-1}$. 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, self-gravity, and drag), demonstrates the validity of our self-gravity implementation in the shearing box setup.

Zoom In Zoom Out Reset image size

All of our streaming instability simulations use a shearing box domain of size ${L}_{x}\times {L}_{y}\times {L}_{z}$ = $0.2H\times 0.2H\times 0.2H$. Here, H is the vertical scale height of the gas, which is initialized with a hydrostatic (Gaussian) vertical profile,

$$\begin{eqnarray}&&{\rho }_{g}={\rho }_{0}{\rm{exp}}\left(\displaystyle \frac{-{z}^{2}}{2{H}^{2}}\right),\end{eqnarray} \tag{ 16 }$$

where ${\rho }_{0}$ is the mid-plane gas density. We choose code units so that the standard gas parameters are unity; ${\rho }_{0}=H={\rm{\Omega }}={c}_{{\rm{s}}}=1$. Given the limited vertical extent of our domain, the gas density does not vary significantly.

The particles are distributed uniformly in x and y and with a Gaussian profile in z. The scale height of this profile is ${H}_{p}=0.02H$. 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 $\eta {v}_{{\rm{K}}}$ as described above, and the lack of a pressure gradient to counteract vertical gravity means that the particles will settle toward the disk mid-plane.

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,

$$\begin{eqnarray}&&\tau \equiv {t}_{{\rm{stop}}}{{\rm{\Omega }}}^{-1},\end{eqnarray} \tag{ 17 }$$

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

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

which is the ratio of the particle mass surface density ${{\rm{\Sigma }}}_{p}$ to the gas surface density ${{\rm{\Sigma }}}_{g}$, and a radial pressure gradient parameter that accounts for the sub-Keplerian gas in real disks

$$\begin{eqnarray}&&{\rm{\Pi }}\equiv \eta {v}_{{\rm{K}}}/{c}_{{\rm{s}}}.\end{eqnarray} \tag{ 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, $\tau =0.3$, metallicity is moderately super-Solar, Z = 0.02, and ${\rm{\Pi }}=0.05$.

With the inclusion of self-gravity, an additional parameter is needed to describe the relative importance of self-gravity and tidal shear. We define a parameter $\tilde{G}$,

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

which describes the strength of self-gravity in the simulation. Physically, varying $\tilde{G}$ is equivalent to changing the gas density (and thus through assuming the same Z, the particle mass density) or the strength of tidal stretching (i.e., changing Ω) or both. Thus, in a real disk, $\tilde{G}$ 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 ${M}_{0}={\rho }_{0}{H}^{3}$ and assuming a radius of 3 au in an MMSN. For $\tilde{G}=0.05$, this equates to ${M}_{0}=6.5\times {10}^{26}\;{\rm{g}}\approx 720{M}_{{\rm{Ceres}}}$. This physical unit conversion depends on $\tilde{G}$, as we discuss further below.

Our fiducial simulation is a relatively low-resolution run with ${N}_{x}\times {N}_{y}\times {N}_{z}$ = $128\times 128\times 128$ gas zones and 2,400,000 particles. For this simulation, we turn on self-gravity at a time ${t}_{{\rm{sg}}}=170{{\rm{\Omega }}}^{-1}$ and set $\tilde{G}=0.05$.

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 $64\times 64\times 64$ with 300,000 particles, and two higher resolutions: $256\times 256\times 256$ with 19,200,000 particles and $512\times 512\times 512$ with 153,600,000 particles. In the next subset, we vary the value of $\tilde{G}$ 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.

Table 1.  Streaming Instability Simulations

Run	Domain Size	Resolution	${N}_{{\rm{par}}}$	τ	Z	$\tilde{G}$	${t}_{{\rm{sg}}}$	Comments
	$({L}_{x}\times {L}_{y}\times {L}_{z})H$	${N}_{x}\times {N}_{y}\times {N}_{z}$					(${{\rm{\Omega }}}^{-1}$)
SI64-G0.05	$0.2\times 0.2\times 0.2$	$64\times 64\times 64$	300,000	0.3	0.02	0.05	400	⋯
SI128-G0.05	$0.2\times 0.2\times 0.2$	$128\times 128\times 128$	2,400,000	0.3	0.02	0.05	170	Fiducial Run
SI256-G0.05	$0.2\times 0.2\times 0.2$	$256\times 256\times 256$	19,200,000	0.3	0.02	0.05	150	⋯
SI512-G0.05	$0.2\times 0.2\times 0.2$	$512\times 512\times 512$	153,600,000	0.3	0.02	0.05	110	⋯
SI128-G0.02	$0.2\times 0.2\times 0.2$	$128\times 128\times 128$	2,400,000	0.3	0.02	0.02	170	⋯
SI128-G0.02_tm20	$0.2\times 0.2\times 0.2$	$128\times 128\times 128$	2,400,000	0.3	0.02	0.02	150	Restarted $20{{\rm{\Omega }}}^{-1}$ earlier
SI128-G0.02_tm10	$0.2\times 0.2\times 0.2$	$128\times 128\times 128$	2,400,000	0.3	0.02	0.02	160	Restarted $10{{\rm{\Omega }}}^{-1}$ earlier
SI128-G0.02_tp10	$0.2\times 0.2\times 0.2$	$128\times 128\times 128$	2,400,000	0.3	0.02	0.02	180	Restarted $10{{\rm{\Omega }}}^{-1}$ later
SI128-G0.02_tp20	$0.2\times 0.2\times 0.2$	$128\times 128\times 128$	2,400,000	0.3	0.02	0.02	190	Restarted $20{{\rm{\Omega }}}^{-1}$ later
SI128-G0.05	$0.2\times 0.2\times 0.2$	$128\times 128\times 128$	2,400,000	0.3	0.02	0.05	170	Fiducial Run
SI128-G0.05_tm20	$0.2\times 0.2\times 0.2$	$128\times 128\times 128$	2,400,000	0.3	0.02	0.05	150	Restarted $20{{\rm{\Omega }}}^{-1}$ earlier
SI128-G0.05_tm10	$0.2\times 0.2\times 0.2$	$128\times 128\times 128$	2,400,000	0.3	0.02	0.05	160	Restarted $10{{\rm{\Omega }}}^{-1}$ earlier
SI128-G0.05_tp10	$0.2\times 0.2\times 0.2$	$128\times 128\times 128$	2,400,000	0.3	0.02	0.05	180	Restarted $10{{\rm{\Omega }}}^{-1}$ later
SI128-G0.05_tp20	$0.2\times 0.2\times 0.2$	$128\times 128\times 128$	2,400,000	0.3	0.02	0.05	190	Restarted $20{{\rm{\Omega }}}^{-1}$ later
SI128-G0.1	$0.2\times 0.2\times 0.2$	$128\times 128\times 128$	2,400,000	0.3	0.02	0.1	170	⋯
SI128-G0.1_tm20	$0.2\times 0.2\times 0.2$	$128\times 128\times 128$	2,400,000	0.3	0.02	0.1	150	Restarted $20{{\rm{\Omega }}}^{-1}$ earlier
SI128-G0.1_tm10	$0.2\times 0.2\times 0.2$	$128\times 128\times 128$	2,400,000	0.3	0.02	0.1	160	Restarted $10{{\rm{\Omega }}}^{-1}$ earlier
SI128-G0.1_tp10	$0.2\times 0.2\times 0.2$	$128\times 128\times 128$	2,400,000	0.3	0.02	0.1	180	Restarted $10{{\rm{\Omega }}}^{-1}$ later
SI128-G0.1_tp20	$0.2\times 0.2\times 0.2$	$128\times 128\times 128$	2,400,000	0.3	0.02	0.1	190	Restarted $20{{\rm{\Omega }}}^{-1}$ later
SI128-G0.05-no_clump	$0.2\times 0.2\times 0.2$	$128\times 128\times 128$	2,400,000	0.3	0.02	0.05	0	⋯
SI128-G0.05-low_clump	$0.2\times 0.2\times 0.2$	$128\times 128\times 128$	2,400,000	0.3	0.02	0.05	40	⋯
SI128-G0.05-med_clump	$0.2\times 0.2\times 0.2$	$128\times 128\times 128$	2,400,000	0.3	0.02	0.05	170	Same as Fiducial Run
SI128-G0.05-high_clump	$0.2\times 0.2\times 0.2$	$128\times 128\times 128$	2,400,000	0.3	0.02	0.05	240	⋯

Download table as:  ASCIITypeset image

Finally, we determine the effect of the state at which self-gravity is initiated by examining two diagnostics in a simulation identical to the fiducial one, but with no self-gravity. The first is the maximum value of the volume-averaged particle mass density, ${\rho }_{p}$, and the second is a weighted version of the volume-averaged particle mass density, defined as $\sqrt{\langle {\rho }_{p}^{2}\rangle }/\langle {\rho }_{p}\rangle$. 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 ${t}_{{\rm{sg}}}=40{{\rm{\Omega }}}^{-1}$ (the "low clumping" case), ${t}_{{\rm{sg}}}=170{{\rm{\Omega }}}^{-1}$ (medium clumping; this is the same as the fiducial simulation), and ${t}_{{\rm{sg}}}=240{{\rm{\Omega }}}^{-1}$ (high clumping). We also run one case with self-gravity turned on from the beginning state ${t}_{{\rm{sg}}}=0$.

All of these runs are labelled by the numerical resolution employed, the strength of self-gravity, and any other modifying information. For example, run SI128-G0.05-low_clump has 128 zones per dimension, $\tilde{G}=0.05$, 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.

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 ${\rho }_{p}$ over the vertical extent of the domain.

Most of our diagnostics require that we locate gravitationally bound mass clumps (after self-gravity has been included) and calculate their mass. Following the basic arguments in Johansen et al. (2011), we calculate ${{\rm{\Sigma }}}_{p}$ at a given time and use a peak-finding 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 self-gravity compared to tidal effects. Mathematically, we seek a radius R such that the tidal term, $3{{\rm{\Omega }}}^{2}R$, and the gravitational acceleration of a test particle at the edge of the clump's Hill sphere, ${{GM}}_{p}/{R}^{2}$ are in balance,

$$\begin{eqnarray}&&{M}_{p}=\displaystyle \frac{3{{\rm{\Omega }}}^{2}{R}^{3}}{G};\end{eqnarray} \tag{ 21 }$$

when the mass enclosed within our test circle equals M_p, we set R as the Hill radius and M_p 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 xy plane. This is alleviated somewhat by setting a minimum distance, equivalent to the Hill radius of a 0.1 ${M}_{{\rm{Ceres}}}$ 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 ${{\rm{\Sigma }}}_{p}$. 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 two-dimensional approach to a fully three-dimensional problem. However, even though clumps exist in three dimensions, the particle layer is vertically thin enough to make the two-dimensional 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 power-law differential mass distribution,

$$\begin{eqnarray}&&\displaystyle \frac{{dN}}{{{dM}}_{p}}={C}_{1}{M}_{p}^{-p},\end{eqnarray} \tag{ 22 }$$

works well. Here, C₁ is a constant that for the purposes of this paper is arbitrary. We also consider the cumulative mass distribution,

$$\begin{eqnarray}&&N(\gt {M}_{p})={C}_{2}{M}^{-c},\end{eqnarray} \tag{ 23 }$$

where C₂ is another arbitrary constant. We can easily translate Equation (22) to a size distribution to find ${dN}/{dr}\propto {r}^{-q}$ where $q=3p-2$.

Under the assumption that the data is drawn from a power-law distribution, the power-law index p of the differential distribution can be determined directly from the set of measured masses ${M}_{p,i}$ 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

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

where ${M}_{p,{\rm{min}}}$ is the minimum value of the planetesimal mass in our data, and n is the total number of planetesimals. The error is given by their Equation (3.2),

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

Using the MLE for p avoids any binning step, which can introduce bias into the estimate. To visually represent the differential distribution, however, we make a local estimate of ${dN}/{{dM}}_{p}$ by taking ${{dM}}_{p}$ to be half the distance in the mass co-ordinate between a given planetesimal and its nearest neighbors. A least squares linear fit in log–log space to ${dN}/{{dM}}_{p}$ versus M_p gives an alternate measure of p, 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 power law 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 power law then our naive fit will over-estimate the value of p. Second, our lower-resolution runs form rather small number of planetesimals, and any method for estimating p in this situation is noisy, possibly biased, and liable to return a non-Gaussian 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.

In this section, we describe some properties of our fiducial run, SI128-G0.05. Self-gravity is switched on at $t=170{{\rm{\Omega }}}^{-1}$, after which the mutual gravitational attraction between particles causes the particle density to increase; in units of the initial mid-plane gas density, the maximum particle density rapidly increases to $\sim 5\times {10}^{3}$, and then slowly increases to $\sim 3\times {10}^{4}$ afterward.

A time progression of the particle surface density is shown by a series of snapshots in Figure 3. The streaming instability produces several azimuthally extended structures at early times, which eventually form into one large clumping structure. After self-gravity 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.

Zoom In Zoom Out Reset image size

The evolution of the mass of these planetesimals is shown in Figure 4; both the total mass and maximum mass are shown. This evolution can be compared to Figure 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 $t-{t}_{{\rm{sg}}}=20{{\rm{\Omega }}}^{-1}$ in our fiducial simulation and the $\tilde{G}=0.05$ 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.

Zoom In Zoom Out Reset image size

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 increases toward 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., Figure 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; $t-{t}_{{\rm{sg}}}=10.9{{\rm{\Omega }}}^{-1}$, $t-{t}_{{\rm{sg}}}=12.3{{\rm{\Omega }}}^{-1}$ and $t-{t}_{{\rm{sg}}}=13.6{{\rm{\Omega }}}^{-1}$. 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 "by-eye" agreement in the values of M_p and ${dN}/{{dM}}_{p}$. 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 p. The best-fit value of p from this simulation is $p=2.0\pm 0.3;$ clearly there is a large error on p. 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" non-self-gravitating simulation at different times. Specifically, we initiated self-gravity at $20{{\rm{\Omega }}}^{-1}$ before, $10{{\rm{\Omega }}}^{-1}$ before, $10{{\rm{\Omega }}}^{-1}$ after, and $20{{\rm{\Omega }}}^{-1}$ after the fiducial run. These runs are included in Table 1; the run appended with "tm20" ("tp20") means self-gravity is initiated at the fiducial case restart time minus (plus) $20{{\rm{\Omega }}}^{-1}$. 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 Figure 5. We fit this data with a power law and now find $p=1.5\pm 0.1;$ the error has been reduced by roughly a factor of three.

Zoom In Zoom Out Reset image size

There is a question of whether or not the properties of the distribution of planetesimal masses depend on the time at which self-gravity 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.

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 self-gravity 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 64³ to 512³, the number of formed planetesimals at this time are 2, 13, 30, and 53, respectively.

Zoom In Zoom Out Reset image size

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 self-gravity was turned on in Figure 7. As both Figures 6 and 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 self-gravity 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.

Zoom In Zoom Out Reset image size

This resolution effect is also demonstrated via the differential mass distribution, shown in Figure 8. We do not include a power-law fit to either the SI64-G0.05 planetesimals since there are only two objects formed at that resolution or to the SI128-G0.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.

Zoom In Zoom Out Reset image size

With the exception of the lowest resolution, the highest mass planetesimals are ${M}_{p}\sim 0.1{M}_{{\rm{Ceres}}}$ for all resolutions. The higher masses (by only a factor of a few) in the 64³ 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 $p=1.5\pm 0.1$ for SI128-G0.05, $p=1.6\pm 0.1$ for SI256-G0.05, and $p=1.6\pm 0.1$ for SI512-G0.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 Figure 8, we plot the power-law index p with the errors given by Equation (25). There does not appear to be any significant trend in the value of p with resolution, and $p\approx 1.4$–1.8.

Our best-fit value of p for the highest-resolution simulation is the same as that found for the 512³ simulation of Johansen et al. (2015); they found a value of p = 1.6 by fitting the cumulative distribution to a power law with an exponential tail. Furthermore, comparing the left panel of Figure 8 to Figure 3 of Johansen et al. (2015) reveals similar scatter about the best-fit line for each resolution.

We vary $\tilde{G}$ to match a similar exploration of this parameter in Johansen et al. (2012); specifically, $\tilde{G}=0.02,0.05,0.1$. Note that with the exception of one of the simulations with $\tilde{G}=0.1$, 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 $\tilde{G}$. 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 $\tilde{G}$, and others are smaller. Given the differences likely present in the clump finding algorithms, we believe that the agreement is pretty strong.

Zoom In Zoom Out Reset image size

In converting the planetesimal mass to physical mass units (to normalize by ${M}_{{\rm{Ceres}}}$), the gravity parameter $\tilde{G}$ 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 $\tilde{G}$. This is shown in Figure 10. While the total and maximum planetesimal masses increase with $\tilde{G}$ as was the case in Figure 9, the differences between the curves are significantly reduced, and the maximum mass of SI128-G0.02 and SI128-G0.05 appear to be roughly the same.

Zoom In Zoom Out Reset image size

As with the previous simulations, we calculate the differential mass distribution. The result is shown in Figure 11. As with the fiducial calculation, we have run an additional four simulations for each value of $\tilde{G}$ corresponding to self-gravity activated at $20{{\rm{\Omega }}}^{-1}$ and $10{{\rm{\Omega }}}^{-1}$ 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.

Zoom In Zoom Out Reset image size

We examine the best-fit power-law slope as a function of $\tilde{G}$ in the right panel of Figure 11. The power-law slope is $p=1.7\pm 0.1$ for $\tilde{G}=0.02$, $p=1.5\pm 0.1$ for $\tilde{G}=0.05$, and $p=1.6\pm 0.1$ for $\tilde{G}=0.1;$ there is some variance in p with $\tilde{G}$, but these values fall within the range of p = 1.4–1.8.

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 $\tilde{G}$ as shown in Figure 12. As with the differential mass distribution, the slope of the distribution remains roughly constant with $\tilde{G}$, but the distribution shifts toward larger mass.

Zoom In Zoom Out Reset image size

In this section, we consider the effect of changing the time at which particle self-gravity is activated. We first consider the maximum particle mass density as a function of time for this series of simulations, which is shown in Figure 13. There is a clear difference in the initial growth of ${\rho }_{p,{\rm{max}}}$ 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.

Zoom In Zoom Out Reset image size

These basic results are corroborated by examination of the total and maximum planetesimal mass evolution, as shown in Figure 14. We should reiterate that the initial stages of growth should be treated with some caution here due to the possible presence of false-positives and overlapping clumps that cause errors in our clump-finding algorithm. However, that the initial evolution is generally consistent with the evolution of the maximum particle density in Figure 13 is encouraging.

Zoom In Zoom Out Reset image size

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 Figure 14, there is a steep growth early on in SI128-SG-weak_clump, followed by a relatively flat evolution, which is then followed by growth again. This second growth period, which happens roughly between $t-{t}_{{\rm{sg}}}=60{{\rm{\Omega }}}^{-1}$ and $t-{t}_{{\rm{sg}}}=100{{\rm{\Omega }}}^{-1}$, 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 SI128-G0.05-weak_clump, we thus chose a time after the second period of planetesimal growth. At the chosen times, SI128-G0.05-no_clump produces 11 clumps, SI128-G0.05-weak_clump produces 10, SI128-G0.05-med_clump produces 13, and SI128-G0.05-strong_clump produces 11.

In Figure 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, $p=2.2\pm 0.4$ for SI128-G0.05-no_clump, $p=2.1\pm 0.4$ for SI128-G0.05-weak_clump, $p=2.0\pm 0.3$ for SI128-G0.05-med_clump, and $p=1.8\pm 0.2$ for SI128-G0.05-strong_clump. These p 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.

Zoom In Zoom Out Reset image size

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 Figure 16.

Zoom In Zoom Out Reset image size

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.