Trans-Neptunian Binaries as Evidence for Planetesimal Formation by the Streaming Instability
David Nesvorny Rixin Li Andrew N. Youdin Jacob B. Simon William M. Grundy

19.pt




^1Department of Space Studies, Southwest Research Institute,
1050 Walnut St., Suite 300, Boulder, CO 80302, USA
^2Steward Observatory & Department of Astronomy, University of Arizona, 933 N. Cherry Avenue, 
Tucson, AZ, 85721, USA
^3JILA, University of Colorado, 440 UCB, Boulder, CO 80309, USA
^4Lowell Observatory, 1400 W. Mars Hill Rd., Flagstaff, AZ 86001, USA
*e-mail:davidn@boulder.swri.edu
Trans-Neptunian Binaries as Evidence for Planetesimal Formation by the Streaming Instability
    David Nesvorný^1,*, Rixin Li^2, Andrew N. Youdin^2, Jacob B. Simon^1,3,
William M. Grundy^4
    present
===============================================================================================



A critical step toward the emergence of planets in a protoplanetary disk consists 
in accretion of planetesimals, bodies 1-1000 km in size, from smaller disk constituents. 
This process is poorly understood partly because we lack good observational constraints on 
the complex physical processes that contribute to planetesimal formation [1]. In the outer solar system, 
the best place to look for clues is the Kuiper belt, where icy planetesimals survived to 
this day. Here we report evidence that Kuiper belt planetesimals formed by the streaming 
instability, a process in which aerodynamically concentrated clumps of pebbles gravitationally 
collapse into ∼100-km-class bodies [2]. Gravitational collapse was previously suggested 
to explain the ubiquity of equal-size binaries in the Kuiper belt [3,4,5]. We analyze 
new hydrodynamical simulations of the streaming instability to determine the model expectations 
for the spatial orientation of binary orbits. The predicted broad inclination distribution 
with ≃80% of prograde binary orbits matches the observations of trans-Neptunian
binaries [6]. The formation models which imply predominantly retrograde binary orbits (e.g., [7]) 
can be ruled out. Given its applicability over a broad range of protoplanetary disk 
conditions [8], it is expected that the streaming instability seeded planetesimal formation 
also elsewhere in the solar system, and beyond.  

The streaming instability (SI) is a mechanism to seed planetesimal formation by aerodynamically 
concentrating particles to high densities [2,9-11]. SI simulations show that particle 
concentration is strong if the (local and height-integrated) solid-to-gas ratio is at least 
modestly enhanced over solar abundances [10]. Global disk 
evolution, including photoevaporation, ice lines, pressure traps and other effects (e.g., [12]) 
can readily produce the required enhancement, suggesting that the SI should commonly operate 
in protoplanetary disks to hatch planetesimals. Alternately, if the first planetesimals formed 
with maximum sizes of ∼1-10 km (by any mechanism), they can subsequently grow by accreting 
mass in mutual collisions, a gradual process known as the collisional coagulation (e.g., [13]).
Previous attempts to discriminate between different formation processes from the size distribution 
of planetesimals have been inconclusive [14,15].

We analyze a suite of vertically stratified 3D simulations of the streaming instability (SI) [16,17]. 
The simulations were performed with the ATHENA code [18], which accounts for the hydrodynamic 
flow of gas, aerodynamic forces on particles, backreaction of particles on the gas flow, and particle
self-gravity. We used the shearing box approximation with at least 512^3 gas cells, more than 
1.3×10^8 particles and appropriate boundary conditions (Methods). Each simulation was 
parametrized by the dimensionless stopping time of participating particles, τ=t_ stopΩ, 
where Ω is the Keplerian frequency, and the local particle-to-gas column density ratio, Z 
(additional parameters are discussed in Methods). We adopted τ=0.3-2, which would correspond 
to sub-cm-size pebbles in the Minimum Mass Solar Nebula (MMSN; [19]) at 45 au if the gas density 
was reduced by photoevaporation [12], and Z=0.02-0.1. Other choices of these 
parameters yield similar results [16,17] as long as the system remains in the SI regime [8].

As the time progresses in our simulations (Figure 1), dense azimuthal filaments 
form, fragment and condense into hundreds of gravitationally-bound clumps. We used an efficient tree-based 
algorithm (PLAN; Methods) to identify all clumps (Figure 1c). 
Unfortunately, the resolution in the ATHENA code does not allow us to follow the gravitational 
collapse of each clump into completion. 
Instead, we measure the total angular momentum, J, and its z-component J_z = J cosθ, 
giving the clump obliquity θ.
The total angular momentum can be compared to that of a critically 
rotating Jacobi ellipsoid: J_ c=0.39 (G M^3 r)^1/2 [20], where G, M and r are the gravitational 
constant, mass and effective radius (obtained from M with ρ=1 g 
cm^-3). We find that J/J_ c≫ 1 (this conclusion is insensitive to the choice of density), 
thus demonstrating that either most of the initial angular momentum must be removed or a typical SI 
clump cannot collapse into an isolated planetesimal.  

The vigorous rotation of the SI clumps established here is conducive to the formation of binary 
planetesimals with properties that closely match observations of the trans-Neptunian binaries.  Specifically, 
in the regime of J/J_ c≫ 1, gravitational collapse is capable of producing a ∼100% binary 
fraction [3] consistent with observations [5]. Binary planetesimals produced by gravitational collapse 
have nearly equal-size components (r_2/r_1>0.5, where r_1 and r_2 are the primary and secondary 
component radii; Figure 2) and large separations (a_ b/(r_1+r_2)>10, where a_ b is the binary 
semimajor axis), just as needed to explain observations [4-6]. Moreover, the matching colors of binary components 
[21] imply that each binary formed with a uniform compositional mix, as expected for gravitational 
collapse (but not random capture). Here we elaborate the prediction of the SI model for the spatial 
orientation of binary orbits.

The clump obliquities obtained in our SI simulations have a broad distribution (Figure 3) with 
≃80% of clumps having θ<90^∘ (prograde rotation relative to the heliocentric orbit) 
and ≃20% having θ>90^∘ (retrograde rotation). This result is relatively insensitive 
to various SI parameters (e.g., τ and Z; Methods and Supplementary Figures 1-4) in the regime 
of strong particle clumping that has been explored. Other concentration mechanisms – e.g., isotropic 
turbulence or secular gravitational instability – lack concrete predictions for rotation, but are not 
expected to produce the same distribution of angular momenta by chance. The obliquity distribution 
shown in Figure 3 is thus a telltale signature of the SI.

To understand how clump obliquities were established, we traced their evolution back in time.
Two stages were identified. To get insights into the initial aerodynamic stage, we computed the 
vorticity of the particle field, ∇× v_ p, where v_ p is the particle 
velocity. We found that the colatitude distribution of the vorticity vectors of dense particle 
clumps is broad (Supplementary Figure 4). This is a consequence of the SI that acts to produce 
the vertical motions needed to tilt the vorticity out of the disk midplane [2,10,17]. The 
preference for prograde rotation is established during the subsequent stage of gravitational 
collapse (because prograde clumps become gravitationally bound more often than the retrograde 
clumps). This preference is consistent, for instance, with gravitational accretion of pebbles 
onto protoplanets [22]. The clump obliquities change little after the initial collapse phase 
(Supplementary Figure 2). Thus the imprint of particle concentration by the SI, modified by 
gravitational collapse, is preserved.

To predict the expected orientation of binary orbits from the SI, we rely on the published results 
of gravitational collapse simulations [3], where the initial value of θ was shown to be a 
good proxy for the binary inclination, i_ b.
The clump obliquity distribution can thus be directly compared to the observed distribution 
of binary inclinations [6]. Here we focus on binaries found in the dynamically cold population 
of the classical Kuiper belt (hereafter cold classicals, or CCs, defined here as members of the 
classical main belt with heliocentric orbit inclinations i<5^∘; see ref. [23] for a 
definition of the classical main belt). The CCs are thought to have formed in-situ at a>40 au and 
survived the epoch of planetary migration relatively unharmed [24].  They probably a relatively 
pristine record of planetesimal formation. 

A comparison of the SI-model predicted binary inclinations with observations (Supplementary Table 1) 
reveals that the two distributions are indistinguishable from each other (Figure 3). Specifically, 
we implemented the Kolmogorov-Smirnov (K-S) test by comparing the cumulative distribution functions 
of observed (20) and model (over 400) binaries. The K-S test indicates that the null hypothesis 
(i.e., the two samples are drawn from the same underlying distribution) cannot be ruled out with more 
than 13% significance. This comparison clinches an argument in favor of the SI. Crucially, there is 
a marked ≃4:1 preference for prograde orbits (i_ b<90^∘). For comparison, ref. 
[7] proposed that the equal-size binaries in the Kuiper belt formed by capture during the coagulation 
growth of planetesimals. Capture in their model presumably occurred as a result of dynamical friction 
from a sea of small planetesimals (the L^2s mechanism in [7]) or via three-body encounters 
(the L^3 mechanism). These capture models can be ruled out based on the observed inclination 
distribution of binaries ([6]; Figure 3), because the L^2s mechanism predicts retrograde binary orbits 
with i_ b∼ 180^∘, whereas the L^3 mechanism implies a ∼3:2 preference for 
retrograde orbits [25].

The SI is expected to occur over a wide range of protoplanetary disk conditions and pebble sizes 
[8], which suggests that the planetesimal formation by the SI was widespread. Previous studies 
explored the SI implications for the initial size distribution (ISD) of planetesimals (e.g., [16,26]),
which can be described by a rolling power-law function with an exponential cut-off at large sizes. 
Attempts to validate the ISD on the size distribution of asteroids and Kuiper belt objects are 
obscured by secondary processes that modified the distributions after formation (e.g., sustained 
impact fragmentation, [27]). Here, we point out that the ISD expected from the SI appears to be 
broadly consistent with the rounded profile of the absolute magnitude distribution of the CCs [28]. 

The distinctive shape of the New Horizons flyby target (486958) 2014 MU_69, which belongs to the 
CC population, provides additional constraints on the planetesimal formation process.
MU_69 is a contact binary consisting of two lenticular lobes, roughly 22× 20 × 7 km 
and 14 × 14 × 10 km in size, connected by a narrow neck [29]. The New Horizons team 
have been interpreting the shape as resulting from gentle gravitational collapse [29]. The CC 
binaries apparently span the full range of component separations and sizes from widely separated 
∼100-km-class binaries to contact/small binaries such as MU_69. They are the key 
to understanding the protoplanetary disk conditions at >30 au during planetesimal formation.




[Youdin and Kenyon(2013)]2013pss3.book....1Y [1] Youdin, A. N., Kenyon, S. J. From 
disks to planets. Planets, Stars and Stellar Systems 3, T. D. Oswalt,  
L. M. French, and P. Kalas (eds.), Dordrecht, 1-62 (2013). 

[Youdin and Goodman(2005)]2005ApJ...620..459Y [2] Youdin, A. N., Goodman, J. 
Streaming instabilities in protoplanetary disks. Astrophys. J. 620, 459-469 (2005). 

[Nesvorný et al.(2010)]2010AJ....140..785N [3] Nesvorný, D., Youdin, A. N., 
Richardson, D. C. Formation of Kuiper belt binaries by gravitational collapse. 
Astron. J. 140, 785-793 (2010). 

[Noll et al.(2008)]2008ssbn.book..345N [4] Noll, K. S., Grundy, W. M., Chiang, E. I., 
Margot, J.-L., Kern, S. D. Binaries in the Kuiper belt. The Solar System Beyond 
Neptune,  M. A. Barucci, H. Boehnhardt, D. P. Cruikshank, and A. Morbidelli (eds.), Tucson, 
345-363 (2008). 

[Fraser et al.(2017)]2017NatAs...1E..88F [5] Fraser, W. C. et al.  
All planetesimals born near the Kuiper belt formed as binaries. Nature Astronomy 1, 
id. 0088 (2017).

[Grundy et al.(2018)]2017NatAs...1E..89F [6] Grundy, W. M. et al. Mutual orbit 
orientations of transneptunian binaries. Icarus, submitted (2019).

[Goldreich et al.(2002)]2002Natur.420..643G [7] Goldreich, P., Lithwick, Y., Sari, 
R. Formation of Kuiper-belt binaries by dynamical friction and three-body 
encounters. Nature 420, 643-646 (2002).

[Yang et al.(2017)]2017A   A...606A..80Y [8] Yang, C.-C., Johansen, A., Carrera, D. 
Concentrating small particles in protoplanetary disks through the streaming instability. 
Astron. Astrophys. 606, A80 (2017).



















[Johansen et al.(2007)]2007Natur.448.1022J [9] Johansen, A., Oishi, J. S., 
Mac Low, M.-M., Klahr, H., Henning, T., Youdin, A. Rapid planetesimal formation 
in turbulent circumstellar disks. Nature 448, 1022-1025 (2007). 

[Johansen et al.(2009)]2009ApJ...704L..75J [10] Johansen, A., Youdin, A., Mac Low, 
M.-M. Particle clumping and planetesimal formation depend strongly on metallicity. 
Astrophys. J. Letters 704, L75-L79 (2009).

[Bai and Stone(2010)]2010ApJ...722.1437B [11] Bai, X.-N., Stone, J. M. Dynamics 
of solids in the midplane of protoplanetary disks: implications for planetesimal formation. 
Astrophys. J. 722, 1437-1459 (2010).

[Carrera et al.(2017)]2017ApJ...839...16C [12] Carrera, D., Gorti, U., Johansen, 
A., Davies, M. B. Planetesimal formation by the streaming instability in a photoevaporating 
disk. Astrophys. J. 839, 16 (2017). 

[Kenyon and Luu(1998)]1998AJ....115.2136K [13] Kenyon, S. J., Luu, J. X.  
Accretion in the early Kuiper belt. I. Coagulation and velocity evolution. 
Astron. J. 115, 2136-2160 (1998).

[Morbidelli et al.(2009)]2009Icar..204..558M [14] Morbidelli, A., Bottke, W. F., 
Nesvorný, D., Levison, H. F. Asteroids were born big. Icarus 204, 
558-573 (2009).

[Weidenschilling(2011)]2011Icar..214..671W [15] Weidenschilling, S. J.
Initial sizes of planetesimals and accretion of the asteroids. Icarus 214, 
671-684 (2011). 

[Simon et al.(2017)]2017ApJ...847L..12S [16] Simon, J. B., Armitage, P. J., Youdin, 
A. N., Li, R. Evidence for universality in the initial planetesimal mass function. 
Astrophys. J. Letters 847, L12-L17 (2017).

[Li et al.(2018)]2018ApJ...862...14L [17] Li, R., Youdin, A. N., Simon, J. B.
On the numerical robustness of the streaming instability: particle concentration and 
gas dynamics in protoplanetary disks. Astrophys. J. 862, 14-29 (2018).

[Stone et al.(2008)]2008ApJS..178..137S [18] Stone, J. M., Gardiner, T. A., Teuben, 
P., Hawley, J. F., Simon, J. B. Athena: a new code for astrophysical MHD. 
Astrophys. J. Supplement Series 178, 137-177 (2008).

[Hayashi(1981)]1981PThPS..70...35H [19] Hayashi, C. Structure of the solar 
nebula, growth and decay of magnetic fields and effects of magnetic and turbulent viscosities 
on the nebula. Progress of Theoretical Physics Supplement 70, 35-53 (1981).

[Poincaré(1885)]1885BuAsI...2..109P [20] Poincaré, H. Mémoires et 
observations. Sur l'équilibre d'une masse fluide animée d'un mouvement de rotation. 
Bulletin Astronomique 2, 109-118 (1885).

[Benecchi et al.(2011)]2011Icar..213..693B [21] Benecchi, S. D., Noll, 
K. S., Stephens, D. C., Grundy, W. M., Rawlins, J. Optical and infrared 
colors of transneptunian objects observed with HST. Icarus 213, 693-709 (2011).

[Johansen and Lacerda(2010)]2010MNRAS.404..475J [22] Johansen, A., Lacerda, P.  
Prograde rotation of protoplanets by accretion of pebbles in a gaseous environment. 
MNRAS 404, 475-485 (2010). 

[Gladman et al.(2008)]2008ssbn.book...43G [23] Gladman, B., Marsden, B. G., Vanlaerhoven, 
C. Nomenclature in the outer solar system. The Solar System Beyond Neptune,  M. A. Barucci, 
H. Boehnhardt, D. P. Cruikshank, and A. Morbidelli (eds.), Tucson, 43-57 (2008).

[Parker and Kavelaars(2010)]2010ApJ...722L.204P [24] Parker, A. H., Kavelaars, J. J. 
Destruction of binary minor planets during Neptune scattering. Astrophys. J. Letters 
722, L204-L208 (2010).

[Schlichting and Sari(2008)]2008ApJ...686..741S [25] Schlichting, H. E., Sari, R. 
The ratio of retrograde to prograde orbits: a test for Kuiper belt binary formation theories. 
Astrophys. J. 686, 741-747 (2008).

[Johansen et al.(2015)]2015SciA....1E0109J [26] Johansen, A., Mac Low, M.-M., Lacerda, P., 
Bizzarro, M. Growth of asteroids, planetary embryos, and Kuiper belt objects by chondrule 
accretion. Science Advances 1, 1500109 (2015).

[Bottke et al.(2005)]2005Icar..175..111B [27] Bottke, W. F. et al. The fossilized size 
distribution of the main asteroid belt. Icarus 175, 111-140 (2005).




[Petit et al.(2016)]2016DPS....4812016P [28] Petit, J.-M. et al. The absolute magnitude 
distribution of cold classical Kuiper belt objects. AAS/Division for Planetary Sciences 
Meeting Abstracts 48, 120.16 (2016).

[Stern et al.(2019)]2019Science [29] Stern, A. et al. Initial results from the 
New Horizons exploration of 2014 MU_69, a small Kuiper belt object. Science,
in press (2019).

[Shannon and Dawson(2018)]2018MNRAS.480.1870S [30] Shannon, A., Dawson, R. 
Limits on the number of primordial Scattered disc objects at Pluto mass and higher 
from the absence of their dynamical signatures on the present-day trans-Neptunian 
populations. MNRAS 480, 1870-1882 (2018). 









Corresponding author

David Nesvorný

Southwest Research Institute

1050 Walnut St., Suite 300

Boulder, Colorado 80302

Phone: (303) 546-0023

Email: davidn@boulder.swri.edu     

Acknowledgments

The work of D.N. was funded by the NASA Emerging Worlds program. R.L. acknowledges support from NASA
grant NNX16AP53H. A.N.Y. acknowledges support from NASA through grant NNX17AK59G and the NSF through 
grant 1616929. The funding sources of J.B.S. are NASA grants NNX13AI58G, NNX16AB42G, 80NSSC18K0640, 
and 80NSSC18K0597. W.M.G contribution was supported in part by NASA Keck PI Data Awards, administrated 
by the NASA Exoplanet Science Institute and in part by data analysis grants from the Space Telescope 
Science Institute (STScI), operated by the Association of Universities for Research in Astronomy, Inc., 
(AURA) under NASA contract NAS 5-26555.

Competing Interest Statement

The authors declare no competing interests.

Author contributions

D.N. suggested a comparison of clump obliquities with binary orbit inclinations and prepared the
manuscript for publication. R.L. ran one of the ATHENA simulations and performed data analyses 
with PLAN. A.N.Y developed scaling relations for planetesimal mass estimates. J.B.S. ran 
two of the ATHENA simulations. W.M.G. provided the data on trans-Neptunian binaries. 
All authors contributed to the interpretation of the results and writing of this paper.

Author information

Reprints and permissions information is available at www.nature.com/reprints.

Correspondence and requests for materials should be addressed to davidn@boulder.swri.edu.

Methods

ATHENA code 
 
Our numerical simulations use the Athena code. 
Athena is a second-order accurate, flux-conservative Godunov code for solving (in this 
particular application) the coupled equations of hydrodynamics and particle dynamics.  
To properly integrate the gas dynamics, we have employed the dimensionally unsplit corner 
transport upwind method [31] and the piecewise parabolic method [32]
to spatially reconstruct gas quantities in a third order accurate manner. The calculation of
numerical fluxes is done via the HLLC Riemann solver 
[33]. An in-depth description of the base Athena algorithm, along with tests of this 
algorithm, can be found in ref. [18].

Particles in the Athena code are treated via the super-particle approach, in which each 
super-particle is a statistical representation of a larger swarm of particles. 
The super-particles' (hereafter just `particles' for brevity) equations of motion 
are integrated via a semi-implicit drift-kick-drift method (see ref. [34]), 
and the momentum exchange between the particles and gas is handled via the triangular 
shaped cloud (TSC) scheme [35,36]. This approach maps the particle 
momentum from the particle locations to the grid-cell centers (where hydrodynamic 
quantities are handled) and conversely the gas velocity to the particle locations.  
A detailed description and tests of the particle integration algorithm can be 
found in ref. [34]. 

Shearing box
 
All of our simulations are carried out in the local, shearing box approximation.  
This approximation treats a co-rotating portion (or “patch") of an accretion disk of 
length scale L ≪ R, where R is the distance from the central star. 
In this limit, the relevant equations can be expanded into a Cartesian basis (x,y,z), 
which is defined relative to disk cylindrical coordinates as x=(R-R_0), y=R_0 ϕ, 
and z = z, where R_0 is a reference distance. 
Furthermore, appropriate source terms must be added to the equations to 
account for the non-inertial reference frame of the domain. See ref. [37] for a detailed 
description of the shearing box and refs. [34,38] for a description of the shearing box 
within the context of particle-gas calculations. 

The co-rotation frequency of the shearing box equals the Keplerian frequency at R_0, 
which we define as Ω.  In our particular setup, we include the vertical component 
of the star's gravity (i.e., the simulations are vertically stratified). 
All of our simulations also employ the orbital advection (or FARGO) algorithm 
[39,40] to analytically integrate quantities along the Keplerian shear flow and integrate 
the perturbed (relative to Keplerian) quantities using the algorithms described above.  
This is done for both the gas and the particles and improves the speed and accuracy of 
the code.  To preserve epicyclic energy to machine precision, Crank-Nicholson differencing 
is used to integrate the non-inertial, shearing box source terms.  Finally, we assume 
an isothermal equation of state for simplicity. 

Boundary conditions
 
The boundary conditions of our domain are shearing-periodic in the radial direction, 
purely periodic in azimuth, and a modified open/outflow in the vertical 
direction.  The radial shearing periodic boundaries are standard in the shearing box 
set-up [37]. Briefly, particles/gas that exit one radial boundary enter at the 
opposite boundary but with a displacement in azimuth and a change in angular momentum 
and energy applied to account for the difference in the orbital position.  The modified 
outflow conditions in the vertical direction consists of exponentially extrapolating 
the gas density into the ghost zones, as this significantly reduces numerical artifacts
[17,38]. Initially, these boundary conditions allow for excellent maintenance of 
hydrostatic equilibrium.  However, as the system evolves, gas will necessarily be lost 
via these vertical boundaries (particularly since our domain is quite small and thus 
the gas density at the vertical boundaries is not significantly different from that 
at the mid-plane).  To ensure mass conservation, we 
renormalize the gas density in every cell at every time step by a factor that ensures 
total mass conservation.

Crucial to driving the streaming instability is the radial pressure gradient of the 
gas.  However, a global pressure gradient is inconsistent with the radial shearing-periodic 
boundaries.  To circumvent this issue, an inward radial force is applied to the particles.  
This force creates an inward drift equivalent to that which would occur as particles 
lose angular momentum to the sub-Keplerian gas in a real disk.  This method is widely 
employed in shearing box calculations of particle-gas mixtures (e.g., [34]) 
and is well validated. 

Particle self-gravity
 
In order to study the formation of planetesimals, the mutual gravitational attraction 
between particles must be included.  When activated (see below), the gravitational force 
is added to the particle equations of motion, and this force is calculated from Poisson's 
equation.  Numerically, Poisson's equation is solved as follows.  The particle mass density 
is mapped to the grid cell centers using the TSC approach.  
This density is then mapped to the nearest time at which 
the radial boundaries are purely periodic (see ref. [37] for a discussion of this mapping 
and these “periodic points"), after which a 3D FFT is employed to transform Poisson's 
equation into Fourier space.  

The gravitational potential is then solved in Fourier space, after which it is transformed 
back to real space using a second 3D FFT and then mapped back to the original non-periodic 
frame.  As our vertical boundaries are open, a Green's function approach is used to solve 
Poisson's equation in the vertical direction (see details in [38,41]). Finally,
the gravitational force is calculated by a second-order, central finite difference of 
the gravitational potential and then mapped from the grid cell centers back to the locations 
of the particles via TSC.  Boundary conditions must also be applied to the gravitational 
potential, and these are essentially the same as for the gas quantities: shearing-periodic 
in x, purely periodic in y, and open in z with the potential extrapolated into the 
ghost zones via a third-order scheme. The algorithm has been tested rigorously in ref. [38]. 

Initial conditions
 
All of our simulations have the following initial conditions. The gas is in
hydrostatic balance with vertical gravity, resulting in an initially Gaussian profile

    ρ_g = ρ_0  exp(-z^2/2H^2) ,
 
where ρ_0 is the mid-plane gas density, and H is the vertical scale height 
of the gas. We set the standard gas parameters equal to unity; ρ_0 = H = 
Ω =  = 1 ( is the sound speed).

The particles are initially distributed with a Gaussian vertical profile, with particle 
scale height H_p = 0.02H, whereas in x and y, the particles are uniformly distributed.  
Initially, all perturbed velocities (gas and particle velocities with the Keplerian shear 
subtracted) are set to zero. Finally, we introduce random noise (sampled from a uniform 
distribution) to the particle locations in all three dimensions to seed the streaming 
instability. 

We first run our numerical simulations without particle self-gravity, 
allowing the streaming instability to fully develop and produce particle clumping.  After 
the streaming instability has saturated (which is generally defined by when the maximum 
particle density is stochastically varying but statistically constant in time), we switch 
on self-gravity.  This method is done largely for numerical convenience, as particle 
self-gravity imposes a performance hit on the integration.  In other works [38,42], we 
have verified that the start time of self-gravity does not largely influence the final 
properties of planetesimals. 

Simulation parameters
 
In general, streaming-induced planetesimal formation calculations can be characterized 
via four dimensionless quantities [16,42]: the dimensionless stopping time,

    τ = t_ stopΩ ,

which is the ratio of the dimensional stopping time t_ stop (i.e., the timescale 
over which particle momentum relative to the gas motion decreases by 1/e) to the dynamical time 
Ω^-1 and characterizes the gas-particle interaction (for a single particle species 
in this case), the particle concentration,

    Z = Σ_p/Σ_g ,
  
which is the ratio of the particle mass surface density Σ_p to the gas surface 
density Σ_g (note that this parameter is related to, but not the same as 
'metallicity'), and a radial pressure gradient parameter that accounts for the 
irradiation from the central star and thus the sub-Keplerian rotating gas in real disks,

    Π = η v_ K/ ,

where v_ K is the Keplerian velocity and  η is related to the degree to which 
gas orbital velocities are sub-Keplerian: 

v_ϕ = v_ K (1-η) in the shearing box approximation.

When particle self-gravity is activated, an additional parameter comes into play, which
is the relative strength of self-gravity to tidal shear.  We quantify this relative 
strength as

    G̃ = 4π Gρ_0/Ω^2.

This parameter can be related to the Toomre Q via G̃ = (4/√(2π)) Q^-1. 
Higher values of G̃ (lower values of Q), which occur at larger distances 
from the central star, equate to weaker tidal shear relative to gravitational binding. 

Simulation setup
 
We have run three simulations in total that span a range of τ, Z and resolution. Run B22 
has τ = 2, Z = 0.1, and a domain size L_x× L_y× L_z = (0.2H)^3 
with resolution N_x× N_y× N_z = 512^3. Run C203 has τ = 0.3, 
Z = 0.02, and a domain size L_x× L_y× L_z = (0.2H)^3 with 
resolution N_x× N_y× N_z = 512^3. Finally, simulation A12 has τ = 2, 
Z = 0.1, and a domain size L_x× L_y× L_z = 0.1H×0.1H×0.2H with resolution 
N_x× N_y× N_z = 512×512×1024. All simulations 
have G̃ = 0.05, which equates to a Toomre Q ≃ 32, Π = 0.05, and are initiated 
with N_ par = 1.34×10^8 (run A12) or 1.54×10^8 particles (B22 and C203). 

The relation between particle size and τ depends on uncertain properties of the gas disk.
If, for reference, we adopt values of the MMSN from ref. [43], a cm-size particle 
at 45 au would have τ = 0.3. Ref. [12] argued that CCs formed when the 
gas density at 45 au was reduced by photoevaporation. If, for example, the gas density was 
reduced ten fold, the range of τ adopted here would correspond to sub-cm-size particles. 
Alternately, ref. [42] considered a gas disk model with surface density Σ(r) ∝ 1/r, 
outer radius at 500 au, and total mass M_ disk=0.05 M_⊙, where M_⊙ is the solar mass. 
This disk model would give τ=0.3 for a 2.5 mm particle and τ=2 for a 1.7 cm 
particle (both at 45 au).

Clump identification
 
To identify and further characterize the properties of planetesimals produced in our simulations, 
we use a newly developed clump-finding tool, PLanetesimal ANalyzer (, [44]).  
It is designed to work with the 3D particle output of Athena and find gravitationally-bound 
clumps robustly and efficiently.  , which is written in 
C⁠.4ex+⁠.4ex+ 
with OpenMP/MPI, is massively parallelized to analyze billions of particles and many 
snapshots simultaneously.  

Briefly, the workflow of our clump-finding algorithm in  is as follows.  The 
approach is based on the halo finder  [45], which allows for fast grouping of 
physically related particles.   first builds a memory-efficient linear Barnes-Hut 
tree representing all the particles [46,47].  Each particle is then 
assigned a density computed from the nearest N_ den particles (N_ den = 64 by default).  
For particles with densities higher than a threshold, δ_ outer = 8ρ_0/G̃, 
 chains them up towards their densest neighbors repetitively until a density 
peak is reached.  All the particle chains linked to the same density peak are combined to 
create a group.  

 then merges those groups by examining their boundaries to construct a list of 
bound clumps.  Based on the total kinematic and gravitational energies, deeply intersected 
groups are merged if bound.  However, two particle groups with a saddle point less dense 
than δ_ saddle=2.5δ_ outer remain separated [45].  Next,  goes 
through each group — or raw clump — to unbind any contamination (i.e., passing-by and 
not bound) particles and gather possibly unidentified member particles within its Hill sphere.  
After discarding those clumps with Hill radii smaller than one hydrodynamic grid cell (Δ x) 
or density peaks less than δ_ peak=3δ_ outer,  outputs the 
final list of clumps with their physical properties derived from particles.

Unlike ,  does not further diagnose sub-structures within clumps.  
Most clumps in our simulations are highly-concentrated, where particles often collapse into 
regions much smaller than the Hill radius but comparable to or smaller than Δ x, 
invalidating the search for sub-cell structures.  The exact internal architecture of 
self-bound clumps is beyond the scope of this work.  Nevertheless,  does a robust 
job of identifying clumps as seen in this work (Figure 1) and other recent work [42].  

Clump properties
 
In order to determine the inclination of bound clumps, we use  to calculate 
their angular momenta in the inertial frame by summing up the contribution from each particle 
[48,49].  Assuming that the i-th particle of a clump is located at (x_i, y_i, z_i) relative 
to the center-of-mass of the clump, the inertial frame angular momentum J = (J_x, J_y, J_z) 
can be written as

    J_x    = ∑_i m_i (y_i ż_̇i̇ - z_i ẏ_̇i̇ - Ω x_i z_i), 
    
    J_y    = ∑_i m_i (z_i ẋ_̇i̇ - x_i ż_̇i̇ - Ω y_i z_i), 
    
    J_z    = ∑_i m_i [x_i ẏ_̇i̇ - y_i ẋ_̇i̇ - Ω (x_i^2 + y_i^2)],

where m_i is the particle mass. The obliquity of such a clump (i.e., the colatitude of 
J) is then calculated from 

    θ = ⟨J, ẑ⟩ = cos^-1(J·ẑ/|J|) .




Code availability

The ATHENA code is available on GitHub (https://github.com/PrincetonUniversity/Athena-Cversion). 
The PLAN code is available on Zenodo (see ref. [44]).

Data availability

The data that support the plots within this paper and other findings of this study are available 
from the corresponding author upon reasonable request.


natexlab#1#1




[Colella(1990)]colella90 [31] Colella, P. Multidimensional upwind methods 
for hyperbolic conservation laws, J. of Computational Physics 87, 171-200 (1990).

[Colella and Woodward(1984)]1984JCoPh..54..174C [32] Colella, P., Woodward, P. R. 
The piecewise parabolic method (PPM) for gas-dynamical simulations. J. of 
Computational Physics 54, 174-201 (1984).

[Toro(1999)]toro99 [33] Toro, E. F. Riemann Solvers and Numerical Models 
for Fluid Dynamics, Springer-Verlag: Berlin, Heidelberg (1999).

[Bai and Stone(2010)]2010ApJS..190..297B [34] Bai, X.-N., Stone, J. M. 
Particle-gas dynamics with Athena: method and convergence. Astrophys. J. 
Supplement Series 190, 297-310 (2010).

[Hockney and Eastwood(1981)]1981csup.book.....H [35] Hockney, R. W., Eastwood, 
J. W. Computer Simulation Using Particles, New York: McGraw-Hill (1981). 

[Youdin and Johansen(2007)]2007ApJ...662..613Y [36] Youdin, A., Johansen, A. 
Protoplanetary disk turbulence driven by the streaming instability: linear evolution 
and numerical methods. Astrophys. J. 662, 613-626 (2007).

[Hawley et al.(1995)]1995ApJ...440..742H [37] Hawley, J. F., Gammie, C. F., Balbus, 
S. A. Local three-dimensional magnetohydrodynamic simulations of accretion disks. 
Astrophys. J. 440, 742-763 (1995).

[Simon et al.(2016)]2016ApJ...822...55S [38] Simon, J. B., Armitage, P. J., Li, R., 
Youdin, A. N. The mass and size distribution of planetesimals formed by the streaming 
instability. I. The role of self-gravity. Astrophys. J. 822, 55-72 (2016). 

[Masset(2000)]2000A   AS..141..165M [39] Masset, F. FARGO: A fast eulerian 
transport algorithm for differentially rotating disks. Astron. Astrophys. Supplement 
Series 141, 165-173 (2000).

[Stone and Gardiner(2010)]2010ApJS..189..142S [40] Stone, J. M., Gardiner, T. A. 
Implementation of the shearing box approximation in Athena. Astrophys. J. Supplement 
Series 189, 142-155 (2010). 

[Koyama and Ostriker(2009)]2009ApJ...693.1346K [41] Koyama, H., Ostriker, E. C. 
Pressure relations and vertical equilibrium in the turbulent, multiphase interstellar medium. 
Astrophys. J. 693, 1346-1359 (2009).

[Abod18]koyama09 [42] Abod, C. P. et al., The mass and size distribution of planetesimals 
formed by the streaming instability. II. The effect of the radial gas pressure gradient. 
ArXiv e-prints arXiv:1810.10018 (2018).

[Chiang and Youdin(2010)]2010AREPS..38..493C [43] Chiang, E., Youdin, A. N. 
Forming planetesimals in solar and extrasolar nebulae. 
AREPS 38, 493-522 (2010).

[Li(2018)]PLAN [44]
Li, R., PLAN: PLanetesimal ANalyzer v0.2 (2018),

<https://doi.org/10.5281/zenodo.1436807>

[Eisenstein and Hut(1998)]1998ApJ...498..137E [45] Eisenstein, D. J., Hut, P. HOP: 
A new group-finding algorithm for N-body simulations. Astrophys. J. 
498, 137-142 (1998).

[Morton(1966)]Morton1966 [46]
Morton, G. M., A Computer Oriented Geodetic Data Base and a New
  Technique in File Sequencing, International Business Machines Co. (1966).

[Barnes and Hut(1986)]1986Natur.324..446B [47] Barnes, J., Hut, P.
A hierarchical O(N log N) force-calculation algorithm. Nature 324, 446-449 (1986). 

[Lissauer and Kary(1991)]1991Icar...94..126L [48] Lissauer, J. J., Kary, D. M. 
The origin of the systematic component of planetary rotation. I - Planet on a circular orbit. 
Icarus 94, 126-159 (1991).

[Dones and Tremaine(1993)]1993Icar..103...67D [49] Dones, L., Tremaine, S.
On the origin of planetary spins. Icarus 103, 67-92 (1993).









