FORMATION OF KUIPER BELT BINARIES BY GRAVITATIONAL COLLAPSE

The existence of binary Kuiper Belt Objects (KBOs) other than Pluto-Charon (Christy & Harrington 1978) has been suspected since the discovery of the Kuiper Belt (KB; Jewitt & Luu 1993), but it was not until 2000 December that the first binary KBO, 1998 WW₃₁, was detected by direct ground-based imaging (Veillet et al. 2001, 2002). Recent observations indicate that ∼30% of 100 km class classical cold KBOs with orbital inclinations i < 5° are binaries (Noll et al. 2008a, 2008b; >0.06 arcsec separation, <2 mag contrast). Binaries with larger primaries, large magnitude differences, and smaller separations may be even more common (Brown et al. 2006; Weaver et al. 2006) and probably require a different formation mechanism (e.g., Canup 2005).

The properties of known binary KBOs differ markedly from those of the main-belt and near-Earth asteroid binaries (Merline et al. 2002; Noll et al. 2008a). The 100 km class binary KBOs identified so far are widely separated and their components are similar in size. These properties defy standard ideas about processes of binary formation involving collisional and rotational disruption, debris re-accretion, and tidal evolution of satellite orbits (Stevenson et al. 1986). They suggest that most binary KBOs may be remnants from the earliest days of the solar system. If so, we can study them to learn about the physical conditions that existed in the trans-Neptunian disk when large KBOs formed.

The KB (Kuiper 1951; Jewitt & Luu 1993) provides an important constraint on planet formation. To explain its present structure, including the large binary fraction among the classical cold KBOs discussed above, we must show how the 100 km size and larger bodies accreted from smaller constituents of the primordial trans-Neptunian disk. Two main possibilities exist: (1) hierarchical coagulation (hereafter HC), where two-body collisions between objects in a dynamically cold planetesimal disk lead to objects' accretion and growth and (2) gravitational instability (hereafter GI), where the gas-particle effects and/or gravitational instabilities produce concentrations of gravitationally bound solids followed by their rapid collapse into large objects. We briefly comment on these theories below.

As for HC, Stern (1996), Stern & Colwell (1997), Kenyon & Luu (1998, 1999), and Kenyon (2002) conducted simulations of the primordial "bottom-up" process involving collisional accumulation of small KBOs into larger ones (also see Kenyon et al. 2008 for a review). They found that two competing physical processes, growth by mergers and erosion by fragmentation, determine the final result. According to these studies, the observed KBOs can only form by HC in ≲10⁸ yr if (1) the orbits in the belt were initially much more circular and planar than they are now (e ∼ i ∼ 10⁻⁴–10⁻² compared to present eccentricities e ∼ 0.1 and inclinations i ∼ 10°) and (2) the initial disk mass was ∼100–1000 times larger than the current KB mass, M_KB ∼ 0.01–0.1 M_Earth (Trujillo et al. 2001a, 2001b; Gladman et al. 2001; Bernstein et al. 2004; Fraser et al. 2008).

The GI hypothesis has been advanced by recent breakthroughs in theory and simulation (see Chiang & Youdin 2010 for a review). The classical GI of a particle-rich nebula mid-plane (Safronov 1969; Goldreich & Ward 1973; Youdin & Shu 2002) can be prevented by even a modest amount of stirring from a turbulent gas disk (Weidenschilling 1980; Cuzzi et al. 1993). However, particles can also clump in a turbulent flow (e.g., Cuzzi et al. 2001, 2008; Johansen et al. 2006). The streaming instability (Youdin & Goodman 2005) is a powerful concentration mechanism by which weak particle clumps perturb the gas flow in a way that increases their amplitude (Youdin & Johansen 2007; Johansen & Youdin 2007). Simulations of rocks in a gas disk find that streaming instability-induced clumping produces gravitationally bound clusters of solids, either with (Johansen et al. 2007) or without (Johansen et al. 2009) large-scale magnetohydrodynamic turbulence. These clumps exceed the mass of compact 100 km radius planetesimals. The local disk metallicity (relevant for the amount of condensed solids) needs to slightly exceed solar abundances in order to counter turbulent stirring and trigger strong clumping (Youdin & Shu 2002; Johansen et al. 2009). Much work remains to determine the relative roles of GI and HC in the solar system and beyond.

1.1. Binary Formation in Hierarchical Coagulation

1.2. New Model for Binary Formation in Gravitational Instability

Benecchi et al. (2009) reported resolved photometric observations of the primary and secondary components of 23 binary KBOs. They found that the primary and secondary components of each observed binary pair have identical colors to within the measurement uncertainties. On the other hand, the wide color range of binary KBOs as a group is apparently indistinguishable from that of the population of single KBOs. These results can be difficult to understand in (some of) the models of binary KBO formation discussed in Section 1.1. Instead, the most natural explanation is that binary KBOs represent snapshots of the local composition mix in a nebula with important temporal and/or spatial gradients.⁵

The observed color distribution of binary KBOs can be easily understood if KBOs formed by GI. The common element invoked by various GI models is the final stage of gravitational collapse when the gravitationally bound pebbles and boulders are brought together, collide, and eventually accrete into large objects. We envision a situation in which the excess of angular momentum in a gravitationally collapsing swarm prevents formation of a solitary object. Instead, a binary with large specific angular momentum forms from local solids, implying identical composition (and colors) of the binary components. Moreover, binaries with similarly sized components are expected to form in this model because similar components maximize the use of the collapsing cloud's angular momentum (Figure 1; Nesvorný 2008).⁶

Zoom In Zoom Out Reset image size

Our model for binary KBO formation is similar to that of binary stars from the collapse of a rotating molecular cloud core (Kratter et al. 2008), and more specifically to binary star formation in fragmenting disks around black holes (e.g., Alexander et al. 2008). It has not been studied in the context of planetary science. For example, while Johansen et al. (2007, 2009) investigated the formation of gravitationally bound concentrations of solids, they did not follow the final stage of gravitational collapse in detail because the spatial resolution of their code was limited by the need to resolve much larger scales of disk dynamics. In a similar context, Tanga et al. (2006) investigated the implications of the gravitational reaccumulation process for the evolution of asteroid shapes.

Here, we conduct N-body numerical simulations of a gravitationally collapsing segment of disk solids to determine whether the observed 100 km class binary KBO could have formed in the GI model. We attempt to "reverse engineer" the conditions that give rise to binary formation by varying the initial set of parameters. This is because precise initial conditions in a bound clump are uncertain due to the complex physics of particles in turbulent accretion disks. We do not attempt to extract our initial data from the Johansen et al. (2007, 2009) simulations because they have low resolution (<10 grid cells) across the densest clumps. We describe our integration method and setup in Section 2. The results are presented in Section 3 and discussed in Section 4.

Our simulations of gravitational collapse were performed with a modified version of the N-body cosmological code PKDGRAV (Stadel 2001), described in Richardson et al. 2000 (also see Leinhardt et al. 2000; Leinhardt & Richardson 2002). PKDGRAV is a scalable, parallel tree code that is the fastest code available to us for this type of simulation. A unique feature of the code is the ability to rapidly detect and treat collisions between particles. We used N = 10⁵ particles per run. Each PKDGRAV particle was given initial mass M = M_tot/N, where M_tot was the assumed total mass of the gravitationally unstable swarm. Initially, the PKDGRAV particles were distributed in a spherical volume with radius R_tot < R_Hill = (GM_tot/3Ω²_Kep)^1/3, in which self-gravity dominates (G is the gravitational constant).

The initial velocities of PKDGRAV particles were set to model the collapse phase that occurs after some GI. Since the exact GI conditions are uncertain due to the modest resolution and uncertainties in the existing instability calculations, we sampled around a range of the initial velocities to see how different assumptions would influence the results. Specifically, we gave the swarm uniform rotation with several different values of Ω ≲ Ω_circ, where Ω_circ = V_circ/R_tot, and $V_{\rm circ}=\sqrt{GM_{\rm tot}/R_{\rm tot}}$ is the speed of a particle in a circular orbit about the cloud at R_tot. In addition, particles were also given random velocities with characteristic speed V_rand < V_circ.

The Keplerian shear was included in the Hill approximation as in Tanga et al. (2004, except that no periodic boundaries were imposed). We also conducted experiments where the Sun was directly included in the simulations as a massive PKDGRAV particle. The results obtained with these two methods were similar. Since $\Omega _{\rm circ}/\Omega _{\rm Kep} = \sqrt{3} (R_{\rm Hill}/R_{\rm tot})^{3/2}$, the shearing effects quickly diminish for R_tot < R_Hill, because the cloud is initially compact and collapses in a fraction of the orbital period.

Given the exploratory nature of our investigation, we neglected certain physical ingredients that should be less significant, but could be added to the next generation of models. Specifically, gas drag was ignored because our estimates show that the effects of gas drag should be small relative to collisional damping inside the gravitationally bound clump (see the Appendix). In addition, "mass loading" (see, e.g., Hogan & Cuzzi 2007) damps turbulence inside dense particle clumps, making it safe to ignore the forcing of particle motions by turbulent gas. The evolution of particle speeds in our simulations is set by gravitational interactions and physical collisions, the dominant effects during the final stage of collapse.

We ignored collisional fragmentation of bodies in the collapsing swarm because the expected impact speeds are low (see Section 3) and we can develop a better understanding of the collapse process with simple models. Note that debris produced by disruptive collisions between bodies in the collapsing swarm are gravitationally bound so that even if fragmentations occur, fragments are not lost. The fragmentation can be included in the next-generation models using scaling laws developed for low-speed collisions between icy bodies (e.g., Leinhardt & Stewart 2009; Stewart & Leinhardt 2009), even though it can be challenging to deal with the full complexity of the collisional cascade.

We divided the integrations into two suites. In the first suite of our "core" simulations, we used a simple physical model of collapse and covered a regular grid of parameter values in M_tot, Ω, and r. Specifically, we used Ω = 0.5, 0.75, 1.0, and 1.25Ω_circ and R_eq = 100, 250, and 750 km, where R_eq is the equivalent radius of a sphere with mass M_tot and ρ = 1 g cm⁻³. The collisions between PKDGRAV particles were treated as ideal mergers.⁷ Also, we used R_tot = 0.6R_Hill and V_rand = 0 in the core runs. The initial rotation vector of the swarm was set to be aligned with the normal to its heliocentric Keplerian orbit.

The initial radius of PKDGRAV particles, R, was set as R = fr, where r is the starting boulder size and f is an inflation factor used here to compensate for the fact that the number of PKDGRAV particles in the simulation is much smaller than the expected number of bodies in the collapsing swarm. Several possible choices of f exist. If PKDGRAV particles are required to mimic the actual collision rate of radius r boulders (case A), then f² = n/N = (R_eq/r)³/N with the initial number of boulders, n, being set by the mass constraint. This choice poses problems during the late simulation stages, however, because f = 3 × 10⁶ with R_eq = 250 km and r = 25 cm. Thus, if f remains constant during the simulation, and a fraction of PKDGRAV particles accrete into a body with mass equivalent to, say, a 50 km radius KBO, the corresponding PKDGRAV particle would have radius R = 50f^1/3 ≈ 7200 km! This is obviously bad because the separation of components in many known binary systems is <10,000 km (Noll et al. 2008a).

A different choice of f would be to use R = R* = R_eq/N^1/3 (case B), which for the above used example implies the initial radius R = 5.4 km and ρ = 1 g cm⁻³ of PKDGRAV particles. This setup severely underestimates the rate of collisions in the collapsing swarm of real submeter boulders, but has the advantage that the late stages of accretion of large objects are treated more realistically, because the corresponding PKDGRAV particles have adequate radii and bulk densities.

We conducted simulations with the two extreme setups A and B discussed above, and also for several intermediate cases. We define these cases by the initial ratio f* = R/R*, where f* = 1 corresponds to case B and f* = (n/N)^1/6 to case A. The intermediate cases with 1 < f* < (n/N)^1/6 are probably more realistic than the two extreme cases. They conservatively use lower-than-realistic collision rates and do not allow the large objects to grow beyond reasonable limits in radius. Specifically, we used f* = 1, 3, 10, 30, and 100.

Thus, with three values of R_eq, four values of Ω, and five values of f*, we have 60 different initial states of the swarm. Four simulations were performed for each state where different random seeds were used to generate the initial positions of PKDGRAV particles in the swarm. We used a 0.3 day time step in the PKDGRAV integrator so that the expected binary orbital periods were resolved by at least ∼100 time steps. We verified that shorter time steps lead to results similar to those obtained with the 0.3 day time step. The integration time was set to T_int = 100 yr, or about 0.6P(30), where P(30) is the orbital period at 30 AU. Together, our core simulations represent 240 jobs each requiring about two weeks on one Opteron 2360 CPU. To increase the statistics in the most interesting cases, 10 simulations with different random seeds were performed for Ω = 0.75Ω_circ, f* = 10, and all R_eq values.

Our second suite of simulations includes a diverse set of jobs in which we tested a broader range of parameters, extended selected integrations over several orbital periods at 30 AU, used different R_tot and V_rand values, included effects of inelastic bouncing of PKDGRAV particles, imposed retrograde rotation of the initial swarm, etc. We describe the results of these simulations in Section 3.

While our core simulations with f*>30 produce massive bodies that are frequently bound in binary systems, the binary separations tend to be very large because the inflated PKDGRAV particles prevent formation of tight binaries. On the other hand, the simulations with f* < 3 show low collision rates and do not produce massive objects in 100 yr. Moreover, as expected, simulations with Ω>Ω_circ lead to the swarm's dispersal due to excess angular momentum. We therefore first discuss the results obtained with intermediate values of f*, which are probably the most realistic ones, and Ω ⩽ Ω_circ. All binary systems produced in these simulations were followed for 10,000 yr to check on their stability and orbital behavior.

The binary systems that form in T_int = 100 yr are usually complex, typically including two or more large objects and hundreds of smaller bodies. Over the next 10,000 yr, these systems clear out by collisions and dynamical instabilities. In all cases analyzed here the final systems are remarkably simple. They typically include a binary with two large objects, and one or two small satellites on outer orbits with separations exceeding, by a factor of a few, the separation of the inner pair. We have not followed these systems for longer time spans. It is likely that most of the small, loosely bound satellites would not survive Gyr of dynamical and collisional evolution in the KB (Petit & Mousis 2004).

Figure 2 shows the primary radius, R₁, and the secondary-over-primary radius ratio, R₂/R₁, obtained for binaries that formed in the runs with intermediate values of f*. Each of these simulations, done for different R_eq, Ω ⩽ Ω_circ, and random seeds, produced at least one binary with similar-size large components. In some cases, more than one separate binary system was found. Values of R₂/R₁ obtained here range between ∼0.3 and 1 with most systems having R₂/R₁ > 0.7. For example, if we limit the statistics to Ω < Ω_circ, about 80% of binary systems have R₂/R₁ > 0.7.

Zoom In Zoom Out Reset image size

We compare our results to observations in Figures 3 and 4. Figure 3 shows the primary magnitude and magnitude difference, $\Delta _{\rm mag}=5\log _{10}(\sqrt{p_1/p_2}R_1/R_2)$, for the simulated binaries and known binary KBOs in the classical KB (p₁ and p₂ are the albedos of the binary components). We assumed p₁ = p₂ = 0.08 and heliocentric distance of 44 AU for the model results. The observed binary KBO parameters were taken from Noll et al. (2008a). Most simulated binaries have Δ_mag < 1, in good agreement with observations. The results that match the present observations the best were obtained with R_eq = 250 km.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The simulated distribution of R₂/R₁ is compared to observations in Figure 4. The match is strikingly good given the various uncertainties and approximations in our core simulations, except for R₂/R₁ < 0.7, where the number of simulated binaries shows a slight excess. Note that the observations are incomplete for small R₂/R₁ values, because it is hard to identify faint satellites near bright primaries.

The binary orbits obtained in our simulations are shown in Figure 5. The semimajor axis values range from ∼10³ to several 10⁵ km. Most eccentricity values are below 0.6 but cases with e > 0.6 do also occur. The observations of binary KBOs show similar trends (Figure 5, top). Notably, the orbits of several binary systems obtained in the simulations with R_eq = 750 km are similar to that of 2001 QW₃₂₂, which has a ≈ 120, 000 km and e ≲ 0.4 (Petit et al. 2008). This suggests that gravitational collapse can provide a plausible explanation for the 2001 QW₃₂₂ system. The large orbit of 2001 QW₃₂₂ is difficult to explain by other formation mechanisms discussed in Section 1.1.

Zoom In Zoom Out Reset image size

The binary inclinations show a wide spread about the plane of the angular momentum of the initial swarm (≲50° with only a few cases having 50° < i < 90°). Only one of the simulated binaries was found to have switched to retrograde rotation with respect to that of the original swarm. The prograde-to-retrograde ratio of binaries produced by GI will therefore mainly depend on the angular momentum vector orientations of the collapsing swarms. The normalized angular momentum of the simulated binary systems, J/J' (see, e.g., Noll et al. 2008a for a definition), ranges between ∼0.4 and ∼5, with larger values occurring for larger separations. For comparison, the known binary KBOs in the classical KB have 0.3 ≲ J/J' ≲ 3.5.

Interestingly, we do not find any strong correlation between the obtained J/J' values of the final binary systems, or equivalently their separation, and the assumed initial rotation Ω of the swarm. Such a correlation would be expected if most of the swarm's angular momentum ends up in J/J'. The lack of it shows how the angular momentum is distributed among the accreting bodies. If there is too much momentum initially (Ω ∼ Ω_circ), only a relatively small fraction of the mass and momentum ends up in the final binary. Indeed, it is clear that much mass is lost in the Ω = Ω_circ case as both R₁ and R₂/R₁ are on the low end of the distribution (Figure 2).

We found that several stable triple systems were produced in the simulations. For example, one of the simulations with Ω = 0.75Ω_circ produced a triple system with R₁ = 126 km, R₂ = 119 km, and R₃ = 77 km, where R₃ denotes the radius of the smallest component on the outer orbit. For comparison, (47171) 1999 TC₃₆ has R₁ = 140 km, R₂ = 129 km, and R₃ = 67 km (Benecchi et al. 2010). The two orbits of the simulated triple are nearly coplanar (Δi = 5°) and have low eccentricities (0.2 and 0.3, respectively). These properties are again reminiscent of (47171) 1999 TC₃₆. The separations of components in the simulated triples, including the one discussed here, tend to be a factor of a few larger than those in (47171) 1999 TC₃₆ (867 and 7411 km, respectively; Benecchi et al. 2010).

Figure 6 illustrates the size distribution of bodies growing in the collapsing swarm for R_eq = 250 km, Ω = 0.75Ω_circ, and f* = 10. Initially, bodies grow by normal accretion for which the growth rate of an object is not a strong function of its mass. Upon reaching a threshold of R ∼ 20 km, however, the largest objects start growing much faster than the smaller ones. This is diagnostic of runaway growth (see, e.g., Kortenkamp et al. 2000 for a review). Runaway growth occurs in the collapsing swarm because the collisional cross section of the largest bodies is strongly enhanced by gravitational focusing.

Zoom In Zoom Out Reset image size

Figure 7 shows the mean dispersion speed, V_disp, of bodies in the collapsing swarm as a function of time. It slowly increases due to dynamical stirring from large bodies but stays relatively low during the whole simulation (V_disp ≲ 2 m s⁻¹). This leads to a situation in which the escape speed, V_esc, of R > 10 km bodies largely exceeds V_disp, and the runaway accretion begins. Note also that the size distribution does not change much after 80 yr, because the large bodies run out of supply. This shows that the integration time span was roughly adequate in this case.

Zoom In Zoom Out Reset image size

We now turn our attention to the results obtained with f* = 1. Our core simulations with f* = 1 show little accretion because the collisional cross section of PKDGRAV particles is small in this case. This suggests that a longer integration time span is needed for f* = 1. We extended several core integrations with f* = 1 to T_int = 1000 yr, or about six orbital periods at 30 AU, and found that large objects accrete in these extended simulations in very much the same way as illustrated in Figure 6. The binary properties obtained in the extended runs with f* = 1 were similar to those discussed above, but better statistics will be needed to compare the results more carefully.

In additional tests, we used the same M_tot as in the core simulations and R_tot = 0.4R_Hill to see how things would work for a very dense initial concentration of solids. With f* = 10 we found that the largest object that grows out of the swarm has R = 150 km (compared to R = 92 km for R_tot = 0.6R_Hill). Notably, large bodies can also rapidly form with f* = 1 in this case, the largest having R = 110 km (compared to R = 22 km for R_tot = 0.6R_Hill). On the other hand, simulations with R_tot = 0.8R_Hill lead to smaller R values, probably because the shearing effects become important when R_tot approaches R_Hill.

This shows that the accretion timescale sensitively depends on the initial concentration of solids in the collapsing cloud. For a reference, with R_tot = R_Hill at 30 AU we obtain a concentration of solids, ρ_solids, about 15 times greater than that of the gas in the standard Minimum Mass Solar Nebula (ρ_gas; Hayashi et al. 1985), while R_tot = 0.4R_Hill leads to ρ_solids/ρ_gas ∼ 230. These values are in the ballpark of the ones produced in the simulations of Johansen et al. (2009) for protoplanetary disks with slightly enhanced metallicity.

We also performed several additional simulations with V_rand ≠ 0 and/or inelastic bouncing⁸ of PKDGRAV particles. These tests showed that binary formation occurs over a broad range of V_rand and restitution coefficient values, so long as the initial V_rand value is significantly smaller than V_circ. Placing a hard upper limit on V_rand as a function of other parameters, however, will require a systematic sampling of parameter space that is beyond the scope of this paper.

We found that the observed propensity for binary KBOs and their properties can be a natural consequence of KBO formation by GI. The binary formation in GI is robust, directly linked to the formation of large KBOs, and does not require finely tuned size distributions invoked by the HC models (see, e.g., Noll et al. 2008a). The common colors of the components of binary KBOs, their orbital parameters, including the wide binary systems such as 2001 QW₃₂₂, and triple systems such as (47171) 1999 TC₃₆, can be readily explained in this context. Moreover, the binary fraction in the KB expected in the GI model is large reaching ∼100% for a broad range of initial parameters. This favorably compares with observations that indicate, when extrapolated to smaller binary separations, that >50% of classical low-i KBOs are binary systems (Noll et al. 2008a).

The inclination distribution of binary orbits can help to constrain KB formation (Schlichting & Sari 2008b). Unfortunately, the binary orbits determined so far typically have a pair of degenerate solutions representing reflections in the sky plane. These solutions have the same a and e, but different inclinations. The very few unique inclination solutions that have been reported up to now seem to indicate that the binary orbits can be prograde (i < 90°, (42355) Typhon/Echidna; Grundy et al. 2008), retrograde (i > 90°, 2001 QW₃₂₂; Petit et al. 2008), or nearly polar (i ∼ 90°, (134860) 2000 OJ₆₇ and 2004 PB₁₀₈; Grundy et al. 2009).

The broad distribution of binary inclinations should be a signature of the formation mechanism rather than that of the later evolution because the long-term dynamical effects should not have a strong impact on the binary orbits with i < 40° and i > 140°, and cannot switch from prograde to retrograde motion (or vice versa; Perets & Naoz 2009). To explain the retrograde orbits in the GI model, we thus probably need to invoke a retrograde rotation of the collapsing clump, while the simulations of Johansen et al. (2007, 2009) seem to generally indicate prograde rotation. This issue needs to be studied in more detail, however, using a better resolution in the dynamical codes. The rotation direction of clumps in the model of Cuzzi et al. (2008) is uncertain.

Our binary formation model could also potentially apply to the Orcus–Vanth and Pluto–Charon systems, although the corresponding large M_tot values were not studied here.

Observations by Brown et al. (2010) imply sizes of Orcus and Vanth of 900 and 280 km, respectively, a mass ratio of 33, if equal densities and albedos are assumed, and the semimajor axis of the binary orbit 8980 ±  20 km. This mass ratio and orbit would be consistent with formation from a giant impact and subsequent outward tidal evolution of the binary orbit. Assuming a factor of two lower albedo for the non-icy Vanth, however, implies sizes of 820 and 640 km and a mass ratio of 2 (Brown et al. 2010). Such parameters could be difficult to reconcile with the impact formation of the Orcus–Vanth system and could rather indicate a different formation mechanism, perhaps akin to that studied in this work. Physical properties of the Orcus–Vanth system need to be determined better to discriminate between different formation models.

Using impact simulations, Canup (2005) was able to explain the main properties of the Pluto–Charon system (e.g., ∼15% mass ratio, J/J' ∼ 0.4) using an oblique, low-speed impact of an object that contained 30%–50% of the current Pluto–Charon mass. It remains to be shown, however, whether such collisions were sufficiently common in the early KB since the relevant timescale could be long (Canup 2005). On the other hand, formation of the Pluto–Charon system by gravitational collapse would require very large M_tot of the collapsing swarm, which can be a challenge for the GI theories. Interestingly, a hybrid formation model (collapse followed by an impact) is also possible, because low-speed collisions between large bodies commonly occur in our simulations.

Note that a precursor binary system similar to Pluto–Charon is needed to explain the capture of Neptune's moon Triton by exchange reaction (Agnor & Hamilton 2006), indicating that these massive binary systems were once common in the outer solar system.

Wide binary systems with similar-size components could have also formed in the inner solar system. Indeed, constraints from the Size Frequency Distribution (SFD) of main-belt asteroids indicate that standard hierarchical coagulation was not the driving force of planetesimal accretion at 2–4 AU (Morbidelli et al. 2009). Instead, asteroids have probably formed by the GI-related processes (Johansen et al. 2007; Cuzzi et al. 2010). The results of gravitational collapse simulations presented here, when scaled to a smaller Hill radius at 2–4 AU, can therefore also be applied to the asteroid belt. If so, it may seem puzzling why wide binaries with similar-size components are not detected in the asteroid belt today.

We speculate that wide asteroidal binaries, if they actually formed, would have been disrupted by collisions and scattering events during the subsequent evolution. For example, even a relatively small impact on one of the two binary components can impart enough momentum into the component's orbit to unbind it from its companion. This can happen when roughly m_iv_i > m_bv_b, where m_i and v_i are the mass and speed of the impactor, respectively, and m_b and v_b are the mass and speed of the binary component (see Petit & Mousis 2004), respectively. For the component radii R₁ = R₂ = 50 km, density ρ = 2 g cm⁻³, v_b ∼ 10 m s⁻¹ (corresponding to separation ∼0.05 R_Hill ≈ 1000 km at 2.5 AU), and v_i = 5.8 km s⁻¹ (Farinella & Davis 1992), this would imply the impactor mass m_i ≳ 10⁻³m_b or, equivalently, impactor radius r_i ≳ 5 km (for ρ = 2 g cm⁻³) for the binary to become unbound.

Since, according to Bottke et al. (2005), there are N_i ∼ 10⁴ asteroids with r_i > 5 km in the present asteroid belt, we can estimate that the present rate of unbinding collisions would be ∼2P_iN_iR²₁ = 2 × 10⁻¹⁰ yr⁻¹, where P_i = 2.8 × 10⁻¹⁸ km⁻² yr⁻¹ is the intrinsic collision probability (see, e.g., Farinella & Davis 1992 for a definition), and the factor of two appears because the impact can happen on any of the two components. This would indicate binary lifetimes comparable to the age of the solar system. The number of relevant impactors N_i, however, was likely much larger in the past, perhaps by a factor of 10–1000 (Weidenschilling 1977; Petit et al. 2001; Levison et al. 2009), than in the present asteroid belt. In addition, gravitational scattering from large planetary embryos, thought to have formed in the main-belt region (Petit et al. 2001), would have also contributed to disruption of wide binaries. It thus seems unlikely that a significant fraction of wide asteroid binaries could have survived to the present time.

A notable exception of an asteroid binary produced by gravitational collapse may be (90) Antiope (Merline et al. 2000), which is the only known asteroid binary with large, equal-size components (R₁ ∼ R₂ ∼ 45 km). We speculate that the small separation of components in the Antiope system (only ∼170 km) could have been a result of the tidal evolution of the original, possibly much wider orbit. Indeed, it has been pointed out that wide binaries with orbits that are significantly inclined (inclinations 392 < i < 1408) undergo Kozai oscillations during which the tidal dissipation is especially effective, and can shrink and circularize the binary orbit (Perets & Naoz 2009). For reference, the current inclination of the Antiope's binary orbit is ∼40° (Descamps et al. 2009). Alternatively, the (90) Antiope system could have formed by impact-induced fission of a 100 km parent asteroid (e.g., Weidenschilling et al. 2001).

The survival of binary KBOs after their formation is an open problem. Petit & Mousis (2004) estimated that several known binary KBOs (e.g., 1998 WW₃₁, 2001 QW₃₂₂, and 2000 CF₁₀₅) have lifetimes against collisional unbinding that are much shorter than the age of the solar system. These estimates were based on an assumed relatively steep SFD extending down to r_i = 5 km, which favors binary disruption, because of the large number of available impactors. When we update Petit & Mousis' estimates with a probably more reasonable SFD of KBOs given by Fraser et al. (2008), which is steep down to 60–95 km and then very shallow (differential power index ∼1.9), we find that a typical 100 km class wide binary KBO is unlikely to be disrupted over 4 Gyr (≲1% probability), except if the KB was much more massive/erosive in the past. This poses important constraints on KB formation as it may indicate that the classical low-i KBOs formed in a relatively quiescent, low-mass environment.

Levison et al. (2008) proposed that most of the complex orbital structure seen in the KB region today (see, e.g., Gladman et al. 2008) can be explained if bodies native to 15–35 AU were scattered to >35 AU by eccentric Neptune (Tsiganis et al. 2005). If these outer solar system events coincided in time with the Late Heavy Bombardment (LHB) in the inner solar system, as argued by Gomes et al. (2005), binaries populating the original planetesimal disk at 15–35 AU would have to withstand ∼700 Myr before being scattered into the KB. Even though their survival during this epoch is difficult to evaluate, due to major uncertainties in the disk's mass, SFD and radial profile, the near absence of binaries among 100 km sized hot classical KBOs (Noll et al. 2008a, 2008b) seems to indicate that the unbinding collisions and scattering events must have been rather damaging. The (617) Patroclus–Menoetius binary system, thought to have been captured into its current Jupiter–Trojan orbit from the 15–35 AU disk (Morbidelli et al. 2005), can be a rare survivor of the pre-LHB epoch, apparently because its relatively tight binary orbit (a = 680 km; Marchis et al. 2006) resisted disruption.

We thank Bill Bottke, Hal Levison, and Alessandro Morbidelli for stimulating discussions, and an anonymous referee for a very helpful report on the manuscript.

While aerodynamic forces are crucial in creating dense clumps, they are less important in the final collapse phase. The aerodynamic stopping time of a rock with radius r, density ρ, and mass m is

$$\begin{equation*} t_{\rm stop} = {\rho r \over \rho _{\rm gas} c_{\rm gas}}{,} \end{equation*}$$

where ρ_gas is the gas density, and c_gas is the sound speed.

For a rough estimate of the collision rate, we assume that the solid mass is distributed in a sphere with fractional radius f_H of the Hill radius, giving a number density n ∼ (M_tot/m)/(f_HR_Hill)³ and a virial speed $v \sim \sqrt{G M_{\rm tot}/(f_{\rm H} R_{\rm Hill})}$. With a geometric cross section, σ ∼ r², the collision time t_coll ∼ 1/(nσv) gives a ratio

$$\begin{equation*} {t_{\rm coll} \over t_{\rm stop}} \sim {\Sigma _{\rm gas} a_\odot ^2 \over M_\odot } \left(M_\odot \over M_{\rm tot} \right)^{1/3} f_{\rm H}^{7/2} \approx 0.05 \sqrt{a_\odot \over 30\ {\rm AU}}{250\,{\rm km} \over R_{\rm eq}}f_{\rm H}^{7/2}{,} \end{equation*}$$

where a_☉ is the distance to the Sun and Σ_gas ∼ ρ_gasc_gas/Ω_Kep is the gas surface density.

We thus estimate that collisions are dominant when collapse begins and f_H ∼ 1. The strong dependence on f_H means that collisions become even more dominant as collapse proceeds.

We also estimate that drag forces do not have a strong effect on a binary that forms by collapse. The KBO size R now exceeds the gas mean free path and turbulent drag applies with a characteristic timescale

$$\begin{equation*} t_{\rm drag} \sim {\rho R \over \rho _{\rm gas} v_{\rm orb}} \approx 8\ {\rm Gyr} \left(a_\odot \over 30\ {\rm AU}\right)^{2.8}\sqrt{{a_{\rm b} \over 10^4\ {\rm km}}{100\,{\rm km} \over R}}{.} \end{equation*}$$

For simplicity, we assumed a binary system with equal mass components, circular binary orbit with separation a_b and orbital speed v_orb. Since t_drag largely exceeds the ∼ Myr lifetime of the gas disk, the effect of gas drag on the binary orbit is negligible.