A Pluto--Charon Concerto II. Formation of a Circumbinary Disk of Debris After the Giant Impact
DD RAFT VERSION F EBRUARY
24, 2021Typeset using L A TEX preprint style in AASTeX63
A Pluto–Charon Concerto II. Formation of a Circumbinary Disk of Debris After the Giant Impact S COTT
J. K
ENYON
Smithsonian Astrophysical Observatory, 60 Garden Street, Cambridge, MA 02138 B ENJAMIN
C. B
ROMLEY
Department of Physics & Astronomy, University of Utah, 201 JFB, Salt Lake City, UT 84112
ABSTRACTUsing a suite of numerical calculations, we consider the long-term evolution of circumbi-nary debris from the Pluto–Charon giant impact. Initially, these solids have large eccentricityand pericenters near Charon’s orbit. On time scales of 100–1000 yr, dynamical interactionswith Pluto and Charon lead to the ejection of most solids from the system. As the dynamicsmoves particles away from the barycenter, collisional damping reduces the orbital eccen-tricity of many particles. These solids populate a circumbinary disk in the Pluto–Charonorbital plane; a large fraction of this material lies within a ‘satellite zone’ that encompassesthe orbits of Styx, Nix, Kerberos, and Hydra. Compared to the narrow rings generated fromthe debris of a collision between a trans-Neptunian object (TNO) and Charon (Bromley &Kenyon 2020), disks produced after the giant impact are much more extended and may be aless promising option for producing small circumbinary satellites.
Keywords: planets and satellites: dynamical evolution — planets and satellites: formation— dwarf planets: Pluto INTRODUCTIONFor nearly a century, the Pluto–Charon system has provided a new window into the physical processesthat shape planetary systems. As the harbinger of the Kuiper belt, the discovery of Pluto contributed a firstglimpse of the architecture of the trans-Neptunian region. Combined with light curve data (Andersson &Fix 1973), the identification of Charon on high quality images revealed a tidally-locked binary planet witha short orbital period (6.4 d, Christy & Harrington 1978). Mutual eclipses later enabled novel extractiontechniques for the characterization of surface properties and their evolution with time (e.g., Buie et al.2010a,b). The discovery of four circumbinary satellites – Styx, Nix, Kerberos, and Hydra – yielded newinsights into circumbinary orbital dynamics and tidal evolution (Showalter et al. 2011, 2012; Cheng et al.2014b; Brozovi´c et al. 2015; Showalter & Hamilton 2015). Finally, the spectacular results of the
NewHorizons mission illustrated complex geological processes on frozen worlds and made new connectionsbetween surface phenomena on Pluto–Charon and the physical properties of objects in the Kuiper belt (e.g., e-mail: [email protected]: [email protected] a r X i v : . [ a s t r o - ph . E P ] F e b Stern et al. 2015; Grundy et al. 2016; McKinnon et al. 2016; Weaver et al. 2016; Stern et al. 2018; Singeret al. 2019; McKinnon et al. 2020; Singer et al. 2021).Together with the Earth–Moon system, Pluto–Charon establishes interesting constraints on the frequencyof ‘giant impacts’, collisions between massive protoplanets often considered as a final stage of planet for-mation (e.g., Agnor et al. 1999; Canup & Asphaug 2001; Canup 2004; Agnor & Asphaug 2004; Kenyon& Bromley 2006; Asphaug 2014). In the current picture, mutual gravitational perturbations place grow-ing protoplanets on crossing orbits (see also Wetherill 1980; Lissauer 1987; Weidenschilling et al. 1997;Goldreich et al. 2004). As protoplanet orbits become chaotic, collision velocities become large enough toenable catastrophic collisions, which may lead to complete mergers or binary planets, with additional ejec-tion of debris (e.g., Canup 2005; Asphaug et al. 2006; Canup 2011; Genda et al. 2012; Stewart & Leinhardt2012; Genda et al. 2015b,a; Kenyon & Bromley 2016; Quintana et al. 2016).Although many studies elucidate the physical state of material following the Moon-forming impact (e.g.,Nakajima & Stevenson 2014; Pahlevan et al. 2016; Lock et al. 2018; Nakajima & Stevenson 2018; Tang& Young 2020), there are few investigations of the evolution of circumbinary debris following the Pluto–Charon collision. Several analyses explain how small satellites on circumbinary orbits become unstable asthe central binary expands from a period of 1–2 d to the current 6.4 d (e.g., Ward & Canup 2006; Lithwick& Wu 2008; Youdin et al. 2012; Cheng et al. 2014b; Bromley & Kenyon 2015b; Woo & Lee 2018). Kenyon& Bromley (2014b) and Walsh & Levison (2015) consider how collisional damping converts eccentriccircumbinary orbits into the more circular orbits required for the growth of solids into 5–20 km satellites.Bromley & Kenyon (2015b) demonstrate how a swarm of cm-sized pebbles could stabilize satellite orbitsthroughout the tidal evolution of the central binary.Here, we examine the evolution of solids orbiting a highly eccentric Pluto–Charon binary. Within theCanup (2005, 2011) smooth particle hydrodynamics (SPH) calculations, there is a broad range of plausibleoutcomes for the semimajor axis, a ≈ R P (where R P is the radius of Pluto) and eccentricity, e ≈ ∼ g up to ∼ × g. The more likely configurations have a mass in debris that is ∼ a = 10 R P , e = 0.6–0.8) retain a largerfraction of the impact debris than more compact binaries ( a = 5 R P , e = 0.2–0.4). Circumbinary systems ofsolids with less mass ( g) preserve more of their initial mass than more massive swarms with ∼ g.All of the 70 calculations maintain enough mass to form several 5–20 km satellites.After establishing the numerical methods for our calculations in §2, we describe several simulations indetail and summarize global results in §3. We conclude with a discussion (§4) and a brief summary (§5). NUMERICAL METHODS2.1.
Basic Parameters
To set the stage for the calculations, we adopt measured properties of the Pluto–Charon system (Stern et al.2015; Nimmo et al. 2017; Stern et al. 2018). For an adopted gravitational constant G = 6 . × − ,Pluto has mass m P = 1 . × g, radius R P = 1188.3 km, and mean density ρ P = 1.854 g cm − . Charonhas mass m C = 1 . × g, radius R C = 606 km, and mean density ρ C = 1.702 g cm − . Combinedwith analyses of HST observations (Brozovi´c et al. 2015; Showalter & Hamilton 2015), detailed n -bodycalculations (Kenyon & Bromley 2019b) provide limits on the mass of Nix ( m N (cid:46) . × g) and Hydra( m H (cid:46) . × g; see also Youdin, Kratter & Kenyon 2012). Although constraints on the masses of Styx( m S ) and Kerberos ( m K ) are more elusive, plausible estimates are m S ≈ × g and m K ≈ g.Thus, the total mass in satellites is m sat ≈ g. These satellites orbit within a ‘satellite zone,’ a ≈ R P , defined as a region that encompasses the orbits of the four satellites (e.g., Kenyon & Bromley2019a; Bromley & Kenyon 2020).The current Pluto–Charon system has a semimajor axis, a ≈ R P , and an orbital period P = 6.387 d(see Stern et al. 2018, and references therein). The orbit is nearly circular (Buie et al. 2012); satellite orbitsare also nearly circular (e.g., Buie et al. 2013; Brozovi´c et al. 2015; Showalter & Hamilton 2015). For thenumerical calculations considered here, we adopt a much more compact binary, a ≈ R P , with a rangeof eccentricities, e ≈ – yr (e.g., Farinella et al. 1979; Dobrovolskis et al. 1997; Peale1999; Cheng et al. 2014a; Correia 2020). This time scale is much longer than the 1000 yr evolution time ofa typical numerical calculation for the evolution of circumbinary debris. Thus, we ignore tidal evolution.2.2. Numerical Code
To calculate the evolution of circumbinary solids in the Pluto–Charon system, we run numerical simula-tions with
Orchestra , an ensemble of computer codes designed to track the accretion, fragmentation, andorbital evolution of solid particles ranging in size from a few microns to thousands of km (Kenyon 2002;Bromley & Kenyon 2006; Kenyon & Bromley 2006, 2008; Bromley & Kenyon 2011a, 2013; Kenyon &Bromley 2016; Kenyon et al. 2016; Bromley & Kenyon 2020). In this application, we employ several
Orchestra tools: a multiannulus coagulation code to derive the dynamical evolution of small solids and agravitational n -body code to enable the solids to respond to the gravitational potential of the central binary.A set of massless tracers enables communication between the coagulation and n -body codes.Using the Orchestra multiannulus coagulation routine, we establish a radial grid, centered on the Pluto–Charon barycenter, of 28 concentric annuli distributed in equal intervals of a / between a in (as listed inTable 1) and a max = 160 R P . Each annulus has 80 mass bins with minimum radius r min = 0 . µ m andmaximum radius r max = 1 m. In a region between a in and a out , we seed the grid with solids having totalmass M , material density 1.5 g cm − , surface density Σ ∝ a − , and a size distribution n ( r ) ∝ r − . . Withthese initial conditions, most of the mass lies in the largest objects; the innermost annulus has somewhatmore mass than the outermost annulus. Solids have initial pericenter distance q and inclination ı as listedin Table 1.Within the coagulation code, solids evolve due to collisional damping from inelastic collisions and elas-tic, gravitational interactions. For inelastic and elastic collisions, we follow the statistical, Fokker-Planckapproaches of Ohtsuki (1992) and Ohtsuki et al. (2002), which treat pairwise interactions (e.g., dynamical Table 1.
Starting Conditions for Coagulation Calculations a Model a ( R P ) e a s ( R P ) a in ( R P ) a out ( R P ) q ı a The columns list the orbital semimajor axis, a in Pluto radii, and eccentricity, e ,of the Pluto–Charon binary; the semimajor axis, a s in Pluto radii, of the inner-most stable circular orbit; the initial inner and outer radii, a in and a out , and thepericenter q of circumbinary solids (all in units of R P ); and the initial inclination ı of circumbinary solids.friction and viscous stirring) between all objects. We assume the mass distribution is fixed in time. Atthe start of each calculation, collisions between m-sized (cm-sized) particles are destructive (accretive).Within a day or two, however, collisional damping reduces relative velocities by factors of 2–3. Althoughcollisions between m-sized objects remain destructive, collisional growth among smaller objects tends toreplenish these objects about as fast as they are destroyed. Several tests allowing collisions to grow andfragment particles suggests this process modifies the overall size distribution little over the course of a typ-ical simulation. In these tests, we used the full calculational approach for two sets of calculations, where(i) collision outcomes are derived and debris is distributed among mass bins and (ii) collision outcomesare ignored (maintaining a constant size distribution) but collision rates are used to calculate collisionaldamping as in set (i). Aside from mass redistribution, the calculations are otherwise identical, using theelastic formalism to treat gravitational interactions that do not result in a physical collision and the inelasticformalism to treat collisional damping. To save resources for the more cpu-intensive n -body calculations,we thus ignore collisional redistribution of mass throughout each calculation.To allow the e and ı of mass bins to react to the gravity of Pluto and Charon, we assign massless tracerparticles to each mass bin. Among all mass bins, we assign 448,000 tracers. Within an annulus, each massbin is allocated ∼
300 tracers ( ∼ a for each tracer within annulus i is a random number between a i − . δa i and a i + 0 . δa i , where a i is the semimajor axis of the annulus in the coagulation grid and δa i is the radial width of the annulus. With this selection, tracers have semimajor axes that span the full rangebetween the inner edge and the outer edge of the coagulation grid. Orbital e i, and ı i, for each tracer followthe e and ı for the assigned mass bin. Tracers with a = a i have pericenters q i = q = a i (1 − e i ) . Other tracersin annulus i have pericenters q i that cluster around q : q − . δ a i (1 − e i ) ≤ q i ≤ q + 0 . δ a i (1 − e i ) .Because δa i is finite, the tracers have a modest spread in pericenters about the adopted q . Tracers initiallyhave random orbital phases; some tracers are initially near orbital pericenter, while others are near apocenter.Although the midplane of the disk lies in the Pluto–Charon orbital plane, some tracers initially lie withinthe orbital plane while others begin their evolution out of the orbital plane. Once assigned to a mass bin,tracers may move to another annulus in response to gravitational interactions with Pluto and Charon and the ( de/dt, dı/dt ) derived in the coagulation code. However, tracers do not move among mass bins.Within the coagulation code, each mass bin k in each annulus i contains a number of particles N ik and atotal mass M ik . The average mass of a particle is then m ik = M ik /N ik . To allow these particles to respondto the motions of the tracers, we assign a number and total mass of particles to each tracer. With N ti,k tracers per mass bin, each tracer is responsible for N (cid:48) i,k = N i,k /N ti,k coagulation particles. To avoid a tracer‘carrying’ fractional coagulation particles, the algorithm assigns integer numbers of coagulation particles toeach tracer. The mass carried by each tracer follows: M (cid:48) ik = m ik N (cid:48) ik . Summed over all tracers, the masscarried by tracers is the total mass in the coagulation grid. The number and mass of coagulation particlesassigned to that tracer effectively follow that tracer as it orbits the system barycenter.Within the n -body code, tracers are treated as massless particles. As the evolution proceeds, these tracersrespond to the gravitational potential of the central binary. In doing so, they generally move to largerdistances from the barycenter. While some are ejected, others shift from an annulus i into another annulus j . This shift generates a change in the mass of bins in annulus i and a corresponding (opposite) change inthe mass of bins in annulus j . At the end of a time step, the combined set of shifts results in new valuesfor the number and total mass of particles in each mass bin. With new numbers of tracers in each mass bin, N (cid:48) i,k and M (cid:48) i,k also change. Thus, the number of particles and the total mass in those particles (but not theaverage mass of a particle) carried by each tracer varies throughout a calculation.To choose appropriate starting conditions, we examine outcomes of the Canup (2011) SPH calculations.Among this suite of simulations, ‘successes’ that produce an intact Charon and surrounding debris yieldbinaries with a ≈ R P and e ≈ a have larger e . To span this rangeefficiently, we consider binaries with ( a, e ) = (5 R P , 0.2), (5 R P , 0.4), (10 R P , 0.6), and (10 R P , 0.8).Although these choices do not include the lowest e and highest e systems, they are sufficient to learn howthe evolution of the debris depends on the properties of the underlying binary.Starting points for the debris surrounding Pluto–Charon follow Canup (2011) and models for the giantimpact that produced the Moon (e.g., Ida et al. 1997). From the Canup (2011) calculations, the mass of thedebris orbiting Pluto–Charon spans the range from ∼ g to ∼ × g. Considering the initial orbitalstate of the debris, we expect considerable losses from dynamical ejections. Thus, the long-term evolutionof a debris field with a mass of g is unlikely to lead to a stable state with a mass comparable to thecurrent mass of the satellites, ∼ g. We therefore consider masses of − g.To select initial values of ( a, e ) for debris particles, we match the final maximum radius of the equivalentcircular orbit derived by Canup (2011), a eq,max (cid:46) R P ; a eq is related to the angular momentum oforbiting particles (e.g., Canup 2004, 2011; Nakajima & Stevenson 2014; Bromley & Kenyon 2020): a eq = a (1 − e ) . (1)Because most of the debris lies near Charon (see Fig. 2 of Canup 2011), we consider calculations withpericenters q (i) near the apocenter of Charon’s orbit around the system barycenter or (ii) within Charon’sorbit. Assigning a range of semimajor axes, a d ≈ R P , for the adopted q yields equivalent circularorbital radii, a eq ≈ Σ( a ) . Independent of the properties of the central binary, we adopt Σ ∝ a − . Although Canup(2011) does not quote Σ( a ) for the debris, Ida et al. (1997) consider similar surface density distributionsin n-body calculations of the formation of the Moon from an extended disk surrounding the Earth. SPHcalculations for the giant impact do not have the resolution to derive a size distribution of small particles.Consistent with the size distributions derived for other impact studies (e.g., Greenberg et al. 1978; Durda1996; Durda & Dermott 1997; Leinhardt et al. 2000; Leinhardt & Stewart 2009), we adopt a size distribu-tion n ∝ r − . for the debris particles. For simplicity, we assume that all particles within annulus i have thesame initial eccentricity e i and inclination ı i . Thus, particles with an initial semimajor axis at the inner edgeof the coagulation grid have a smaller e i than those with a larger initial semimajor axis.Compared to the central binary, orbiting solids with the adopted initial surface density profile have rel-atively little total angular momentum. Following Canup (2005), the Pluto–Charon system has an angularmomentum L P C ≈ q Ω( m P + m C ) a / (1 + q ) , where q is the mass ratio, Ω is the angular velocity, and a is the semimajor axis (see also Canup 2011). For the adopted masses, L P C ≈ . × g cm s − ( × g cm s − ) for a = 5 R P (10 R P ). The initial configurations of circumbinary solids we considerhave L s ≈ − × g cm s − ( M = 10 g) to L s ≈ − × g cm s − ( M = 10 g). For themost massive systems, the angular momentum in orbiting solids is only ∼
5% to 10% of L P C .Table 1 summarizes the starting points for each calculation. For each combination of a and e of thecentral binary, we choose two values for the initial inclination, ı = 0.025 and 0.25. The first choice allowsfor maximum collisional damping throughout the evolution; the second choice provides limits on the abilityof collisions to evolve a thick disk into the thin disk required for the growth of satellites (Kenyon & Bromley2014b; Walsh & Levison 2015).Among calculations with e = 0.4, we consider two choices for a in : (i) a in = 8 R P places the solids close tothe binary; (ii) a in = 18 R P places solids farther out in the system. Both of these choices match the a eq,max limits derived in Canup (2011). Aside from testing whether the evolution of solids depends on their initialangular momentum, these calculations allow us to connect results for compact Pluto–Charon binaries with a = 5 R P to those for wide binaries with a = 10 R P . In the calculations described below, the mass distributionat 1000 yr depends on the ability of collisional damping to prevent ejections by the central binary. Whenmaterial is close to the binary, collisions and damping are more numerous, but solids spend more time nearthe binary. More distant solids have lower particle densities and less damping, but these solids spend lesstime near the binary. Following the evolution for two sets of starting conditions provides insight into therelative importance of collisional damping and gravitational perturbations from the binary.To give the initial conditions of the models additional context, Table 1 lists a s the semimajor axis of theinnermost stable circular orbit for each model. For the adopted masses of Pluto and Charon, the numericalresults of (Holman & Wiegert 1999) yield a s a = 2 .
175 + 4 . e − . e . (2)In the compact binaries we consider, a s ≈ R P (18 R P ) for e = 0.2 (0.4). Most of the initial mass thenlies outside the innermost stable orbit. The wider binaries have a s ≈ R P (45 R P ) for e = 0.6 ( e = 0.8).At the start of each calculation, most of the mass is on unstable orbits with a < a s .Our solutions to the evolution equations conserve mass and energy to machine accuracy. Timesteps forthe coagulation calculations range from a few seconds at t = 0 to sec near t = 1000 yr. Within the n -bodycode, the adaptive integrator divides coagulation steps into smaller steps as needed to maintain the integrityof the calculation. With 112 cpus, typical calculations require 30000–50000 timesteps and 2–3 cpu sec perstep. Nearly all of this time is spent evolving the dynamics of tracer particles. Over the course of a 1000 yrrun, calculations conserve mass and energy to better than one part in .2.3. Operations
All calculations proceed as follows. For a time step ∆ t , the coagulation code derives the changes in the e and ı of each mass bin from collisions and orbital interactions with all other mass bins. Collision rates arederived from the particle-in-a-box algorithm (e.g., Kenyon & Bromley 2002, 2004, 2008). For each massbin k in annulus i , the collision rate with solids in mass bin l in annulus j depends on the number density ofparticles ( n ik and n jl ), the collision cross-section, the relative velocity, and the overlap of orbits (see sec. 2of Kenyon & Bromley 2008). When i = j , the overlap is 1; otherwise, the overlap is approximately theratio of the volume of annulus j that lies within the volume of annulus i .At the end of the coagulation step, each tracer is assigned a target eccentricity e t and inclination ı t based onits current e i and ı i and the changes in e and ı for its mass bin from the coagulation calculation. These targetsresult in time derivatives for e and ı during a timestep ∆ t , de/dt = ( e t − e i ) / ∆ t and dı/dt = ( ı t − ı i ) / ∆ t .Before the n -body step, tracers are also assigned new numbers and masses of particles based on the numberand mass within each mass bin. Although the n -body code does not use this information, each tracer carriesa changing mass of solids based on the evolution of solids in the coagulation grid.Within the n -body step, we solve the equations of motion for a set massless tracers orbiting the Pluto–Charon binary, using an adaptive algorithm with sixth-order time steps, based on Richardson extrapolation(Bromley & Kenyon 2006, 2011a). The code follows objects in the center-of-mass frame, calculating forcesby direct summation. Aside from the extensive tests reported in Bromley & Kenyon (2006, 2011a), thecode has been used for investigations of migration (Bromley & Kenyon 2011b), Saturn’s rings (Bromley& Kenyon 2013), and the formation and stability of the Pluto–Charon system system (Kenyon & Bromley2014a; Bromley & Kenyon 2015a; Kenyon & Bromley 2019a,b,c; Bromley & Kenyon 2020).Throughout the n -body calculations, tracers evolve with their derived de/dt and dı/dt and respond to thegravitational potential of Pluto and Charon. Orbits evolve with ∼
200 steps per binary orbit. The tracers’eccentricity e and inclination ı are shifted incrementally at rates set by the coagulation code. These changesto e and ı are implemented by small adjustments to the direction of travel and (if necessary) incrementalshifts in position, without affecting the other osculating orbital elements.At the end of the n -body step, tracers have new positions, velocities, and orbital elements a , e , and ı .Tracers with new a that places them outside their ’old’ annulus are placed in new annuli; the number andmass of particles assigned to that tracer move out of the old mass bin into a new mass bin within the newannulus. Some tracers collide with Pluto or Charon; others are ejected beyond the outer limits of the n -bodycalculation space, a (cid:38) R P . After these tracers are deactivated for the remainder of the calculation, thecoagulation particles associated with these tracers are removed from the coagulation grid. When an activetracer has a semimajor axis outside of the coagulation grid, its mass is removed from the old annulus andis not placed in a new annulus. If that tracer returns to the grid before a collision with Pluto–Charon orejection, the mass that it carries also returns to the grid.Although assigning tracers to annuli based on their current position ( x, y, z ) seems reasonable, placementbased on a is more in the spirit of the coagulation code. The collision and Fokker-Planck algorithms deriverates based on particle volumes, V = 4 πa ∆ aH , where ∆ a ≈ δa + ea , δa is the physical width of theannulus, and H is the vertical scale height which depends on the orbital inclination ı (Kenyon & Bromley2008). For circumbinary solids, the physical extent of particle orbits is much larger than the physical widthof each annulus. Thus, placing tracers in annuli according to a recovers the correct volume for calculationsof collision and stirring rates.New orbital elements for tracers inform the e and ı of mass bins in the coagulation code. Within eachmass bin, the median e and ı and their inter-quartile ranges for tracers assigned to that mass bin establishthe new e and ı for the mass bin. The inter-quartile ranges allow for quantitative monitoring of the accuracyof the median in measuring the typical e and ı for a set of tracers. Typically, the inter-quartile ranges aresmall.This sequence provides a closed-loop algorithm that allows the mass bins and the tracers to respond tothe gravity of Pluto–Charon and the orbiting solids. The coagulation particles tell the tracers how to reactto solid material orbiting Pluto–Charon. In turn, the tracers tell the solid material how to react to Pluto–Charon. As long as time steps are not too long, the lag between the coagulation and n -body steps does notintroduce significant offsets in the evolution of the mass bins and the tracer particles.To maximize the accuracy of this approach, the codes have a set of control parameters. During the co-agulation step, changes in the ( e, i ) of particles in each mass bin never exceed 1%. When changes exceedthis limit, the time step is reduced and repeated. During the n -body step, the adaptive integrator adjusts thenumber of adaptive steps per Pluto–Charon orbit to track collisions between tracers and Pluto and Charon.Several examples in Bromley & Kenyon (2006) illustrate the accuracy of the n -body code. After completionof the n -body step, algorithms within the coagulation code verify that changes in the median ( e, i ) amongtracers within a mass bin are still less than 1% and that the changes in the inter-quartile ranges among thetracers is less than 2%. Larger changes trigger a reduction in the time step. As an independent check onthis process, we compared results with the maximum level of changes set at 0.5% and could not measure asignificant difference between outcomes over 1000 yr of evolution. RESULTSFor each of the starting conditions in Table 1, we performed simulations with five initial masses, M =10 g, × g, g, × g, and g. All calculations follow the same evolutionary sequence.During their first pass around the barycenter, tracers get a push from the central binary. Concurrently,collisions damp the orbital e and ı of every tracer. For the rest of a calculation, the behavior of tracersdepends on the relative importance of damping and gravitational accelerations that try to remove solidsfrom the region where circular orbits are unstable. For some tracers, damping circularizes their orbits atlarge a and reduces their inclinations before the binary ejects them from the system. Among other tracers,collisions fail to circularize orbits. Although tracers with large initial ı damp into the orbital midplane ofthe central binary, many are ejected through the midplane beyond the system’s Hill sphere.After 1000 yr of evolution, most surviving tracers lie within a disk surrounding the binary. These tracershave e (cid:46) − . The disk has a sharp inner edge, where solids follow orbits just outside the unstable region, a (cid:38) a s . The peak in the disk surface density ( Σ ) lies at somewhat larger distances from the binary and is Figure 1.
Snapshots of the time evolution of material orbiting a central binary with a = 5 R P and e = 0.4. Each panelshows the positions of solids in the ( a eq , e ) plane for q = 7 R P and the times (in yr) indicated in the legend. Cyanpoints show positions of individual tracer particles. Magenta points denote regions of high density; the intensity ofthe color provides a measure of the density. Initially, most of the solids lie on high e orbits with pericenter close tothe apocenter of Charon’s orbit. As the calculation proceeds, collisional damping places more and more solids into adisk. Within the disk, the position of peak surface density moves outward with time. within or close to the satellite zone. Beyond the peak surface density at a = a m , Σ declines rather steeply, Σ ∝ a − p with p ≈ ı ≈
0; a few have ı ≈ a eq ≈ e ≈
1. Tracers with larger a eq have arange of e ; this range grows with increasing a eq . At a eq ≈ R P , some tracers have e ≈ a = 5 R P , e = 0.4, a in = 8 R P , a out = 32 R P , q = 7 R P , and ı = 0.25) with an initial mass of g. To place these results in the context of previous studies, we plot e (later, ı ) as a function of a eq instead of the semimajor axis. Initially, solids have large e ; a eq is then muchsmaller than the semimajor axis. As collisional damping reduces e , a eq ≈ a .In this calculation, solids initially have a broad range of e , from e ≈ R P to e ≈ R P .Without collisional damping, all of these orbits are unstable. However, collisional damping acts rapidly tocircularize these orbits. After 0.01 yr, some solids lie in a thin ring with a eq ≈ R P (upper left panelof Fig. 1). Concurrently, the central binary attempts to eject particles from the system, placing them withinthe ensemble of tracers at a eq ≈ R P and e (cid:46) a eq (cid:46) R P and e (cid:46) a eq move to larger a eq (Fig. 1, upper right panel).Some of these lie within the narrow ring at small a eq ; others lie within a fledgling disk extending from ∼ R P to ∼ R P . Almost all of the rest make up the swath of particles with a range of e at a eq ≈ R P . However, there are a few particles with large a eq and large e . These solids are on their way outof the Pluto–Charon system.From 1 yr to 1000 yr, the evolution turns into a competition between collisional damping and dynamicalexcitation by the central binary. Throughout each orbit, collisional damping tries to circularize particle or-bits. Solids with large e receive a kick from the binary near the pericenters of their orbits. This combinationof processes tends to push particles onto more circular orbits at larger a . As the system evolves, solidsseparate into the two groups illustrated in the middle and lower panels of Fig. 1. A group of particles with e ≈ R P . This peak gradually moves to larger a eq . The second group of particles lies within the broad swath at the upper right of each panel, with large e at small a eq and small e at large a eq . Within this swarm, collisional damping acts to prevent ejections;combined with kicks from the binary, particles evolve to larger and large a eq within the disk. For roughly2/3 of the particles, collisional damping is too slow; these particles are eventually ejected from the system.Of the remaining solids, nearly all lie within the disk. With a final mass of × g, this disk has sufficientmass to produce a few larger objects with the combined mass of the current satellites ( ∼ g, Kenyon &Bromley 2019b).In the lower right panel of Fig. 1, tracer particles outside the disk constitute only 6% of the remainingmass in the system. Based on the evolution from 100 yr to 1000 yr, we estimate that 1/3 to 1/2 of thismaterial will be placed into the disk at a (cid:38) R P . Because these solids will contribute negligibly to thedisk in the satellite zone, we halted the calculation at t ≈ ı = 0.25. In only a few days, a substantial fraction of particles has small inclination and lieswithin the narrow ring at 10–20 R P or a swarm of high inclination particles with a eq ≈ R P . Over thenext year, many particles damp into the disk at radii a eq ≈ R P . Others remain in a high inclinationgroup that moves out to 40–50 R P . Aside from a few outliers, the high inclination group disappears from10 yr to 1000 yr. All other tracers lie within the disk.At 1000 yr, the disk in Fig. 1–2 is well-defined and fairly stable. Inside 100 R P , the disk contains a littlemore than twice the mass of the Pluto–Charon satellite system, M d ( a ≤ R P ) ≈ . × g. With loweccentricity, e (cid:46) ı (cid:46) R P from the system barycenter. For a binary with a = 5 R P and e =0.4, nearly all orbits in this region are unstable (Eq. 2, Table 1). As the evolution proceeds, the binarypushes material to larger distances; collisional damping reduces orbital e and ı to keep this material in thesystem. Over the next 1000 yr, the surface density at 10–20 R P drops from close to 100 g cm − to less than0.1 g cm − .While the binary ejects solids from the unstable region, damped solids fill the disk beyond 50 R P . From0.1 yr to 1 yr, the surface density at 40–100 R P grows by a factor of 100, from − – − g cm − to − – − g cm − . During the next 1000 yr, the surface density in the disk increases by another factor of 10.With little change in Σ during the last 100–200 years of evolution, the surface density profile of the disk isstable at the end of this calculation.1 Figure 2.
As in Fig. 1 for the inclination. As time proceeds, particles with ı ≈ a eq . Some solidsrapidly damp into a thin disk with a small scale height. Others are ejected out of the system with a range of ı relativeto the orbital plane of the central binary. At the end of the calculation, most of the remaining solids lie within thevertically thin disk; a few with large ı are on their way out of the system. Throughout the evolution, material in the inner part of the disk at 10–20 R P is more likely to be ejectedthan to remain in the disk. At 0.1–1 yr, the inner disk contains ∼ g of solids. After 1000 yr, only ∼ × g remains at 20–35 R P . Although some inner disk material finds stability in the outer disk at50–150 R P , these solids are a small fraction of the material ejected from 10–25 R P . At the end of thiscalculation, ∼ g ( ∼ g) lies at 40–60 R P (60–120 R P ). Overall, the mass of the solids falls fromthe initial g to a final . × g. Nearly all of the mass remaining at 1000 yr is in the disk.In addition to a redistribution of mass, the surface density profile gradually becomes smoother with time.After a few weeks of evolution, the profile is irregular and varies chaotically beyond 50 R P . Within ayear, the profile is smoother, with a wavy appearance beyond 30 R P . These waves probably result fromthe action of the central binary, which perturbs solids as collisional damping circularizes their orbits. By1000 yr, the surface density profile is fairly featureless, with only a few small-scale wiggles superimposedon a monotonically decreasing profile.The model 6 calculations follow nearly the same evolution shown in Figs. 1–3. With a smaller q , solidshave a larger initial e but the same initial ı . As the binary attempts to kick solids out of the system, col-lisional damping tries to circularize their orbits. Solids that cross Charon’s orbit feel larger gravitationalperturbations than those that always stay outside of the binary orbit. Larger kicks make it harder for damp-ing to reduce e and ı and to keep these particles in the system. Despite the actions of the central binary,damping reduces inclinations fairly rapidly, generating a thin disk in the Pluto–Charon orbital plane within2 Figure 3.
Surface density profiles at 0.1 yr (orange points), 1 yr (green points), and 1000 yr (purple points) forthe calculation shown in Figs. 1–2. The vertical dashed line marks the position of the innermost stable circular orbitaround the central binary. Vertical solid lines indicate the extent of the satellite zone. M = 1 − × g), collisional damping tends to lag gravitational perturbations from the binary. Although e and ı drop with time, perturbations from the binary push particles to larger and larger a . At larger a ,the solids occupy a somewhat smaller volume due to the smaller e (which limits their radial excursions)and the smaller ı (which limits their vertical excursions), which increase rates and hence the effectivenessof damping. Combined with smaller perturbations from the binary at large a , enhanced damping allowsformation of a disk at a (cid:38) R P . After 1000 yr, remaining solids have formed a disk with a size similarto the solids in model 5. However, these disks have less mass than the model 5 disks. Their shallowersurface density distributions have a larger fraction of mass at larger distances, a characteristic of the initiallylarger kicks from the central binary and the ineffectiveness of collision damping at small semimajor axis, alesssima s .In contrast, calculations with the model 6 initial conditions and large initial masses are able to retain alarger fraction of their initial mass. When the initial mass is larger, collisional damping – which scalesas the square of the mass – is more effective. The vertical extent of the swarm decays more rapidly intothe midplane; more particles lie within the dense ring at 20–30 R P and within the nascent disk at largerdistances. As the calculations proceed, the binary still manages to eject most of the initial mass from thesystem. However, the final disk masses are remarkably similar to those in model 5. Although there are3 Figure 4.
Distribution of surviving particles in ( a eq , e ) space at 1000 yr for calculations with the model 3–6 initialconditions and M = g. Cyan points indicate the positions of individual tracers. Magenta pixels denote regionsof high density; the intensity of the color provides a measure of the particle density in each pixel. Labels list e for thecentral binary, q , and ı . In all panels, the cyan and magenta points delineate an extended disk with an inner edge at ∼ R P , a maximum density at 30–40 R P , and an outer edge at 200–300 R P . The morphology of the points in theupper and the lower middle panels for q = 7 demonstrates that the outcome is fairly independent of ı . Comparisonof the upper two and lower two panels illustrate the impact of q . When ı = 0.25, the outcome at 1000 yr is fairlyindependent of q . At 1000 yr, systems with larger ı have more particles at higher eccentricity than those with smaller ı . Above the disk, particles in limbo have e ≈ a eq and a broad range of e at large a eq . Over time, some ofthese particles will be incorporated into the disk; the central binary will eject others from the system. modest differences among the calculations, massive disks with the model 5 starting conditions have diskswith a physical structure close to the structure of disks with the model 6 parameters (see Table 2).Reducing the initial inclination of the particles ( ı = 0.025, models 3–4) has a modest impact on theresults. With a smaller inclination, the solids occupy a smaller volume. Collisional damping is then moreeffective. At the same time, particles on eccentric orbits see larger perturbations in the gravitational potentialdue to the central binary. As particles experience larger accelerations from the binary, collisional dampingcircularizes their orbits somewhat more rapidly compared to solids in the model 5–6 calculations. Overall,these two processes approximately cancel: the model 3 and 4 calculations retain roughly the same amountof mass as their model 5 and 6 counterparts; the overall extent and mass distribution within the disks arealso fairly similar.Fig. 4 compares the distribution of e as a function of a eq for surviving particles at 1000 yr in calculationswith the model 3–6 initial conditions and M = g. The upper two panels show results for models 3and 4; the lower two panels plot outcomes for models 5 and 6. All four panels have a similar morphology.4 Figure 5.
As in Fig. 4 for the inclination. At 1000 yr, all calculations show a clear disk in the midplane of the Pluto–Charon binary; the disk begins at an inner edge of 20 R P , reaches peak surface density at 30–40 R P , and extends to200–300 R P beyond the right edge of the panels. Above the disk, particles in limbo have ı ≈ a eq andsmaller ı at large a eq . Over time, some of these particles will be incorporated into the disk; the central binary willeject others from the system. A large swarm of particles lies in the orbital plane of the binary. The inner edge of this distribution has a eq ≈ R P ; the outer edge is beyond the right limits of the plot. Each swarm has a density maximum at30–40 R P . The position of the density maximum varies little from one panel to the next; there is no obviousdependence of the disk morphology on q or e .In all four panels, there is a second swath of particles extending from ( a eq , e ) = (30, 1) to larger a eq . Atlarge a eq , the particles have a broad range of eccentricity, e ≈ q = 7 R P clearly have fewerparticles in this swarm than those with q = 5 R P . At a eq ≈ R P , there is also a clear gap betweenthe disk particles and those outside the disk. In the panels for q = 5 R P , this gap is occupied by otherparticles in limbo.Based on the evolution up to 1000 yr, we anticipate that many of the particles in limbo will eventuallybecome part of the disk. In all of these calculations, the lower edge of the swarm in ( a eq , e ) space formsa nearly straight line from ( a eq , e ) ≈ (35, 1) to ( a eq , e ) ≈ (140–180, 0.1–0.2). Surviving solids above thislocus will have final disk radii beyond 140–180 R P and contribute little mass to the disk.Fig. 5 repeats Fig. 4 for the inclination. As in Fig. 4, all of the panels display a similar morphology. A largeswarm of particles with ı ≈ ∼ R P to beyond 150 R P . With ı ≈
0, these disk particles lie5in the Pluto–Charon orbital plane. Peak surface density within this disk is at 30–40 R P . Unlike Fig. 4, theupper three panels in Fig. 5 show few particles outside of the Pluto–Charon orbital plane . Although someof these survivors might eventually become part of the disk, the binary will probably eject most of them.The lowest panel in Fig. 5 has a large group of particles from ( a eq , ı ) = (40, 0.25) to ( a eq , ı ) = (150, 0.).These particles are also in limbo; the competition between collisional damping and excitation from thebinary is more balanced than for the stable particles already in the midplane and for the less fortunate solidsejected from the system. From an analysis of the evolution at 10–1000 yr, particles within this swath willeventually land in the disk at large a eq or will be ejected from the Pluto–Charon system. The trend in ı with a eq indicates how evolution into the disk proceeds: as damping reduces ı , the binary pushes a particleto larger a eq . The location where the swath crosses the midplane is the location where particles damp toroughly zero inclination. As the evolution proceeds beyond 1000 yr, we expect that some of particles withhigher inclination will damp into the disk at a eq larger than 150 R P . Because our main interest is on thelocation of solids near the satellite zone, we chose not to follow this evolution in detail.Fig. 6 compares the surface density distributions for models 3–6 and M = g. All exhibit the samegeneral trend: a sharp edge at 20 R P , a clear maximum at 30–35 R P , and a gradual decline out past 150 R P .Aside from this trend, there are several clear differences. Models with q = 7 R P ( q = 5 R P ) have a broad(narrow) peak. The shape and location of the peak does not depend on ı . At large a eq , the decline of thesurface density is fairly independent of q and ı ; in all calculations Σ declines by a factor of 300–500 from30–35 R P to 150 R P .Although three of the surface density distributions have a monotonic decline from 35 R P to 150 R P ,results for model 6 show a clear wave at 80–150 R P . In this example, the surface density profile evolves asshown in Fig. 3. After one year, disk material lies close to the innermost stable orbit. As the system evolves,most of these solids are ejected. However, some solids in the inner disk move onto stable orbits well outsidethe innermost stable orbit. As the outer disk stabilizes, the surface density has a wavy appearance generatedby the interplay between the binary trying to excite the orbits and collisional damping working to maintainsolids on circular orbits. During these interactions, solids find a stable a and e . The time scale for the wavesto damp (and for the surface density to become smooth) depends on the properties of the inner binary andthe initial q of circumbinary solids. In less eccentric binaries, the waves damp more rapidly compared tothe waves in disks surrounding more eccentric binaries. Similarly, waves generated by solids with initiallysmaller q take longer to damp than those with initially larger q .In the example shown in Fig. 6, the wave results from particles with initially large q . At later times, theparticles have modest eccentricity, e (cid:46) ı (cid:46) − ; collisional damping hasnot quite circularized their orbits. These particles are visible as a thickening of the e vs a eq distributionat large a eq in the lowest panel of Fig. 4. In the other calculations shown in Fig. 4, collisional dampingoperates somewhat faster and circularizes nearly all of the orbits in 1000 yr. In model 6, another 1000 yrof evolution would allow collisional damping to finish circularizing the orbits of these particles; the surfacedensity distribution will then follow the other results more closely.Although mean-motion resonances (MMRs) may eventually appear in these surface density distributions,none are present after 1000 yr of evolution. For the calculations illustrated in Fig. 6, low order MMRs suchas the 3:1 to 7:1 lie inside the innermost stable orbit for an e = 0.4 binary. In these four examples, the14:1, 15:1, 16:1, and 17:1 MMRs at 29–33 R P are within the region of peak surface density at ∼ R P .However, the time scale for MMRs to develop in this region is (cid:38) Figure 6.
Surface density distributions for model 3–6 calculations with M = g. The legend indicates ( a, e ) forthe central binary and ( q , ı ) for the solids in each calculation. Vertical grey lines show the boundaries of the satellitezone. The surface density rises from a sharp inner edge at ∼ R P to a sharp peak at 35 R P ; the surface densitythan falls over two orders of magnitude from the peak to 150 R P . Ripples in the surface density at 80–150 R P for themodel with ( q , ı ) = (7, 0.25) are typical of disks that have not quite settled into an equilibrium state. et al. 1997; Peale 1999; Cheng et al. 2014a; Correia 2020), lower order MMRs will sweep through thisregion (Ward & Canup 2006; Lithwick & Wu 2008; Cheng et al. 2014b). Collisional damping can stabilizesolids against the perturbations induced by these MMRs (Bromley & Kenyon 2015b).To illustrate the impact of the initial a i and q on the evolution of solids, Fig. 7 plots the ”survivor fraction”as a function of the initial semimajor axis of solids. For solids with initial semimajor axis at 8–10 R P , from2% ( M = 10 g) to 10% ( M = 10 g) remains in the system after 1000 yr. These fractions are amazinglyindependent of the initial q . At 10–20 R P , the survivor fraction remain fairly independent of initial a and q for a given initial mass. At 20–25 R P , however, systems with larger initial q retain a larger and largerfraction of solids as the initial a increases. In contrast, systems with an initial q at 5 R P have a low survivorfraction independent of initial a . Clearly, systems with larger q retain much more of their initial mass thanthose with smaller q .The variation of survivor fraction with initial a and q is a result of the interplay between gravitationalforcing from the binary and collisional damping. In the restricted three-body problem, perturbations in thepotential due to the binary are a strong function of the semimajor axis, with leading terms ∝ a − insideCharon’s orbit and ∝ a − outside the binary (Lee & Peale 2006; Leung & Lee 2013; Bromley & Kenyon2015a, 2021). Near the innermost stable orbit, the binary pumps up the eccentricity of circumbinary solids,or ejects then directly during strong, close encounters. For a Pluto–Charon binary with a = 5 R P and e =7 Figure 7.
Fraction of tracers that survive 1000 yr of evolution as a function of their initial semimajor axis for abinary with a = 5 R P , e = 0.4, and various initial masses as indicated in the legend. Solid (dashed) curves showresults for systems of solids with pericenters near the apocenter of Charon’s orbit (with pericenters midway betweenthe apocenter and pericenter of Charon’s orbit), indicated by the vertical grey lines). Although the ability of the centralbinary to solids at 8–25 R P is fairly independent of q , the retention of solids on more distant, more eccentric orbitsdepends on q . a = 18 R P ; inside this distance, objects on completely circular orbitscannot survive. Solids on eccentric orbits are rapidly ejected. At larger distances, near the peak of thesurface density distribution ( a ≈ R P ), damping is fast enough to preserve solids even when their initialeccentricity is high.Collisional damping depends on the particle density n and the relative velocity v as n v . Initially, particleshave large e and occupy a large volume; collisional damping rates are then low. The binary proceeds to ejectparticles by increasing their a and e . As a grows, however, particles spend less time near the binary. Atlarge a , collisional damping dominates perturbations from the central binary. Damping lifts the pericentersof particle orbits, reducing binary perturbations near pericenter. As orbits circularize, n grows faster than v falls; damping rates grow rapidly. Throughout circularization, particles orbit at larger a due to binaryforcing, but with smaller e and larger q . Although moving to larger a reduces perturbations from the binary and damping rates, damping tends to dominate when particles have larger a and smaller e . Thus, orbitscircularize at semimajor axes much larger than the original semimajor axis.Although damping is effective at large a when solids have pericenters near Charon’s apocenter (e.g., q ≈ Q c ), the binary dominates when q (cid:28) Q c (e.g., Lee & Peale 2006; Leung & Lee 2013; Bromley &Kenyon 2015a, 2021). Inside Charon’s orbit, most orbits are chaotic except for a few stable regions nearPluto (Winter et al. 2010; Giuliatti Winter et al. 2013, 2014, 2015). When solids with a (cid:38) R P travel8inside Charon’s orbit, strong, repeated interactions with the binary carry a high risk of collision or ejections.For some tracers, damping raises the pericenter sufficiently to escape these outcomes. For most tracers,however, close encounters dominate; the binary ejects the tracer from the system before damping acts. Inthis way, solids with initial q (cid:28) Q c have a small survivor ratio independent of initial a .Tracer survival also depends on the initial mass of the swarm. Because damping depends on N ik , systemswith more mass damp more effectively. With all orbits passing close to (or within) the binary, very efficientdamping has the drawback of circularizing many orbits at semimajor axes a (cid:46) a s . These orbits are unstable;thus, these particles are ejected rapidly. More massive systems of solids damp a larger fraction of their massinside the region where orbits are unstable. In this way, more massive systems lose a larger fraction of theirinitial mass to dynamical ejections.Calculations with the model 7–8 parameters begin with the same binary parameters and the same initialmasses and inclinations as models 3–4, but have a different radial distribution of solids. With an inner radiusof 19 R P instead of 8 R P , all of the solids begin with large eccentricities, e ≈ e ≈ e allow the binary more chances to eject particlesbefore collisions circularize orbits. Compared to the results for model 5 (model 6), the model 7 (model 8)calculations retain less mass in the circumbinary disk. Aside from the lower mass, the disks have roughlythe same extent and radial distribution of mass in all of the models (see Table 2).Despite differences in the amount of retained mass, all of the model 7–8 calculations evolve in a similarway. Solids near the center-of-mass receive an immediate kick from the binary. For solids farther away, col-lisional damping has a few days to reduce e and ı before the gravity of the binary impacts their trajectories.After the first few orbits, collisional damping eases some solids into a circumbinary ring at 20–30 R P anda few others into a disk at 30–100 R P . As the ring and disk form, the binary ejects many other solids fromthe system. Within the first year or two, the ring and disk are fairly well-defined; as times goes by, the ringexpands and damping places more solids into the disk.Fig. 8 compares the remarkably similar surface density distributions at t = 1000 yr for models 3, 4, 7, and8 and M = 10 g. All have negligible mass inside 20 R P , a sharp peak at 30–40 R P , and a steep declinein Σ from the peak to 200–300 R P . Distributions for models 7 and 8 are somewhat wavier at 100–200 R P than for models 3 and 4. Another 1000 yr of damping would smooth out any waves in all of the surfacedensity distributions. Despite the wave, it is clear that these systems have much of their mass inside thesatellite zone.Calculations with larger initial mass in solids follow the same evolution path as those with the lowestmass. When M is larger, collisional damping is more effective; orbits circularize more rapidly. Becausethe solids evolve towards orbits with smaller a and shorter periods, the binary has more opportunities toeject solids from the system. Over the first few years of evolution, systems with larger M suffer moreejections than those with smaller M . Survivors are placed on orbits with larger a .Fig. 9 shows density distributions at 1000 yr for model 5 calculations with different M . At a eq = 10–25 R P , the five calculations have a sharp edge, where Σ grows by 3–4 orders of magnitude. Each disk has aclear maximum in Σ at 30–60 R P ; the position of the peak moves to larger a eq when M is larger. Beyondthe peak, Σ falls rapidly with a eq . For several calculations, the drops in σ is somewhat wavy. Over time,these waves decline in amplitude.In these calculations, the slope of the surface density distribution is a strong function of the initial mass.In systems with the smallest M ( g), Σ roughly follows a power-law with slope p ≈ Figure 8.
Surface density distributions for model 5–8 calculations with M = g. The legend indicates ( a, e ) for the central binary and ( q , ı ) for the solids in each calculation. The vertical grey lines show the boundaries of thesatellite zone. The surface density rises from a sharp inner edge at ∼ R P to a sharp peak at 35 R P ; the surfacedensity than falls over two orders of magnitude from the peak to 150 R P . Ripples in the surface density at 80–150 R P for the model with ( q , ı ) = (7, 0.25) are typical of disks that have not quite settled into an equilibrium state. mass grows, the slope of the power-law drops; for M = g, p ≈ Σ( a ) relative to systems withless initial mass.Despite the sharp peaks in Σ , surviving solids beyond the satellite zone contain a large fraction of the finalmass. In the most massive models with M = g, roughly equal amounts of mass lie within the satellitezone and at a eq (cid:38) R P . In model 5, for example, the mass in the satellite zone is roughly 100 times themass of the Pluto–Charon satellites. In the least massive models, the sharper drop in Σ( a eq ) leaves a smallerfraction of the final mass in solids beyond the satellite zone.Outcomes of calculations with different combinations of ( a, e ) for the Pluto–Charon binary and q for thesolids lead to similar results. For each model and M , Table 2 summarizes the mass remaining in activetracers at 1000 yr ( M f ), the mass in disk material with e ≤ M d ), the inner and outer edge of the disk( a , a ), the location of peak surface density ( a m ), and the maximum surface density ( Σ m ). Disk edgesare defined as the point where the surface density first rises above − g cm − (inner edge) or where thesurface density first drops below − g cm − outside the peak surface density.The Table collates several interesting features summarized in the discussions of models 3–8. The ratio of M f and M d to the initial mass M declines with M . Systems with a small initial mass, M = 10 g, place0 Figure 9.
As in Fig. 8 for model 5 and various M as indicated in the legend. Systems with larger initial mass havesurface density maxima at larger semimajor axes.
10% to 30% of small solids into a disk. When the initial mass is 100 times larger, only ∼
5% of the initialmass remains in the system. In all cases, the inner edge of the disk is 3–4 times the separation of the innerbinary, with a ≈ R P when a = 5 R P and a ≈ R P when a ≈ R P . Peak surface densityoccurs at roughly twice a . Although a and a m are independent of M , the outer edge of the disk correlateswell with M , ranging from a ≈ R P for M = 10 g to a ≈ R P for M = 10 g.Systems with large disk radii have significant amount of solids outside the satellite zone.To highlight one outcome from Table 2, Fig. 10 summarizes differences in the final surface density dis-tributions with different central binaries. When the Pluto–Charon binary is compact and nearly circular,peak surface density is close to the inner edge of the satellite zone. As e grows, peak Σ moves outward.In Pluto–Charon binaries with ( a, e ) = (5 R P , 0.6), the surface density reaches a broad maximum at 50–60 Pluto–Charon, roughly coincident with the orbits of Kerberos and Hydra. For e = 0.8, the edge of thedisk rests in the middle of the satellite zone; Σ has a sharp peak at ∼
85 Pluto–Charon, well outside thesatellite zone.In each calculation at t = 1000 yr, very little of the disk lies within regions where circumbinary orbits areunstable (e.g., Holman & Wiegert 1999; Doolin & Blundell 2011). In numerical simulations of circumbinarytest particles, the innermost stable orbits for systems with mass ratios similar to Pluto–Charon lie at 14 R P ( a, e = 5 R P , 0.2), 17–18 R P ( a, e = 5 R P , 0.4), 38–42 R P ( a, e = 10 R P , 0.6), and 46 R P ( a, e = 10 R P ,0.8). Most of the results in Table 2 place the inner edge of the disk in the stable region of the system (seealso Table 1). In the few examples where disk material lies within an unstable region, the binary will ejectthese solids within another few thousand years (Kenyon & Bromley 2019a).1 Table 2.
Disk Properties at 1000 yr a Model M M f M d a a a m Σ m Model M M f M d a a a m Σ m a For each model number, the columns list the initial mass in solids M , the final mass in solids M f , and the final mass in the disk, M d , inunits of g; the inner and outer disk radius, a and a (in units of R P ); the radius of maximum surface density, a m , in R P , and themaximum surface density Σ m in units of g cm − . Although the binary clearly sets the inner edge of the disk, collisional damping and the gravitationalpotential of the binary combine with the initial conditions to establish the overall disk morphology. Incompact binaries with a = 5 R P , some circumbinary solids initially have orbits with semimajor axes that arestable for particles on circular orbits. Others have semimajor axes within the unstable zone. As collisionaldamping circularizes the solids, the binary potential encourages them to move onto orbits with larger a .Over 1000 yr, the binary evacuates the volume where orbits are unstable, leaving behind solids on stableorbits. Because solids were nearly on stable orbits at t = 0, the final orbits for many solids lie close to theinner edge of the stable region.2 Figure 10.
Comparison of surface density distributions at 1000 yr for calculations with M = g and thecombinations of e , a in , and q indicated in the legend. The vertical grey lines indicate the boundaries of the satellitezone. In compact Pluto–Charon binaries with a = 5 R P and e = 0.2–0.4, the surface density reaches a maximum at theinner edge of the satellite zone, ∼ R P . In wider binaries with a = 10 R P and e = 0.6–0.8, peak surface densitylies at larger a eq , 50–60 R P when e = 0.6 and 80–90 R P when e = 0.8. Systems with e = 0.8 have peak Σ outside thesatellite zone. In contrast, solids orbiting less compact, more eccentric binaries all initially lie within the unstable region.Although collisional damping has a somewhat bigger challenge in circularizing the orbits of these solids,wider binaries with longer orbital periods push on the solids more gently. As the solids move farther andfarther away from the binary, damping has time to circularize the orbits. Because the unstable region aroundthe binary is large, solids have to find stable orbits well outside the binary.In these examples, particles on circular orbits at the inner edge of the satellite zone are stable (unstable) inbinaries with a = 5 R P and e = 0.2–0.4 ( a = 10 R P and e = 0.6–0.8). This difference results in very differentdisk morphologies: circumbinary disks in compact binaries have significant mass within the satellite zone;disks in wider binaries have most of their mass outside the satellite zone.To conclude this section, Fig. 11 compares results for Σ from model 5 (green symbols, M = g) andmodel 11 (orange symbols, M = g) with the surface density distribution for debris from the impact ofa TNO with Charon (purple symbols; Bromley & Kenyon 2020). In the TNO collision, 60% to 70% of the g ejected by the impact survives to form a circumbinary ring. Although this ring contains less materialthan required to build the known circumbinary satellites, the surface density maximum coincides with themaximum of disks surrounding a Pluto–Charon binary with e = 0.6. The inner edge of this ring also roughlyoverlays the inner edge of the disk around the more eccentric binary.3 Figure 11.
Comparison of surface density distributions generated from the impact of a TNO on Charon (BK2020,purple symbols Bromley & Kenyon 2020) with the final states of debris surrounding the giant impact that formsthe Pluto–Charon binary. Green (orange) symbols indicate results for M , e, a in , q = g , . , R P , R P ( g , . , R P , R P ). The vertical grey lines indicate the boundaries of the satellite zone. The sharp outer edge of rings produced from TNO collisions with Charon provides a sharp contrast withthe disks generated from giant impacts. From the peak at 80 R P to 100 R P , the surface density withinthe ring drops by nearly three orders of magnitude. Within disks from the giant impact, the decline in Σ is a more modest factor of ten. The gradual decline in Σ is a hallmark of all of the extended disks in ourcalculations. DISCUSSIONOver the past decade, the giant impact hypothesis has become the most popular model for the formationof the Pluto–Charon binary planet (e.g., McKinnon et al. 2017; Stern et al. 2018; McKinnon et al. 2020, andreferences therein). In the commonly accepted picture, two planets orbiting the Sun within the protosolarnebula suffer a glancing collision (Canup 2005, 2011). The lower mass secondary loses momentum and iscaptured by the more massive primary; debris from the collision orbits the newly-formed, compact binary.In a more recent option, the planets first form a very wide binary (Rozner et al. 2020). Over time, the orbitdecays, which facilitates a collision with a geometry similar to that in a giant impact.Immediately after the collision, the circumbinary debris lies on very eccentric orbits that are intrinsicallyunstable. Within this debris, solid particles follow one of two dynamical paths (Kenyon & Bromley 2014b;Walsh & Levison 2015; Kenyon & Bromley 2019c; Bromley & Kenyon 2020). Solids that do not physically4collide with other solids cannot circularize their orbits. Within a month, the binary either ejects these solidsor accretes them. When an orbiting solid collides with another solid, the large relative collision velocitiesguarantee catastrophic disruption; particles within the debris are then much smaller than the original solids.After only a few such collisions, 1–10 km solids can be ground down into 1–10 m solids (Kenyon &Bromley 2014b; Bromley & Kenyon 2020). Because small solids damp effectively, this process allowssolids to circularize their orbits on short time scales.In a first attempt to chart the evolution of solids following the giant impact, Walsh & Levison (2015)performed a set of numerical calculations of high eccentricity particles with initial semimajor axes betweenthe 3:1 and 5:1 orbital resonances with Charon . Although some solids eventually cross Charon’s orbitand are ejected from the binary, initial orbits have pericenters well outside the apocenter of Charon’s orbit.Within this suite of calculations, systems with sufficiently large numbers of inter-particle collisions perorbit develop circumbinary disks with relatively few losses from dynamical ejections. The combination ofcollisional damping and dynamical interactions with Charon generate extended disks on stable circumbinaryorbits. For binaries with e = 0.0–0.3, disk material orbits outside the 4:1 orbital resonance with Charon andinside the 7:1 resonance. Disks orbiting more eccentric binaries are more extended than those orbitingmore circular binaries. With little loss in mass from dynamical ejections, all disks have enough mass in thesatellite zone to grow several 10–20 km moons. The lack of mass outside the orbit of Hydra eliminates thepossibility of satellite growth outside the satellite zone.The calculations described here complement and extend these results. In models 1–14, the solids beginon eccentric orbits that cross Charon’s orbit (models 2, 4, 6, 8, 10, 12, and 14) or have pericenters at theapocenter of Charon’s orbit (models 1, 3, 5, 7, 9, 11, and 13). Compared to Walsh & Levison (2015), thecentral binary has a broader range of eccentricity, e = 0.2–0.8 instead of e = 0.0–0.3. The range of initialsolid mass considered here is also somewhat larger than in Walsh & Levison (2015).Despite differences in initial conditions, outcomes are similar. In all of the 70 calculations discussed here,the binary ejects a substantial fraction of solids from the system. Compared to Walsh & Levison (2015), theinitial orbital states of the solids are much more extreme; thus, collisional damping is unable to circularizethe orbit of every solid in the system. Despite the ejections, disk formation is ubiquitous. The inner structureof the disk is fairly similar to results described in Walsh & Levison (2015), with an inner edge just outsidethe region where circumbinary orbits are unstable and a maximum Σ at somewhat larger radii.However, the disks considered here are much more extensive than those in Walsh & Levison (2015), withouter radii of 150–400 R P . The initial conditions are responsible for this difference. By starting on orbitsthat cross or nearly cross the orbit of Charon, the solids considered here receive much larger kicks fromthe central binary and are pushed much farther away from the system barycenter. Collisional damping stilleffectively circularizes orbits, but generates much larger disks. Within these disks, the mass in solids rangesfrom ∼ e =0.2–0.4, the surface density peak lies at the inner edge of the satellite zone (Fig. 10). These disks havea factor of 5–6 less mass near the orbit of Hydra than near the orbits of Styx and Nix. In contrast, disksorbiting Pluto–Charon binaries with e = 0.6 have most of their mass near the orbits of Kerberos and Hydra For a binary semimajor axis a = 10 R P , the 3:1, 5:1, and 7:1 resonances have semimajor axes of ∼ R P , 30 R P , and 36 R P . e = 0.8 form disks with nearly all of theirmass outside the satellite zone at a (cid:38) R P .Clearly, understanding which disks form satellites with masses and orbital architectures similar to thePluto–Charon satellites requires a more extensive set of calculations, However, results in Kenyon & Bromley(2014b) and Walsh & Levison (2015) provide useful guidance. Among the systems considered here, thosewith e ≈ e = 0.2(0.8), newly-formed moons would need to migrate outward (inward) to reach the same locations as theknown satellites. For all binaries, disks have a fraction of their mass outside the satellite zone. Formationof massive satellites in this region is plausible (Kenyon & Bromley 2014b). If these disks form satellites,the moons must either be small enough to avoid detection by New Horizons (Weaver et al. 2016) or migratefrom large distances into the satellite zone (Kenyon & Bromley 2014b).For any circumbinary disk in the Pluto–Charon system, several evolutionary processes impact satellites asthey grow within the disk. Right after the giant impact, tidal forces begin to circularize and then to expandthe Pluto–Charon orbit (e.g., Farinella et al. 1979; Dobrovolskis et al. 1997; Peale 1999; Cheng et al. 2014a;Correia 2020). If satellites grow rapidly as in Walsh & Levison (2015), orbital resonances sweep throughthe satellite zone (Ward & Canup 2006; Lithwick & Wu 2008; Cheng et al. 2014b). Satellite orbits arethen rapidly destabilized. More slowly growing satellites as in Kenyon & Bromley (2014b) may retainsufficient mass in small particles to stabilize the orbits of large moons (Kenyon & Bromley 2014b; Bromley& Kenyon 2015b). Dynamical calculations of massless tracer particles suggest a ∼ R P . Although gravity of theSun does not affect solids within the satellite zone, Lidov–Kozai oscillations can place particles in the outerdisk on high e orbits. These particles might then collide with solids in the inner disk, adding their massto the satellite zone and enhancing the formation of large moons. We have not included the Sun in thecalculations discussed here, but it is straightforward to do so. Adding the Sun might reduce the extent ofthe circumbinary disks and enable the formation of more massive satellites in the inner disk.Despite the popularity of the giant impact picture for satellite formation in Pluto–Charon, it is plausiblethat satellites form in the debris from an impact between Charon and a TNO (Bromley & Kenyon 2020).Because debris from such a collision naturally forms a circumbinary ring, this model is attractive. Satellitesform only near or within the satellite zone, with a much smaller chance of moon formation outside thesatellite zone. If the TNO collision occurs after the completion of tidal evolution, newly-formed satellitesneed not dodge expanding orbital resonances from the binary. The circumbinary rings of Bromley & Kenyon(2020) are also safe from Lidov–Kozai oscillations.Choosing among possible models requires calculations that include the physical processes outlined above.It seems plausible that at least one configuration of a circumbinary disk or ring will yield satellites similarto the known satellites. Aside from explaining the origin of the satellites, these calculations may constrainthe properties of the giant impact (by establishing a preference for a small range of initial binary a and e ) orthe space density of TNOs after Pluto–Charon reaches its current tidally locked state.6 SUMMARYWe consider the dynamical evolution of circumbinary solids following the giant impact that formed Pluto–Charon. Starting from initial orbits that cross or nearly cross Charon’s orbit, particles respond to the time-varying binary potential and damping from particle-particle collisions. On time scales of 100–1000 yr,the central binary ejects a large fraction of the solids into the Solar System. At the same time, collisionaldamping circularizes the orbits of many other particles. The combined action of damping and gravitationalpushes from Pluto–Charon yields a vertically thin circumbinary disk (see also Walsh & Levison 2015). Thedisk extends from close to the innermost stable circumbinary orbit at 20–40 R P to 200–400 R P . Withinthe disk, the maximum surface density is often within the ‘satellite zone,’ the circumbinary volume thatcontains the orbits of the four small satellites orbiting Pluto–Charon.Although all initial conditions lead to disk formation, outcomes depend on the properties of the Pluto–Charon binary and the initial orbital architecture of the circumbinary debris from the giant impact. Widerbinaries with a = 10 R P and e = 0.6–0.8 retain a larger fraction of the impact debris than more compactbinaries with a = 5 R P and e = 0.2–0.4. Debris fields with less mass, ∼ g, and a smaller fraction oforbits with pericenters inside Charon’s orbit preserve more of their initial mass than those with more massand Charon-crossing orbits. Despite these differences, all calculations retain enough mass to form several5–20 km satellites. However, some configurations may enable the growth of large satellites well outside theorbit of Hydra, which is precluded by New Horizons imaging data.With the calculations of Bromley & Kenyon (2020), there are now two plausible paths for the formationof the Pluto–Charon circumbinary satellites. In the giant impact model, the conditions within circumbinarydisks are favorable for growing small moons in less than 1 Myr. As these moons form, they must survivetidal expansion of the central binary. In a small impact model, debris from the collision between Charonand a TNO evolves into a circumbinary ring of material roughly coincident with the satellite zone. As thebinary has already evolved into its current configuration, satellites generated within this ring need not worryabout tidal expansion.Choosing between these two scenarios requires additional calculations that include tidal evolution andinteraction with the Sun and allow for orbital migration of the satellites. If a robust choice is possible,these calculations could provide additional constraints on the formation of Pluto–Charon by isolating arange of initial a and e that is more conducive to the growth and survival of Styx, Nix, Kerberos, andHydra. Outcomes of impacts between Charon and large TNOs could yield better estimates of the crateringfrequency on Pluto–Charon and the early evolution of large solids in the Kuiper belt (e.g., Singer et al. 2019;Kenyon & Bromley 2020).Some of the data (ascii or binary output files and C programs capable of read-ing them) generated from previous numerical studies of the Pluto–Charon systemare available at a publicly-accessible repository (https://hive.utah.edu/) with these urlshttps://doi.org/10.7278/s50d-w273-1gg0 (Kenyon & Bromley 2019a), https://doi.org/10.7278/S50D-HAJT-E0G0 (Kenyon & Bromley 2019b), https://doi.org/10.7278/S50D-EFCY-ZC00 (Kenyon& Bromley 2019c), https://doi.org/10.7278/S50D4AKFQZFC (Bromley & Kenyon 2020), andhttps://doi.org/10.7278/S50D5Q2MFDBT (Kenyon & Bromley 2020). For this paper, data are available atthe same repository, with the url https://doi.org/10.7278/S50DSSMBHHXN.7We acknowledge generous allotments of computer time on the NASA ‘discover’ cluster. Advice andcomments from M. Geller and an anonymous referee improved the presentation of these results. Portionsof this project were supported by the NASA Outer Planets and
Emerging Worlds programs through grantsNNX11AM37G and NNX17AE24G. REFERENCES
Agnor, C., & Asphaug, E. 2004, ApJL, 613, L157Agnor, C. B., Canup, R. M., & Levison, H. F. 1999,Icarus, 142, 219Andersson, L. E., & Fix, J. D. 1973, Icarus, 20, 279Asphaug, E. 2014, Annual Review of Earth andPlanetary Sciences, 42, 551Asphaug, E., Agnor, C. B., & Williams, Q. 2006,Nature, 439, 155Bromley, B. C., & Kenyon, S. J. 2006, AJ, 131, 2737—. 2011a, ApJ, 731, 101—. 2011b, ApJ, 735, 29—. 2013, ApJ, 764, 192—. 2015a, ApJ, 806, 98—. 2015b, ApJ, 809, 88—. 2020, arXiv e-prints, arXiv:2006.13901—. 2021, AJ, 161, 25Brozovi´c, M., Showalter, M. R., Jacobson, R. A., &Buie, M. W. 2015, Icarus, 246, 317Buie, M. W., Grundy, W. M., & Tholen, D. J. 2013,AJ, 146, 152Buie, M. W., Grundy, W. M., Young, E. F., Young,L. A., & Stern, S. A. 2010a, AJ, 139, 1117—. 2010b, AJ, 139, 1128Buie, M. W., Tholen, D. J., & Grundy, W. M. 2012,AJ, 144, 15Canup, R. M. 2004, ARA&A, 42, 441—. 2005, Science, 307, 546—. 2011, AJ, 141, 35Canup, R. M., & Asphaug, E. 2001, Nature, 412, 708Cheng, W. H., Lee, M. H., & Peale, S. J. 2014a, Icarus,233, 242Cheng, W. H., Peale, S. J., & Lee, M. H. 2014b, Icarus,241, 180Christy, J. W., & Harrington, R. S. 1978, AJ, 83, 1005Correia, A. C. M. 2020, A&A, 644, A94Dobrovolskis, A. R., Peale, S. J., & Harris, A. W. 1997,in Pluto and Charon, ed. S. A. Stern & D. J. Tholen(Tucson: University of Arizona Press), 159Doolin, S., & Blundell, K. M. 2011, MNRAS, 418,2656Durda, D. D. 1996, Icarus, 120, 212Durda, D. D., & Dermott, S. F. 1997, Icarus, 130, 140 Farinella, P., Milani, A., Nobili, A. M., & Valsecchi,G. B. 1979, Moon and Planets, 20, 415Genda, H., Fujita, T., Kobayashi, H., Tanaka, H., &Abe, Y. 2015a, Icarus, 262, 58Genda, H., Kobayashi, H., & Kokubo, E. 2015b, ApJ,810, 136Genda, H., Kokubo, E., & Ida, S. 2012, ApJ, 744, 137Giuliatti Winter, S. M., Winter, O. C., Vieira Neto, E.,& Sfair, R. 2013, MNRAS, 430, 1892—. 2014, MNRAS, 439, 3300—. 2015, Icarus, 246, 339Goldreich, P., Lithwick, Y., & Sari, R. 2004, ARA&A,42, 549Greenberg, R., Hartmann, W. K., Chapman, C. R., &Wacker, J. F. 1978, Icarus, 35, 1Grundy, W. M., Binzel, R. P., Buratti, B. J., et al. 2016,Science, 351, aad9189Holman, M. J., & Wiegert, P. A. 1999, AJ, 117, 621Ida, S., Canup, R. M., & Stewart, G. R. 1997, Nature,389, 353Kenyon, S. J. 2002, PASP, 114, 265Kenyon, S. J., & Bromley, B. C. 2002, AJ, 123, 1757—. 2004, AJ, 127, 513—. 2006, AJ, 131, 1837—. 2008, ApJS, 179, 451—. 2014a, ApJ, 780, 4—. 2014b, AJ, 147, 8—. 2016, ApJ, 817, 51—. 2019a, AJ, 157, 79—. 2019b, AJ, 158, 69—. 2019c, AJ, 158, 142—. 2020, The Planetary Science Journal, 1, 40Kenyon, S. J., Najita, J. R., & Bromley, B. C. 2016,ApJ, 831, 8Lee, M. H., & Peale, S. J. 2006, Icarus, 184, 573Leinhardt, Z. M., Richardson, D. C., & Quinn, T. 2000,Icarus, 146, 133Leinhardt, Z. M., & Stewart, S. T. 2009, Icarus, 199,542Leung, G. C. K., & Lee, M. H. 2013, ApJ, 763, 107Lissauer, J. J. 1987, Icarus, 69, 2498