THE FORMATION OF PLUTO'S LOW-MASS SATELLITES

Scott J. Kenyon; Benjamin C. Bromley

doi:10.1088/0004-6256/147/1/8

1. INTRODUCTION

The Pluto–Charon binary is an icy jewel of the solar system. Discovered in 1978 (Christy & Harrington 1978; Noll et al. 2008, and references therein), the central binary has a mass ratio q_PC ≡ M_C /(M_P + M_C) = 0.10 (where M_P ≈1.3 × 10²⁵ g is the mass of Pluto and M_C ≈1.5 × 10²⁴ g is the mass of Charon; Buie et al. 2006), an orbital semimajor axis a_PC ≈17 R_P (Pluto radii; 1 R_P ≈ 1160 km; Young & Binzel 1994; Young et al. 2007), and an orbital period T_PC ≈ 6.4 d (e.g., Sicardy et al. 2006; Buie et al. 2006). Four satellites—Styx, Nix, Kerberos, and Hydra—orbit at semimajor axes, a_Styx ≈ 37 R_P, a_Nix ≈ 43 R_P, a_Kerberos ≈ 50 R_P, and a_Hydra ≈ 57 R_P (Weaver et al. 2006; Showalter et al. 2011, 2012, 2013). Curiously, the orbital periods of the satellites are roughly 3 (Styx), 4 (Nix), 5 (Kerberos), and 6 (Hydra) times the orbital period of Charon. Although the satellite masses are uncertain by at least an order of magnitude (e.g., Buie et al. 2006; Tholen et al. 2008), they are much less than the mass of the Pluto–Charon binary. For an albedo A and a mass density ρ = 1 g cm⁻³, the masses are M_S ≈ 9 × 10¹⁶A^−3/2 g, M_N ≈ 8.4 × 10¹⁸A^−3/2 g, M_K ≈ 2.6 × 10¹⁷A^−3/2 g, and M_H ≈ 1.4 × 10¹⁹A^−3/2 g (e.g., Buie et al. 2006; Showalter et al. 2011, 2012; Youdin et al. 2012). For A ≈ 0.04–1 (e.g., Marcialis et al. 1992; Roush et al. 1996; Stansberry et al. 2008; Brucker et al. 2009), the combined mass of the four small satellites is M_s ≈ 3–300 × 10¹⁹ g ≈2.3–230 × 10⁻⁶ M_P.

The architecture of Pluto–Charon and the four smaller satellites provides fascinating tests for the dynamical stability of multiple systems (e.g., Lee & Peale 2006; Tholen et al. 2008; Süli & Zsigmond 2009; Winter et al. 2010; Peale et al. 2011; Youdin et al. 2012). Dynamical fits of orbits to observations yield constraints on the masses, radii, and albedos of the satellites. Analyses using the restricted three-body problem or n-body simulations identify stable orbits surrounding Pluto–Charon and place additional constraints on the masses of the satellites. Taken together, current theoretical results for Nix, Kerberos, and Hydra suggest masses (albedos) close to the lower (upper) limits derived from dynamical fits to observations, with M_s ≲ 2.5 × 10²⁰ g and A ≳ 0.2 (Youdin et al. 2012).

The system also challenges planet formation theories. Currently, most ideas center on a giant impact scenario, where the Pluto–Charon binary forms during a collision between two icy objects with masses M₁ ≈ M_P and M₂ ≈ qM_P/(1 − q), where q ≈ 0.25–0.33 (McKinnon 1989; Canup 2005, 2011). Immediately after the impact, the binary has an initial separation of 5–10 R_P. Over 1–10 Myr, tidal evolution synchronizes the rotation of Pluto–Charon and then circularizes and expands the binary to its present separation (Farinella et al. 1979; Dobrovolskis et al. 1997).

In some scenarios, the small satellites also form during the giant impact (e.g., Stern et al. 2006; Canup 2011). Prior to the discovery of Kerberos/Styx, Ward & Canup (2006) developed an elegant model where the satellites move outward during the expansion of the Pluto–Charon orbit. If the Pluto–Charon orbit is eccentric throughout the expansion, resonant interactions with Charon can maintain the 1:4:6 period ratio of Charon, Nix, and Hydra. Ward & Canup (2006) comment that a lower-mass satellite might orbit with a 5:1 period ratio. However, they also note that most of the impact debris probably developed high-eccentricity orbits. High-velocity collisions of debris particles complicate the ability of Charon to maintain a 3:4:5:6 period ratio for four satellites.

Lithwick & Wu (2008) consider another path for satellite formation. Their analytic and numerical calculations suggest that maintaining the 4:1 (Nix) and 6:1 (Hydra) period ratios in an expanding binary requires mutually exclusive eccentricity evolution for the Pluto–Charon binary (see also Peale et al. 2011). Tidal interactions between Nix and Charon might establish the 4:1 period ratio, while interactions between Nix and Hydra establish the 4:6 period ratio. However, they consider this possibility unlikely. Instead, they propose that the Pluto–Charon binary captured many small planetesimals into a disk. Collisional evolution then led to the formation of Styx, Nix, Kerberos, and Hydra; orbital migration later drove the orbits into their current configuration.

Here, we consider whether the Pluto–Charon small satellites grow approximately in situ from material ejected during the impact. We begin by demonstrating that a giant impact capable of forming Pluto–Charon is a natural outcome of scenarios for planet formation at 15–40 AU. After the impact, collisions among smaller objects orbiting Pluto-Charon produce a disk or a ring of debris. For typical debris masses M_d ≈ 10¹⁹–10²¹ g, collisions produce several satellites with masses comparable to those of Styx–Hydra at semimajor axes reasonably close to the current orbits of the satellites. Once small satellites form, they may migrate through the leftover debris. These results support the idea that satellites with masses comparable to those of Styx, Nix, Kerberos, and Hydra grow out of material in a circumbinary disk.

We outline the steps involved in this picture in Section 2. In Section 3, we derive quantitative constraints on (1) the probability of giant impacts (Section 3.1), (2) the likelihood that Pluto was surrounded by a cloud of small objects prior to impact (Section 3.2), (3) the spreading time for the circumbinary ring (Section 3.3), (4) timescales and possible outcomes for satellite growth within the ring (Section 3.4), and (5) plausible migration scenarios (Section 3.5). In Section 4, we consider the implications of our results for the New Horizons mission (Stern 2008). We conclude with a brief summary in Section 5. To isolate basic conclusions from descriptions of the numerical simulations, we use short summary paragraphs throughout the text.

2. PHYSICAL MODEL

To construct a predictive theory for the formation and evolution of satellites around the Pluto–Charon binary, we focus on the giant impact scenario. In this picture, two massive protoplanets within the protoplanetary³ debris disk collide and produce a binary planet embedded in a ring of debris. Aside from the basic idea of a giant impact, our discussion diverges from Earth–Moon models. Unlike the Earth–Moon event, (1) the impact of icy protoplanets does not form a disk of vapor and molten material, (2) a single object with roughly the mass of Charon survives the impact, and (3) debris from the impact lies well outside the Roche limit of Pluto (Canup 2005, 2011, and references therein).

Our approach has some features in common with recent scenarios where coagulation in material just outside Saturn's ring system leads to the formation of several of the innermost moons (e.g., Charnoz et al. 2010, 2011). Unlike Saturn, the evolution of the Pluto–Charon binary has a large impact on the formation and dynamical evolution of the satellites (e.g., Ward & Canup 2006; Lithwick & Wu 2008). The large velocity dispersion of debris from the Pluto–Charon impact also complicates coagulation calculations compared to Saturn's rings, where the low velocity dispersion allows rapid growth once material escapes the Roche limit⁴ (Salmon et al. 2010; Charnoz et al. 2011).

Our calculations follow some aspects of the Ward & Canup (2010) method for deriving the formation and evolution of satellites within the circumplanetary disks of gas giant planets (see also Lunine & Stevenson 1982; Canup & Ward 2002; Mosqueira & Estrada 2003a, 2003b; Mosqueira et al. 2010; Sasaki et al. 2010; Ogihara & Ida 2012). Like Ward & Canup (2010), we consider the evolution of a circumplanetary disk where small particles grow into satellites. For satellites around a gas giant, Ward & Canup (2010) derive the evolution of a circumplanetary disk of gas and dust fed from the surrounding protoplanetary disk. They include the growth and migration of satellites within the gaseous disk. For Pluto–Charon, there is no gaseous component of the disk. Thus, we consider the evolution of a pure particle disk (see also Ruskol 1961, 1963, 1972; Estrada & Mosqueira 2006). Instead of deriving the radial disk structure and the growth of satellites analytically, we use a hybrid coagulation/n-body code to follow the growth and migration of satellites as a function of radial distance from the central binary (see also Ogihara & Ida 2012).

Although satellite formation is a continuous process within an evolving protoplanetary disk, it is convenient to list the main steps along the path from the formation of proto-Pluto and proto-Charon to the growth of stable satellites orbiting the Pluto–Charon binary.

1.
Proto-Pluto formation—within the protoplanetary disk, objects with R ≈ 500–1000 km grow from much smaller planetesimals at semimajor axes a ≈ 20–50 AU on timescales of 10–100 Myr (e.g., Kenyon & Bromley 2012). Typical calculations produce several Pluto-mass objects per AU (e.g., Kenyon & Bromley 2008, 2010).
2.
Long-term accretion—once 500–1000 km protoplanets form, they continue to attract material from the protoplanetary disk. Protoplanets accrete some small objects (e.g., Kenyon & Bromley 2008, 2010) and capture others (e.g., Ruskol 1972; Weidenschilling 2002; Goldreich et al. 2002). Once a few objects lie in bound or temporary orbits, they may interact with others passing close to the planet (e.g., Ruskol 1963; Durda & Stern 2000; Stern 2009; Pires dos Santos et al. 2012). Although the amount of material acquired by these processes is uncertain, they probably add little to the overall mass of the planet.
3.
Impact—forming a binary planet similar to Pluto-Charon requires collisions with an impact parameter b = sin θ ≈ 0.75–0.95, where θ is the impact angle (Canup 2005, 2011). After the collision shears off some material from the more massive planet (Pluto), the secondary (Charon) has an orbit with a_PC ≈ 5–10R_P and e_PC ≈ 0.1–0.8 (Canup 2011). Additional icy debris with a total mass M_d ≲ 10²³ g lies at distances of 1–30 R_P from Pluto.
4.
Circumbinary ring formation—after the binary forms, Pluto and Charon attempt to accrete and to eject the debris. Stable orbits around the binary have minimum semimajor axes ranging from a ≈ 2.5 a_PC for e_PC = 0.1 to a ≈ 4.3 a_PC for e_PC = 0.8 (see Holman & Wiegert 1999). To achieve these orbits, debris particles must gain angular momentum and lose orbital energy. If a low-mass cloud of captured objects surrounds proto-Pluto prior to impact, interactions between this cloud and debris from the impact produces a disk or ring of debris surrounding Pluto–Charon. Without this cloud, interactions among the debris particles and between the debris and Pluto–Charon leads to a somewhat smaller disk or ring. In either case, debris particles initially have large orbital e; high-velocity collisions are destructive and remove material from the binary system.
5.
Satellite growth—as Pluto–Charon clear debris from the vicinity of their orbit, destructive collisions and radiation pressure begin to deplete material from the circumbinary disk or ring. Gravitational scattering spreads material to larger semimajor axes. Eventually, collisional damping overcomes secular perturbations from the central binary; growth by mergers begins. Growth leads to an ensemble of small satellites with masses and orbits set by the initial properties of the circumbinary disk.
6.
Satellite migration—once small satellites form, they try to clear their orbits by scattering smaller objects. Gravitational scattering among smaller objects tries to fill the orbits of the larger objects. When the large objects dominate, they reduce the surface density along their orbits and increase the surface density at smaller and larger semimajor axes. These surface density enhancements are not smooth; thus, the large objects feel a torque from the small objects. If the sum of all of the torques does not vanish, the large objects migrate radially inward or outward (e.g., Goldreich & Tremaine 1982; Ward 1997; Kirsh et al. 2009). In the right circumstances, migration rates can be as large as 3–10 R_P every 1000 yr (Bromley & Kenyon 2013).
7.
Tidal expansion—over long timescales of 0.1–10 Myr, tidal interactions between Pluto and Charon modify the orbital period and the rotational periods (Dobrovolskis et al. 1997). On similar timescales, interactions between the binary and the surrounding disk can also modify the orbital period (Lin & Papaloizou 1979a, 1979b). Tidal forces rapidly synchronize Charon's rotational period with the orbital period on a timescale t_sync ≈ 10²–10⁴ yr for a_PC ≈ 5–10 R_P (see also Peale 1999). From an initial semimajor axis a_PC ≈ 5–10 R_P, tidal forces expand the orbit to the current a_PC = 17 R_P in 0.2–20 Myr (see also Farinella et al. 1979; Dobrovolskis et al. 1997).

Each of these steps involves complex physical processes that influence satellite growth. Accurate modeling of the full sequence is well beyond analytic theory and current numerical simulation. Here, we adopt the results of previous calculations of planet formation in the protoplanetary disk (e.g., Kenyon 2002; Kenyon & Bromley 2004a, 2004c, 2008, 2010, 2012). This information allows us to estimate (1) the frequency of giant impacts (Section 3.1), and (2) the amount of protoplanetary material orbiting Pluto prior to impact (Section 3.2).

Understanding the next step in the evolution—the formation and evolution of a circumbinary ring of icy material—is challenging. Using constraints from detailed smoothed particle hydrodynamics (SPH) simulations (Canup 2005, 2011), we consider several physical processes that shape the physical extent of a circumbinary ring and develop constraints on the likely properties of the ensemble of small particles in the ring (Section 3.3). We use these results as input to detailed numerical simulations of satellite growth (Section 3.4) and migration (Section 3.5).

With this discussion, we generate realistic physical models piecewise to trace a reasonable path from the formation of proto-Pluto and proto-Charon to an ensemble of satellites orbiting the Pluto–Charon binary. These models provide quantitative, testable predictions for the upcoming New Horizons flyby.

3. CALCULATIONS

To develop clear constraints on this picture, we consider satellite formation in the context of a standard model for planet formation in the protoplanetary disk (e.g., Safronov 1969; Greenberg et al. 1984; Wetherill & Stewart 1989; Kokubo & Ida 1996; Weidenschilling et al. 1997; Kenyon & Luu 1998; Nagasawa et al. 2005). In this model, the surface density of solids follows a power law

$\begin{equation} \Sigma (a) = \Sigma _0 x_m a^{-k} , \end{equation} \tag{ 1 }$

where Σ₀ is a normalization factor designed to yield a surface density equivalent to the "minimum mass solar nebula" at 1 AU (Weidenschilling 1977; Hayashi 1981), x_m is a scaling factor to allow a range of initial disk masses, and k = 1–2. For a measured in AU, Σ₀ = 30 g cm⁻² (see Kenyon & Bromley 2010, and references therein).

When the protoplanetary disk first forms, the solids consist of small 0.1–10 μm particles within a much more massive gaseous disk (e.g., Youdin & Kenyon 2013). Eventually, the solids aggregate into km-sized or larger planetesimals (e.g., Chiang & Youdin 2010; Youdin 2011). Binary collisions among planetesimals produce larger and larger protoplanets throughout the disk. At 20–50 AU, this process yields a large ensemble of Pluto-mass objects within 50–100 Myr after the formation of the Sun (e.g., Stern 1995; Stern & Colwell 1997; Kenyon & Luu 1999; Kenyon & Bromley 2004c; Kenyon et al. 2008).

3.1. The Setup: Frequency of Giant Impacts

When Pluto-mass protoplanets first form, they contain a small fraction of the mass in solid material. As long as the mass in planetesimals is larger than the mass in protoplanets, dynamical friction between the two sets of objects dominates viscous stirring among protoplanets (Goldreich et al. 2004). Thus, protoplanet orbital eccentricities remain small. Collisions among protoplanets are then very rare. When protoplanets contain more than half of the total mass, viscous stirring dominates dynamical friction. Protoplanet eccentricities increase; chaotic growth begins. During chaotic growth, protoplanets grow through giant impacts with other protoplanets (e.g., Chambers 2001, 2013; Kominami & Ida 2002; Kenyon & Bromley 2006). Giant impacts remain common until protoplanets clear their orbits of other protoplanets and leftover planetesimals.

Outcomes of giant collisions among protoplanets depend on the impact parameter, b, and the impact velocity v_imp (Asphaug et al. 2006; Leinhardt et al. 2010; Canup 2011). Usually, v_imp ≳ v_esc, where v_esc is the mutual escape velocity of a merged pair of protoplanets. Head-on impacts with b ≈ 0 then produce mergers with little or no debris. Grazing impacts with b ≈ 1 yield two icy planets on separate orbits and some debris. When b ≈ 0.75–0.95, collisions often produce a binary planet—similar to the Earth–Moon or Pluto–Charon system—embedded in a disk or ring of debris.

The frequency of collisions that produce a binary planet depends on the distribution of possible impact parameters. During the late stages of evolution, the system of protoplanets is flattened, with a vertical scale height larger than the Hill radius of a Pluto-mass planet (Wetherill & Stewart 1993; Weidenschilling et al. 1997; Chambers & Wetherill 1998; Kenyon & Bromley 2006, 2010). Assuming all impact orientations in the plane of the disk are equally likely (e.g., Agnor et al. 1999; Chambers 2013), the probability that a collision will produce a binary planet is roughly the probability of an impact with b ≈ 0.75–0.95, ∼33%.

To estimate the frequency of these impacts among icy objects during the early history of the solar system, we examine results from numerical calculations of planet formation at 15–150 AU around a solar-type star (Kenyon & Bromley 2008, 2010, 2012). In these calculations, we follow the evolution of protoplanets and planetesimals as a function of semimajor axis a and time t. Calculations begin with an initial ensemble of planetesimals with maximum radius R₀ and roughly circular orbits in a disk with a surface density Σ(a) (see Equation (1)).

Figure 1 shows the frequency (number per target per Myr) of impacts with b = 0–1 as a function of time for one set of calculations⁵ from Kenyon & Bromley (2010). The curves plot the frequency of impacts between each target with mass M_t ≳ 10²⁵ g and a projectile with a mass M_p ≳ q_p(M_t + M_p) where q_p = 0.25. As summarized in the caption, color indicates the semimajor axis of a disk annulus, ranging from 30–36 AU (violet) to 54–66 AU (magenta). Solid (dashed) lines show results for disks with k = 1 and x_m = 0.33 (0.10).

Throughout the evolution of the solar system, giant impacts are common (see also McKinnon 1989; Stern 1991; Canup 2005, 2011). Initially, all planetesimals are smaller than the target mass; the frequency of giant impacts is zero. As planetesimals grow into large protoplanets, the frequency of giant impacts increases. Because protoplanets grow faster in the inner disk than in the outer disk (Kenyon & Bromley 2004a, 2004b), giant impacts occur earlier in the inner disk (violet and turquoise curves) than in the outer disk (orange and magenta; see also Lissauer 1987; Stern 1995; Kenyon & Bromley 2005). In each disk annulus, the giant impact frequency per target rapidly rises to a clear peak when the surface density of targets and projectiles is largest. As large objects collide and merge into larger protoplanets, the giant impact frequency declines. The rate of decline, p_i∝t⁻¹, is the rate expected for collisional evolution in a protoplanetary disk (Kenyon & Bromley 2002; Dominik & Decin 2003; Wyatt 2008).

Several factors contribute to the strong variation of p_i with initial disk mass. In these calculations, the number of Pluto-mass objects N_P grows with increasing disk mass, N_P∝x_m (Kenyon & Bromley 2008, 2012). For giant impacts requiring two Pluto-mass objects, the number of possible pairs of impactors scales with the square of the disk mass. In more massive disks, the larger mass of leftover planetesimals circularizes the orbits of the largest planets more effectively. Thus, gravitational focusing factors also scale with disk mass. Combining these two factors, the impact probability grows rapidly with disk mass, $p_i \propto x_m^k$ , with k ≈ 3–4.

To put these results in perspective, we define t_N as the time to the next binary-producing giant impact on a single target. If p_i is the probability of a giant impact and if the distribution of impact orientations is random, the probability of an impact capable of producing a binary planet is p_bp ≈ 0.33p_i. The time $t_N \approx p_{bp}^{-1}$ . In a massive disk at 30–36 AU, t_N ≈ 5–10 Myr at an evolution time of 50–100 Myr. Binary planets then occur often, but they may be destroyed by another giant impact 5–10 Myr later. After an evolution time of 500 Myr to 1 Gyr (1–3 Gyr), t_N ≈ 150–300 Myr (500 Myr to 1 Gyr). Thus, binary icy planets form throughout the evolution of the planetary system.

3.2. Before the Giant Impact: Long-term Accretion

After Pluto-mass objects form, they continue to interact with material from the protoplanetary disk. Although protoplanets directly accrete some of this material, they can also capture material onto temporary or bound orbits (e.g., Ruskol 1972). To estimate the amount of material Charon or Pluto capture in the 0.5–1 Gyr in between giant impacts, we consider a simple model. Objects from the protoplanetary disk enter the planet's Hill sphere at a rate $\dot{M}_H$ . This material collides with other objects from the protoplanetary disk at a rate $\dot{M}_H \tau _{{\rm pd}}$ , where τ_pd is the optical depth of the protoplanetary disk. Collisions produce a single large remnant and some debris. The planet captures a fraction f_c of this material into bound orbits. Although other interactions are possible (e.g., Suetsugu et al. 2011; Pires dos Santos et al. 2012; Suetsugu & Ohtsuki 2013), this approach yields a reasonable initial estimate for captured material (e.g., Ruskol 1972; Weidenschilling 2002; Estrada & Mosqueira 2006).

To derive the capture rate $\dot{M}_c$ , we make several plausible assumptions. The optical depth of the protoplanetary disk is τ_pd ≈ 3Σ_pd/4ρ_pdR_pd, where Σ_pd is the surface density and (R_pd, ρ_pd) are the radius and mass density of a typical object. During the late stages of planet formation, R_pd ≈ 100 km (Kenyon & Bromley 2012). For simplicity, we set ρ_pd = 1 g cm⁻³. The protoplanetary disk has vertical scale height H_z > R_H and surface density Σ_pd = d_pdΣ₀x_m(a/1 AU)^−k (see Equation (1)), where Σ₀ = 30 g cm⁻², x_m = 1, k = 1, and d_pd < 1 is a depletion factor which measures the fraction of disk material lost to collisional grinding (e.g., Kenyon et al. 2008) and dynamical interactions (e.g., Morbidelli et al. 2008). For epochs when giant impacts are common, d_pd ≈ 0.1. Within the Hill sphere, collisions at a ≳ γR_H (with γ ≈ 0.3–0.4) produce no bound material (e.g., Toth 1999; Martin & Lubow 2011). For collisions at smaller a, f_c ≈ 0.001 (e.g., Ruskol 1972; Weidenschilling 2002).

With these assumptions, disk objects enter the Hill sphere of the planet at a rate $\dot{M}_H \approx 3 \pi \Sigma _{{\rm pd}} \Omega _{{\rm pd}} (\gamma R_{H})^2$ , where Ω_pd is the angular velocity of the planet around the Sun and the factor of three accounts for a small degree of gravitational focusing. The capture rate is then $\dot{M}_c \approx f_c \dot{M}_H \tau _{{\rm pd}}$ :

$\begin{eqnarray} \dot{M_c} &\approx& 3 \times 10^{10} \left({ d_{{\rm pd}} \over 10^{-1} } \right)^2 \left({ f_c \over 10^{-3} } \right) \left({ \gamma _c \over 0.333 } \right)^2\nonumber\\ &&\times \left({ \Sigma _0 \over {\rm 30\ g\, cm^{-2}} } \right)^2 \left({ R \over R_P} \right)^2 \left({ R_{{\rm pd}} \over {\rm 100\ km}} \right)^{-1}\nonumber\\ &&\times \left({ a \over {\rm 20\ AU} } \right)^{-3/2} {\rm g\ yr^{-1}}. \end{eqnarray} \tag{ 2 }$

After ∼0.1–1 Gyr, modest reservoirs of material orbit Charon–Pluto mass planets. Pluto captures 3–30 × 10¹⁸ g of orbiting debris. With a smaller radius and gravity, Charon attracts 1–10 × 10¹⁸ g.

To estimate the collision frequency of captured material, we assume each fragment has radius R_f ∼ 1 km, mass M_f ∼ 4 × 10¹⁵ g, orbital semimajor axis a_f ∼ 1000 R_P (∼0.25R_H), eccentricity e_f ∼ 0.1, and inclination i_f ∼ π/20. With 10³–10⁴ fragments on randomly oriented orbits, the collision rate is n_fσ_fv_f ≈ 1–10 (R_f/1km)⁻¹ Gyr⁻¹, where n_f is the number density, σ_f is the cross-section, and v_f is the relative velocity. In a cloud of 1 km objects, collisions are rare. Smaller fragments collide more frequently.

To constrain the total specific angular momentum L of the cloud, we set the specific angular momentum per fragment $L_f \approx (G M_Pa_f (1 - e_f^2))^{1/2}$ . If each vector $\bf L_f$ is randomly oriented, the cloud has a typical $\bf L \approx$ 0 with a standard deviation of roughly $\sqrt{N} L_f$ . Thus, there is a reasonable probability that a cloud of fragments has a net prograde or retrograde rotation around Pluto.

Prior to impact, a low-mass cloud of fragments probably surrounds Pluto. This mass depends on the typical particle size and the depletion of the protoplanetary disk. If collisions within Pluto's Hill sphere yield large fragments with R_f ≳ 1 km, the fragments occupy a large volume, rarely collide, and have a small net prograde or retrograde angular momentum. Despite a potentially significant mass, this material makes little or no contribution to the formation of low-mass satellites.

If collisions within Pluto's Hill sphere yield very small fragments, they collide more frequently and may produce a collisional cascade. This process slowly reduces the circumplanetary mass. With a small net prograde or retrograde $\bf L$ , collisions may produce a disk of material with a ≪ a_f (e.g., Brahic 1976, 1977). Although the properties of this disk are uncertain, debris from the Pluto–Charon impact may interact with pre-existing disk material. Captured material may then contribute to the formation of low-mass satellites.

For New Horizons, captured material at large a is challenging to detect and an unlikely hazard. Although 1 km fragments are detectable, they have a small filling factor and are not worth a detailed search. Within Pluto–Charon's Hill sphere, the total collision probability between a captured fragment and New Horizons is ≲ 10⁻¹³ throughout the flyby.

3.3. After the Giant Impact: Dynamical Spreading of the Circumbinary Ring

Shortly after the formation of the Pluto–Charon binary, icy material with a total mass of ≲ 10²³ g lies at distances of 1–30 R_P (Canup 2011). As the evolution proceeds, most of this material resides in a ring just outside the binary system. To understand whether a circumbinary ring of debris might expand into a disk, we evaluate timescales for interactions among the ensemble of debris particles and between the binary and the debris. Particles of debris initially orbit with a ≈ 1–30 R_P (Canup 2011). The particles have radius R, mass density ρ = 1 g cm⁻³, and total mass M_d. SPH simulations of the Pluto–Charon impact suggest R ≲ 10 km and M_d ≈ 10²¹–10²³ g (Canup 2011). From more general simulations of impacts, the debris probably has a range of sizes from ≲0.1 km up to 10 km (e.g., Leinhardt & Stewart 2012).

Immediately after impact, collisions among debris particles produce a collisional cascade which tends to circularize particle orbits, decrease particle radii, and remove mass from the system. Although the details of this process are uncertain, it probably leads to a ring of material with surface density Σ at a ≈ 20R_P. The Pluto–Charon binary then has semimajor axis a_PC = 5–10 R_P, orbital period T_PC = 1–3 d, and eccentricity e_PC ≈ 0.1–0.8.

After the ring forms, two opposing mechanisms—tides from the binary and collisional damping among ring particles—drive the evolution. Tidal forces set the inner edge of the ring (e.g., Pringle 1991; Holman & Wiegert 1999) and impose a forced eccentricity, e_f ≈ 1.1(a_PC/a)e_PC, on ring particles (e.g., Heppenheimer 1974; Murray & Dermott 1999; Wyatt et al. 1999; Moriwaki & Nakagawa 2004; Lee & Peale 2006; Tsukamoto & Makino 2007). Angular momentum transfer from the binary to the ring leads to a maximum eccentricity, e_m ≈ 2e_f (e.g., Wyatt et al. 1999; Moriwaki & Nakagawa 2004). The timescale for secular perturbations from tides to raise e is

$\begin{equation} t_s \approx 300 \left({ a \over 20 R_P} \right)^{7/2} \left(1 + {3 \over 2} e_{PC}^2 \right)^{-1} T_{PC}. \end{equation} \tag{ 3 }$

Aside from a, this timescale is independent of the properties of the ring. For T_PC = 1 d, t_s ≈ 1 yr.

Physical collisions transfer angular momentum among ring particles (Goldreich & Tremaine 1978; Hornung et al. 1985), which tends to circularize their orbits (e.g., Lynden-Bell & Pringle 1974; Goldreich & Tremaine 1978). If Ω is the angular velocity of ring particles with number density n, cross-section σ, and relative velocity v, the timescale for physical collisions is t_c ≈ (dn/dt)⁻¹ ≈ (nσv)⁻¹ ≈ ρR/ΣΩ. Thus,

$\begin{equation} t_c \approx 450 \left({R \over {\rm 1\ km} } \right) \left({ a \over 20 R_P} \right)^{7/2} \left({10^{22}\ {\rm g} \over M_d} \right) T_{PC}. \end{equation} \tag{ 4 }$

In an optically thick ring with M_d/Ra² ≳ 1.5, the collision time is similar to or shorter than the secular time.

When collisions are frequent, angular momentum transfer among ring particles slows the precession, prevents high-velocity collisions, and lessens the impact of destructive collisions. To conserve the angular momentum absorbed from the central binary, the ring must expand to larger a (Lin & Papaloizou 1979a, 1979b). Tides input angular momentum into the ring at a rate $\dot{L} \approx f q_{PC}^2 a_{PC}^2 \Omega _{PC}^2$ , where f is a constant of order unity (Lin & Papaloizou 1986a). Setting the rate of angular momentum input equal to $\dot{L} = \dot{a} \Omega a / 2$ , the expansion rate $\dot{a}$ for the ring is

$\begin{equation} \dot{a} \approx 5 f \left({a \over 20 R_P} \right)^{-1/2} R_P{\rm yr^{-1}}. \end{equation} \tag{ 5 }$

The ring spreads from 20 R_P to 60 R_P in 5–10 yr, ∼20–40 collision times.

Although understanding the early evolution of the ring requires a detailed analytic (e.g., Rafikov 2013) or numerical (e.g., Scholl et al. 2007; Paardekooper et al. 2012) calculation, these estimates suggest that a narrow ring of material plausibly expands into a disk at least as fast as satellites grow from the debris (Equation (4); see Section 3.4 below). For the very large collision velocities of particles with eccentricity e_m, destructive collisions can reduce typical particle sizes from ∼1 km to 0.1 km over 20–40 collision times. Significant mass removal on this timescale also seems likely.

3.4. Growth of Satellites in the Particle Disk

As the circumbinary ring of small particles expands, collisional damping reduces particle random velocities. When random velocities are low enough, collisions produce merged objects instead of a cloud of fragments. Dynamical friction lowers the random velocities of the largest objects, promoting additional growth. Eventually satellites reach sizes of 1–10 km, when they may begin to migrate through the disk.

To perform numerical calculations of planet and satellite formation, we use Orchestra, an ensemble of computer codes for the formation and evolution of planetary systems. Currently, Orchestra consists of a radial diffusion code that computes the time evolution of a gaseous or particulate disk (Bromley & Kenyon 2011a), a multiannulus coagulation code that computes the time evolution of a swarm of planetesimals (Kenyon & Bromley 2004a, 2008), and an n-body code that computes the orbits of gravitationally interacting protoplanets (Bromley & Kenyon 2006, 2011a, 2013). The coagulation code uses statistical algorithms to evolve the mass and velocity distributions of low-mass objects with time; the n-body algorithm follows the individual trajectories of massive objects. Within the coagulation code, Orchestra includes algorithms for treating interactions between coagulation mass bins and the n-bodies (Bromley & Kenyon 2006). To treat interactions between small particles and the n-bodies more rigorously, the n-body calculations include tracer particles. Massless tracer particles allow us to calculate the evolution of the orbits of small particles in response to the motions of massive protoplanets. Massive tracer particles enable calculations of the response of n-bodies to the changing gravitational potential of small particles (Bromley & Kenyon 2011a, 2011b, 2013).

We perform calculations on a cylindrical grid with inner radius a_in and outer radius a_out. The model grid contains K concentric annuli with widths δa_i = 0.035a_i centered at semimajor axes a_i. Calculations begin with a cumulative mass distribution $N_c(m_{ik}) \propto m_{ik}^{-q^{\prime }}$ of particles with mass density ρ = 1 g cm⁻³ and maximum initial mass m₀. For comparison with investigators that quote differential size distributions with N∝m^−q, q' = q + 1. Here, we adopt an initial $q^{\prime }_0$ = 0.17; thus most of the initial mass is in the largest objects.

Planetesimals have horizontal and vertical velocities h_ik(t) and v_ik(t) relative to a circular orbit. The horizontal velocity is related to the orbital eccentricity, $e_{ik}^2(t)$ = 1.6 (h_ik(t)/V_{K, i})², where V_{K, i} is the circular orbital velocity in annulus i. The orbital inclination depends on the vertical velocity, sin ²i_ik(t) = (2(v_ik(t)/V_{K, i})²).

The mass and velocity distributions of the planetesimals evolve in time due to inelastic collisions, drag forces, and gravitational encounters. As summarized in Kenyon & Bromley (2004a, 2008), we solve a coupled set of coagulation equations that treats the outcomes of mutual collisions between all particles with mass m_j in annuli a_i. We adopt the particle-in-a-box algorithm, where the physical collision rate is nσvf_g, n is the number density of objects, σ is the geometric cross-section, v is the relative velocity, and f_g is the gravitational focusing factor. Depending on physical conditions in the disk, we derive f_g in the dispersion or the shear regime (Kenyon & Bromley 2012). For a specific mass bin, the solutions include terms for (1) loss of mass from mergers with other objects and (2) gain of mass from collisional debris and mergers of smaller objects.

Collision outcomes depend on the ratio $Q_c/Q_D^*$ , where $Q_D^*$ is the collision energy needed to eject half the mass of a pair of colliding planetesimals to infinity and Q_c is the center of mass collision energy. From detailed n-body and smooth particle hydrodynamics calculations (Benz & Asphaug 1999; Leinhardt et al. 2008; Leinhardt & Stewart 2009),

$\begin{equation} Q_D^* = Q_b R^{\beta _b} + Q_g \rho R^{\beta _g}, \end{equation} \tag{ 6 }$

where $Q_b R^{\beta _b}$ is the bulk component of the binding energy, $Q_g \rho R^{\beta _g}$ is the gravity component of the binding energy, and R is the radius of a planetesimal. We adopt Q_b = 2 × 10⁵ erg g⁻¹ cm^0.4, β_b = −0.40, Q_g = 0.22 erg g⁻² cm^1.7, and β_g = 1.30 (Leinhardt et al. 2008; Leinhardt & Stewart 2009).

For two colliding planetesimals with masses m₁ and m₂, the mass of the merged planetesimal is

$\begin{equation} m = m_1 + m_2 - m_{ej} , \end{equation} \tag{ 7 }$

where the mass of debris ejected in a collision is

$\begin{equation} m_{ej} = 0.5 (m_1 + m_2) \left(\frac{Q_c}{Q_D^*} \right)^{9/8}. \end{equation} \tag{ 8 }$

To place the debris in our grid of mass bins, we set the mass of the largest collision fragment as m_L = 0.2m_esc and adopt a cumulative mass distribution $N_c \propto m^{-q_d}$ with q_d = 0.833, roughly the value expected for a system in collisional equilibrium (Dohnanyi 1969; Williams & Wetherill 1994; Tanaka & Ida 1996; O'Brien & Greenberg 2003; Kobayashi & Tanaka 2010). This approach allows us to derive ejected masses for catastrophic collisions with $Q_c \sim Q_D^*$ and for cratering collisions with $Q_c \ll Q_D^*$ (see also Wetherill & Stewart 1993; Williams & Wetherill 1994; Tanaka & Ida 1996; Stern & Colwell 1997; Kenyon & Luu 1999; O'Brien & Greenberg 2003; Kobayashi & Tanaka 2010).

To compute the evolution of the velocity distribution, we include collisional damping from inelastic collisions and gravitational interactions. For inelastic and elastic collisions, we follow the statistical, Fokker–Planck approaches of Ohtsuki (1992) and Ohtsuki et al. (2002), which treat pairwise interactions (e.g., dynamical friction and viscous stirring) between all objects in all annuli. As in Kenyon & Bromley (2001), we add terms to treat the probability that objects in annulus i_k interact with objects in annulus i_l where k, l = 1–K (Kenyon & Bromley 2004a, 2008). We also compute long-range stirring from distant oligarchs (Weidenschilling 1989).

In the n-body code, we directly integrate the orbits of objects with masses larger than a pre-set "promotion mass" m_pro. The calculations allow for mergers among the n-bodies. Additional algorithms treat mass accretion from the coagulation grid and mutual gravitational stirring of n-bodies and mass batches in the coagulation grid. To treat interactions among proto-satellites, we set m_pro = 10¹⁷–10¹⁹ g.

Our calculations follow the time evolution of the mass and velocity distributions of objects with a range of radii, r_ij = r_min to r_ij = r_max. To track the amount of material ejected by radiation pressure, we set r_min = 20 μm. The upper limit r_max is always larger than the largest object in each annulus.

For this initial exploration of satellite growth in a circumbinary ring, we derive results for (1) pure coagulation calculations with fragmentation and no n-body component and (2) hybrid calculations with the n-body component but no fragmentation. Aside from enabling a broader set of calculations per unit cpu time, these two sets of calculations probably bracket the likely timescales for satellite growth around the Pluto–Charon binary. With no stirring from the central binary, the pure coagulation calculations set a firm lower limit on the growth time. With no collisional damping or dynamical friction from small objects and no stirring from the central binary, hybrid calculations without fragmentation likely establish an upper limit on the growth time. Together, the two sets of calculations show that the likely growth time is comparable to or longer than the time for the ring to spread (Section 3.3).

We perform both sets of calculations on a grid of 32 concentric annuli extending from a_in = 15 R_P to a_out = 50 R_P. This cylindrical grid provides a reasonably accurate representation of plausible orbits of debris particles from the giant impact around a circular binary system with an initial separation of 5 R_P (Canup 2011). The inner edge of this grid lies roughly at the innermost stable orbit around a newly formed Pluto–Charon binary. The outer edge allows the grid to cover a radial extent roughly two times larger than the current radial extent (39–57 R_P) of the satellite system. For a ring of fixed mass M_d and aspect ratio a_out/a_in ≈ 3.33, we expect satellites to form at roughly the same relative positions with roughly the same mass on a formation time which scales as $a_{{\rm in}}^{2.5}$ . Thus, choosing a_in = 30 R_P and a_out = 100 R_P results in formation times that are roughly 5.5 times longer than those described below.

3.4.1. Pure Coagulation Calculations

These calculations begin with a ring of material surrounding a single central object with mass equivalent to the Pluto–Charon binary. The total mass in the ring ranges between M_d = 3 × 10¹⁹ g and M_d = 3 × 10²¹ g, with an initial surface density distribution Σ∝a⁻¹. Within each annulus of the ring, the solid particles range in size from 20 μm to 0.1 km. Radiation pressure probably ejects particles with R ≲ 20 μm on short timescales (e.g., Burns et al. 1979; Poppe & Horányi 2011). Although particles with R ≳ 0.1 km might lie within the ring, adopting a maximum particle size of 0.1 km provides a rough upper limit on the amount of mass lost to a collisional cascade. We consider calculations starting with several 1 km objects in the next section. To allow the system to find an equilibrium between the large orbital eccentricity expected of material ejected from the central binary and the small eccentricity expected from collisional damping, we adopt an initial eccentricity e₀ = 0.1 and an initial inclination i₀/e₀ = 0.5.

All calculations follow a similar pattern. With the large initial (e, i) of all particles, collisions among large particles are destructive and produce copious amounts of smaller particles. Among the smaller particles, collisional damping is effective. Through dynamical friction, smaller particles damp the orbits of the larger particles. Thus, the (e, i) for all particles declines with time. Collisions become less destructive. Eventually, collisional velocities become small enough to promote growth over destruction. During this phase, roughly 25% of the initial mass is ground down into particles smaller than 20 μm. Radiation pressure ejects this material from the binary (e.g., Burns et al. 1979; Kennedy & Wyatt 2011).

After the system reaches a low velocity equilibrium, the largest particles begin to grow rapidly. During a very short phase of runaway growth, objects grow from roughly 0.1 km to 3–30 km. As the largest objects grow, they gravitationally stir the smaller particles. Once the largest objects contain most of the mass, stirring dominates collisional damping. The orbital (e, i) of the smaller objects rapidly increases, reducing gravitational focusing factors and halting runaway growth. At the end of runaway growth, each annulus contains a handful of proto-satellites with radii of 3–30 km.

Following runaway growth, the evolution of the system stalls. The proto-satellites rapidly accrete any small particles left in their annuli. Because the small particles contain so little mass, this growth of the large particles is very small. In the particle-in-a-box approximation adopted for these pure coagulation calculations, the probability of collisions among large objects is small. Thus, the sizes of the proto-satellites reach a rough steady state set approximately by the initial mass in each annulus of the grid.

Figure 2 illustrates the growth of satellites for calculations with M_d = 3 × 10¹⁹ g (thin lines) and M_d = 3 × 10²¹ g (thick lines). As summarized in the caption, each line represents the evolution in the radius of the largest object in annuli ranging in distance from 15 R_P to 50 R_P from the center-of-mass. In these disks, satellite growth is very rapid. After an initial period of destructive collisions, it takes only 1–100 yr for an ensemble of 0.1 km fragments to grow into 10 km satellites. The growth time is inversely proportional to the initial disk mass, t_g ≈ 100(3 × 10²⁰ g/M_d) yr. Once satellites reach a maximum radius of 5–50 km, growth stalls. Satellites remain at a roughly constant radius for thousands of years.

These calculations demonstrate the inevitable growth of 5–30 km satellites from an ensemble of much smaller objects. Independent of the initial orbital eccentricity, collisional damping reduces the internal velocity dispersion to levels that promote mergers during collisions. Although the collisional cascade removes roughly 25% of the initial mass, small objects grow rapidly into satellites with radii of 5–30 km.

Despite the efficiency of small satellite formation, coagulation calculations produce too many satellites. In these models, each of the 32 annuli in the calculation yields 0–3 satellites with radii of 5–30 km. Even with small number statistics, this result is much larger than the current number of known satellites. However, the coagulation calculations do not allow large-scale dynamical interactions among satellites. For the models shown in Figure 2, satellites in adjacent annuli have orbital separations of 5–15 R_P. Thus, satellites are close enough to perturb the orbits of their nearest neighbors, leading to chaotic interactions as in simulations of terrestrial planet formation (Chambers 2001; Kokubo & Ida 2002; Kominami & Ida 2004; Kenyon & Bromley 2006).

3.4.2. Hybrid Calculations

To investigate dynamical interactions among newly formed satellites, we now consider a suite of hybrid simulations with Orchestra. Objects in the coagulation code that reach a preset mass, m_pro, are promoted into the n-body code. When a satellite with (m_ij, e_ij, i_ij) in an annulus with semimajor axis a_i is promoted, we assign (e_n, i_n) = (e_ij, i_ij), a_n = a_i + gΔa_i, and a random orbital phase, where g is a random number in the range −0.5 to +0.5. The n-body code converts these orbital elements into an initial (x, y, z) position and an initial (v_x, v_y, v_z) velocity vector. The n-body code follows the trajectories of all promoted satellites. Algorithms within Orchestra allow n-bodies to interact with coagulation particles (Bromley & Kenyon 2006, 2011a).

The starting conditions for these calculations are identical to those for the pure coagulation models. We consider three cases for the initial disk mass, m_d = 3 × 10¹⁹ g (low-mass disk), m_d = 3 × 10²⁰ g (intermediate-mass disk), and m_d = 3 × 10²¹ g (massive disk). For each case, we adopt a different promotion mass: m_pro = 3 × 10¹⁷ g (low mass), m_pro = 10¹⁸ g (intermediate mass), and m_pro = 3 × 10¹⁸ g (high mass). These promotion masses are a compromise between the "optimal" promotion mass, m_pro = 3–30 × 10¹⁶ g, required to allow newly formed n-bodies to adjust their positions and velocities to existing n-bodies and the practical limits required to complete calculations in a reasonable amount of time. For each combination of disk mass and promotion mass, we ran 20–25 simulations.

To check the accuracy of our results for the intermediate- and high-mass disks, we ran an additional 12 simulations for each of two cases with promotion masses a factor of three larger and a factor of three smaller than the promotion masses listed above. When the promotion mass is a factor of three larger, the chaotic growth phase begins later and lasts longer. Often, chaotic growth in the outer disk is well-separated in time from chaotic growth in the inner disk. This unphysical behavior contrasts with calculations with smaller promotion masses, where the timing of chaotic growth changes smoothly from the inner disk to the outer disk. As a result, these calculations have less scattering and radial mixing among proto-satellites. When the promotion mass is a factor of three smaller, the onset, character, and duration of chaotic growth changes very little. In all cases, the final number of satellites is fairly independent of promotion mass.

Because fragmentation has a small impact on the results of pure coagulation calculations per unit cpu time, these hybrid calculations do not include fragmentation. Neglecting fragmentation increases the mass available for satellite growth by roughly 25%. Without fragmentation, collisions fail to produce copious amounts of 1 m to 100 m particles. Collisional damping among these debris particles and dynamical friction between the debris and larger survivors reduces the orbital e and i, aiding runaway growth. Thus, these hybrid calculations artificially increase the evolution time required for satellites to reach R ≳ 10 km.

Figure 3 shows the evolution of the semimajor axes for promoted satellites orbiting a single central object with the mass of the Pluto–Charon binary. All objects begin with r ≲ 0.1 km (m ≲ 4 × 10¹² g). The evolution time required for these objects to reach the promotion mass scales with the initial disk mass, t_pro ≈ 10³(3 × 10¹⁹ g/m_d) yr. This timescale is roughly a factor of ten longer than derived for pure coagulation calculations with fragmentation. In the inner disk, objects have shorter orbital periods and shorter collision times. Thus, the first promoted objects appear near the inner edge of the disk. As the calculation proceeds, satellites are promoted farther and farther out in the disk. Eventually, a few promoted objects appear in the outer disk.

Every calculation experiences chaotic growth (Goldreich et al. 2004; Kenyon & Bromley 2006). During chaotic growth, the satellites scatter one another throughout the grid. Because objects grow fastest in the most massive disks, chaotic growth begins earlier—∼100 yr—in the high-mass disk and much later—∼10⁴ yr—in the low-mass disk. Chaotic growth is also "stronger" in more massive disks. In more massive disks, massive large satellites are more numerous and generate larger radial excursions of smaller satellites than in lower-mass disks.

Throughout the evolution shown in Figure 3, the n-bodies slowly accrete the very small planetesimals remaining in the coagulation grid. The number of leftovers correlates well with the number of n-bodies. Thus, the total mass in leftover planetesimals within the coagulation grid approaches zero as the number of n-bodies reaches a minimum.

In all of these calculations, the final number of satellites correlates inversely with the initial disk mass. In massive disks with M_d = 3 × 10²¹ g, there are usually 2 or 3 massive satellites with r ≈ 50–80 km at the end of the calculation. As we reduce the initial disk mass, the calculations produce more satellites with smaller masses. In intermediate-mass disks with M_d = 3 × 10²⁰ km, chaotic growth leaves all of the initial mass in 4–5 satellites with radii of 20–30 km. In low-mass disks, there are very few mergers during chaotic growth. After ∼10⁷ yr, there are typically 7–9 satellites with radii of 7–12 km.

To test how outcomes change with the initial size distribution, we consider a suite of 15 simulations with two additional 10 km planetesimals placed randomly in the grid. These large planetesimals stir their surroundings and thus counteract collisional damping by nearby small planetesimals. With much larger masses than any other planetesimals in the grid, they have large collisional cross-sections and can grow very rapidly.

Figure 4 compares the time evolution of the semimajor axes for n-bodies in calculations with (lower panel) and without (upper panel) two additional massive planetesimals. In standard calculations with M_d = 3 × 10¹⁹ g, there is a short period of chaotic growth at ∼10⁴–10⁵ yr followed by a long quiescent period with an occasional merger. When there are two large planetesimals at the onset of the calculation, these planetesimals slowly accumulate small planetesimals. Meanwhile, small planetesimals in the rest of the disk rapidly grow to sizes (4–5 km) that allow promotion into the n-body code. From Figure 4, the timescale for these objects to reach the promotion mass is roughly 10³–10⁴ yr independent of the two large satellites. However, stirring by the two large satellites promotes an earlier oligarchic growth phase which enables more uniform growth of the largest planetesimals. Thus, there are many more promoted objects and a more intense chaotic growth phase. In addition to the two initial massive planetesimals, one or two other planetesimals accrete all of the other promoted objects. After ∼10⁵ yr, there are only a few small n-bodies left. The 3–4 massive satellites accrete these objects in a few hundred thousand years.

This result is typical of all calculations in low- or intermediate-mass disks that begin with a few large (10 km) planetesimals placed randomly throughout the grid. Stirring by the large objects slows runaway growth and allows a larger group of oligarchs to grow rapidly. Although these oligarchs reach the promotion mass on similar timescales, many more of them reach the promotion mass. Typically, calculations with initial large planetesimals produce 2–4 times as many n-bodies as calculations without large planetesimals. Stirring by the 10 km planetesimals and interactions among the many n-bodies leads to a very active phase of chaotic growth. During chaotic growth, the radial excursions of the n-bodies are 2–3 times larger in the semimajor axis. Collisions among n-bodies are more frequent, leading to a system with fewer, but more massive, satellites. With a set of two initially large planetesimals, low- (intermediate-) mass disks produce 3–4 (2–3) massive satellites instead of 7–9 (4–5).

For all of these hybrid calculations, the radius evolution follows a standard pattern (Figure 5). As in Figure 2, the largest objects in the coagulation grid grow from 0.1 km to 1–10 km. Once they are promoted into the n-body grid, they interact with objects in the coagulation grid and all of the promoted objects in the n-body grid. Unlike the pure coagulation calculations, growth does not stall at 5–30 km. The extra interactions between promoted objects and the larger volume sampled by scattered objects enables a few large objects to accumulate nearly all of the remaining mass in the grid. Thus, these objects reach radii that are 30% to 100% larger (and 2–8 times more massive) than the largest objects in the pure coagulation calculations.

3.4.3. Summary of Coagulation and Hybrid Calculations

The suite of pure coagulation and hybrid calculations demonstrates that the Pluto–Charon satellites grow rapidly from a disk of small planetesimals. In calculations with fragmentation, destructive collisions grind down roughly 25% of the initial disk mass into 20 μm particles which are ejected from the binary system. Collisional damping among larger debris particles reduces orbital eccentricity of 0.1–10 m objects. Dynamical friction between these small planetesimals and the surviving large planetesimals damps the orbital eccentricities of the largest objects, greatly reducing the frequency of destructive collisions and enabling rapid growth of surviving planetesimals. Thus, fragmentation removes mass and aids the eventual rapid growth of satellites.

In hybrid calculations without fragmentation, less collisional damping and dynamical friction slow growth considerably. However, satellites still grow on timescales similar to the expansion time for the central binary. The number of satellites N_s, typical satellite radius R_s, and the time for satellites to reach their final mass t_f scales with the initial mass of the planetesimal disk. For calculations without any large (10 km) fragments at t = 0 and M_d ≈ 3–300 × 10¹⁹ g, we infer

$\begin{equation} \left. \begin{array}{@{}l@{}}N_s \\ R_s \\ t_f \end{array} \right\rbrace = \left\lbrace \begin{array}{@{}ll}2^{22.5 - {\rm log} M_d} \\ 27.5 \left({ 3 \times 10^{20}\ {\rm g} \over M_d } \right)^{0.425}&{\rm km}\ {\rm no \ initial \ large \ fragments} \\ 10^4 \left({ 3 \times 10^{20}\ {\rm g} \over M_d } \right) \ {\rm yr}. \end{array} \right. \end{equation} \tag{ 9 }$

In low- and intermediate-mass disks, calculations with several 10 km fragments produce fewer satellites with larger masses. For M_d ≈ 3–30 × 10¹⁹ g, N_s ≈ 3–5 and R_s ≈ 15–60 km instead of the N_s ≈ 4–10 and R_s ≈ 7–30 km for calculations without large fragments. Because large fragments tend to sweep up any small planetesimals, calculations with fragments yield a smaller spread in the number and masses of planetesimals than calculations without fragments.

3.4.4. Observations: Comparisons and Predictions

To conclude this section, Figure 6 plots pairs of (e, i) for satellites that survive for 10⁷ yr. For most satellites, the orbital inclination is small: the average inclination is 〈i〉 = 0 fdg 1 ± 0 fdg 2. Despite the apparent excess of high i objects at small e for the low- and intermediate-mass disks, the average i and the dispersion in the average do not vary with initial disk mass. The average e grows with disk mass, however, with 〈e〉 = 0.02 ± 0.02 for M_d = 3 × 10¹⁹ g, 〈e〉 = 0.03 ± 0.02 for M_d = 3 × 10²⁰ g, and 〈e〉 = 0.06 ± 0.06 for M_d = 3 × 10²¹ g. For small and intermediate initial disk masses, the average eccentricity is very low. For large disk masses, the 〈e〉 is a factor of 2–3 larger.

Although the final orbital elements of satellites are independent of the promotion mass, e depends on the initial masses of the largest planetesimals. When we add several large planetesimals to the initial distribution of 0.1 km and smaller planetesimals, the calculations yield fewer satellites with larger masses (Figure 4). As in calculations with larger total disk masses, a few large satellites stir up their surroundings more than many smaller satellites. These satellites then have larger e.

The orbital elements of model satellites agree reasonably well with the observed elements of the Pluto–Charon satellites. Roughly 10% of the model satellites have e ≲ e_H, where e_H is the orbital eccentricity of Hydra. Nearly all (83%) satellites have orbital inclinations i ≲ i_H, where i_H is the inclination of Hydra. Satellites with integer period ratios are also common. For calculations with M_d = 3 × 10²¹ g, 60% of the outer satellites have period ratios reasonably close to the 5:3, 2:1, or 3:1 commensurabilities with the inner satellite. Disks with smaller initial masses produce more closely-packed, lower-mass satellites which often lie close to the 3:2 and 4:3 commensurabilities. For these lower-mass disks, 50% to 70% of satellites have period ratios within a few percent of the 4:3, 3:2, 5:3, 2:1, or 3:1 commensurabilities.

These results are encouraging. Model satellites often lie close to the observed orbital commensurabilities in the Pluto–Charon satellites. The derived inclinations of model satellites also agree very well with observations. Although lower-mass disks produce satellites with the smallest e, derived eccentricities for all calculations are a factor of 4–20 larger than observed. Thus, the calculations have two successes—the small inclinations and the likelihood of close orbital commensurabilities—and one failure—the lack of very small eccentricities for the satellites.

In addition to satellite formation, pure and hybrid coagulation calculations leave behind remnant disks of small particles. When satellites are close to their final masses, remnant disk particles have large orbital eccentricities and large vertical scale heights. At this epoch, it is more likely for small disk particles to collide with one another than to collide with a 10–15 km satellite. A collisional cascade gradually depletes the remnant disk (e.g., Kenyon & Bromley 2008, 2010, and references therein). For conditions in the Pluto–Charon system, the remnant disk mass at ∼10⁴ yr is roughly 10¹⁵–10¹⁷ g and declines with time as t⁻¹ (see also Dominik & Decin 2003; Kenyon & Bromley 2004a). Thus, the disk has a fairly substantial mass for 1–10 Myr. At later times, the loss of material from collisional grinding reaches an equilibrium with material captured from collisions between small Kuiper Belt objects and the satellites (e.g., Stern et al. 2006; Poppe & Horányi 2011). Based on current analyses of the capture rate, this equilibrium mass is roughly 10¹³–10¹⁴ g.

If satellites have much larger radii, ∼40–70 km, they accrete small particles faster than the collisional cascade removes them. The remnant disk mass then declines more rapidly, on timescales of roughly 10⁵ yr instead of 1–10 Myr. Direct albedo measurements from New Horizons will test this possibility.

To explore the outcomes of formation models in more detail, we now consider numerical calculations of satellite migration through the circumbinary disk surrounding Pluto–Charon. Aside from testing analytic estimates for essential migration parameters, these calculations allow us to estimate the frequency of migration-driven mergers and to examine the evolution of e and i for an ensemble of migrating satellites.

3.5. Migration of Satellites Through the Circumbinary Disk

3.5.1. Background

Throughout the growth process outlined in Section 3.4, gravitational scattering plays an important role in the local structure of the disk. The Hill radius,

$\begin{equation} R_H\approx 10 \left({R \over 1\ {\rm km} } \right) \left({a \over {20\, R_P} } \right) {\rm km} , \end{equation} \tag{ 10 }$

defines the region where the gravity of a satellite overcomes the gravity of Pluto–Charon. If the satellite can clear its orbit of small objects, it can open a gap in the radial distribution of solids. For convenience, we define the Hill radius necessary for a satellite to open up a gap in a dynamically cold disk surrounding the Pluto–Charon binary (e.g., Rafikov 2001; Bromley & Kenyon 2013):

$\begin{eqnarray} R_{{\rm gap}}&\approx& 19 \left({ v_r \over {75\ {\rm cm\ s^{-1}}} } \right)^{2/3} \left({ M_d \over 10^{20}\ {\rm g} } \right)^{1/3} \left({ R \over {\rm 1\ km} } \right)^{-1/3}\nonumber\\ &&\times \left({ a_c \over 20 R_P\ } \right)^{1/3} \ {\rm km}. \end{eqnarray} \tag{ 11 }$

Satellites with R_H ≈ 20 km have physical radii R ≈ 2 km (Equation (10)). For circumbinary disks with M_d ≈ 10²⁰ g, satellites with R ≳ 2 km open gaps in the disk. All four Pluto–Charon satellites have R ≳ 3 km. Thus, each can open a gap in a cold circumbinary disk.

For satellites with R ≳ R_gap, there are two possible modes of migration (e.g., Bromley & Kenyon 2011b, 2013). If the satellite can (1) clear a gap in its corotation zone and (2) migrate across this gap in one synodic period, the satellite undergoes fast migration. Satellites with R_H > R_gap (R_H = R_gap) satisfy the first (second) condition. Although the size of the gap grows with R, the synodic period decreases with R. Thus, there is a maximum Hill radius—defined as R_fast—which allows a satellite to satisfy the second condition. With this definition, satellites with R_gap ≲ R_H ≲ R_fast undergo fast migration. For the Pluto–Charon binary,

$\begin{equation} R_{{\rm fast}}\approx 20 \left({ M_d \over 10^{20}\ {\rm g} } \right)^{1/2} \left({ a \over 20 R_P\ } \right)^{3/2} \ {\rm km}. \end{equation} \tag{ 12 }$

Depending on the disk properties and satellite location, this Hill radius corresponds to objects with physical radii of a few kilometers. Coupled with the limits on R_gap, satellites with physical radii of 2–5 km undergo fast migration. These objects drift radially inward or outward at a rate

$\begin{equation} \dot{a}_{{\rm fast}} \approx \pm 5 \left({ M_d \over 10^{20}\ {\rm g} } \right) \left({ a \over 20 R_P\ } \right)^{3/2} {\rm km\ yr^{-1}}. \end{equation} \tag{ 13 }$

On timescales comparable to the accretion time of 200 yr, 2–5 km objects drift roughly 1 R_P. Thus, growth and migration are simultaneous.

Massive satellites with R_H > R_fast migrate through the disk at a rate that depends on the disk viscosity. In gaseous disks with a large viscosity, this "gap" migration can transport gas giants from 5 to 10 AU to within a few stellar radii of a solar-type star in ∼1 Myr (Lin & Papaloizou 1986b; Ward 1997; Nelson & Papaloizou 2004). In a disk of particles, the smaller viscosity derived from gravitational scattering leads to a smaller gap migration rate (Ida et al. 2000; Cionco & Brunini 2002; Kirsh et al. 2009; Bromley & Kenyon 2011b; Ormel et al. 2012). For a circumbinary disk in Pluto–Charon, the expected migration rate is (Bromley & Kenyon 2013)

$\begin{equation} \dot{a}_{II} \approx -0.15 \left({ M_d \over 10^{20}\ {\rm g} } \right) \left({ R \over {\rm 20\ km} } \right) \left({ a \over 20 R_P\ } \right)^{1/2} {\rm km\ yr^{-1}}. \end{equation} \tag{ 14 }$

Satellites with R ≈ 10–20 km migrate 1 R_P every 10⁴ yr. This rate is slower than the growth rate. Thus, large satellites grow faster than they migrate.

3.5.2. Migration of a Single Satellite

To explore satellite migration in a disk surrounding the Pluto–Charon binary, we use the N-body component of the Orchestra code. In this mode, the calculations directly track interactions between satellites and massive "super-particles" that represent small particles in the disk. For massive objects within a disk around a single central object, our calculations (Bromley & Kenyon 2011a, 2011b) reproduce published results of previous investigators (e.g., Malhotra 1993; Hahn & Malhotra 1999; Kirsh et al. 2009). Our investigation of migration within Saturn's rings includes extensive tests with gaps in the disk, orbital resonances, and the gravitational perturbations of distant massive moons outside resonance (Bromley & Kenyon 2013). To test the theoretical limits on migration rates in a disk surrounding Pluto–Charon, we first consider calculations of a lone satellite with R ≈ R_gap − R_fast within a particle disk around a single central object with the combined mass of Pluto–Charon. We then consider how a binary central object modifies the mode and rate of migration. Because the hybrid calculations often produce many small satellites, we conclude this section with simulations of multiple satellites orbiting a central binary.

For this suite of simulations, we adopt a disk surface density distribution, Σ∝a⁻¹, with a fixed mass of 3 × 10²⁰ g. The disk extends from a = 20 R_P to a = 70 R_P around a single object with a radius of 1 R_P or a binary with a separation of 5 R_P. We follow an ensemble of super-particles and satellites in an annulus of full width of 50 R_P. Super-particles have masses of 1/2000^th the mass of the satellite. These objects interact with the satellite and the central mass, but not with each other. Satellites have fixed bulk mass density, ρ = 1 g cm⁻³.

At the start of each simulation, super-particles are dynamically cold. Thus, the disk is geometrically thin, with H_z/R_H ≲ 1. Particle trajectories evolve solely by interactions with the central (single or binary) object and massive satellites. Unlike our simulations of Saturn's A ring (Bromley & Kenyon 2013), there is no collisional damping among super-particles. These initial conditions are ideal to assess a satellite's ability to migrate through the disk and lead to fairly robust upper limits on the migration rate. In a more realistic disk, collisional damping, dynamical friction, and viscous stirring generally produce increases in the velocity dispersion and vertical scale height of disk particles on timescales comparable to the growth time (see Section 3.4). Because migration rates fall off as the inverse cube of the mean disk particle eccentricity (Ida et al. 2000; Kirsh et al. 2009; Bromley & Kenyon 2011b), we expect that migration timescales in a real disk are somewhat longer than our estimates for idealized disks.

These simulations assume a constant mass of small particles in the disk. In reality, collisions and interactions with the central binary deplete the small particles which drive the migration of proto-satellites. Shortly after the Pluto–Charon impact, however, collisions are destructive. At this epoch, 1–10 km objects that survive the impact may migrate through the expanding disk. While collisional damping calms disk material, large survivors first migrate slowly through the disk (Equation (14)); later, they may migrate more rapidly (Equation (13)). Once collisions yield larger merged objects, small particles deplete on the collision time, which is similar to the timescale for fast migration. During this phase, newly formed 1–10 km objects may migrate rapidly through the disk. As these objects continue to grow, they gradually open up gaps and deplete the disk of small particles. Because our simulations do not currently provide an accurate picture of gap migration with accretion, we restrict our calculations to conditions appropriate for fast migration during the early evolution of the disk.

Figure 7 illustrates results for single satellites around a single central mass (left panels) and a central binary system (right panels). At a = 20 R_P (Figure 7; upper left panel), satellites with R ≈ 1–4 km should open a gap and undergo fast migration (Equations (11) and (12)). Satellites with R = 3 km migrate at rates $\dot{a} \approx$ 20 km yr⁻¹, close to the rate predicted in Equation (13). As the satellite radius increases from 3 km to 5–10 km, the migration rate drops considerably. These larger satellites have too much inertia to maintain the fast radial drift rate. Larger satellites with R ≳ 10 km migrate at rates close to predicted rates for gap migration (Equation (14)).

Migration rates also depend on the semimajor axis of the satellite (Figure 7; lower left panel). When the velocity dispersion of the disk is fixed, it is harder for smaller satellites with R = 3–4 km to open up a gap in the disk at larger a than at smaller a (Equation (11)). Thus, small satellites migrate much more slowly at large a. However, R_fast also increases with a. Larger satellites migrate more rapidly at larger a than at smaller a. For a larger satellite with R = 8 km, the migration rate slows considerably from a = 60 R_P (roughly 20 km yr⁻¹) to a = 40 R_P (roughly 10 km yr⁻¹) to a = 20 R_P (≲ 1 km yr⁻¹). The relative change in the migration rate follows theoretical predictions (Equations (13) and (14)).

The presence of a central binary changes these results modestly (Figure 7; right columns). The most obvious impact of the Pluto–Charon binary is the increased velocity dispersion of disk particles. Larger velocity dispersion reduces the effectiveness of torque exchange between the satellite and disk particles in the satellite's corotation zone, making it more difficult for a satellite to open up a gap in the disk (Equation (11)). Thus, smaller satellites migrate more slowly around a binary. There is some evidence from the simulations that larger satellites may migrate more rapidly around a binary. In these cases, larger satellites first clear their corotation zone and halt their radial drift. As the simulation proceeds, stirring from the binary feeds some material into the corotation zone, commencing a kind of "attenuated" migration (Bromley & Kenyon 2011b). We speculate that this behavior may allow larger objects to drift more quickly through a circumbinary disk.

Simulations with a variety of disk masses confirm these general trends. Under the right conditions in disks with masses larger than ∼10¹⁸ g, satellites with R ≈ 1–20 km can undergo fast, gap, or attenuated migration with rates ranging up to 20 km yr⁻¹. In lower-mass disks, small satellites cannot clear the gap required to initiate fast migration. Although large satellites can form gaps in disks with M_d ≲ 10¹⁸ g, the minimum radius required to open a gap is larger than the maximum radius for fast migration. Thus, large satellites cannot migrate in the fast mode. Instead, these objects migrate a factor of ten more slowly in the gap migration mode.

To test whether the disk can circularize the orbits of small satellites, we derive predicted rates for 4–10 km objects in disks with M_d = 1–10 × 10¹⁷ g surrounding a single central object. For satellites in this size range, damping rates are independent of radius. The damping timescales with the disk mass:

$\begin{equation} t_{{\rm damp}} \approx {e \over \dot{e}} \approx 10^5 \left({10^{17}\ {\rm g} \over M_d} \right) \ {\rm yr}. \end{equation} \tag{ 15 }$

Even in a low-mass disk, damping times are comparable to the formation time for 4–10 km satellites.

Estimating damping rates for satellites orbiting a central binary is complicated by precession, orbital resonances, and other dynamical issues. For the small rates indicated by simulations with a single central object, accurate estimates require a substantial investment of cpu time. When the Pluto–Charon binary has a circular orbit, dynamical friction between the satellite and disk particles will still circularize the satellite's orbit. In an eccentric Pluto–Charon binary, the damping time probably depends on the timescale for the disk and tides to circularize the binary. Because these processes act on similar timescales, deriving the circularization time for a satellite orbit is very cpu intensive and beyond the scope of this study.

3.5.3. Migration of Multiple Satellites

To investigate whether ensembles of proto-satellites can migrate, we examine a second suite of simulations (Figure 8). Here, we maintain the same disk surface density distribution and focus on satellites with R = 4 km or 10 km orbiting within a wide annulus spanning orbital distances of 35R_P to 55R_P. Satellites and disk material orbit a compact Pluto–Charon binary with orbital separation a = 5 R_P. For material at 35–55 R_P, stirring by the binary has a smaller impact than at 20 R_P.

The top panel of Figure 8 follows a simulation of five 10 km satellites initially evenly spaced in orbital separation. With typical Hill radii of 0.2 R_P, these satellites interact weakly among themselves. All have radii larger than R_fast; migration rates are slow. Modest scattering of disk particles yields some inward or outward motion and occasional stronger gravitational interactions among the satellites. However, these satellites are fairly stable over 1000 yr and longer timescales.

Simulations with ensembles of five 4 km objects produce more complicated outcomes. With much smaller Hill radii, these satellites interact very weakly among themselves. At the onset of each simulation, however, all begin to migrate in the fast mode. Some migrate inward; others migrate outward. When neighboring satellites migrate in the same direction, migration is limited by the stirred up wake of disk particles left behind by the first satellite to pass through that portion of the disk (Bromley & Kenyon 2011b). Satellites do not cross these wakes. Thus, groups of small satellites migrating in the same direction drift radially inward/outward a small distance before the migration stalls.

Small satellites migrating in opposite directions lead to more interesting outcomes. Often, these satellites scatter one another, sometimes exchanging orbits (Figure 8; middle panel). Depending on the sense of migration and the state of the disk after scattering, satellites may undergo multiple scattering events before settling down in roughly stable configurations within a stirred up disk. Sometimes, satellites merge. Because merged satellites are massive, they migrate more slowly and end up in a stable configuration. In our suite of simulations, mergers are less common (∼10%) than scattering events.

Simulations with a mix of 4 km and 10 km satellites also lead to interesting outcomes (Figure 8; lower panel). If a small satellite undergoing fast migration lies between two larger satellites, its larger migration rate guarantees that it will catch up to one of its slowly moving neighbors. When it does catch up, it may "bounce off" the wake of its neighbor and begin migrating in the opposite direction. In roughly half of all interactions, the larger satellite either scatters or merges with the smaller satellite. If scattering places the small satellite on an orbit with a pericenter inside a_min, the central binary can eject the small satellite out of the system.

3.5.4. Summary of Migration Calculations

For satellites with R ≈ 4–10 km in a disk surrounding the Pluto–Charon system, migration is ubiquitous. In a cold, geometrically thin disk, satellites in this size range can undergo fast mode or slower, gap migration at rates very close to those predicted by analytic theory. These rates lead to significant radial motion, ∼1–10 R_P, on timescales comparable to the growth time of 1–10 km satellites, ∼10²–10³ yr.

Interactions between the disk and the central binary can augment or reduce migration rates. Although changes to the rates are modest, smaller (larger) satellites generally migrate more slowly (rapidly). Changes are more significant closer to the binary, where the jitter of particle orbits in the disk modifies the number of particles in the corotation zone of a satellite. Thus, small (large) satellites that form closer to the binary are less (more) likely to migrate than satellites that form farther away from the binary.

Ensembles of migrating satellites produce a variety of interesting dynamical phenomena. In systems with several small satellites or a mix of large and small satellites, differential migration enables large-scale scattering events and satellite mergers. When scattering places a satellite on an orbit with a pericenter smaller than a_min, the central binary ejects the satellite from the system. Mergers among migrating satellites allow lower-mass disks to produce more massive satellites.

The ability of satellites to migrate through an evolving disk depends on the ratio of the vertical scale height H_z to the Hill radius R_H. When collisional damping dominates particle stirring, H_z/R_H ≲ 1. Satellites can open gaps in the disk and migrate at the nominal rates. In a disk dominated by particle stirring, H_z/R_H ≈ 3(a/20R_P)^1/2. When H_z/R_H ≳ 1, satellites cannot open gaps in the disk. Migration is then "attenuated," with a rate smaller by a factor of roughly (H_z/R_H)³ ≈ 30–100 at a ≈ 20–60 R_P (Ida et al. 2000; Kirsh et al. 2009; Bromley & Kenyon 2011b). Despite this reduction, the attenuated gap migration rate is still substantial, ∼1 R_P every 10 Myr, similar to the timescale for tides and dynamical interactions to modify the semimajor axis of the binary orbit.

The ability of satellites to migrate efficiently also depends on the depletion rate of small particles in the disk. Small particles drive the migration of large particles. As particles grow from 0.1 km and smaller objects into 1 km and larger objects, the mass of the disk capable of driving migration declines. When the largest objects have R ≈ 1–5 km, most of the disk mass is in small objects. Thus, the large objects migrate at the nominal rate. Once the largest objects reach sizes of 10–30 km, small objects contain less than half of the initial disk mass. Migration is then slower than the nominal rate.

Scaling the migration rates derived in the simulations to the hotter disks expected from the coagulation and hybrid calculations in Section 3.4, satellites with radii similar to Kerberos and Styx may migrate significantly as they grow. Although larger satellites such as Nix and Hydra migrate more slowly, radial motion of a few Pluto radii seems likely after they reach their final sizes.

Although remnant particle disks with low masses, ∼10¹⁷ g, cannot drive migration on short timescales, they are effective in circularizing the orbits of satellites with sizes of 1–10 km. The timescale to circularize the orbit of a satellite is a factor of 10–100 times shorter than the timescale for tides to expand the orbit of the Pluto–Charon binary. Large final eccentricities, e ≈ 0.02–0.1, are a common feature in hybrid simulations of satellite formation (Section 3.4). Once satellites reach their final masses on timescales of 10⁴–10⁵ yr, remnant disks with masses of 10¹⁶–10¹⁷ g are ubiquitous. Thus, the migration simulations demonstrate that satellites embedded in a low-mass disk will have the small eccentricities observed in the Pluto–Charon system.

This discussion suggests two plausible phases for migration in a circumbinary disk surrounding Pluto–Charon. During the early stages of the evolution, rapid migration of small satellites is possible when collisional damping dominates particle stirring. During this phase, growing satellites can migrate at least several Pluto radii. Later on, when stirring by newly formed satellites dominates collisional damping, attenuated migration in a low-mass disk slowly modifies satellite orbits as the inner binary expands. During this period, satellites may migrate much less than a Pluto radius, especially if the collisional cascade removes most of the small particles remaining in the disk.

4. DISCUSSION

Our analysis in Section 3 suggests that the steps outlined in Section 2 provide a reasonable path to the formation of small satellites close to their current positions within debris from the giant impact. This process begins (step 1) with the growth of Pluto-mass planets in a circumstellar disk of solid material. Once Plutos form, giant impacts are common (step three). During intervals between impacts (step two), Plutos capture ∼10¹⁹ g of debris. The debris has large a and probably does not interact with material from the impact. Following impact (step 4), material ejected from the impact lies at a ≈ 20 R_P. In a reasonably massive disk composed of 0.1–1 km particles, angular momentum transfer from the binary to the ring and among ring particles spreads the ring to the current positions of the satellites on fairly short timescales. As the ring spreads, high-velocity collisions diminish the sizes of the largest particles and remove mass from the ring. Eventually, collisional damping overcomes secular perturbations from the binary, enabling the growth of satellites within a broad ring. As satellites grow, they migrate through the ring. Although our calculations do not yield satellites at the same positions as Styx–Hydra, they yield satellites at similar positions. Often, stable satellites lie in approximate resonances.

Despite the plausibility of this picture, there are definite uncertainties. During the interval prior to impact, Pluto and Charon probably capture a modest amount of material. However, it is not clear whether any of this material can interact with debris from the impact. After impact, there is a broad range of likely values for the binary eccentricity and the properties of the ejecta. Within this range, it is not certain that collisions can damp velocities and redistribute angular momentum fast enough to overcome secular stirring, spread the ring into a disk, and enable the growth of smaller particles. Although it is plausible that a few large particles with R ≈ 1–10 km escape destruction during ring formation and serve as the seeds for large satellites (Section 3.4), they may be unable to flow out in an expanding ring of smaller objects. Even if satellites form within material in an expanded disk of material, they must survive the expansion of the Pluto–Charon binary. For satellites at a ≲ 30 R_P, survival seems unlikely (Peale et al. 2011). To survive at larger a, satellites must avoid interactions that destabilize their orbits.

Once a stable disk forms at a ≈ 15–70 R_P, our calculations yield testable predictions for the masses and orbital configurations of the Pluto–Charon satellites. The timescale for satellite growth depends on the properties of the binary and the ring of debris. In any configuration, there is a period where collisional damping must overcome tidal forcing from the central binary. Although we have not calculated this evolution in detail, our estimates in Section 3.3 imply a timescale of at least 1–10 yr. As this evolution proceeds, the ring probably loses a significant fraction of its initial mass. Once collisions calm the disk, the formation time for 1–10 km satellites is

$\begin{equation} t_f \approx 10^3 \left({3 \times 10^{19} \ {\rm g} \over M_d } \right) \left({ a \over 30 R_P} \right)^{5/2} \ {\rm yr}. \end{equation} \tag{ 16 }$

During oligarchic and chaotic growth, proto-satellites radially migrate through the leftover debris. Although proto-satellites cannot fall into the Pluto–Charon binary, satellites can traverse a substantial fraction of the circumbinary disk throughout chaotic growth. Migration can induce mergers among lower-mass satellites, enabling more massive satellites to form in lower-mass disks. Ensembles of satellites with similar masses migrate in step with one another, enabling satellites to find stable equilibrium orbits around Pluto–Charon.

Migration rates depend on the nature of the central object. Stirring by a central binary can heat the particle disk and slow migration of satellites with radii of O(1) km. The binary may also increase migration rates of larger satellites. Migration is driven by the asymmetry of inward and outward satellite-disk scattering events. Compared to a single central mass, a binary preferentially disturbs a larger fraction of the orbits of the inwardly scattered disk particles and prevents their return to the satellite. We plan additional simulations to explore whether this possibility is sufficient to produce the derived changes in migration rates.

Migration rates are also sensitive to the mass of small objects remaining in the disk. When satellites reach sizes of 1–10 km, small particles still contain most of the mass. If collisional damping can maintain a small vertical scale height, satellites migrate 1–10 R_P through the disk. As these satellites grow, they migrate more and more slowly due to the larger vertical scale height and smaller total mass in small particles. Eventually, migration stalls completely. Throughout this process, interactions between the satellites and the remnant disk of small particles circularize the orbits of the satellites.

Together with previous results from Youdin et al. (2012), our calculations indicate several conclusions about each satellite.

1.
The semimajor axis of Styx is almost as close as possible to the minimum stable circumbinary semimajor axis a_min. In our framework, Styx has a low mass because (1) the debris ring had less mass close to Pluto–Charon, (2) it was scattered out of a more massive region of the ring by Nix or Hydra, or (3) after formation at large a, it migrated into an innermost stable orbit.
2.
The orbit of Kerberos lies within a small stable region between the orbits of Nix and Hydra (Youdin et al. 2012). More massive satellites in this region are unstable. From the coagulation and migration calculations, Kerberos is a smaller remnant of the accretion process.
3.
The masses of Nix and Hydra are consistent with the coagulation models. The central cores of either satellite might have survived the giant impact and then accreted debris leftover from the impact. Alternatively, the satellites could form in a disk of small particles orbiting Pluto–Charon.

Aside from informing the next generation of numerical simulations for satellite formation, the New Horizons mission can test several predictions of our calculations.

1.
In our simulations, more massive disks generally produce fewer small satellites. Forming four (or more) satellites orbiting Pluto–Charon implies circumbinary disks with relatively small masses, M_d ≲ 3 × 10²⁰ g. Together with numerical simulations of the orbital stability of the satellites (Youdin et al. 2012), these results require satellite albedos, A ≈ 0.4–1. This range for A is similar to the measured A ≈ 0.5 for Pluto and Charon (see Marcialis et al. 1992; Roush et al. 1996; Brown 2002; Buie et al. 2010a, 2010b, and references therein) and much larger than the A for many Kuiper belt objects (e.g., Stansberry et al. 2008; Brucker et al. 2009). Direct measurements of A by New Horizons will test these constraints.
2.
All of the hybrid calculations leave several very small satellites orbiting within an extended disk of debris beyond the orbit of Hydra. The radii of typical leftovers—R ≲ 1–3 km—imply optical magnitudes of 28 or fainter. Testing this prediction is a challenge for Hubble Space Telescope (Showalter et al. 2011, 2012), but satellites this size or smaller are easily visible during the flyby of New Horizons. For objects with R ≈ 1–3 km, the most likely range of semimajor axes is a ≈ 70–90 R_P. Smaller satellites might have a ≈ 100–200 R_P. Both sets of small satellites should have small inclinations relative to the Pluto–Charon orbital plane, similar to the observed i for the known satellites. If smaller particles in a leftover debris disk (see the next item) have sufficient mass, these outer satellites should also have circular orbits. Smaller masses in leftovers imply larger e for small outer satellites. We plan numerical calculations to test the longevity of these satellites in the tidal field of the Sun–Neptune system.
3.
Estimates for the spreading time in Section 3.3 suggest scattering will produce ensembles of smaller particles, R ≲ 0.01–10 m, in an extended disk well beyond the orbits of the known satellites. For a mass of ∼10¹⁵ g in leftover small particles with typical radius R, the predicted optical depth of the disk is τ_l ∼ 10⁻¹⁰R⁻¹ at semimajor axes a ≈ 70–200 R_P. Larger (smaller) disk masses imply proportionately larger (smaller) optical depths. Although this mass is significant, the predicted optical depth is much smaller than the τ ≲ 10⁻⁷ estimated for clouds of 1–100 μm particles produced by impacts of small Kuiper belt objects on the known satellites (Poppe & Horányi 2011). New Horizons measurements of the optical depths of small particles surrounding the known satellites and at larger distances should place interesting constraints on the long-term evolution of the satellite system and its interactions with the Kuiper belt.
4.
Measuring the amount of material at a ≳ 100 R_P tests formalisms for the capture frequency (Section 3.2). In our estimates, most of the mass lies at large semimajor axes, a ≈ 1000 R_P. Failure to identify any material at a ≈ 100–1000 R_P suggests that capture is less frequent than our estimates. In principle, the amount of material at these distances provides a measure of the capture efficiency.
5.
Composition data from New Horizons should also test the giant impact scenario. If the relative contamination from captured material (Section 3.2) or impacting Kuiper belt objects (e.g., Stern et al. 2006; Poppe & Horányi 2011) is small, satellites should have roughly the same composition as Charon. The broad variety of physical characteristics among Kuiper belt objects (see Barucci et al. 2008) suggests that satellites with significant contributions of material from the Kuiper belt should have very different compositions from satellites grown from the debris of a giant impact (see also Stern 2009). Thus, direct measurement of the compositions of individual satellites and gradients in the composition among any smaller satellites and debris should constrain scenarios for satellite formation.
6.
Observed shapes of satellites might also test formation models. In general, large fragments from the impact should be more irregular than satellites grown from a circumbinary disk/ring of much smaller particles (see Tanga et al. 2009a, 2009b, and references therein).

To develop predicted configurations for the Pluto–Charon system during a New Horizons encounter, we ran n-body simulations of Pluto–Charon, the known small satellites, and several predicted satellites embedded within a disk of 0.5–2 × 10⁶ massless tracer particles. On 10² yr timescales (roughly 5000 orbits of the inner binary), the larger satellites clear many tracers from their orbits (Figure 9, left panel). On 10⁴ yr timescales, the four inner satellites clear most of the tracers from the inner disk (Figure 9, right panel). Although the three outer satellites can clear tracers along their orbits, most tracers remain between satellite orbits in the outer disk.

**Figure 9.** Predicted configuration for the Pluto–Charon system. The Pluto–Charon binary (represented by the two largest white disks), the four small satellites—Styx, Nix, Kerberos, and Hydra (represented by the four small white disks), and three smaller satellites (represented by the green disks) lie within an extended ensemble of solid particles shown as small blue dots. Left panel: system configuration after 10² yr of a computer simulation with two million massless tracer particles surrounding the known and predicted moons. Right panel: system configuration after 10⁴ yr of a computer simulation with 0.5 million massless tracer particles. On short timescales, satellites remove tracer particles along their orbits. On long timescales, satellites remove nearly all tracers in the inner disk; low-mass satellites in the outer disk begin to clear their orbits. On timescales of ≳ 10⁶ yr, satellites will clear the inner disk of small particles (see also Youdin et al. 2012); collisional erosion and scattering will probably remove most of the particles in the outer disk.
Download figure:
Standard image High-resolution image

On much longer (1–10 Myr) timescales, the inner satellites clear their orbits completely (see also Youdin et al. 2012). In the outer disk, a few small satellites with radii of a few km cannot clear more than a few hundred km around their orbits. Thus, many tracers remain. The predicted optical depth and extent of the disk at late times depend on (1) the number of small satellites beyond Hydra, (2) the outcomes of collisions between small disk particles, and (3) the amount of mass in small particles produced from recent impacts on the satellites (e.g., Stern 2009; Poppe & Horányi 2011) and capture from the Kuiper belt. If (1) there are few small satellites, (2) collisions damp and do not destroy small disk particles, and (3) captures and impacts are common, we expect New Horizons to detect a modest disk of small particles with rings and gaps similar to those shown in the right panel of Figure 9. If, however, (1) there are many small satellites beyond Hydra's orbit, (2) collisions destroy small disk particles, and (3) captures and impacts are relatively rare, the New Horizons encounter may reveal little or no disk material. Once New Horizons has mapped the Pluto–Charon system completely, improved numerical simulations will enable a more complete examination of the formation and history of the satellite system.

5. SUMMARY

We describe a theoretical framework for the origin of Pluto's small satellites in the context of a giant impact that produces the Pluto–Charon binary (e.g., Canup 2011). Our main results follow.

1.
Throughout the early history of the solar system, giant impacts capable of producing the Pluto–Charon binary have a frequency of 1 per 100–300 Myr. These impacts leave behind enough debris to enable satellite formation (Canup 2011). Once the debris lies within a ring, angular momentum transport among ring particles may enable the ring to spread into a disk.
2.
In a ring or a disk of small particles, R ≲ 0.1 km, surrounding Pluto–Charon, satellites with R ≈ 10–80 km form rapidly on a timescale of 10³–10⁴ yr. The formation time is similar to the migration timescale.
3.
After the impact, objects with R ≈ 1–10 km might survive collisional disruption during the ejection process. Calculations of disks with small particles and a few seeds produce fewer satellites with larger masses than calculations with small particles alone.
4.
Radial migration encourages mergers of small satellites within an evolving disk surrounding Pluto–Charon.
5.
Scattering events among differentially migrating satellites sometimes place objects on orbits with pericenters close to the central binary. Unless interactions with disk particles can damp the orbit, the Pluto–Charon binary ejects these satellites from the system.
6.
For identical masses of the central planet(s), the surrounding disk, and a small satellite embedded in the disk, migration rates around a binary planet differ slightly from those around a single planet. For Pluto–Charon, the drift rate of a satellite the size of Styx may be as high as 5 km yr⁻¹.
7.
The current number of small Pluto–Charon satellites strongly favors low-mass disks (see also Youdin et al. 2012). In our calculations, disks with masses M_d ≈ 3–10 × 10¹⁹ g have the best chance at producing 3–5 satellites with orbital semimajor axes a ≈ 40–60 R_P. Matching observations requires Pluto's satellites to have albedos A ≈ 0.4–1. New Horizons data will test this prediction.
8.
Model satellites produced in formation calculations have orbital i close to those observed, but the typical e is generally too large. Orbital migration can reduce e to acceptable levels.
9.
The calculations predict a few very small, R ≲ 1–3 km, satellites and an extended disk of even smaller particles, R ≈ 1–100 cm, beyond the current orbit of Hydra. Detecting these satellites and the disk from the ground is very challenging. If they are present, New Horizons should detect them easily.

Our results demonstrate that numerical calculations can produce simulated satellite systems with properties reasonably close to those observed. Prior to the New Horizons flyby of Pluto–Charon, we expect to refine the predictions considerably. Once New Horizons probes the masses, orbital architecture, and composition of the Pluto–Charon binary, a rich interplay between the data and the numerical simulations will enable a much more robust theory for the formation of satellites (planets) in binary planetary (stellar) systems.

We acknowledge generous allotments of computer time on the NASA "discover" cluster and the SI "hydra" cluster. Advice and comments from M. Geller and A. Youdin greatly improved our presentation. We thank two anonymous referees for thoughtful comments that honed our analysis. Portions of this project were supported by the NASA Astrophysics Theory and Origins of Solar Systems programs through grant NNX10AF35G and the NASAOuter Planets Program through grant NNX11AM37G.

THE FORMATION OF PLUTO'S LOW-MASS SATELLITES

Article metrics

Permissions

Author e-mails

Author affiliations