A Pluto-Charon Sonata III. Growth of Charon from a Circum-Pluto Ring of Debris

TL;DR

This paper investigates the formation of the Pluto-Charon system through numerical simulations of Charon's growth within a circum-Pluto debris ring, challenging the graze-and-merge impact model and exploring satellite formation scenarios.

Contribution

The study provides new insights into Charon's rapid formation and the ejection of large debris, offering an alternative to existing impact models for Pluto-Charon formation.

Findings

Charon analogs form within 30-100 days in simulations.

Large debris (145-230 km) are ejected by the binary system.

Small particles may survive ejection of larger remnants.

Abstract

Current theory considers two options for the formation of the Pluto-Charon binary (Canup 2005, 2011; Desch 2015). In the `hit-and-run' model, a lower mass projectile barely hits the more massive Pluto, kicks up some debris, and remains bound to Pluto (see also Asphaug et al. 2006). In a `graze-and-merge' scenario, the projectile ejects substantial debris as it merges with Pluto (see also Canup 2001). To investigate the graze-and-merge idea in more detail, we consider the growth of Charon-mass objects within a circum-Pluto ring of solids. Numerical calculations demonstrate that Charon analogs form rapidly within a swarm of planetesimals with initial radii of 145-230 km. On time scales of roughly 30-100 days, newly-formed Charon analogs have semimajor axes, a = 5-6 Pluto radii, and orbital eccentricities, e = 0.1-0.3, similar to Charon analogs that remain bound after hit-and-run collisions with Pluto. Although the early growth of Charon analogs generates rings of small particles at a = 50-275 Pluto radii, ejection of several 145-230 km leftovers by the central Pluto-Charon binary removes these small solids in 10-100 yr. Simple estimates suggest small particles might survive the passage of 10-20 km objects ejected by the central binary. Our results indicate that the Pluto-Charon circumbinary satellite system was not formed by a graze-and-merge impact when the formation of Charon within a circum-Pluto disk leads to the ejection of several 100-200 km particles through the orbital plane of the Pluto-Charon binary. If a growing Charon ejects only much smaller particles, however, graze-and-merge impacts are a plausible formation channel for the Pluto-Charon binary and an ensemble of small, circumbinary satellites.

Tables

Table 1. Table 1: Results for n 𝑛 n -body calculations a a The columns list the ratio of the initial mass of the swarm relative to the mass of Charon, M 0 / m C subscript 𝑀 0 subscript 𝑚 𝐶 M_{0}/m_{C} , the initial radius of solids orbiting Pluto r 0 subscript 𝑟 0 r_{0} , the ratio of the final to initial mass for Pluto ( μ 1 = m 1 / m P subscript 𝜇 1 subscript 𝑚 1 subscript 𝑚 𝑃 \mu_{1}=m_{1}/m_{P} ), the ratio of the mass of the Charon analog to Charon’s mass ( μ 2 = m 2 / m C subscript 𝜇 2 subscript 𝑚 2 subscript 𝑚 𝐶 \mu_{2}=m_{2}/m_{C} ), and the average semimajor axis and eccentricity of the Charon analog and their dispersions.

$M_0/m_C$	$r_0$ (km)	${\mu }_1$	${\mu }_2$	$a_2$ ($r_P$)	$e_2$
0.85	230	1.02 $\pm$ 0.02	0.61 $\pm$ 0.07	6.25 $\pm$ 1.75	0.26 $\pm$ 0.12
0.85	185	1.02 $\pm$ 0.02	0.54 $\pm$ 0.07	4.90 $\pm$ 0.97	0.20 $\pm$ 0.06
0.85	145	1.02 $\pm$ 0.02	0.48 $\pm$ 0.10	5.02 $\pm$ 1.46	0.19 $\pm$ 0.07
1.25	230	1.02 $\pm$ 0.01	0.79 $\pm$ 0.11	5.39 $\pm$ 1.40	0.20 $\pm$ 0.08
1.25	185	1.03 $\pm$ 0.02	0.75 $\pm$ 0.10	6.74 $\pm$ 2.64	0.20 $\pm$ 0.06
1.25	145	1.02 $\pm$ 0.02	0.72 $\pm$ 0.14	5.39 $\pm$ 1.62	0.21 $\pm$ 0.08
1.75	230	1.04 $\pm$ 0.03	1.04 $\pm$ 0.31	6.68 $\pm$ 1.29	0.28 $\pm$ 0.10
1.75	185	1.05 $\pm$ 0.04	1.02 $\pm$ 0.20	6.52 $\pm$ 1.62	0.26 $\pm$ 0.10
1.75	145	1.04 $\pm$ 0.02	0.98 $\pm$ 0.14	6.04 $\pm$ 1.83	0.23 $\pm$ 0.06

Table 2. Table 2: Tracer statistics at 150–300 d a a The columns list the physical variables and the fraction of calculations with the indicated ranges of f t subscript 𝑓 𝑡 f_{t} for listed physical variable at an evolution time of 150–300 d.

Variable	$f\left(f_t=0\right)$
$m_2\le$ 0.7 $m_C$	0.14	0.15	0.15	0.09	0.48
$m_2>$ 0.7 $m_C$	0.52	0.16	0.06	0.05	0.21
$a_2\le 6r_P$	0.36	0.19	0.10	0.04	0.32
$a_2>6r_P$	0.31	0.12	0.10	0.10	0.38
$e_2\le 0.22$	0.32	0.19	0.07	0.05	0.37
$e_2>0.22$	0.35	0.12	0.13	0.09	0.32
$e_0$ = 0.1	0.33	0.16	0.08	0.07	0.35
$e_0$ = 0.2	0.32	0.12	0.12	0.07	0.36
$e_0$ = 0.3	0.31	0.17	0.13	0.05	0.34
$e_0$ = 0.4	0.37	0.16	0.07	0.07	0.33
$r_0$ = 145 km	0.28	0.16	0.08	0.06	0.42
$r_0$ = 185 km	0.32	0.17	0.11	0.08	0.33
$r_0$ = 230 km	0.41	0.15	0.11	0.06	0.28
$M_0$ = 0.85 $m_C$	0.19	0.14	0.13	0.07	0.47
$M_0$ = 1.25 $m_C$	0.16	0.16	0.09	0.07	0.52
$M_0$ = 1.75 $m_C$	0.66	0.17	0.08	0.06	0.04

Table 3. Table 3: Correlations of Physical Variables with f t subscript 𝑓 𝑡 f_{t} at 150–300 d a a The columns list the physical variables, the correlation coefficients for the Pearson ( r 𝑟 r ), Spearman rank-order ( r S subscript 𝑟 𝑆 r_{S} ) and Kendall ( τ 𝜏 \tau ) tests along with the probabilities for a lack of correlation between the listed variable and f t subscript 𝑓 𝑡 f_{t} at 150–300 d.

Variable	$r$	$r_S$	$\tau$	$p_r$	$p_S$	$p_{\tau }$
$m_2$	$-$0.35	$-$0.59	$-$0.35	$4.0\times {10}^{-15}$	0.0	$4.1\times {10}^{-30}$
$a_2$	$-$0.03	$+$0.06	$+$0.10	$5.1\times {10}^{-1}$	$1.8\times {10}^{-1}$	$9.3\times {10}^{-4}$
$e_2$	$-$0.02	$-$0.09	$+$0.00	$6.0\times {10}^{-1}$	$4.7\times {10}^{-2}$	$9.0\times {10}^{-1}$
$e_0$	$+$0.01	$-$0.31	$-$0.02	$9.1\times {10}^{-1}$	$8.8\times {10}^{-12}$	$4.8\times {10}^{-1}$
$r_0$	$-$0.14	$-$0.63	$-$0.11	$1.7\times {10}^{-3}$	0.0	$2.7\times {10}^{-4}$
$M_0$	$-$0.31	$-$0.91	$-$0.35	$2.1\times {10}^{-12}$	0.0	$5.9\times {10}^{-30}$

Equations

$$r_H={\left(\frac{m_1+m_2}{3m_P}\right)}^{1/3}\bar{a}=\frac{r_1}{r_P}{\left(\frac{1+q}{3}\right)}^{1/3}\bar{a},$$

$$\alpha \gtrsim {\left(\frac{3{\rho }_P}{{\rho }_0}\right)}^{1/3}\frac{1+q^{1/3}}{(1+q{)}^{1/3}}.$$

$$\frac{de}{dt}=\frac{e}{t_{damp}}$$

$$\frac{dı}{dt}=\frac{ı}{t_{damp}},$$

$$t_{damp}=1250{\left(\frac{a}{r_P}\right)}^{3/2}\left(\frac{r}{1\mathrm{cm}}\right)\left(\frac{10^{22}\mathrm{g}}{M_t}\right)\sec .$$

$$E\approx \frac{3G{m}_C^2}{5r_C},$$

$$T=\frac{E}{3Nk}$$

Full text

planets and satellites: dynamical evolution — planets and satellites: formation — dwarf planets: Pluto

1 INTRODUCTION

Inside the protosolar nebula, Pluto and other trans-Neptunian objects (TNOs) begin their lives as disparate small particles with a broad range in sizes and (probably) semimajor axes (e.g., Armitage, 2013; Birnstiel et al., 2016; Lammer & Blanc, 2018; Raymond et al., 2018). Within a few thousand yr, these particles grow to cm and larger sizes (e.g., Dullemond & Dominik, 2005; Brauer et al., 2008; Birnstiel et al., 2010; Drażkowska & Dullemond, 2018; Lenz et al., 2019). During the next 0.1–0.3 Myr (e.g., Najita & Kenyon, 2014), agglomeration and concentration processes generate km-sized or larger planetesimals (e.g., Youdin & Goodman, 2005; Johansen et al., 2006, 2007; Lambrechts & Johansen, 2012; Birnstiel et al., 2016; Arakawa & Nakamoto, 2016; Blum, 2018). In some theories, planetesimals become planets by accreting cm-sized pebbles (e.g., Kenyon & Bromley, 2009; Bromley & Kenyon, 2011; Chambers, 2016; Johansen & Lambrechts, 2017; Lenz et al., 2019, and references therein). In others, successive collisional mergers of planetesimals and larger objects eventually produce several planets (e.g., Stern & Colwell, 1997; Kenyon & Bromley, 2008, 2010; Raymond et al., 2011; Kenyon & Bromley, 2012; Hansen & Murray, 2012; Schlichting et al., 2013, and references therein).

Numerical calculations suggest the transformation of a swarm of planetesimals into a planetary system is often a chaotic process, with numerous ‘giant impacts’ between planet-mass objects (e.g., Agnor et al., 1999; Chambers, 2001; Asphaug et al., 2006; Genda et al., 2012; Chambers, 2013; Asphaug, 2014; Quintana et al., 2016). Many of these binary collisions simply add to the ever-growing mass of the newly-formed planet. Others may completely disrupt the more massive object of the pair. Sometimes, the geometry allows the impactor to eject a cloud of debris or to survive, leading to the formation of a binary planet similar to the Earth–Moon (e.g., Hartmann & Davis, 1975) or the Pluto–Charon system (e.g., McKinnon, 1989).

Although rocky and icy planets form in similar ways, collision outcomes depend on the geometry and the velocity of the impact (e.g., Benz & Asphaug, 1999; Leinhardt & Richardson, 2002; Asphaug et al., 2006; Leinhardt & Stewart, 2012; Emsenhuber & Asphaug, 2019; Arakawa et al., 2019). Among planetary embryos with $a\lesssim$ 2 au and $e\gtrsim$ 0.1, collisions have typical relative velocities $v_{rel}\gtrsim$ 3–10 $\rm km~{}s^{-1}$. Despite their large relative velocities, rocky planetary embryos in the terrestrial zone also have large binding energies. Thus, terrestrial planet formation is fairly robust; Earth-mass or larger planets grow from much less massive objects on time scales of 10–100 Myr (Kenyon & Bromley, 2006; Chambers, 2008; Raymond et al., 2011; Hansen & Murray, 2012; Genda et al., 2015). At $a\gtrsim$ 20 au, the relative velocities of icy planetary embryos are smaller, $v_{rel}\lesssim$ 1–3 $\rm km~{}s^{-1}$. With much smaller binding energies, however, it is much easier for collisions to prevent the growth of massive planets (e.g., Kenyon & Bromley, 2008, 2010, 2012). Dynamical interactions with gas giants at smaller $a$ can also halt growth and limit the production of Earth-mass or larger planets beyond 20 au (e.g., Levison & Morbidelli, 2003; Levison et al., 2008; Morbidelli et al., 2008).

To explore the physical conditions required for a giant impact to produce the Pluto–Charon system, Canup (2005, 2011) performed an extensive suite of SPH calculations with a range of compositions, masses, spins, and impact angles for a collision between proto-Pluto and proto-Charon. In ‘graze-and-merge’ collisions with impact velocities $v_{i}$ comparable to the escape velocity $v_{esc}$ ($v_{c}\approx$ 1.0–1.2 $v_{esc}$), the impactor merges with Pluto and ejects a massive disk of debris (for other applications, see Leinhardt et al., 2010; Canup et al., 2013; Nakajima & Stevenson, 2014; Barr, 2016; Barr & Bruck Syal, 2017, and references therein). Collisions within the disk eventually generate a lower density satellite with the mass of Charon. ‘Hit-and-run’ collisions with $v_{c}\approx$ 1.0–1.2 $v_{esc}$ allow a proto-Charon to survive and remain bound to Pluto (for other examples of hit-and-run collisions, including those where the impactor escapes the target, see Asphaug et al., 2006; Chambers, 2013; Emsenhuber & Asphaug, 2019). Many hit-and-run collisions generate a disk of debris with a mass sufficient to produce the four known circumbinary satellites. From simulations of planet formation within the protosolar nebula, both types of collisions are probably common (e.g., Kenyon & Bromley, 2012, 2014).

Here, we describe a suite of numerical simulations designed to establish whether a Charon-mass satellite can form out of a disk of debris orbiting a Pluto-mass planet. Although several analyses of the Pluto–Charon system currently appear to preclude the possibility of Charon formation from a graze-and-merge collision (e.g., Canup, 2005, 2011; McKinnon et al., 2017; Bierson et al., 2018; Stern et al., 2018, and references therein), it is useful to explore the evolution of a massive, circum-Pluto disk of debris to test alternate formation models (e.g., Desch, 2015; Desch & Neveu, 2017). Our goal is to learn whether the physical conditions required for the growth of a Charon-mass satellite are consistent with the structures generated by SPH calculations of giant impacts (Canup, 2005, 2011).

Tracking the growth of Charon within a circum-Pluto disk also places limits on plausible hit-and-run models. As we show below, Charon analogs reach their final mass on short times scales, $\sim$ weeks. At this point, the Charon analog lies within an expanding cloud of debris on an orbit similar to the published end states of hit-and-run impacts. The evolution of this debris as a function of initial conditions helps us to understand (i) the amount of material Pluto and Charon might accrete well after the impact and (ii) outcomes that are more (or less) favorable to the formation of the circumbinary satellites as the Pluto–Charon binary expands from tidal interactions (see also Walsh & Levison, 2015).

Following a description of numerical methods (§2), we summarize the growth of Charon analogs from a disk of massive planetesimals (§3) and the dynamical evolution of massless tracers responding to the growth of a massive satellite (§4). After discussing how the results impact our understanding of the formation and evolution of the Pluto–Charon system (§5), we conclude with a brief summary.

2 NUMERICAL METHODS

To track the evolution of circumplanetary debris in the Pluto system, we run a series of numerical simulations with Orchestra, a parallel C++/MPI hybrid coagulation + $n$-body code that follows the accretion, fragmentation, and orbital evolution of solid particles ranging in size from a few microns to thousands of km (Kenyon, 2002; Bromley & Kenyon, 2006; Kenyon & Bromley, 2008; Bromley & Kenyon, 2011, 2013; Kenyon & Bromley, 2016; Kenyon et al., 2016). The ensemble of codes within Orchestra includes a multi-annulus coagulation code, an $n$-body code, and radial diffusion codes for solids and gas. Several algorithms link the codes together, enabling each component to react to the evolution of other components.

Here, we use the $n$-body code to track the time evolution of the orbits of massive particles and massless tracers orbiting Pluto. In the calculations, the tracers serve as analogs of small particles of debris with negligible mass compared to the mass of Charon, $m_{C}$ = $1.586\times 10^{24}$ g. Pluto has initial mass $m_{P}$ = $1.303\times 10^{25}$ g and radius $r_{P}$ = 1183 km (Stern et al., 2015; Nimmo et al., 2017). Calculations begin with an ensemble of massive particles with total mass $M_{0}$ and $N_{t}$ = 14000 tracers in orbit around Pluto. Each massive particle has initial mass $m_{0}$, radius $r_{0}$, mass density $\rho_{0}$ = 1.75 $\rm g~{}cm^{-3}$, and a semimajor axis, $a_{0}\approx$ = 3–11 $r_{P}$. Roughly half of the tracers are randomly placed in semimajor axis among the massive particles. The rest have a random $a$ = 10–55 $r_{P}$.

For the adopted physical properties of Pluto and the orbiting solids, the fluid Roche limit lies at $a_{R}\approx$ 2.5 $r_{P}$ (e.g., Aggarwal & Oberbeck, 1974; Weidenschilling et al., 1984; Canup & Esposito, 1995; Hyodo & Ohtsuki, 2014). Material at the Roche limit has an orbital period $P\approx$ 10 hr. Solid particles with $a\lesssim a_{R}$ relax to a non-spherical equilibrium shape that depends on $a$ and the material strength (e.g., Holsapple & Michel, 2006, 2008). Although collision rates depend on geometric cross-sections, we ignore modest differences in rates for ellipsoidal particles in this pilot study. For the 150–250 km solid particles considered here, self-gravity sets the strength against tidal stresses. These particles are probably stable against tides at $a\gtrsim 0.7a_{R}\gtrsim$ 2 $r_{P}$. Although orbital motion can place solids inside this limit, 150–250 km solids with large gravity are more likely to experience surface and perhaps internal fractures than tidal shredding (e.g., Holsapple & Michel, 2006; Sharma et al., 2006; Holsapple & Michel, 2008; Sharma, 2009; Quillen et al., 2016). During the much longer period of time spent outside the Roche limit, they probably recover their original structure. Thus, we ignore the possibility of tidal shredding when solids lie inside the Roche limit.

Collision outcomes also depend on tidal stresses (e.g., Canup & Esposito, 1995; Ohtsuki et al., 2013; Hyodo & Ohtsuki, 2014). Two colliding particles with masses $m_{1}$ and $m_{2}=qm_{1}$ cannot merge if the sum of their radii is larger than their mutual Hill sphere,

[TABLE]

where $\bar{a}=0.5(a_{1}+a_{2})$. Setting $\bar{a}=\alpha~{}r_{P}$ and $\tilde{r}=(r_{1}+r_{2})/r_{H}\lesssim 1$, accretion is possible when (Hyodo & Ohtsuki, 2014):

[TABLE]

With $\rho_{P}\approx\rho_{0}$, $\alpha\approx$ 2.3 for $q$ = 1; $\alpha\approx$ 1.4 for $q$ = 0 (see also Weidenschilling et al., 1984). Canup & Esposito (1995) derive a similar criterion. In our calculations, collisions between two $n$-bodies inside this limit are statistically unlikely. After the conclusion of each calculation, we identify $n$-bodies with positions inside 2.5 $r_{P}$ and verify that they do not collide with other $n$-bodies while at a distance $r\lesssim$ 2.5 $r_{P}$. For the suite of simulations discussed here, no $n$-body violates these constraints.

Compared to the final states of SPH calculations for graze-and-merge collisions (for example, Fig. 1 of Canup, 2005), our initial conditions are somewhat more radially extended. Roughly 24 hr after the collision, the SPH calculations show solids extending from close to the surface of the central planet, 1–2 $r_{P}$, to $\sim$ 10–15 $r_{P}$; 70% of this material initially lies outside the Roche limit. Because the central planet has a short rotational period, $\sim$ 2.5 hr, solids that orbit synchronously with the planet lie well inside the Roche limit at $a_{c}\approx$ 1.05 $r_{P}$. Over time, material orbiting the central planet eventually evolves into a disk (Canup, 2005). To avoid calculating the more complicated evolution of solids inside the Roche limit (e.g., Canup & Esposito, 1995; Ohtsuki et al., 2013; Hyodo & Ohtsuki, 2014), the disks considered here have an inner radius $a_{in}\approx$ 3 $r_{P}$, an outer radius $a_{out}\approx$ 11 $r_{P}$, and a surface density distribution $\Sigma\propto a^{-n}$ with $n\approx$ 1.7–1.8. Although Canup (2005) does not quote $\Sigma(a)$ for circum-Pluto disks, Ida et al. (1997) consider similarly shallow surface density distributions in $n$-body calculations of the formation of the Moon from an extended disk surrounding the Earth. Consistent with energy equipartition through gravitational scattering as circum-Pluto material evolves into a disk (Brahic, 1976; Stewart & Wetherill, 1988; Ohtsuki et al., 2002), massive particles and tracers have the same initial eccentricity $e_{0}$ and inclination $e_{0}$/2 (see also Ida et al., 1997).

Throughout a calculation, the adaptive integrator adjusts time steps to resolve collisions (i) among all massive particles (including Pluto) and (ii) between tracers and other massive particles. Although tracers may merge with massive particles, the massive particles do not respond to the tracers. When a tracer or massive particle has $a\gtrsim a_{max}$ or $e>1$, the code flags it as ejected and freezes its position and velocity for the remainder of the calculation. In these calculations, $a_{max}$ = 5400 $r_{P}$, slightly smaller than the Hill radius of the Pluto–Charon binary planet, $R_{H,PC}\approx$ 6750 $r_{P}$.

To avoid starting with a gravitationally unstable ring of massive particles, calculations begin with a sufficiently large $e_{0}$. Defining the initial velocity dispersion $v_{0}$, angular velocity $\Omega_{0}$, and surface density $\Sigma_{0}$, a ring is unstable when $v_{0}\Omega_{0}<\pi G\Sigma_{0}$ (e.g., Chiang & Youdin, 2010). To cast this condition in terms of $m_{P}$, $e_{0}$, and $M_{0}$, we set (i) $M_{0}\ =\pi\Sigma_{0}a_{0}^{2}$, where $a_{0}$ is a characteristic semimajor axis for the swarm, and (ii) $v_{0}\ =ev_{K}$, where $v_{K}$ is the local orbital velocity. The stability condition is then $e_{0}\gtrsim M_{0}/m_{P}$. Requiring that the mass reservoir orbiting Pluto be sufficient to produce a Charon ($M_{0}\gtrsim 0.1m_{P}$), $e_{0}\gtrsim$ 0.1.

In this set of calculations, we ignore fragmentation when two large objects physically collide and merge. For a collision involving two icy objects with $r$ = 100 km on orbits with $a\sim$ 5–6 $r_{P}$ and $e\sim 0.2$, the typical mass lost to debris is less than 1% (Kenyon & Bromley, 2014, 2015, 2017, and references therein). To simplify this initial set of calculations, we ignore the impact of this debris on the evolution of massive and massless particles in the $n$-body code. After each calculation, we verify that collisions between $n$-bodies are at velocities insufficient to eject more than 1% of the combined mass of colliding objects in debris that escapes the merged object. This condition is more stringent than the velocity limits required for mergers inside the Roche limit (e.g., Hyodo & Ohtsuki, 2014).

To survey outcomes as a function of initial conditions, we consider a range of $m_{0}$, $e_{0}$, and $M_{0}$. The range of total mass $M_{0}$ – 0.85 $m_{C}$ to 1.75 $m_{C}$ – brackets the current mass of Charon. We choose $e_{0}$ = 0.1, 0.2, 0.3 or 0.4 to cover many outcomes from detailed SPH calculations of giant impacts (Canup, 2005, 2011). Although rings of particles with $e_{0}\ \approx$ 0.1 are formally unstable for $M_{0}$ = 1.75 $m_{C}$, these calculations allow us to examine trends in outcomes as a function of $e_{0}$ in more detail. The three choices for $r_{0}$ – 145 km, 185 km, and 230 km – yield a manageable number of initial massive particles and allow us to learn whether outcomes depend on $r_{0}$. With 56 cores, each calculation requires a few hours to several days of computer time. We perform 12–15 calculations with each setup to test the sensitivity of outcomes to shot noise.

Comparing the final states of massless tracers with the Pluto–Charon circumbinary satellites requires (i) knowledge of the physical properties of the satellites and (ii) an understanding of the processes required to convert debris left over from the formation of Charon into a few small satellites. From observations with HST (Brozović et al., 2015; Showalter & Hamilton, 2015), the satellites have nearly circular orbits in the orbital plane of the Pluto–Charon binary with semimajor axes $a_{S}\approx$ 36.06 $r_{P}$ (Styx), $a_{N}\approx$ 41.16 $r_{P}$ (Nix), $a_{K}\approx$ 48.84 $r_{P}$ (Kerberos), and $a_{H}\approx$ 54.72 $r_{P}$ (Hydra). Direct $n$-body calculations place robust constraints on the masses of Nix and Hydra (Kenyon & Bromley, 2019b): $m_{N}\lesssim 4.7\times 10^{19}$ g and $m_{H}\lesssim 5.5\times 10^{19}$ g. Although the masses of Styx and Kerberos have larger uncertainties, the $n$-body calculations suggest total masses of the circumbinary satellite system $M_{SNKH}\lesssim 10^{20}$ g (Kenyon & Bromley, 2019b). Together with size measurements from New Horizons (Weaver et al., 2016), these results suggest the circumbinary satellites have mass densities smaller than Pluto and Charon, $\rho_{SNKH}\lesssim$ 1.5 $\rm g~{}cm^{-3}$.

3 RESULTS: EVOLUTION OF MASSIVE PARTICLES

In all calculations, the sequence of merger events follows one of several well-defined paths. Often, several early collisions generate 3–4 objects with masses equal to 2 $m_{0}$. These larger objects then accrete most of the other smaller planetesimals. Sometimes, binary collisions among the $N$ initial objects generate $\sim N/2$ larger objects; mergers among this swarm eventually produce several large object(s) orbiting Pluto. Rarely, one object rapidly merges with two other objects in the swarm and then accumulates material from the rest of the swarm. As the largest objects grow in mass, they often eject several small planetesimals from the system. Once only 3–4 $n$-bodies remain in any of these scenarios, they either merge to form a single massive object, collide with Pluto, or are ejected from the system.

As the largest $n$-bodies merge and grow, the dynamical evolution of the full ensemble is rather simple. Within a few orbits, short-range scattering events raise (lower) the typical $e$ ($\imath$) of most $n$-bodies. After several mergers, dynamical friction tends to reduce the $e$ of the more massive $n$-bodies and raise the $e$ of the less massive $n$-bodies. Additional mergers accelerate this process, until several $n$-bodies have $e\approx$ 1. During the next 1–2 orbits, these $n$-bodies either merge with another $n$-body or are ejected. Once a single $n$-body has more than 90% of the mass outside Pluto, it rapidly clears its orbit of any remaining $n$-bodies. Across all of the calculations, Pluto and the most massive $n$-bodies eject up to nine low mass particles from the system; 3–4 ejections is typical.

Fig. 1 illustrates the time evolution of the number of $n$-bodies in typical calculations. After the first merger of two $n$-bodies, it takes only a day or two for collisions to reduce the number of $n$-bodies by half. After 10–50 days, only one (sometimes two) large object(s) orbit(s) Pluto. Among all of the calculations, the time scale for the number of $n$-bodies to fall to 0, 1, or 2 ranges from a few days to $\sim$ 50–100 days.

The outcomes of the evolution are fairly insensitive to the initial conditions. Although shot noise is important in every calculation, the time scale to reduce the number of $n$-bodies from $N$ to 1–2 scales inversely with $M_{0}$: more massive collections of $n$-bodies evolve more rapidly than less massive swarms. However, systems with $e_{0}$ = 0.1 evolve as fast as systems with $e_{0}$ = 0.4. Calculations with smaller planetesimals initially evolve somewhat more rapidly than those with larger planetesimals, but this difference is lost after several mergers generate an ensemble of objects with masses 2-3 times larger than the initial mass.

Fig. 2 shows the time evolution of the mass of each $n$-body for the set of calculations in Fig. 1. Initially, all of the $n$-bodies have the same mass. Within a few hours, several have merged with other $n$-bodies and doubled, tripled, or quadrupled in mass. In some cases, one $n$-body gains mass more rapidly than the others and remains the most massive object until the end of the calculation (Fig. 2, middle panel). In others, several $n$-bodies have comparable masses throughout the evolution, until one emerges as most massive (Fig. 2, upper and lower panels). However, the evolution proceeds, after 10–50 days only 0, 1, or 2 $n$-bodies orbit(s) Pluto.

Sometimes, 1–2 objects escape collisions for 10–100 days (Fig. 2, green object in the lower panel and the cyan object in the upper panel). Throughout their evolution, larger objects tend to stir their eccentricities to larger and larger values. Usually, the largest object in the swarm accretes this wayward planetesimal. Sometimes, stirring leads to an ejection or to accretion by Pluto.

Despite the chaotic evolution in mass, the most massive $n$-body has a final orbit with $a\approx$ 5–10 $r_{P}$ (Fig. 3). During the chaotic growth phase when the number of $n$-bodies declines dramatically, the semimajor axes (eccentricity) of each $n$-body vary from $\sim$ 3 $r_{P}$ to 15–25 $r_{P}$ ($\lesssim$ 0.1 to $\gtrsim$ 0.5). Sometimes, interactions between the last 2–3 $n$-bodies leads to impacts on Pluto, ejections from the system, and no surviving companions to Pluto. Usually, however, a single massive companion remains in orbit at 5–10 $r_{P}$ from Pluto.

Fig. 3 also illustrates the different types of semimajor axis evolution for the swarm of $n$-bodies. In the top panel, the larger initial masses of individual objects stir up the orbits of their nearest neighbors, expanding the extent of the swarm in $a$. Although the orbits then begin to cross, few $n$-bodies move from the inner edge of the swarm to the outer edge of the swarm over the course of the evolution. For the calculation shown lower panel, there is less stirring among lower mass objects; orbit crossing is more common. The semimajor axis evolution is more chaotic, with several small objects moving from small $a$ to large $a$ (and sometimes back to small $a$) as the largest objects grow.

In all of the calculations, several objects reach $a\gtrsim$ 20–25 $r_{P}$ with large $e$ (e.g., the green object in the top panel of Fig. 3 at 0.3–0.4 days). These objects are rarely accreted by Pluto or another object in the swarm. Usually, they are ejected from the system.

Statistics for the full ensemble of calculations reveals several characteristic outcomes (Table 1). Throughout the evolution, Pluto usually suffers one or more large impacts from the orbiting swarm. The typical gain in Pluto’s mass ranges from $\delta m_{1}\approx 0.02~{}m_{P}$ when $M_{0}=0.85~{}m_{C}$ to $\delta m_{1}\approx 0.05~{}m_{P}$ when $M_{0}=1.75~{}m_{C}$. The dispersion in this mass gain is comparable to the mass gain. In some calculations, Pluto accretes material roughly equivalent to the mass of Charon.

The final mass $m_{2}$ of the largest $n$-body orbiting Pluto correlates with the initial mass in solids. When $M_{0}=0.85~{}m_{C}$, $m_{2}$ is roughly half the mass of Charon. As the initial mass in solids grows, $m_{2}$ also grows, reaching 5% larger than the mass of Charon when $M_{0}=1.75~{}m_{C}$. The ratio $m_{2}\ /M_{0}$ also correlates with $M_{0}$, declining from 0.65 when $M_{0}$ is small to 0.57 when $M_{0}$ is large. Calculations with small $M_{0}$ are more efficient in concentrating the initial mass into a single large $n$-body orbiting Pluto.

The final mass also correlates with the initial sizes of the solids. Ensembles of small solids are less efficient at producing a massive $n$-body than ensembles of large solids (Table 1). When $r_{0}$ is smaller, it is easier for the largest $n$-body in the swarm to scatter the smallest $n$-body away from the rest of the swarm. Repeated dynamical interactions between the largest and smallest $n$-body often result in the ejection of the smallest $n$-body. When $r_{0}$ is larger, there is less of a contrast between the masses of the smallest and largest $n$-bodies, limiting the number of ejections. Fewer dynamical ejections allow the largest $n$-body to reach a larger $m_{2}$.

Comparisons among the mass distributions suggest correlations with $r_{0}$ but not $e_{0}$ (Fig. 4). When $M_{0}$ = 0.85 $m_{C}$ (Fig. 4, lower panel), the mass distributions for $r_{0}$ = 145 km (orange curve), $r_{0}$ = 185 km (green curve), and $r_{0}$ = 230 km (purple curve) clearly differ from one another. Using a K–S test (Press et al., 1992), the probability that these distributions are drawn from the same parent population is negligible, $\lesssim 10^{-5}$ to $10^{-4}$. However, mass distributions for calculations with $M_{0}$ = 1.25 $m_{C}$ are nearly identical, with large K–S probabilities, 30% to 90%, of being drawn from the same parent population. At $M_{0}$ = 1.75 $m_{C}$, the mass distributions for the three $r_{0}$ have a 5% to 10% probability of being drawn from the same parent population.

Despite the variation of $m_{2}$ with initial conditions, there is remarkably little variation of the final semimajor axis $a_{2}$ or eccentricity $e_{2}$. Typically, $a_{2}\ \approx$ 5–7 $r_{P}$ and $e_{2}\approx$ 0.2. With a dispersion of 1.0–2.5 $r_{P}$, the overall range in $a_{2}$ is large, with a minimum of 3 $r_{P}$ and a maximum of 11.5 $r_{P}$. In some calculations, the final orbit of Pluto’s companion is nearly circular, with $e_{2}\approx$ 0.05–0.10. In others, the final orbit is highly eccentric $e_{2}\approx$ 0.4–0.5.

Analyses of the variation of $a_{2}$ or $e_{2}$ with initial conditions yields no strong correlations. There is a weak correlation between $e_{2}$ and $e_{0}$, where systems with smaller $e_{0}$ produce a massive satellite with smaller $e_{2}$. However, the trend is significant at less than the 1.5-$\sigma$ level. The final eccentricity of a massive satellite is insensitive to $M_{0}$ or $r_{0}$; $e_{2}$ ranges from 0.05 to 0.40 with a typical $e_{2}$ $\approx$ 0.2. Among all of the massive satellites produced in the calculations, $a_{2}$ is remarkably uncorrelated with $e_{0}$, $r_{0}$, or $M_{0}$. For any combination of $e_{0}$, $r_{0}$, and $M_{0}$, $a_{2}$ ranges from 3 $r_{P}$ to 10–11 $r_{P}$, with a typical $a_{2}$ $\approx$ 5–6 $r_{P}$.

4 RESULTS: EVOLUTION OF MASSLESS TRACER PARTICLES

The swarm of massless tracer particles evolves in step with the growth of the largest $n$-bodies. At the start of each calculation, the tracers adjust their orbits to the gravitational field of all of the massive objects. Within a few orbits, some tracers impact Pluto or one of the massive $n$-bodies. Strong gravitational interactions lead to the ejection of several others from the system. The rest of the tracers settle down into eccentric orbits around Pluto.

As each system evolves, the orbits of the tracers gradually expand away from the orbits of the $n$-bodies. Early on, many tracers collide with Pluto or another $n$-body. As the evolution proceeds, the small dynamical kicks each tracer receives from the $N$ remaining $n$-bodies slowly scatters them to larger and larger semimajor axes. With no collisional damping, the orbital eccentricities of many tracers gradually grow with time. Eventually, these tracers are ejected. Other tracers are able to maintain a fairly small $e$ and end up on wide orbits around the central Pluto–Charon binary.

Throughout the radial expansion and ejection of tracers, there is little excitation of tracer orbits in $z$. During the first few days of evolution, the vertical scale height of the tracers is roughly constant. As they pass close to massive $n$-bodies, their eccentricities grow; however, their inclinations are largely unchanged. In systems that develop structure in the $x-y$ plane (see below), there is little evidence for structure with $z$. Tracers with large $e$ are ejected at inclinations fairly similar to their starting inclinations.

Fig. 5 illustrates the positions of swarms of tracers in the $x-y$ plane at the end of several calculations. When the Charon analog has a small mass, $m_{2}\ \approx$ 0.6–0.8 $m_{C}$ (upper panel), the ensembles of tracers occupy fairly circular rings with radius of 100–175 $r_{P}$. Aside from a dense set of tracers in the well-defined rings, there is often a halo of tracers on orbits with larger semimajor axes, $a\approx$ 200–275 $r_{P}$. Although the rings of tracers look nearly circular, tracer orbits are eccentric, with typical $e\approx$ 0.2–0.5.

In systems that produce a more massive Charon analog ($m_{2}\ \approx$ 0.9–1.1 $m_{C}$; Fig 5, lower panel), tracers suffer more energetic interactions with the $n$-bodies. Aside from having more frequent impacts with Pluto or other $n$-bodies, these systems eject more tracers. In some cases, the few remaining tracers have a rather chaotic distribution in the $x-y$ plane (lower right panel). In others, snapshots contain partial rings of tracers at distance of 100–275 $r_{P}$ from the central binary (lower left and lower middle panels). Because the tracers from the empty portion of the ring have been ejected, the partial ring pattern appears to rotate as the evolution proceeds.

Throughout many calculations, the distribution of tracers in the $x-y$ plane shows a clear spiral pattern (Fig. 6). These structures are often more compact than the rings displayed in Fig. 5, with typical radii of 50–100 $r_{P}$. Often, the tracers follow a clean one-armed spiral (lower left panel). Sometimes, the tracers lie in two distinct concentrations, with a weak trailing spiral from each structure (upper middle panel). Other ensembles of tracers show very tightly wound spirals (lower middle panel) or more open systems resembling a spiral galaxy (upper right panel).

As with the ring systems in Fig. 5, systems with spiral structure tend to persist when the mass of the Charon analog is small, $m_{2}\ \lesssim$ 0.8 $m_{C}$. When the final mass of the Charon analog is larger, short-lived spiral patterns are present as the mass of this analog grows with time. Once the mass exceeds $\sim$ 0.9 $m_{C}$, however, the Charon analog tends to eject small leftover massive particles out through the orbital plane. These massive particles take many tracers with them, destroying the spiral pattern.

Despite the persistent ring or spiral structure among the tracers in many calculations, the typical positions of tracers are well outside the orbits of the current circumbinary satellites. In roughly half of the calculations, the semimajor axes for tracers at 150–300 d are $a_{t}\sim$ 125–175 $r_{P}$, much larger than the semimajor axis of the outermost satellite Hydra, $a_{H}\approx$ 55 $r_{P}$. Only $\sim$ 10% of systems retain tracers close to the orbit of Hydra.

Among all of the calculations, the fraction of tracers remaining in the system after 150–300 days, $f_{t}$, correlates with the final mass of the Charon analog (Fig. 7). In the upper part of the plot, there is a clear progression from $f_{t}\approx$ 0.5–0.6 when $m_{2}\ \approx$ 0.4 $m_{C}$ to $f_{t}\approx$ 0.1 when $m_{2}\ \approx$ $m_{C}$. For systems with fewer tracers, $f_{t}\lesssim$ 0.1, there is a secondary trend where $f_{t}$ drops from $\sim$ 0.1 at $m_{2}\ \approx$ 0.4 $m_{C}$ to less than $10^{-3}$ at $m_{2}\approx$ 0.6–0.7 $m_{C}$. From the color-coding of the points, these trends appear to be independent of $e_{0}$, the initial eccentricity of the $n$-bodies in the swarm.

Aside from the strong correlation of $f_{t}$ with $m_{2}$, the number of systems with few if any tracers is another striking feature of Fig. 7. Among the $\sim$ 475 calculations shown, $\sim$ 15% have only 1–14 tracers left after $\sim$ 300 days of evolution. A much larger fraction, $\sim$ 33%, have no tracers. These percentages correlate well with the final mass of the Charon analog. Of the $\sim$ 235 calculations where the final mass of the Charon analog is $m_{2}\ \lesssim$ 0.7 $m_{C}$, $\sim$ 15% have no tracers after 150–300 d of evolution. In contrast, roughly 50% of the $\sim$ 240 calculations where $m_{2}\ \gtrsim$ 0.7 $m_{C}$  have no tracers at the 150–300 d mark (Table 2, first two rows).

Survival of tracers depends less on the final semimajor axis or eccentricity of the Charon analog than on its mass (Table 2, rows 3–6). Although there is a slight tendency for systems with larger $a_{2}$ to retain more tracers, Charon analogs with $a_{2}\lesssim$ 6 $r_{P}$ and $a_{2}\gtrsim$ 6 $r_{P}$ are roughly equally likely to have $f_{t}$ = 0 or $f_{t}\gtrsim$ 0.1. The calculations display a similar trend with $e_{2}$: Charon analogs with $e_{2}\lesssim$ 0.22 are somewhat more likely to retain a larger fraction of tracers than those with $e_{2}\gtrsim$ 0.22. However, the differences do not seem significant.

For any final orbit of the Charon analog, systems are much more likely to have a large fraction of their initial complement of tracers or no tracers at all than to have some intermediate fraction of tracers (Table 2). Roughly equal numbers of systems (33%) either have no tracers or most of their tracers. Another 15% have less than 0.1% of their initial set of tracers; 10% (7%) have 0.1% to 1% (1% to 10%). The relative lack of systems with 0.1% to 10% of their tracers provides a quantitative measure of the ability of dynamical ejections of massive $n$-bodies to remove tracers efficiently on short time scales.

To quantify correlations between $f_{t}$, the initial conditions for the planetesimals, and the final properties of the Charon analog in more detail, we rely on a parametric measure (Pearson’s linear correlation coefficient $r$) and two non-parametric measures (the Spearman rank-order correlation coefficient, $r_{S}$, and Kendall’s $\tau$; Press et al., 1992). Table 3 lists the results. Confirming the visual impression from Fig. 7, all three tests measure a strong correlation between $f_{t}$ and the final mass of the Charon analog $m_{2}$. Probabilities of no correlation range from $p_{r}\approx 4\times 10^{-15}$ for the Pearson linear correlation coefficient to $p_{\tau}\approx 4\times 10^{-30}$ for Kendall’s $\tau$ to zero for the Spearman rank-order test.

Although Kendall’s $\tau$ suggests a 3$\sigma$ correlation between $f_{t}$ and $a_{2}$ (in the sense that Charon analogs with larger $a_{2}$ retain more tracers), the other two tests measure no significant correlation ($p_{r}$ = 0.51 and $p_{S}$ = 0.18; Table 3). The Spearman rank-order test indicates that the correlation between $f_{t}$ and $e_{2}$ barely has a 2$\sigma$ significance; the other tests measure no correlation ($p_{r}$ = 0.60 and $p_{\tau}$ = 0.90). Thus, the fraction of remaining tracers depends on the final mass of the Charon analog, but not its semimajor axis or eccentricity.

The correlation coefficients also generally confirm the visual impression that the number of tracers remaining after 150–300 d is independent of $e_{0}$. The Pearson ($p_{r}$ = 0.91) and Kendall ($p_{\tau}$ = 0.48) tests suggest very high probabilities that $e_{0}$ and $f_{t}$ are uncorrelated. Although the Spearman rank-order test indicates a very low probability that $e_{0}$ and $f_{t}$ are uncorrelated ($p_{S}\approx 10^{-11}$, with $r_{S}$ = $-0.31$), a plot of $f_{t}$ as a function of $e_{0}$ (not shown) looks completely uncorrelated. We suspect that the low probability from the Spearman rank-order test results from the marginally larger than average number of systems with $e_{0}$ = 0.4 that retain fewer than 10 tracers after 150–300 d and the marginally smaller than average number of systems with $f_{t}\approx$ 0.001-0.1 (Table 2).

In Fig. 7, there are several examples where the $n$-bodies have ejected nearly all of the initial complement of tracers ($f_{t}\lesssim 10^{-2}$) despite the production of a very low mass Charon analog ($m_{2}\ \lesssim$ 0.1 $m_{C}$). As the collisional and dynamical evolution proceeds, these systems initially produce a massive Charon analog ($m_{2}\ \gtrsim$ 0.75 $m_{C}$) with a small semimajor axis ($a_{2}\lesssim$ 4 $r_{P}$) and a very eccentric orbit ($e_{2}\gtrsim$ 0.4–0.5). Interactions with several smaller left over $n$-bodies generates an orbit with larger and larger $e_{2}$; eventually, the Charon analog collides and merges with Pluto. Although 1–2 low mass $n$-bodies may survive this evolution on high $e$ orbits at large $a$, most of the tracers are ejected. Given the small numbers, it is impossible to judge whether these systems are associated with any particular $e_{0}$.

Although the ability of a system to retain tracers is independent of $e_{0}$, swarms that start the evolutionary sequence with larger $n$-bodies lose their tracers more rapidly than swarms with initially smaller $n$-bodies (Table 2). For this set of calculations, the fraction of systems with no remaining tracers after 150–300 days is 28% ($r_{0}$ = 145 km), 32% ($r_{0}$ = 185 km), and 41% ($r_{0}$ = 230 km). Coupled with our result that the final mass of a Charon analog correlates with $r_{0}$, systems with a massive Charon analog are less likely to retain their tracers after 150–300 days.

Fig. 8 illustrates the trend of $f_{t}$ as a function of $m_{2}$ and $r_{0}$. The vast majority of the systems with $m_{2}\ \gtrsim$ 0.7 $m_{C}$ and $f_{t}\gtrsim$ 0.1 have $r_{0}$ = 145 km. Only a few have $r_{0}$ = 185 km or 230 km. Among systems with $m_{2}\ \gtrsim$ 0.7 $m_{C}$ and $f_{t}\lesssim$ 0.01, most have $r_{0}$ = 185 km or $r_{0}$ = 230 km (Table 2).

Curiously, the fraction of systems with $f_{t}\approx 10^{-4}-10^{-1}$ is fairly independent of $r_{0}$ (Table 2). Despite the large number of systems with $f_{t}$ = 0 and $f_{t}\gtrsim$ 0.1, only $\sim$ 7% of systems have $f_{t}$ = 0.01–0.1; another 10% have $f_{t}$ = 0.001–0.01 and 15% have $f_{t}=10^{-4}$ to $10^{-3}$. As with tracer retention as a function of the physical properties of the Charon analog, this behavior – where the systems tend to have $f_{t}\gtrsim$ 0.1 or $f_{t}$ = 0 – illustrates the effectiveness of scattered $n$-bodies in removing tracers from the system.

The statistical tests generally confirm these conclusions. All of the tests measure a small probability for the lack of a correlation between $r_{0}$ and the fraction of tracers remaining after 150–300 d of evolution (Table 3). For the Pearson and Kendall tests, the correlation between $r_{0}$ and $f_{t}$ is significant at the 3$\sigma$ or 4$\sigma$ level. The Spearman rank-order test measures a probability of zero that $r_{0}$ and $f_{t}$ are uncorrelated. These results are mainly due to the relation between $m_{2}$ and $r_{0}$, where systems with larger planetesimals generally produce a more massive Charon analog (§3).

Among calculations where the evolution leads to the ejection of most tracers and the merger of the Charon analog with Pluto, there is no obvious preference for the initial radius of planetesimals. This conclusion is based on small-number statistics. Of the 13 calculations with $m_{2}\ \lesssim$ 0.2 $m_{C}$, seven have $f_{t}\lesssim 10^{-3}$ at 150–30 d. Although four of the seven have $r_{0}$ = 145 km, one has $r_{0}$ = 185 km, and two have $r_{0}$ = 230 km.

Tracer retention also depends on the initial mass of the swarm (Fig. 9). Swarms with $M_{0}$ = 0.85 $m_{C}$ cannot produce a Charon analog with the mass of Charon. Despite having a lower mass Charon analog, these systems may still lose nearly all of their tracers after 150–300 d. Although 47% of the evolutionary calculations retain at least 10% of their tracers, 13% lose 99% to 99.9% and another 14% lose 99.9% to 99.99% (Table 2). Roughly 20% lose all of their tracers. The fraction of lost tracers loosely correlates with the final mass of the Charon analog: systems with $m_{2}\ \lesssim$ 0.7 $m_{C}$ are more likely to retain a substantial population of tracers than those with more massive Charon analogs.

When $M_{0}$ = 1.25 $m_{C}$ (Fig. 9, green points), dynamical evolution generates a massive $n$-body with $m_{2}\ \gtrsim$ 0.95 $m_{C}$ ($m_{C}$) roughly 20% (10%) of the time. Compared to calculations where $m_{2}\ \lesssim$ 0.95 $m_{C}$, these systems retain fewer tracers ($f_{t}\lesssim$ 0.1) than those with a smaller massive $n$-body ($f_{t}\lesssim$ 0.8–0.9). Independent of the final mass of the Charon analog, these calculations retain fewer tracers overall: 41% lose at least 99% of their initial complement of tracers and 16% lose all of their tracers (Table 2).

Systems with $M_{0}$ = 1.75 $m_{C}$ have an even harder time retaining tracers for 150–300 d (Fig. 9, blue points). Very few of these systems retain their tracers (Table 2); 66% lose all of them. Unlike calculations with lower masses, these almost always generate a Charon-mass (or larger) satellite. In contrast with the lower $M_{0}$ trials studied here, the lack of tracers is more a function of the final mass of the satellite than the initial mass of the swarm.

The statistical tests indicate that the initial mass of the swarm is more important in setting $f_{t}$ than the initial radius of the planetesimals (Table 3). Probabilities that $M_{0}$ and $f_{t}$ are uncorrelated range from $2\times 10^{-12}$ for the Pearson test to $6\times 10^{-30}$ for Kendall’s $\tau$ to zero for the Spearman rank-order test. These probabilities are comparable to those derived for the correlation between $f_{t}$ and $m_{2}$. This result is not surprising: the initial mass of the swarm is much more important in setting the final mass of the Charon analog than the initial radius of planetesimals.

To explore the long-term behavior of the Charon analog and any surrounding tracers, we followed the evolution of selected systems for 10–100 yr. During this period, the physical properties of the Charon analogs change little. Their masses remain constant. Their orbital semimajor axes and eccentricities settle onto constant values that differ little from those at 1 yr. Thus, the binary planet composed of Pluto and a Charon analog is stable.

As the Charon analog reaches a stable state, the evolution in the tracer population is dramatic (Fig. 10). All systems with a massive Charon analog, $m_{2}\ \gtrsim$ 0.9 $m_{C}$, lose 99.9% of their tracers. When $m_{2}\ \gtrsim$ 0.95 $m_{C}$, only two of 63 calculations retain any tracers, one with four and another with one. Both of these calculations have $M_{0}$ = 1.25 $m_{C}$. None of the calculations with $M_{0}$ = 1.75 $m_{C}$ that produce a Charon analog with $m_{2}\ \gtrsim$ 0.9 $m_{C}$ retain any tracers after 10–100 yr of dynamical evolution.

The ability to retain tracers after 10–100 yr clearly correlates with the mass of the Charon analog. In systems with $m_{2}\ \approx 0.5~{}m_{C}$, the number of tracers remaining after 10–100 yr is similar to the number at 150–300 d. Low mass Charon analogs are not as efficient at ejecting leftover $n$-bodies; most leftovers simply collide and merge with Pluto.

From the formation of a massive Charon analog to 10–100 yr, the timing of tracer ejection is random. In some systems, a few remaining low mass $n$-bodies make close approaches to either Pluto or the Charon analog early in the evolution. These are rapidly ejected and remove any remaining tracers on their way out of the system. In other calculations, it may take many months for Pluto or Charon to eject several left over $n$-bodies. Although all systems with a massive Charon analog have no remaining tracers at 10–100 yr, the epoch when the last tracers are ejected depends more on the chaotic nature of a particular system than the mass of the Charon analog or the time when it reaches its final mass.

Other properties of the system – initial conditions ($e_{0}$, $r_{0}$ and $M_{0}$) and the final orbit of the Charon analog ($a_{2}$ and $e_{2}$) – have little impact on the tracer population at 10–100 yr. Once the mass of the Charon analog is set, the tracers ‘forget’ the initial conditions and respond only to the dynamics created by the Charon analog. When the Charon analog is massive, (i) it creates a tracer population with somewhat larger semimajor axes $a_{t}$ and eccentricities $e_{t}$ than lower mass Charon analogs, (ii) it ejects more left over $n$-bodies through the population of tracers, and (iii) its stronger gravity more strongly destabilizes tracers with large $e_{t}$. As a results, tracers do not survive in a system with a massive Charon.

To understand whether this dynamical evolution of the tracers is characteristic, we performed several calculations with a simple collisional damping model. For each tracer with semimajor axis $a$, the change in $e$ and $\imath$ is

[TABLE]

and

[TABLE]

where $t_{damp}$ is a damping time appropriate for collisions among particles with radii $r$ and total mass $M_{t}$ (Bromley & Kenyon, 2015):

[TABLE]

If this damping can overcome excitation by massive $n$-bodies, tracers might remain on circular orbits far from the central binary.

Experiments with $r\approx$ 1–100 cm and $M_{t}=10^{22}-10^{23}$ g yield similar results. Although some tracers manage to achieve fairly circular orbits at $a\approx$ 50–60 $r_{P}$, the set of lower mass $n$-bodies sweeping through the tracers every orbit tends to scatter a few of them to larger $a$. The longer damping times at larger $a$ preclude efficient damping; tracers at larger $a$ are then ejected. Over the course of 10–300 days, the cumulative impact of the $n$-bodies removes nearly all of the tracers orbiting a central Pluto–Charon binary.

5 DISCUSSION

The calculations described here are the first to explore the growth of a Charon mass satellite following a graze-and-merge collision between a planetary embryo and Pluto. In systems where the debris has an initial mass $M_{0}$ = 0.85–1.75 $m_{C}$ at semimajor axes $a\approx$ 3–11 $r_{P}$, collisional evolution almost always leads to a massive satellite with a significant fraction of Charon’s mass. When $M_{0}$ = 1.25 $m_{C}$ (1.75 $m_{C}$), $\sim$ 10% (55%) of the calculations yield a satellite with $m_{2}\gtrsim$ $m_{C}$. Overall, 12% ($M_{0}$ = 1.25 $m_{C}$) to 20% ($M_{0}$ = 1.75 $m_{C}$) of the calculations produce a satellite with $m_{2}\approx$ 0.95–1.05 $m_{C}$. Given sufficient mass in a circum-Pluto ring of planetesimals at 3–11 $r_{P}$, massive satellite formation is a fairly robust outcome of a graze-and-merge impact (see also Canup & Asphaug, 2001; Leinhardt et al., 2010; Asphaug, 2014).

The approach adopted here neglects the impact of tidal stresses on the shapes and collision outcomes of massive solids inside the Roche limit. When the solids in these calculations venture inside the Roche limit, they probably do not have time to develop the expected equilibrium shapes (e.g., Quillen et al., 2016). Thus, a more accurate treatment of shapes probably has little impact on the growth of Charon mass objects. Although we confirm that no collisions occur inside the Roche limit, an improved treatment would begin with material inside the Roche limit (as suggested by the SPH calculations; Canup, 2005, 2011) and follow the evolution of solids inside and outside the Roche limit. Although it seems unlikely that including material inside the Roche limit would change results significantly, we plan to examine this possibility in a future study.

In these calculations, we have assumed that the mass density of individual solids within the circum-Pluto ring of planetesimals nearly matches the mass density of Charon. This assumption runs counter to the SPH results of Canup (2005, 2011), who predicted a much smaller mass density for solids within the debris from a graze-and-merge collision between two partially differentiated embryos. In this scenario, most of the rocky material within the impactor merges with Pluto. Any circum-Pluto material is then mostly ice. However, Desch (2015) and Desch & Neveu (2017) suggest the colliding embryos probably have undifferentiated crusts. Graze-and-merge impacts then leave behind substantial rocky material in a circum-Pluto disk.

The degree of differentiation of colliding embryos depends on their formation and evolution histories. In a protosolar disk where embryos grow rapidly from pebble accretion (e.g., Chambers, 2016; Johansen & Lambrechts, 2017), large solids may retain significant amounts of short-lived radioactivities, which serve as a major heat source following formation. Analysis of data from the New Horizons mission appears to rule out this possibility (McKinnon et al., 2017). The slower growth of embryos from 1–10 km planetesimals (e.g. Kenyon & Bromley, 2008, 2010, 2012) allows heating only from longer-lived radioactivities (e.g., Malamud & Prialnik, 2015) and may be more consistent with the partially differentiated objects required in hit-and-run or graze-and-merge impacts. Given the uncertainty in formation histories, it is worthwhile to consider in more detail how the outcomes of giant impacts depend on the formation mechanism.

5.1 Growth Time and Long-term Dynamical Evolution of Charon Analogs

Starting with swarms of 145–230 km planetesimals orbiting Pluto at 3–11 $r_{P}$, Charon analogs grow rapidly. Within a day, more than half of the planetesimals have merged into a larger object. After 10 days, only a Charon analog and several low mass objects remain. By 100 days, the Charon analog has usually accreted or ejected all of the leftovers.

Although the growth of Charon analogs seems rapid, the 100 day time scale is in line with theoretical expectations (see also Arakawa et al., 2019, for a discussion of SPH calculations of satellite evolution after a giant impact). At $a\approx$ 50–60 $r_{P}$, it takes $\sim$ 0.1 $(M_{0}/{\rm 3~{}\times~{}10^{20}~{}g)}$ Myr to convert an initial mass $M_{0}$ of small planetesimals orbiting Pluto into several large objects that contain most of the initial mass of the swarm (Kenyon & Bromley, 2014). For the initial masses considered here, $M_{0}\approx 3\times 10^{24}$ g, the growth time is $10^{4}$ times smaller, only 10 yr. In addition to the scaling between growth time and the initial mass of the swarm, ensembles of solids evolve on a time scale proportional to the orbital period (e.g., Lissauer, 1987, 1993; Goldreich et al., 2004; Kenyon & Bromley, 2006). Scaling results at 50–60 $r_{P}$ to 5–6 $r_{P}$ yields an expected growth time of $\sim$ 100 days, which agrees with the typical growth time of the Charon analogs described here.

The final orbits of Charon analogs in these calculations are similar to those of the Charon-mass survivors of hit-and-run collisions (Canup, 2005, 2011). Thus, the subsequent tidal evolution of Pluto and Charon should be fairly independent of the formation mechanism. With $a_{C}\approx$ 5–6 $r_{P}$ and $e_{C}\approx$ 0.1–0.4, it takes $\sim$ 1 Myr for tidal forces between Pluto and Charon to circularize the orbit and to migrate Charon to its current semimajor axis (Farinella et al., 1979; Dobrovolskis et al., 1997; Peale, 1999; Cheng et al., 2014a; Barr & Collins, 2015). These analyses assume that the immediate environment of Pluto–Charon has no other large objects with significant gravity. For the outcomes of graze-and-merge scenarios considered here, the Pluto–Charon binary clears out the inner 25–50 $r_{P}$ of the binary within 1–100 yr. In systems where the leftover solids are much smaller than 150–200 km, mass removal might take somewhat longer, but the clearing time should still be much smaller than 1 Myr (see also Winter et al., 2010; Giuliatti Winter et al., 2013, 2014, 2015; Kenyon & Bromley, 2019a).

5.2 Long-term Thermal Evolution of Charon Analogs

Although examining the long-term internal evolution of Charon requires detailed structural calculations with state-of-the-art equations-of-state (e.g., Robuchon & Nimmo, 2011; Malamud & Prialnik, 2013; Desch, 2015; Malamud & Prialnik, 2015; McKinnon et al., 2016; Hammond et al., 2016; Desch & Neveu, 2017; McKinnon et al., 2017; Bierson et al., 2018), it is worth considering whether the rapid formation time might establish somewhat different initial conditions for this evolution. Once it reaches its final mass, Charon has a gravitational binding energy

[TABLE]

where we assume for simplicity a uniform sphere of constant density. With $m_{C}=1.6\times 10^{24}$ g and $r_{C}$ = 606 km, $E\approx 2\times 10^{33}$ erg. Setting the internal energy $U=NkT=E/2$, where $N$ is the number of molecules in Charon, $k$ is Boltzman’s constant, and $T$ is the internal temperature:

[TABLE]

Making the simple assumption that Charon is 100% water ice, $N=m_{C}/18m_{H}$, where $m_{H}$ is the mass of a hydrogen atom. Then, $T\approx$ 90 K. A Charon composed of 50% water ice and 50% silicon has $T\approx$ 70 K.

It is unlikely that Charon analogs formed during a graze-and-merge impact radiate the gravitational binding energy efficiently. Adopting an effective temperature $T_{e}$ = 90 K (70 K), Charon has a luminosity $L_{C}\approx 2\times 10^{20}$ ${\rm erg~{}s^{-1}}$ ($L_{C}\approx 6\times 10^{19}$ ${\rm erg~{}s^{-1}}$). The time scale required to radiate $\sim$ 50% of the binding energy is $\sim$ 0.2 Myr (0.5 Myr). Thus, Charon retains the heat of formation over a time comparable to the tidal expansion time scale.

Despite the relatively slow cooling compared to planetary embryos grown from km-sized or larger planetesimals, the $T$ = 70–90 K derived above is similar to the initial temperature adopted for several detailed evolutionary calculations, $T\approx$ 100 K (Malamud & Prialnik, 2015) and $T$ = 60–100 K (Desch & Neveu, 2017). In these studies, radioactive decay of 235U and other long-lived unstable elements powers a gradual rise in the core temperature $T_{c}$, which reaches $T_{c}\approx$ 500–1000 K $\sim$ 1 Gyr after the impact. This core temperature is sufficient to melt water ice and produce a differentiated planet (see also McKinnon et al., 2017). Although the core temperature is sensitive to composition and the details of the equation-of-state, it is rather insensitive to $E$.

Aside from radioactive decay, other sources of heating are also negligible. The energy in the orbit, $E_{o}\approx Gm_{P}m_{C}/a\sim 2\times 10^{33}$ erg, is comparable to the gravitational binding energy of Charon. Unless the rotational periods of Pluto and Charon are much less than 1 day, the available rotational energy is roughly equal to the orbital energy. Thus, tidal expansion has modest impact on the thermal energy of either Pluto or Charon (Dobrovolskis et al., 1997; Barr & Collins, 2015). Collisions with 100–200 km Kuiper belt objects only add significant energy if the impacts occur at velocities much larger than the orbital velocity of Pluto–Charon around the Sun, which is very unlikely. Overall, the rapid formation of Charon in the graze-and-merge scenario does not change expectations for its thermal evolution.

5.3 Future Prospects: Survival of the Circumbinary Satellite System

However Charon forms, the dynamical evolution of massless tracers described in §4 places strong constraints on the aftermath of hit-and-run and graze-and-merge impacts. In either scenario, some debris orbits the system barycenter well outside Charon’s orbit. For typical surface densities within the debris at 30–50 $r_{P}$, collisional damping rapidly circularizes the orbits of small particles with radii $r\lesssim$ 0.1 km (Bromley & Kenyon, 2015). On time scales of $10^{4}-10^{5}$ yr, collisional evolution can then transform swarms of small particles into an ensemble of 5–20 km satellites (Kenyon & Bromley, 2014). However, if the debris contains several massive particles with radii $r\gtrsim$ 150–200 km, dynamical ejection of only one or two of these objects by the Pluto–Charon binary is sufficient to disrupt the entire swarm of circumbinary particles.

In the calculations described here, the massless tracers serve as proxies for small particles produced in the giant impact or from the collisions of massive $n$-bodies as they merge to produce a Charon analog orbiting Pluto. If the small particles constitute a fraction $f_{s}$ of the initial mass in solids, the total mass in small particles at the end of a calculation is $M_{s}\approx f_{t}f_{s}M_{0}$, where $f_{t}$ is the fraction of tracers remaining after 10–100 yr of evolution. Recalling $M_{SNKH}$ as the total mass of the known circumbinary satellites, setting $M_{s}\gtrsim M_{SNKH}$ allows the surviving tracers to have sufficient mass to form Styx, Nix, Kerberos, and Hydra on a reasonable time scale. With $M_{0}\approx$ 1–2 m_{C}\ $\lesssim 3\times 10^{24}$ g and $M_{SNKH}\lesssim 10^{20}$ g, we require $f_{t}f_{s}\gtrsim 10^{-4}$.

Our calculations suggest this constraint on $f_{t}f_{s}$ is not achievable with the initial conditions considered here. From SPH simulations of giant impacts, having more than 1% of the debris in the form of small particles seems unlikely (e.g., Canup, 2005, 2011; Arakawa et al., 2019). Debris from the collisions of massive icy $n$-bodies is also unlikely to exceed 1% (§2). For swarms of massive $n$-bodies capable of producing a massive Charon analog, $f_{t}\lesssim 10^{-4}$. Thus, the mass remaining in small particles is not large enough to produce the current circumbinary satellites.

Aside from the limited available mass, any surviving tracers are well beyond the orbits of the circumbinary satellites. Scattering by massive $n$-bodies places tracers on eccentric orbits with typical semimajor axes $a_{t}\approx$ 125–175 $r_{P}$, outside the current orbit of Hydra, $a_{H}\approx$ 55 $r_{P}$. On time scales of $10^{1}-10^{6}$ yr, massive collections of small particles with radii $r\approx$ 1 m to 1 km can spread inward from $\sim$ 150 $r_{P}$ to the orbit of Hydra (Bromley & Kenyon, 2015). If sufficient numbers of tracers survive the formation of a Charon analog, the current ensemble of circumbinary satellites might form in a spreading disk of debris.

Growing Charon within a swarm of lower mass planetesimals might allow enough massless tracers to survive the ejections of any leftovers. For material orbiting Pluto at 50-200 $r_{P}$, the orbital velocity is $\sim$ 0.07–0.14 $\rm km~{}s^{-1}$. With escape velocities $v_{e}\approx$ 0.14–0.23 $\rm km~{}s^{-1}$, massive planetesimals with $r$ = 145–230 km stir up nearby tracers to velocities larger than the local escape velocity. Thus, tracers are ejected. Lower mass planetesimals with $r$ = 10–20 km have a factor of ten smaller escape velocities and therefore stir tracers to velocities $\sim$ 10% of the local escape velocity. Small particles with these velocities might survive collisional evolution within circumbinary rings (Kenyon & Bromley, 2014; Walsh & Levison, 2015; Bromley & Kenyon, 2015).

Following the collisional evolution of swarms of particles with $r\lesssim$ 10–20 km at semimajor axes $a\approx$ 5–20 $r_{P}$ requires an accurate treatment of fragmentation (e.g., Kenyon & Bromley, 2014). For typical velocities of 0.05–0.10 $\rm km~{}s^{-1}$ at $a\approx$ 5–10 $r_{P}$, collisions of equal mass 10–20 km objects eject 5% to 10% of the combined mass in small particles. Unlike the calculations of 145–230 km objects here, the large amount of debris generated in collisions of much smaller objects may change the growth of Charon analogs and the ejection of small particles by surviving $n$-bodies. In calculations of icy objects with fragmentation, leftover planetesimals are often small (e.g., Kenyon & Bromley, 2008, 2010, 2014). Dynamical friction between large objects and collisional debris maintains large objects on fairly circular orbits. Thus, there is some chance that starting with a swarm of much lower mass planetesimals might allow the survival of a circumbinary swarm of debris at 30–60 $r_{P}$. We plan to consider these kinds of calculations in a future study.

If small particles at 25–100 $r_{P}$ can survive the immediate aftermath of a graze-and-merge collision, this material must somehow endure the tidal evolution of the binary (e.g., Ward & Canup, 2006; Lithwick & Wu, 2008; Cheng et al., 2014b; Bromley & Kenyon, 2015; Smullen & Kratter, 2017; Woo & Lee, 2018). As the binary circularizes and expand, orbits with small integer ratios of the binary orbital period become unstable. As an example, the 5:1 resonance passes through material orbiting at 30 $r_{P}$ (just inside the current orbit of Styx) when the central binary has a period of 4.5 days. When the expansion is rapid, some material survives on high eccentricity orbits. However, current models suggest circularization and expansion are rather slow. Few solids survive this expansion (see also Walsh & Levison, 2015).

Collisional damping allows solids orbiting at 30–60 $r_{P}$ to survive tidal evolution of the central binary (Bromley & Kenyon, 2015). When a resonance encounters a ring of m-sized or somewhat larger particles, the eccentricities of the particles begin to grow. However, collisions among the swarm damp the eccentricity as fast as the resonance excites it. Damped particles are then transported out with the resonance. Although some particles are lost, most survive the expansion and lie within resonances as the expansion ends.

This mechanism offers a way to produce circumbinary satellites after tidal expansion following either a graze-and-merge or a hit-and-run impact. If the debris at 30–60 $r_{P}$ is composed of km-sized or smaller particles, collisions will circularize the debris on time scales short compared to the tidal expansion time scale (Kenyon & Bromley, 2014; Bromley & Kenyon, 2015). As the binary expands, these particles become trapped in resonance. Once tidal expansion ceases, collisional processes convert the survivors into a few small satellites near resonance. The main challenge is for the solids to avoid growing into satellites before tidal expansion ends. We plan to describe numerical simulations of this process in a future paper.

6 SUMMARY

We consider the evolution of massive disks of planetesimals orbiting a planet with mass similar to Pluto. Within a few weeks, 145–230 km survivors of a graze-and-merge impact grow into a Charon analog with an orbit – $a_{C}\approx$ 5–6 $r_{P}$ and $e_{C}\approx$ 0.1–0.3 – similar to the orbits of Charons that survive a hit-and-run collision and remain bound to Pluto. In a typical calculation, the Charon analog contains $\sim$ 60% of the initial mass in large objects. Pluto typically accretes roughly 25% of the initial mass and ejects the rest.

When a swarm of planetesimals produces a Charon analog with $m_{2}\gtrsim$ 0.9 $m_{C}$, nearly all massless tracers are ejected on time scales of months to decades. Initially, the growth of the Charon analog places tracers on eccentric orbits with $a\approx$ 50–250 $r_{P}$. Although collisional damping would probably circularize the tracers on long time scales (e.g., Kenyon & Bromley, 2014; Bromley & Kenyon, 2015), the central binary ejects several leftover massive planetesimals through the tracers. These ejections disrupt tracer orbits; eventually all of the tracers are also ejected. Thus, there is no circumbinary disk of solids in which to grow satellites with properties similar to the known small satellites.

A simple dynamical analysis suggests that solids at $a\gtrsim$ 30 $r_{P}$ would survive the passage and subsequent ejection of less massive objects with $r\lesssim$ 10–20 km. If Charon grows efficiently within a circum-Pluto debris disk and ejects leftovers with radii no larger than 10–20 km, it might then be possible to retain a circumbinary disk of solids and form small satellites at 30–60 $r_{P}$. We plan to test this scenario in a future study.

We acknowledge generous allotments of computer time on the NASA ‘discover’ cluster. Advice and comments from M. Geller greatly improved our presentation. We thank the referee for a timely and useful report. Portions of this project were supported by the *NASA * Outer Planets and Emerging Worlds programs through grants NNX11AM37G and NNX17AE24G.

Binary output files from the simulations and C programs capable of reading the binary files are available at a publicly accessible repository (https://hive.utah.edu/) with digital object identifier https://doi.org/10.7278/S50D-EFCY-ZC00.