Trans-Neptunian Binaries as Evidence for Planetesimal Formation by the Streaming Instability

TL;DR

This paper provides evidence that Kuiper belt planetesimals formed via the streaming instability, supported by simulations matching observed binary orbit inclinations, thus offering insights into planetesimal formation processes.

Contribution

It demonstrates that the streaming instability can explain the formation and orbital characteristics of trans-Neptunian binaries, ruling out retrograde formation models.

Findings

80% of binaries have prograde orbits

Simulations match observed inclination distribution

Retrograde models are inconsistent with observations

Abstract

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. 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. Gravitational collapse was previously suggested to explain the ubiquity of equal-size binaries in the Kuiper belt. 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. The formation models which imply predominantly retrograde binary orbits can be ruled out. Given its applicability over a broad range of protoplanetary disk conditions, it is expected that the streaming instability seeded planetesimal formation also elsewhere in the solar system, and beyond.

Equations

$${\rho }_g={\rho }_0\exp \left(\frac{-z^2}{2H^2}\right),$$

$$\tau =t_{\mathrm{stop}}\Omega ,$$

$$Z=\frac{{\Sigma }_p}{{\Sigma }_g},$$

$$\Pi =\eta v_{\mathrm{K}}/c_{\mathrm{s}},$$

$$\tilde{G}=\frac{4\pi G{\rho }_0}{{\Omega }^2}.$$

$$J_x$$

$$J_y$$

$$J_z$$

$$\theta=\langle\mbox{\boldmath{$J$}},\hat{z}\rangle=\cos^{-1}\left(\frac{\mbox{\boldmath{$J$}}\cdot\hat{z}}{|\mbox{\boldmath{$J$}}|}\right)\ .$$

Full text

Trans-Neptunian Binaries as Evidence for Planetesimal Formation by the Streaming Instability

David Nesvorný1,∗, Rixin Li2, Andrew N. Youdin2, Jacob B. Simon1,3, William M. Grundy4

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

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 $\sim$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 $\simeq$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 $\sim$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 5123 gas cells, more than $1.3\times 10^{8}$ particles and appropriate boundary conditions (Methods). Each simulation was parametrized by the dimensionless stopping time of participating particles, $\tau=t_{\rm stop}\Omega$, where $\Omega$ is the Keplerian frequency, and the local particle-to-gas column density ratio, $Z$ (additional parameters are discussed in Methods). We adopted $\tau=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\theta$, giving the clump obliquity $\theta$. The total angular momentum can be compared to that of a critically rotating Jacobi ellipsoid: $J_{\rm c}=0.39(GM^{3}r)^{1/2}$ [20], where $G$, $M$ and $r$ are the gravitational constant, mass and effective radius (obtained from $M$ with $\rho=1$ g cm*-3*). We find that $J/J_{\rm c}\gg 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_{\rm c}\gg 1$, gravitational collapse is capable of producing a $\sim$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_{\rm b}/(r_{1}+r_{2})>10$, where $a_{\rm 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 $\simeq$80% of clumps having $\theta<90^{\circ}$ (prograde rotation relative to the heliocentric orbit) and $\simeq$20% having $\theta>90^{\circ}$ (retrograde rotation). This result is relatively insensitive to various SI parameters (e.g., $\tau$ 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, $\nabla\times v_{\rm p}$, where $v_{\rm 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 $\theta$ was shown to be a good proxy for the binary inclination, $i_{\rm 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^{\circ}$; 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 $\simeq$4:1 preference for prograde orbits ($i_{\rm b}<90^{\circ}$). 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 L2s mechanism in [7]) or via three-body encounters (the L3 mechanism). These capture models can be ruled out based on the observed inclination distribution of binaries ([6]; Figure 3), because the L2s mechanism predicts retrograde binary orbits with $i_{\rm b}\sim 180^{\circ}$, whereas the L3 mechanism implies a $\sim$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 MU69, which belongs to the CC population, provides additional constraints on the planetesimal formation process. MU69 is a contact binary consisting of two lenticular lobes, roughly $22\times 20\times 7$ km and $14\times 14\times 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 $\sim$100-km-class binaries to contact/small binaries such as MU69. They are the key to understanding the protoplanetary disk conditions at $>$30 au during planetesimal formation.