A Pluto–Charon Sonata: The Dynamical Architecture of the Circumbinary Satellite System

To explore the extent of the plausible discovery space for new satellites in the Pluto–Charon system, we consider 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 kilometers (Kenyon 2002; Bromley & Kenyon 2006; Kenyon & Bromley 2008, 2016; Bromley & Kenyon 2011a, 2013; Kenyon et al. 2016). The ensemble of codes within Orchestra includes a multiannulus 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 orbits of massive satellites and many massless tracers around the Pluto–Charon binary. For the satellites, we adopt the masses, radii, and orbital elements listed in Table 1 (Brozović et al. 2015; Showalter & Hamilton 2015; Weaver et al. 2016) and initial state vectors from Table 8 of Brozović et al. (2015). To derive the position and velocity of Pluto relative to the system barycenter, we adopt an orbital period of 6.387 day and ${m}_{P}$ = 1.303 × 10²⁵ g, ${r}_{P}$ = 1183 km, and f_P ≲ 0.006 for the mass, radius, and oblateness of Pluto; Charon has mass ${m}_{C}$ = 1.587 × 10²⁴ g = 0.12 ${m}_{P}$, radius ${r}_{C}$ = 606 km = 0.51 ${r}_{P}$, and oblateness f_C ≲ 0.005 (e.g., Young & Binzel 1994; Person et al. 2006; Brozović et al. 2015; Stern et al. 2015; McKinnon et al. 2017; Nimmo et al. 2017, and references therein). To simplify assigning orbits for massless tracer particles, we rotate the cartesian coordinate system of Brozović et al. (2015) to place the angular momentum vector L in the z-direction.

In this study, we do not consider how errors in the measured positions and velocities of Pluto–Charon and the smaller satellites might impact outcomes of the calculations. For a nominal distance of 40 au from the Earth, a 0001 error in the centroid of the point-spread function for a small satellite on an HST image corresponds to ∼30 km. Typical uncertainties in the measured positions of Nix/Hydra (Kerberos/Styx) on a single HST image is ∼0002–0003 (001). Assuming Poisson statistics for the 3–10 HST images acquired per epoch, the likely error in the nominal state vector of Brozović et al. (2015) is ≲30–100 km for each of the satellites. Compared to the observed semimajor axes (4.2 × 10⁴–6.5 × 10⁵ km) or the derived radii of their Hill spheres (200–700 km), this uncertainty is small. Several test calculations with initial positions differing by 20–50 km from the nominal positions for the small satellites yield identical results to those starting from the nominal initial state vector. Thus, we consider results only for one initial state vector.

To evolve the ensemble of massless tracers, we place them on orbits with initial semimajor axis a, e, and ${\imath }$ relative to the binary center of mass. These orbits are similar to a set of "most circular" circumbinary orbits with identical orbital elements (e.g., Lee & Peale 2006; Lithwick & Wu 2008b; Youdin et al. 2012; Leung & Lee 2013; Bromley & Kenyon 2015a, 2015b, 2017). Adding circumbinary satellites to the Pluto–Charon binary limits the set of stable orbits for massless tracers (Lithwick & Wu 2008b; Youdin et al. 2012). Our goal is to find the set of stable tracer orbits for an adopted set of properties for the known satellites.

For each suite of simulations, we derive results using a symplectic integrator that divides the Pluto–Charon orbit into N steps and maintains a constant time step throughout the integration. As outlined in the Appendix, several tests demonstrate that a minimum N = 40 enables calculations with negligible drift over 100–500 Myr in a and e for the Pluto–Charon binary and for small satellites with their nominal masses. The algorithm has also been verified with test simulations in previous papers (e.g., Duncan et al. 1998; Bromley & Kenyon 2006, 2011a, 2011b, 2013, 2017).

On the NASA "discover" computer system, we perform calculations on either one processor (six-bodies or seven-bodies) or 56 processors (Pluto–Charon binary with massless tracers, with or without the small satellites). With one processor (cpu), a system with Pluto–Charon and the four known satellites evolves 4.2 Myr per cpu-day. Adding another small satellite reduces the evolution time to 3.0 Myr per cpu-day. At these rates, completing a typical 100 Myr calculation requires 25–34 days. Multiprocessor calculations with the central binary, 56 × 56 = 3136 tracers, and the four small satellites complete 1.9 Myr of evolution per day. As the Pluto–Charon system ejects tracers, the calculations move somewhat faster. Typical 10–15 Myr calculations for these systems finish in 4–8 days.

For computational convenience, we perform calculations of tracers over a small range in semimajor axis a. When the small satellites are included, the range in a covers regions from (i) well inside to just outside the orbit of Styx, (ii) just inside the orbit of one of the small satellites to just outside the orbit of the next satellite with larger a, or (iii) just inside to well outside the orbit of Hydra. Instead of having a uniform density of tracers with a, this procedure generates an overlap of tracers corotating with each small satellite. Aside from allowing us to compare results for tracers with similar a from different calculations, these starting conditions provide a better measure of the survival rate for corotating tracers.

To derive physical properties for the four small satellites (Table 1), we rely on published observations and n-body simulations. Orbital elements—the semimajor axis a, eccentricity e, inclination i, and orbital period P_orb—are from detailed analyses of HST imaging data (Brozović et al. 2015; Showalter & Hamilton 2015).

Comprehensive imaging data from HST and the New Horizons flyby demonstrate that all of the satellites have irregular, oblong shapes (Showalter & Hamilton 2015; Weaver et al. 2016). Aspect ratios are roughly 2:1:1 (Styx and Kerberos), 1.5:1:1 (Nix), or 3:2:1 (Hydra). For numerical simulations, we adopt equivalent spherical radii ${r}_{i}=\sqrt[3]{{x}_{i}{y}_{i}{z}_{i}}$, where x_i, y_i, and z_i are the three dimensions quoted in Weaver et al. (2016). The 1σ errors in these radii are ±2.5 km for Styx, Nix, and Kerberos and ±8.5 km for Hydra.

Robust analyses of the HST imaging data yield satellite masses and 1σ errors (in units of 10¹⁸ g) m_S ≲ 15 (Styx), m_N = 45 ± 40 (Nix), m_K = 16.5 ± 9 (Kerberos), and m_H = 48 ± 42 (Hydra). Fits to the HST astrometry are rather insensitive to the masses of the small satellites (Brozović et al. 2015). Using the adopted radii, the satellites have mass density ρ_S ≲ 25 ${\rm{g}}\,{\mathrm{cm}}^{-3}$, ρ_N = 1.3 ± 1.3 ${\rm{g}}\,{\mathrm{cm}}^{-3}$, ρ_K = 18 ± 10 ${\rm{g}}\,{\mathrm{cm}}^{-3}$, and ρ_H = 1.2 ± 1.2 ${\rm{g}}\,{\mathrm{cm}}^{-3}$. Although the nominal densities for Nix and Hydra are reasonably close to the density of either Charon (ρ_C = 1.70 ${\rm{g}}\,{\mathrm{cm}}^{-3}$) or Pluto (ρ_P = 1.85 ${\rm{g}}\,{\mathrm{cm}}^{-3}$), results for Styx and Kerberos are unphysical.

For the calculations in this paper, we adopt the HST masses for Nix and Hydra. Revising the mass density of Styx and Kerberos to the physically plausible ρ_S ≈ ρ_K ≈ 1 ${\rm{g}}\,{\mathrm{cm}}^{-3}$ yields masses, m_S = 0.6 and m_K = 1, that are consistent with the detailed analyses of HST imaging data. Long-term n-body simulations with these masses yield stable systems over evolution times of 500 Myr. Test calculations show that doubling the adopted masses for Styx and Kerberos does not change the results significantly. We plan to report a detailed study of satellite masses in a separate publication.

For massless tracers orbiting the Pluto–Charon binary with semimajor axes a_S ≲ a ≲ a_H, plausible regions of stability are set by the Hill radius

$$\begin{eqnarray}&&{r}_{H}={\left(\displaystyle \frac{m}{3{m}_{\mathrm{PC}}}\right)}^{1/3}\,a,\end{eqnarray} \tag{ 1 }$$

where m is the mass of a satellite, ${m}_{\mathrm{PC}}$ is the combined mass of Pluto and Charon, and a is the semimajor axis of a nearby satellite. The fifth column of Table 1 lists the Hill radius of each satellite. In systems with several massive planets or satellites, massive objects typically clear out a zone with a half-width $\delta a\approx {{Kr}}_{H}$ around their orbits, with K = 8–10 (e.g., Wisdom 1980; Petit & Henon 1986; Gladman 1993; Chambers et al. 1996; Deck et al. 2013; Fang & Margot 2013; Fabrycky et al. 2014; Kratter & Shannon 2014; Mahajan & Wu 2014; Pu & Wu 2015; Morrison & Kratter 2016; Obertas et al. 2017; Weiss et al. 2018, and references therein).

Naive application of these constraints leaves little space for satellites with a_S ≲ a ≲ a_H. Setting K ≈ 10 precludes other stable satellites between the orbits of Styx–Nix and Kerberos–Hydra, but allows moons in a small region a ≈ 54,000–56,000 km between the orbits of Nix and Kerberos. In terms of the binary separation ${a}_{\mathrm{PC}}$, this stable region has a ≈ 2.75–2.9 ${a}_{\mathrm{PC}}$. The weaker constraint K ≈ 8 expands the stable region between Nix and Kerberos and enables a second stable region at $a\approx 2.22-2.27\,{a}_{\mathrm{PC}}$ (43,500–44,500 km) between the orbits of Styx and Nix. Stable satellites between the orbits of Kerberos and Hydra are still prohibited.

The Hill condition also allows stable satellites with $a\approx 1.84-2.15\,{a}_{\mathrm{PC}}$ inside the orbit of Styx and $a\gtrsim 3.75\,{a}_{\mathrm{PC}}$ beyond the orbit of Hydra. For $a\lesssim 2.15\,{a}_{\mathrm{PC}}$, it seems likely that the Pluto–Charon binary, Nix, and Hydra will eventually drive out massive satellites. Far outside the orbit of Hydra ($a\gtrsim 50\,{a}_{\mathrm{PC}}$), the gravity of Hydra, the Sun, and the major planets combine to preclude stable satellites for some range of a, e, and i (Michaely et al. 2017). Inside this region, orbits are generally stable. Some orbits with $a\approx 1.84\mbox{--}2.15\,{a}_{\mathrm{PC}}$ might be stable; however, we expect satellites with $a\gtrsim 3.75\,{a}_{\mathrm{PC}}$ are stable on longer timescales.

There have been many studies of the stability of circumbinary orbits (e.g., Dvorak et al. 1989; Holman & Wiegert 1999; Pilat-Lohinger et al. 2003; Musielak et al. 2005; Pichardo et al. 2005, 2008; Verrier & Evans 2009; Farago & Laskar 2010; Doolin & Blundell 2011; Kratter & Shannon 2014; Li et al. 2016; Smullen & Kratter 2017), but only a few results are generally applicable to the Pluto–Charon binary (see also Youdin et al. 2012, and references therein). Nearly all orbits with a ≳ 3 ${a}_{\mathrm{PC}}$ are stable. For coplanar, prograde circumbinary satellites with inclination ı ≈ 0° orbiting a Pluto–Charon binary with ${e}_{\mathrm{PC}}\lesssim 5\times {10}^{-5}$ (Brozović et al. 2015), the innermost stable orbit has a semimajor axis a₀ ≈ 1.7–2 ${a}_{\mathrm{PC}}$. Some retrograde orbits (${\imath }\approx \pi$) with a ≈ 1–2 ${a}_{\mathrm{PC}}$ are stable (Doolin & Blundell 2011). Although orbits with a ≈ 1–3 ${a}_{\mathrm{PC}}$ and ${\imath }$ ≫ 0° are generally less stable than their ${\imath }$ ≈ 0° counterparts, there are islands of stability with a ≈ 1–2 ${a}_{\mathrm{PC}}$ and some ${\imath }$.

For very small particles, radiation pressure may destabilize circumbinary orbits with a ≳ 3 ${a}_{\mathrm{PC}}$ (Burns et al. 1979; Hamilton & Burns 1992; Thiessenhusen et al. 2002; Poppe & Horányi 2011; Pires dos Santos et al. 2013; Porter & Grundy 2015). In the Pluto–Charon system, solar radiation has little impact on large particles with radii r ≳ 1 mm. However, particles with r ≲ 10 μm are ejected on short timescales, ≲200–300 yr. At intermediate radii, the ejection time depends on the optical properties of the grains. With a primary focus on the stability of more massive satellites, we forgo computationally intensive calculations of the solar radiation force over the course of 10–300 Myr integrations. Thus, our analysis of "massless" tracers probes the stability of 1 mm and larger particles with negligible self-gravity. We return to the impact of radiation pressure in Section 4.

To identify plausible locations for small circumbinary satellites in the current Pluto–Charon system, we follow sets of massless test particles orbiting the Pluto–Charon binary. In these initial tests, we do not include the four small satellites. Rather than attempt to duplicate previous efforts in complete detail (e.g., Doolin & Blundell 2011), our goal is to locate stable regions for prograde/retrograde orbits with small inclination to the orbital plane of the binary and for polar orbits with initial ${\imath }$ ≈ 90°. To facilitate this goal, we define a "survivor fraction" as the fraction of tracers in orbit after 1–2 Myr of evolution, f_s = N_s/N₀, where N_s is the number of survivors and N₀ is the initial number of tracers. Table 2 summarizes f_s for these calculations.

For the first tests, we examine the orbital evolution of massless tracers with initial eccentricity e₀ = 10⁻⁵ and a range of semimajor axes, a₀ = 1.0–1.5 ${a}_{\mathrm{PC}}$ and a₀ = 1.45–2.1 ${a}_{\mathrm{PC}}$, on prograde (${\imath }$ ≈ 0), polar (${\imath }\approx \pi /2$), and retrograde (${\imath }\approx \pi$) orbits. Because we do not mix tracers with different inclinations, these tests require six distinct calculations with 3136 tracers in each calculation. Within this set, only one tracer on polar orbits survives. A comparison of the orbital evolution of this tracer with others on polar orbits suggests it will be ejected in ≲1 Myr. Thus f_s = 0 for polar orbits with a₀ = 1.0–2.1 ${a}_{\mathrm{PC}}$. Among tracers on prograde orbits, none with a₀ = 1.0–1.5 ${a}_{\mathrm{PC}}$ survive; however, 29% of those with a₀ ≈ 1.45–2.1 ${a}_{\mathrm{PC}}$ remain. Tracers on retrograde orbits are more stable: 20% (98%) of the survivors have a₀ = 1.0–1.5 ${a}_{\mathrm{PC}}$ (1.45–2.1 ${a}_{\mathrm{PC}}$).

Figure 1 shows the distribution of a and e for particles that survive for 1–2 Myr. The lone polar survivor has a ≈ 2.05 ${a}_{\mathrm{PC}}$ and e ≈ 0.04. Among prograde tracers, survivors have a ≳ 1.70–1.75 ${a}_{\mathrm{PC}}$ and e ≈ 0.01–0.1. Most are concentrated in a cloud with a ≈ 2 ${a}_{\mathrm{PC}}$ and e ≈ 0.03. Retrograde survivors have a ≳ 1.3 ${a}_{\mathrm{PC}}$; the typical e ranges from roughly 0.1 at a ≈ 1.4 ${a}_{\mathrm{PC}}$ to 0.01 at a ≈ 2.1 ${a}_{\mathrm{PC}}$.

Zoom In Zoom Out Reset image size

Tracers with larger a₀ are more likely to survive. When a₀ = 1.95–2.65 ${a}_{\mathrm{PC}}$ (a₀ = 2.60–3.25 ${a}_{\mathrm{PC}}$), roughly 2% (10%) on polar orbits survive for 1–2 Myr (Table 2). The survivor fraction is much larger for tracers on prograde orbits—90% for a₀ = 1.95–2.65 ${a}_{\mathrm{PC}}$ and 100% for a₀ = 2.60–3.25 ${a}_{\mathrm{PC}}$. Some of the prograde survivors have orbital elements similar to the known circumbinary satellites. All of the retrograde tracers survive,

Figure 2 shows (a, e) for survivors with a₀ = 1.95–3.25 ${a}_{\mathrm{PC}}$. Prograde and retrograde tracers have a clear trend in e(a), with e ≈ 0.01–0.1 at a ≈ 2 ${a}_{\mathrm{PC}}$ and e ≈ 0.001–0.01 at a ≈ 3.2 ${a}_{\mathrm{PC}}$. Sets of polar tracers have little or no trend in e with a. However, there is an abrupt inner edge to the distribution of polar tracers at a ≈ 2.2 apc. For all sets of tracers, there is a broad range of e at every stable a.

Zoom In Zoom Out Reset image size

The minimum semimajor axes for circumbinary particles— ${a}_{c}/{a}_{\mathrm{PC}}$ ≈ 1.7 for prograde tracers, 2.2 for polar tracers, and 1.30–1.35 for retrograde tracers—agree with previous results. For example, Doolin & Blundell (2011) infer ${a}_{c}/{a}_{\mathrm{PC}}$ ≈ 1.75 (prograde), 2.2 (polar), and 1.3 (retrograde). The Holman & Wiegert (1999) fit to a suite of simulations yields ${a}_{c}/{a}_{\mathrm{PC}}$ ≈ 2 for prograde orbits; retrograde orbits have a smaller a_c (Wiegert & Holman 1997). The somewhat larger a_c for prograde orbits from Holman & Wiegert (1999) matches the location of the high-density cloud of survivors in Figure 1. Given the smaller number of simulations with shorter duration performed by Holman & Wiegert (1999), it is possible that their calculations identified the most likely a_c rather than the true a_c.

We now consider the stability of tracers in systems with Pluto–Charon and the four small satellites. The state vector of Brozović et al. (2015) establishes initial positions and velocities of Pluto–Charon, Styx, Nix, Kerberos, and Hydra. For simplicity, we transform the coordinates to a cartesian system where the orbital plane of Pluto–Charon defines z = 0. Tracers begin on nearly circular (e₀ ≈ 10⁻⁵) prograde orbits with low (${\imath }\approx {10}^{-5}$) or high (${\imath }\approx \pi /2$) inclination. For three sets of calculations, the initial semimajor axis a₀ of each tracer is randomly distributed between 0.975 a_i and 1.025 a_i+1 where i = 1, 2, or 3, a₁ = a_S, a₂ = a_N, a₃ = a_K, and a₄ = a_H. As a shorthand, we identify (i) Styx–Nix tracers (i = 1; a₀ = 0.975 a_S − 1.025 a_N), (ii) Nix–Kerberos tracers (i = 2; a₀ = 0.975 a_N − 1.025 a_K), and (iii) Kerberos–Hydra tracers (i = 3; a₀ = 0.975 a_K − 1.025 a_H). In a fourth (fifth) calculation, tracers have initial semimajor axes ranging from 0.7 a_S to 1.025 a_S (0.975 a_H to 1.125 a_H). The n-body code follows the orbits of tracers until all have been ejected or for 10 Myr.

Among massless tracers on prograde orbits inside the orbit of Styx, roughly 20% survive 10 Myr of dynamical evolution. In the first year, nearly half are ejected; less than one-third remain after 0.5 Myr. At the end of the evolutionary sequence, the system ejects a few tracers per megayear. If this rate is maintained indefinitely, all would be ejected in 100–200 Myr. It seems likely that the rate will slow; thus, some are likely to survive for several gigayears.

The survival rate for tracers placed on prograde orbits outside the orbit of Styx and inside the orbit of Hydra is much smaller (Table 2). During the first Myr, nearly 90% are ejected. Over the next 9 Myr, the ejection rate slows considerably and eventually falls to roughly 1–2 tracers per Myr. After 10 Myr, the survival fractions are f_s = 0.02 (Styx–Nix tracers), 0.004 (Nix–Kerberos tracers), and 0.02 (Kerberos–Hydra tracers). Over the next 100–200 Myr, the small but finite removal rate during the first 10 Myr implies the removal of all tracers initially placed inside the orbit of Hydra.

Tracers starting from just inside the orbit of Hydra to roughly 1.2 a_H have a larger survival rate. During the first 1–2 Myr of evolution, Nix and Hydra eject nearly half of the tracers. After this initial flurry, the ejection rate slows to a trickle. At 10 Myr, f_s = 0.45.

For prograde tracers, the distributions of (a, e) show several clear features. Inside the orbit of Styx, the density of survivors peaks at 1.8–2.0 ${a}_{\mathrm{PC}}$ (Figure 3). Within the corotation zones of Nix and Hydra, a few tracers orbit with e = 0.005–0.05 (Figures 4–5). There are essentially no tracers outside any of the corotation zones or in the corotation zones of Styx and Kerberos. Survivors also have a strong concentration outside the orbit of Hydra, extending over a region bounded by $a\approx 3.6\mbox{--}4.0\,{a}_{\mathrm{PC}}$ and e ≈ 10⁻⁴–10⁻¹, with a strong density maximum at $(a,e)\approx (3.8\,{a}_{\mathrm{PC}},0.005)$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Among tracers in the corotation zones of Nix and Hydra, the evolution of e follows a standard pattern. For several megayears, e gradually grows from ≲0.001 to 0.05. At some point, perturbations by the central binary, Nix, and Hydra generate e ≳ 0.05 and the tracer begins to cross the orbit of either Nix or Hydra. After another few thousand years, the tracer leaves the system.

Within the strong concentration of tracers inside the orbit of Styx and outside the orbit of Hydra, e traces a similar evolution. Tracers that achieve e ≳ 0.02 suffer small oscillations in a and an increasingly larger e until the pericenter of their orbit approaches the orbit of Styx/Nix (tracers originally inside the orbit of Styx) or Hydra (tracers originally outside the orbit of Hydra). The gravity of either Nix or Hydra then ejects them from the system. The upper envelope of the cyan points in the upper right corner of Figure 5 show the clear growth of e with increasing a that characterizes dynamical ejections (see also Duncan & Levison 1997; Gladman et al. 2002).

The evolution of polar tracers is somewhat different. Unlike prograde tracers, polar tracers with orbits that cross within the corotation zone of a massive satellite are rapidly ejected from the system. Tracers with orbits that cross between the satellites are harder to remove. These tracers spend most of their time well outside the "clearing zone" of the massive satellites and take somewhat longer to eject. Still, some regions are cleared rapidly: it takes only 50 kyr to eject 100% (90%) of polar tracers with initial semimajor axes inside the orbit of Styx (between the orbits of Styx and Nix). Despite the slow removal rate for tracers inside the orbit of Nix after 0.1 Myr, all are lost after 6 Myr.

Polar tracers survive more easily outside the orbit of Nix (Table 2). For orbits between Nix–Kerberos (Kerberos–Hydra), it takes 0.1 Myr (0.3 Myr) to reduce the initial number of tracers by 50%. The removal rate then slows considerably. After 10 Myr, the survivor fractions are f_s = 0.21 (Nix–Kerberos tracers) and f_s = 0.14 (Kerberos–Hydra tracers). To examine the removal rate at later times, we extended these calculations to 20 Myr. At this point the survivor fractions have dropped to f_s = 0.18 (Nix–Kerberos) and f_s = 0.08 (Kerberos–Hydra). In both cases, the removal rate at 20 Myr suggests that the Pluto–Charon satellite system will eventually eject all of the tracers with polar orbits between Nix and Hydra, on a timescale of ∼50 Myr for Nix–Kerberos tracers and ∼30 Myr for Kerberos–Hydra tracers.

The longer lifetime for polar tracers between the orbits of Nix and Kerberos is due to the larger Hill spacing factor (K = 16) relative to the K = 10 factor for the Kerberos–Hydra pair. With Hydra's larger nominal mass and its closer orbit to Kerberos, the volume available for extra satellites is much larger between Nix and Kerberos than it is between Kerberos and Hydra. The larger (smaller) volume results in a slower (faster) removal process for tracers on unstable orbits between Nix and Kerberos (Kerberos and Hdyra).

Outside the orbit of Hydra, tracers on polar orbits are much more stable. After 1 Myr (3 Myr), the massive satellites have ejected only 30% (36%) of tracers on polar orbits with a₀ = 3.2–4.0 ${a}_{\mathrm{PC}}$. After 20 Myr, the survivor fraction exceeds 0.5 (Table 2). On the basis of the slow rate of removal for this set of polar tracers, many will survive for ≳500 Myr.

Despite differences in the dynamical evolution between tracers on polar orbits and tracers on low-inclination prograde orbits, the distributions of (a, e) have some common features (Figures 6–7). At any time, the range in e for survivors is very large, with a maximum e of roughly 0.05 for prograde orbits and roughly 0.01 for polar orbits. Tracers excited to larger e are rapidly ejected. In both sets of calculations, dynamical evolution generates a dense cloud of tracers with $(a,e)\approx (3.8\,{a}_{\mathrm{PC}},0.05)$. Within this cloud, the maximum stable e is 0.01–0.02 instead of 0.05–0.10.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To conclude this section, Figures 8 and 9 display the configuration of tracers within the Pluto–Charon satellite system after 10 Myr. Animations associated with each figure illustrate the time evolution of the complete population. Because we calculate the evolution of tracers in bands that overlap each satellite, the overall population of tracers is somewhat larger along the orbit of each satellite than in the volume between the satellites.

For the ensemble of prograde tracers, it takes 100–300 yr for Nix and Hydra to begin clearing out particles along their orbits. Slight density maxima appear along the orbits of the lower mass satellites. Over roughly 3000 yr, Nix nearly clears out material in its Hill volume, except for a narrow corotation zone. With a longer orbital period, Hydra clears out a similar region in 10–30 kyr. During this period, the satellites eject tracers in the density maxima along the orbits of Styx and Kerberos. As the evolution proceeds, the satellites clear out the Hill volumes of Styx and Kerberos (including their corotation zones), while Nix and Hydra continue to work on removing tracers inside their corotation zones (0.1–0.3 Myr). Once the evolutionary sequence is complete, there are large sets of tracers just inside the orbit of Styx and just outside the orbit of Hydra, along with a few tracers in the corotation zones of Nix and Hydra. Otherwise, the system is fairly empty (Figure 8).

The evolution of polar tracers shows several contrasting features (Figure 9 and associated animation). Without the four circumbinary satellites, polar tracers inside the orbit of Styx are unstable (Figure 1). Nearly all of these tracers disappear within 1000 yr. As the binary ejects these tracers, Nix and then Hydra begin to clear polar tracers that intersect their orbits. Unlike systems of prograde tracers in the orbital plane of the binary, polar tracers orbiting in the corotation zones of each satellite are ejected as rapidly as other tracers within the Hill volume. After 10 kyr (100 kyr), the Hill volume of Nix (Hydra) is nearly empty. As Nix completes clearing its Hill volume, the central binary ejects the last few tracers inside the orbit of Styx; Kerberos also begins to eliminate tracers from its orbit. At 1 Myr, the volume inside the orbit of Nix and the Hill volumes of Kerberos and Hydra are nearly devoid of tracers. After 10–20 Myr, many tracers remain between the orbits of Nix–Kerberos and Kerberos–Hydra and well outside the orbit of Hydra.

As the Pluto–Charon binary and small satellites clear polar tracers out of the corotation zones and the region inside the orbit of Nix, they also clear out material near the 9:2 resonance with the orbit of the central binary. Visible as a narrow dark band in the positions of tracers between the orbits of Nix and Kerberos in Figure 9, the density of tracers in this region is roughly 60% of the density in the rest of this group. The outer edge of the distribution of tracers between the orbits of Kerberos and Hydra is close to the 11:2 resonance. However, there are no tracers between this resonance and the orbit of Hydra. Outside the orbit of Hydra, there is another drop in density at the 7:1 resonance with the binary. This drop is not visible as a dark band in Figure 9; yet, the density of particles is roughly 70% of the density in the rest of the group outside the orbit of Hydra. Although our calculations did not examine the evolution of polar tracers with more distant orbits, the central binary and the small satellites remove polar tracers at the n:1 resonances (with n = 2–7) or the n:2 resonances (with n = 3, 5, 7, 9, 11, and 13) on timescales of ≲1–10 Myr.

Although the calculations with massless tracers establish likely locations for small particles in the Pluto–Charon system, they do not allow us to constrain sets of stable orbits for small satellites with mass. To test whether massive satellites can have as long lifetimes as massless tracers, we select a group of massless, prograde test particles that have survived 10 Myr of circumbinary evolution (Figures 3–5), assign each a radius of 2 km and a mass of 4 × 10¹⁶ g (ρ_s = 1.2 ${\rm{g}}\,{\mathrm{cm}}^{-3}$), and evolve them with Pluto–Charon, Styx, Nix, Kerberos, and Hydra. To avoid gravitational interactions among the survivors, we follow the evolution of only one survivor in addition to the central binary and the four known satellites of Pluto–Charon.

For each set of survivors, the lifetime as a massive satellite is a strong function of initial conditions (Table 2; Figure 10). Among those with orbits inside the orbit of Styx (near the corotation zone of Nix), survival times range from 1 yr to ≳300 Myr (30 yr to ≳300 Myr). Satellites starting within the corotation zone of Hydra last somewhat longer, from roughly 20 kyr to ≳300 Myr. Despite the large range in lifetimes, 90% of the satellites initially inside the orbit of Styx and near the corotation zone of Nix are ejected within 100 Myr; 79% of satellites initially near the corotation zone of Hydra are ejected after 100 Myr. To examine the survival rate at later times, we extended calculations for survivors to 300 Myr. At this point, the survivor fractions are f_s = 0.03 (satellites inside the orbit of Styx), 0.07 (satellites corotating with Nix), and 0.10 (satellites corotating with Hydra). On the basis of the gradual increase in e and ${\imath }$ for these survivors at 200–300 Myr, the likelihood that 2 km satellites inside the orbit of Hydra survive for the age of the solar system is small.

Zoom In Zoom Out Reset image size

Massive moons on prograde orbits outside the orbit of Hydra are stable. After 100 Myr of circumbinary evolution, ∼3.5% of the satellites are ejected. Another 7% are ejected after 150 Myr. For each surviving satellite, there are no obvious trends in the evolution of a, e, or ${\imath }$. Although it is possible that the orbits of some satellites might suffer significant perturbations at later times, the long lifetimes of these moons suggests that many are stable on 4–5 Gyr timescales.

To investigate the potential for small moons on polar orbits, we select a group of tracers that survive 10–20 Myr of evolution on polar orbits among the known Pluto–Charon satellites (Figures 6–7). As with prograde moons, we assign each a radius of 2 km and evolve each one together with the known satellites. Because polar tracers do not survive inside the orbit of Nix, we only consider the evolution of small moons between the orbits of Nix–Kerberos and Kerberos–Hydra and outside the orbit of Hydra.

Between the orbits of Nix–Kerberos and Kerberos–Hydra, several small moons survive 300 Myr of dynamical evolution (Table 2; Figure 11). Moons orbiting between Nix and Kerberos often have short lifetimes of 1–10 kyr; others survive for less than 1 Myr. Roughly 21% (14%) complete a 100 Myr (300 Myr) calculation on a stable orbit. With a minimum lifetime of roughly 1 Myr, moons on polar orbits between Kerberos and Hydra generally last somewhat longer. However, few survive for 10 Myr. After 100 Myr (300 Myr), only ∼3.5% (0%) are still on fairly stable orbits.

Zoom In Zoom Out Reset image size

Outside the orbit of Hydra, small moons on polar orbits generally survive for 300 Myr. For these satellites, perturbations from the known moons are fairly small. Orbiting well inside the Pluto–Charon Hill sphere, these satellites are also fairly immune to jostling from the giant planets and other passersby. As with small moons on prograde orbits outside the orbit of Hydra, we expect that many of these can survive for ≳1 Gyr.

Aside from the ability of the symplectic integrator to maintain the orbits of Pluto–Charon and surrounding satellites (see the Appendix), these calculations have several limitations. In simulations with N_t = 3136 massless tracers, the accuracy of the derived survival rates (f_s; Table 2) from Poisson statistics ranges from ∼0.1% (when f_s ≲ 0.05) to ∼2% (when f_s ∼ 1). Although adding more tracers reduces the uncertainty in f_s, significant improvement requires a factor of 10 larger N_t, which is computationally expensive.

To complete a large set of calculations in a modest period of real time, we limit the duration of each simulation to a time that is much shorter than the age of the solar system. Shorter simulations limit the accuracy of f_s. For systems with massless tracers, the ejection rate for tracers at 10–30 Myr is small and dominated by shot noise; extending the calculations is feasible, but computationally expensive. Addressing the stability of r = 2 km satellites with calculations that cover the age of the solar system is also time-consuming, requiring ∼5 real-time years per system.

To estimate the uncertainty from this limitation, we analyze a set of calculations where a massless tracer or a massive satellite was ejected at some time t_e ≳ 0.1 Myr. For each of these calculations, we predict the likelihood of ejection according to the orbital elements at some earlier time, t < t_e. In systems likely to experience an ejection, oscillations in a, e, and (sometimes) ${\imath }$ grow slowly with time. In nearly all cases where t_e ≳ 0.1 Myr, this growth begins at $t\ll {t}_{e}$. Thus, we can predict ejections with a fair degree of certainty. For the examples we considered, the accuracy of the predictions was ∼95%. In turn, the uncertainty in f_s is ≲5%.

Overall, uncertainties of ∼2% to 5% in the survival rates allow reasonably robust estimates of the likelihood of identifying new satellites in various regions of the Pluto–Charon system. Improving these estimates significantly requires order-of-magnitude larger investments in cpu time and is beyond the scope of the present effort.

As another cpu-saving measure, our symplectic integrator does not attempt to resolve collisions between massless tracers and massive satellites or among massive satellites. Thus, we do not attempt to predict impact rates on Pluto, Charon, or the smaller satellites. Previous studies suggest that large particles with initial semimajor axes beyond the orbit of Styx rarely impact Charon (Smullen & Kratter 2017). None impact Pluto. Although radiation pressure allows a larger fraction of small particles to impact Pluto or Charon (e.g., Poppe & Horányi 2011; Pires dos Santos et al. 2013; Porter & Grundy 2015), more than 75% of small grains are ejected into the solar system, where they become part of the Kuiper Belt (e.g., Smullen & Kratter 2017).

Direct n-body simulations of circumbinary satellites confirm previous conclusions for the semimajor axis of the innermost stable orbit (e.g., Doolin & Blundell 2011). Massless prograde (retrograde) satellites with small ${\imath }$ relative to the plane of the binary orbit are unstable when the initial semimajor axis a₀ ≲ 2.1 ${a}_{\mathrm{PC}}$ (a₀ ≲ 1.7 ${a}_{\mathrm{PC}}$). Some prograde survivors have orbital elements similar to those of the known small satellites in the Pluto–Charon system. Stable satellites on polar orbits must have a₀ ≳ 2.2 ${a}_{\mathrm{PC}}$.

For the adopted masses of Styx, Nix, Kerberos, and Hydra, possible orbits for other stable satellites on circumbinary orbits are tightly constrained. On timescales ranging from a few years to 10 Myr, nearly all massless tracers with prograde or polar orbits and $0.95\,{a}_{S}\lesssim a\lesssim 1.1\,{a}_{H}$ are ejected. Survivors on prograde orbits lie well inside the orbit of Styx or within the corotation zones of Nix and Hydra. Despite the lack of polar survivors inside the orbit of Nix, some tracers remain on orbits between Nix–Kerberos or Kerberos–Hydra. On the basis of the time evolution of e and the loss rate at 5–20 Myr, we suspect that nearly all will be ejected within 100 Myr. For prograde and polar orbits, many tracers outside the orbit of Hydra (a ≈ 3.6–4.0 ${a}_{\mathrm{PC}}$) are stable on ≳10 Myr timescales. The orbital evolution of these tracers suggests they will remain stable over gigayear timescales.

Experiments with massive satellites confirm these conclusions. Moons with radius r = 2 km and mass m = 4 × 10¹⁶ g on prograde orbits with initial semimajor axis a₀ ≈ 1.7–2.1 ${a}_{\mathrm{PC}}$ are unstable on timescales ranging from a few years to 100–300 Myr. Prograde satellites corotating in the orbits of Nix or Hydra are ejected on timescales of 30 yr to 300 Myr. Small moons outside the orbit of Hydra (a ≈ 75,000 km to 80,000 km) survive for 150 Myr and are likely on stable orbits.

Some small moons on polar orbits survive 300 Myr of dynamical evolution. Among those with a₀ between the orbits of Nix–Kerberos (Kerberos–Hydra), a few remain on stable orbits. However, nearly all are ejected. Satellites with polar orbits beyond Hydra are generally stable.

These results agree with expectations based on the Hill radius of each satellite. From previous calculations of multiplanet or multisatellite systems orbiting single or binary central objects, stability requires K ≳ 8–10 (e.g., Wisdom 1980; Petit & Henon 1986; Gladman 1993; Chambers et al. 1996; Deck et al. 2013; Fang & Margot 2013; Fabrycky et al. 2014; Kratter & Shannon 2014; Mahajan & Wu 2014; Pu & Wu 2015; Morrison & Kratter 2016; Obertas et al. 2017). With no stable tracers between the orbits of Styx and Hydra, our calculations support the more conservative K ≳ 10.

The conservative limit for satellite stability disagrees with results from Porter & Stern (2015), who predicted stable locations for satellites before the New Horizons flyby. After reanalyzing HST data and suggesting a somewhat smaller (larger) mass for Nix (Hydra), they describe a suite of numerical calculations that indicate a broad range of stable orbits for particles between the orbits of each pair of satellites. Formally, these results imply K ≲ 2–3. With typical durations of 1700 yr, however, their integrations are too short to infer robust lifetimes for extra satellites in the Pluto–Charon system (e.g., Youdin et al. 2012). Although Orchestra calculations confirm a large survivor fraction for tracers with a_S ≲ a₀ ≲ a_H at 1000–3000 yr (see the animations associated with Figures 8–9), longer-term calculations demonstrate that particles with 0.95 a_S ≲ a₀ ≲ 1.1 a_H are ejected on a broad range of timescales that are often much larger than a few thousand years (see also Youdin et al. 2012). On the basis of these longer integrations, we conclude that stable orbits with a_S ≲ a ≲ a_H for additional satellites are rare.

Aside from placing robust limits on satellite stability, the n-body calculations place strong constraints on the fate of debris in the Pluto–Charon system. In every simulation with massless tracers, all of the lost tracers are ejected from the system. None of the ejecta impact Pluto, Charon, or other massive satellites.

To explore options for detecting Pluto–Charon satellites with JWST, we rely on published descriptions of the instruments (e.g., Beichman et al. 2012; Doyon et al. 2012; Rieke et al. 2015) and sensitivity estimates from the JWST Pocket Guide.³ With fields of view ranging from 74'' × 113'' (MIRI; the Mid-Infrared Instrument) to 22 × 22 (NIRISS; Near-Infrared Imager and Slitless Spectrograph) to 2 × 22 × 22 (NIRCam; Near-Infrared Camera), all of the JWST instruments can probe satellites well outside Hydra's orbit. Although each instrument has tools to mask out part of the field, MIRI is not sensitive enough to detect Nix or Hydra at 10 μm (see below). NIRCam and NIRISS have complementary masking options and similar sensitivity for K-band imaging, with 10σ detections of 10 nJy sources expected in exposure times of 10⁴ s.

Among the known satellites, Nix and Hydra have V = 23–24 and equivalent spherical radii of roughly 20 km. Adopting optical/IR colors of the Sun for these high-albedo objects (V–K = 1.5; Weaver et al. 2016), K ≈ 22. For point sources having a flux of roughly 10³ nJy, a 100 s integration with a JWST imager yields a 10σ detection. With V = 26–27 and K ≈ 25, the smaller Styx and Kerberos require a factor of 10–20 longer integration time for a 10σ detection.

If Nix and Hydra have the mid-IR colors of the Sun, we expect 5 μm (10 μm) fluxes of 200 nJy (60 nJy). At their long-wavelength cutoffs of roughly 4.5 μm, NIRCam and NIRISS can achieve signal-to-noise ratio (S/N) = 10 detections in 200 s integrations. Using MIRI, we expect S/N = 10 (3) detections of Nix or Hydra with 10⁴ s integration times at 5 μm (10 μm). MIRI's coronagraph is unavailable at 5 μm; high-quality detections at 10 μm require much longer integration times than shorter wavelength observations with either NIRCam or NIRISS.

If there are any 1–2 km high-albedo satellites in the Pluto–Charon system, they are 100–200 times fainter than Nix/Hydra. Integrations that yield 10σ detections of Styx and Kerberos result in 3σ detections for these putative satellites. Longer integrations improve the likelihood of identifying any 1–2 km satellites in the system.

Overall, detecting the known satellites of Pluto–Charon with JWST is straightforward. A program that acquires the equivalent of 5–10 10³ s integrations yields excellent S/N for Nix/Hydra, very good S/N for Kerberos/Styx, and sufficient S/N to detect 1–2 km objects from the orbit of Styx to well outside the orbit of Hydra.

The Pluto–Charon system has a long and distinguished history of stellar occultation observations (e.g., Person et al. 2013; Boissel et al. 2014; Bosh et al. 2015; Gulbis et al. 2015; Pasachoff et al. 2016, 2017, and references therein). Aside from demonstrating the presence of Pluto's atmosphere (Elliot et al. 1989) and measuring physical conditions within the atmosphere (e.g., Hubbard et al. 1988; Elliot et al. 2003; Pasachoff et al. 2005), data from occultations provide detailed astrometry on the Pluto–Charon orbit, the ephemerides of the Pluto–Charon system, and information on dust and other small objects orbiting Pluto–Charon. Although Charon has been probed with occultations (e.g., Gulbis et al. 2015), the lone published attempt to monitor an occultation by Nix did not yield the expected drop in stellar flux (Pasachoff et al. 2016).

For the Pluto–Charon system with a semimajor axis a ≈ 40 au, occultations of small satellites with r ≈ 1–2 km are in the Fresnel limit (e.g., Bailey 1976; Cooray & Farmer 2003; Roques et al. 2006; Nihei et al. 2007; Bickerton et al. 2008; Schlichting et al. 2009, 2012; Wang et al. 2009, and references therein). On this scale, the duration of an occultation is

$$\begin{eqnarray}&&{\rm{\Delta }}{\text{}}t\,=\,2{\text{}}D{\phi }_{\star }\ (1+{A}^{1/2})/v\end{eqnarray} \tag{ 2 }$$

where D is the distance from the Earth to Pluto, ${\phi }_{\star }$ is the angular size of the occulted star, v is the relative velocity, and $A={(r/{R}_{\star })}^{2}$ is the amplitude of the occultation. In the amplitude, r is the radius of the satellite and ${R}_{\star }=D{\phi }_{\star }$ is the projected size of the star at Pluto.

Achieving a large amplitude requires that the projected size of the star is comparable to the radius of the satellite. With ${R}_{\star }$ = 1–10 km and D = 40 au, ${\phi }_{\star }\approx (2\mbox{--}20)\times {10}^{-10}$. If the physical radius of the star is similar to the solar radius, the distance to the star is d ≈ 10–100 pc.

Although the amplitude of the occultation can be large, the duration is very short. For a 1 km (10 km) satellite with ${R}_{\star }$ = 1 km, Pluto's orbital motion of v = 5 $\mathrm{km}\,{{\rm{s}}}^{-1}$ implies Δt = 0.8 s (8 s). Resolving the event requires integration times of 0.04 s (0.4 s). The relatively long duration of occultations of 10–20 km satellites makes the event much easier to resolve than occultations by 1–2 km satellites.

Relaxing our assumption of spherical satellites makes occultation observations somewhat more challenging. From New Horizons data, the aspect ratios range from roughly 1.5:1:1 for Nix to roughly 2:1:1 for Styx and Kerberos to roughly 2.6:1.8:1 for Hydra. For a worst-case scenario where the long axis of the satellite lies perpendicular to its path across the star, it is prudent to reduce the integration times by a factor of 2 to ≲0.02 s for Styx/Kerberos and by a factor of 2–3 to ≲0.2 s for Nix/Hydra.

Although ground-based occultations and JWST IR imaging are the only current options for detecting smaller satellites in the Pluto–Charon system, we briefly consider whether planned missions are also capable of discovering new satellites.

The Wide Field Infrared Survey Telescope (WFIRST) is a 2.4 m space telescope concept with a wide-field imager and a coronagraphic instrument (Green et al. 2012; Spergel et al. 2013). Despite having a sensitivity comparable to that of HST, the wide-field imager is not designed for the high-contrast high-resolution observations required to detect existing or additional satellites in the Pluto–Charon system. The coronagraphic instrument is expected to achieve contrast ratios of 10⁻⁹ or larger on scales of roughly 1'' (see also Burrows 2014; Lacy et al. 2018). Detecting very small satellites ($r$ ≪ 1 km) and any faint debris in the system would be straightforward. The likely field of view has a radius of only 2''; thus, it would not be possible to center the coronagraph on Pluto and search for faint satellites beyond the orbit of Hydra. Still, it is worth exploring whether moving Pluto–Charon off the edge of the field would allow deep searches for faint satellites.

The Origins Space Telescope is a 6–9 m telescope concept with four or five proposed instruments (e.g., Meixner et al. 2017; Battersby et al. 2018). Designed to operate at wavelengths λ ≳ 5 μm, current plans do not include an imaging instrument similar to those on JWST or WFIRST. Unless designs change, this facility is unlikely to enable discovery of faint satellites in the Pluto–Charon system.

The Large Ultraviolet/Optical/Infrared telescope (LUVIOR; formerly known as the Advanced Technology Large-Aperture Space Telescope or ATLAST) is a 8–16 m ultraviolet–optical–IR telescope concept with a variety of proposed instruments (e.g., Postman et al. 2009; Feinberg et al. 2014; Thronson et al. 2016; Arney et al. 2017; Bolcar et al. 2017). Scaling from quoted sensitivity estimates for instruments on the 6.5 m JWST, IR imagers on LUVOIR could detect small satellites beyond the orbit of Hydra with radii r ≈ 0.5–1 km. Improvements in instrument efficiency would enable detection of smaller satellites. If the proposed coronagraph could reach a contrast of 10⁻⁹ over a larger field of view, detection of meter-sized objects might be possible.

The Habitable Exoplanet Observatory (HabEx) is a 4–8 m telescope concept designed to image exoplanets and to detect biosignatures in the spectra of exoplanets around nearby stars (e.g., Mennesson et al. 2016). With an aperture comparable to that of JWST and wavelength coverage similar to that of LUVOIR, discovering satellites with r ≲ 1 km requires a coronagraph with new technology (e.g., Ruane et al. 2017, 2018). As with LUVOIR, reaching proposed contrasts of ∼10⁻⁹ enables detection of much smaller satellites.

The n-body calculations discussed here place new constraints on the dynamics of the Pluto–Charon satellite system. For the adopted masses of the four small satellites, the system is as tightly packed as possible. There are few stable orbits where satellites with negligible mass can exist between the orbit of Styx and the orbit of Hydra. This conclusion is independent of inclination angle: satellites on polar orbits are about as unstable as prograde satellites orbiting in the plane of the binary system.

There is also limited space for stable satellites on circumbinary orbits inside the orbit of Styx⁴  (see also Holman & Wiegert 1999; Pichardo et al. 2008; Doolin & Blundell 2011; Youdin et al. 2012). From simulations without any small satellites, stable orbits could exist at 1.3–2.1 ${a}_{\mathrm{PC}}$ (retrograde) and 1.7–2.1 ${a}_{\mathrm{PC}}$ (prograde). Polar orbits inside the orbit of Styx are unstable. The new calculations yield a set of massless (r ≳ 1 mm) survivors on prograde orbits with initial semimajor axes a₀ ≈ 1.8–2.1 ${a}_{\mathrm{PC}}$.

Well outside the orbit of Hydra, massless satellites on prograde and polar orbits are stable. With the nominal masses, Nix and Hydra clear out tracers with a ≲ 1.08 a_H on timescales of 5–10 Myr. On longer timescales, orbits with a ≈ 1.08–1.12 a_H are probably also unstable. Tracers with a ≳ 1.15 a_H appear to be stable.

Tests with massive satellites confirm these results. Although some massless tracers on prograde orbits survive for as long as 10 Myr near the corotation zones of Nix and Hydra, small satellites with r ≈ 2 km and m ≈ 4 × 10¹⁶ g on identical orbits have lifetimes as short as 10–1000 yr. Others survive for 300 Myr. Inside the orbit of Styx, small moons are ejected on timescales ranging from 1 yr to nearly 300 Myr. Small moons on polar orbits between Nix and Hydra are similarly unstable, with lifetimes ranging from 10³ yr to 300 Myr. Given the loss rate in our calculations, it seems unlikely that additional satellites with r ≳ 2 km inside the orbit of Hydra can survive for the age of the solar system.

Massive satellites outside the orbit of Hydra fare much better. After 150 Myr of evolution, roughly 90% of those on prograde orbits survive. Survivors tend to have larger initial semimajor axes than ejected moons. Small moons on polar orbits outside the orbit of Hydra are just as unlikely to be ejected after 300 Myr of evolution. Thus, searches for new satellites should concentrate on regions outside the orbit of Hydra.

Aside from clarifying our understanding of circumbinary dynamics, these results help to interpret the lack of new satellite detections from New Horizons imaging data. With additional satellites generally precluded at a ≲ 1.1 a_H, the volume for discovering new, stable satellites within the overall New Horizons footprint a ≲ 1.6 a_H is restricted. Although New Horizons should have failed to detect new satellites inside the orbit of Hydra, the lack of 2–3 km satellites outside the orbit of Hydra is surprising. In Kenyon & Bromley (2014b), nearly every calculation left several 1–3 km satellites outside the orbit of Hydra. Although it seems unlikely that these satellites might have an albedo significantly smaller than the 55% of Kerberos, overestimates in their predicted sizes by a factor of two are plausible. Gravitational interactions with the known satellites may have led to complete ejection from the Pluto–Charon system instead of simple scattering beyond the orbit of Hydra.

If satellites scattered beyond the orbit of Hydra are simply smaller than the New Horizons threshold of 2 km, they await detection with a new generation of instruments. With typical Hill radii r_H ∼ 100 km at semimajor axes a ≈ 1.5 a_H, several satellites with r ∼ 1 km could orbit stably outside the orbit of Hydra for many gigayears. Thus, it is still possible to discover new moons of Pluto and test numerical models for satellite formation.

The calculations of massless tracers orbiting between Styx and Hydra also provide new insights into the lack of emission from small particles detected from New Horizons data (Lauer et al. 2018). Prior to the New Horizons flyby, several studies derived upper limits on small moons and dust emission from direct imaging (e.g., Steffl et al. 2006; Marton et al. 2015) and occultations (Boissel et al. 2014; Throop et al. 2015). Theoretical studies based on n-body simulations predicted steady-state optical depths from a balance between dust production from impacts on Pluto–Charon and the smaller satellites and losses from radiation pressure and dynamical ejections (Stern et al. 2006; Pires Dos Santos et al. 2011; Poppe & Horányi 2011; Pires dos Santos et al. 2013; Porter & Grundy 2015). The new results described here demonstrate that all orbits with a ≈ 1.7–3.7 ${a}_{\mathrm{PC}}$ are dynamically unstable on timescales ranging from several decades to 10 Myr. For particles with r ≲ 100 μm, radiation pressure shortens these lifetimes considerably (Burns et al. 1979; Hamilton & Burns 1992; Poppe & Horányi 2011; Pires dos Santos et al. 2013; Porter & Grundy 2015). Without additional dust production from impacts, our results imply that there should be no dust or larger particles inside the orbit of Hydra.

If small particles or dust emission are ever detected in the Pluto–Charon system, the mass and location of small particles place interesting constraints on the masses of the four small satellites. If limits on the production rate for small particles from impacts can be established from New Horizons data, the number of survivors between the orbits of Styx and Hydra is sensitive to the mass of Nix and Hydra: smaller masses for these satellites allow larger masses in dust and larger particles (see also Stern et al. 2006; Pires Dos Santos et al. 2011; Poppe & Horányi 2011; Pires dos Santos et al. 2013; Porter & Grundy 2015).

Identifying new satellites or small particles outside the orbit of Hydra would also improve estimates for the mass of Hydra. From the n-body calculations, the innermost stable orbit outside Hydra is much more sensitive to Hydra's mass than to the mass of Nix or the other satellites (see also Michaely et al. 2017). Accurate measurements of the orbital elements for any new satellite would thus provide new limits on Hydra's mass.

Although the current launch date for JWST is not until 2021 March, the observations proposed here are similar to the suite of HST imaging data collected prior to the New Horizons flyby (e.g., Buie et al. 2006; Steffl et al. 2006; Weaver et al. 2006; Steffl & Stern 2007; Tholen et al. 2008; Showalter et al. 2011, 2012; Brozović et al. 2015; Showalter & Hamilton 2015). Starting in Cycle 2, various HST imaging programs sought to constrain the properties of dusty structures and satellites in the Pluto–Charon system. These programs acquired many images per HST orbit, with exposure times ranging from a few seconds to several minutes. As outlined in Brozović et al. (2015), multiple images per orbit enable a robust analysis procedure that eliminates light from background stars and provides the best possible S/N for images of very faint satellites.

We anticipate that a JWST observing program to detect faint satellites in the Pluto–Charon system would be similar to a typical Hubble program. Within a single visit, the likely total exposure time with JWST imagers would probably be several times longer than an HST orbit. As outlined in Gordon et al. (2015) for MIRI, multiple short and moderate exposures would allow the same type of analysis procedure as performed by Brozović et al. (2015).

Aside from the ability of JWST instruments to perform to specifications, the main uncertainty in any JWST program is the overhead involved in acquiring a target, maintaining a fix on the target, conducting the observations, and performing the housekeeping needed for the health of the satellite. Various performance analyses suggest the typical overhead in an observing program is roughly 30% (Gordon et al. 2012a, 2012b). Thus, the program outlined here seems feasible.

Over the past few decades, various groups have observed Pluto and Charon occult fairly bright stars to infer the extent and physical properties of their atmospheres and to plan for the New Horizons flyby (e.g., Person et al. 2006, 2013; Sicardy et al. 2011; Bosh et al. 2015; Dias-Oliveira et al. 2015; Gulbis et al. 2015; Throop et al. 2015; Pasachoff et al. 2016, 2017; Sicardy et al. 2016, and references therein). With exposure times of 0.25 to 5 s on 1 to 2.5 m telescopes, the typical S/N ranges from roughly 10 to better than 100.

Detailed astrometric analyses demonstrate that occultations by the Pluto–Charon system are fairly frequent. For 2008–2015, predictions by Assafin et al. (2010) indicate 30–300 possible occultations per year for each of Pluto, Charon, Nix, and Hydra. Although we are unaware of similar published predictions for Styx and Kerberos, it seems likely from Assafin et al. (2010) that the frequency of occultations for both of these small moons is similar to that for Nix and Hydra. The uncertain orbital path of Pluto across the sky is the main limitation in these predictions (e.g., Benedetti-Rossi et al. 2014; Holman & Payne 2016, and references therein).

Many studies have explored the possibility of using occultations to detect small Kuiper Belt objects beyond 40–50 au (KBOs; e.g., Bailey 1976; Brown & Webster 1997; Alcock et al. 2003; Cooray & Farmer 2003; Cooray 2003; Chang et al. 2006, 2011; Roques et al. 2006; Bickerton et al. 2008; Bianco et al. 2009, 2010; Schlichting et al. 2009; Wang et al. 2009, 2010; Schlichting et al. 2012; Zhang et al. 2013, and references therein). Aside from ground-based optical observations, HST and Rossi X-Ray Timing Explorer have yielded promising sets of data to search for serendipitous occultations of stellar sources by KBOs. On the ground and with HST, exposure times of 0.02 s yield light curves with sufficient S/N to detect KBOs with radii of roughly 1 km.

On the basis of this discussion, it seems straightforward to detect occultations of stars by Pluto's known small satellites with modest aperture ground-based telescopes. However, the likelihood of detecting a serendipitous occultation event from an unknown satellite is small. Assuming a satellite with a diameter D₆ ≈ 2 km is positioned randomly in an orbit with semimajor axis a₆ ≈ 75,000 km (outside Hydra), the chance of having it along the same occultation path as any of the other small satellites is less than 2D₆/2πa₆ ≲ 10⁻⁵. With N of these satellites, the probability is still small, 10⁻⁵ N, unless N is very large.

Overall, the best strategy to detect and to characterize new satellites involves initial imaging observations with JWST followed by ground-based occultation observations. Although JWST data are fairly expensive and limited by spacecraft constraints, observations with any of the imaging instruments sample broad swaths of available discovery space. If new satellites are detected, occultations can then provide additional information on the orbits and shapes.

Over the past few decades, various techniques have revealed ∼15 circumbinary planetary systems (Thorsett et al. 1993; Doyle et al. 2011; Orosz et al. 2012a, 2012b; Welsh et al. 2012; Kostov et al. 2013, 2014, 2016; Schwamb et al. 2013; Bailey et al. 2014; Hinse et al. 2015; Getley et al. 2017; Jain et al. 2017). Aside from systems with main-sequence stars, the binaries include a binary pulsar and a low-mass X-ray binary. Planet masses range from approximately Neptune up to approximately seven times Jupiter. Often, the circumbinary planets orbit in the plane of the inner binary and are reasonably close to the innermost stable orbit. In several, the orbit of the planet is somewhat tilted with respect to the inner binary. Sometimes, the planet is well outside the innermost orbit.

The formation and evolution of circumbinary planets have attracted intense theoretical interest (e.g., Pierens & Nelson 2007, 2008a, 2008b; Paardekooper et al. 2012; Pierens & Nelson 2013; Rafikov 2013; Bromley & Kenyon 2015a; Kley & Haghighipour 2015; Rafikov & Silsbee 2015a, 2015b; Silsbee & Rafikov 2015a, 2015b; Hamers et al. 2016; Li et al. 2016; Vartanyan et al. 2016; Mutter et al. 2017a, 2017b; Fleming et al. 2018; Quarles et al. 2018; Pierens & Nelson 2018; Thun & Kley 2018; Zanazzi & Lai 2018, and references therein). In addition to changing the structure and evolution of a planet-forming gaseous disk, the central binary provides a challenging environment for the growth of Earth-mass and larger planets from kilometer-sized and larger planetesimals. Once planets form, the central binary efficiently removes them from resonant orbits. Multiplanet systems are particularly prone to disruption.

Although Kepler 47 is the only circumbinary planetary system with more than one planet (e.g., Orosz et al. 2012a; Kostov et al. 2013; Hinse et al. 2015), there are numerous multiplanet systems orbiting single stars (e.g., Lissauer et al. 2011, 2012, 2014; Fabrycky et al. 2014; Rowe et al. 2014; Winn & Fabrycky 2015; Ballard & Johnson 2016; Sinukoff et al. 2016; Udry et al. 2017; Weiss et al. 2018, and references therein). Many of these are closely packed, with little or no space for additional planets between the innermost and outermost planets. In each multiplanet system, the planets are often of similar size, with the outermost planet usually the largest (see also Ciardi et al. 2013; Millholland et al. 2017). Sometimes, the planets have very different sizes.

Theoretical studies of closely packed multiplanet systems focus on their origin and stability (e.g., Hansen & Murray 2012, 2013; Rein 2012; Kratter & Shannon 2014; Najita & Kenyon 2014; Raymond & Cossou 2014; Schlaufman 2014; Batygin & Laughlin 2015; Malhotra 2015; Pu & Wu 2015; Steffen & Hwang 2015; Morrison & Kratter 2016; Mustill et al. 2017; Pan & Schlichting 2017, and references therein). Formation at several astronomical unit distances followed by migration through a circumstellar gaseous disk is a favored explanation for many close-in multiplanet systems. However, in situ growth is also a possibility. In both mechanisms, it is unclear how the most closely packed planetary systems form and maintain their stability on gigayear timescales.

Together with previous theoretical studies of the Pluto–Charon system (Canup 2005; Lee & Peale 2006; Ward & Canup 2006; Lithwick & Wu 2008a, 2008b; Canup 2011; Youdin et al. 2012; Cheng et al. 2014; Kenyon & Bromley 2014b; Bromley & Kenyon 2015b; Desch 2015; Pires et al. 2015; Walsh & Levison 2015; Michaely et al. 2017; McKinnon et al. 2017; Smullen & Kratter 2017; Woo & Lee 2018), the results described here inform our understanding of circumbinary planets and closely packed planetary systems orbiting single or binary stars. Despite forming in a relatively gas-free environment, the Pluto–Charon system has issues with in situ formation, migration, orbital resonances, and long-term stability similar to those in exoplanetary systems. Nevertheless, the small satellites probably (i) grew from a ring of debris, (ii) found stable orbits close to resonances with the central binary, and (iii) maintained these orbits for ∼4 Gyr (Weaver et al. 2016; Robbins et al. 2017). Working out the details of this history for the Pluto–Charon satellites and exoplanetary systems will enrich theories of planet formation.