CIRCUMBINARY CHAOS: USING PLUTO'S NEWEST MOON TO CONSTRAIN THE MASSES OF NIX AND HYDRA

Nix, Hydra, and P4 are detected in reflected visible light. Their diameters depend on their unknown albedos, A, as

$$\begin{equation} D_{\rm Nix} \approx 25 A^{-1/2}\,{\rm km} \end{equation} \tag{ 1a }$$

$$\begin{equation} D_{\rm Hyd} \approx 30 A^{-1/2}\,{\rm km} \end{equation} \tag{ 1b }$$

$$\begin{equation} D_{\rm P4} \approx 8 A^{-1/2}\,{\rm km}\,. \end{equation} \tag{ 1c }$$

To relate apparent magnitudes to size, we adopt a Charon-like phase coefficient (Buie et al. 1997; T08).

Mass estimates from photometry rely on an assumed density. With ρ₁ = ρ/(1 g cm⁻³), the mass ratios relative to Pluto are

$$\begin{equation} \mu _{\rm Nix} \approx 6.4 \times 10^{-7} \rho _1A^{-3/2} \end{equation} \tag{ 2a }$$

$$\begin{equation} \mu _{\rm Hyd} \approx 1.1 \times 10^{-6} \rho _1A^{-3/2} \end{equation} \tag{ 2b }$$

$$\begin{equation} \mu _{\rm P4} \approx 2.0 \times 10^{-8}\rho _1A^{-3/2} \,, \end{equation} \tag{ 2c }$$

for Nix, Hydra, and P4. To express (without loss of generality) masses in terms of albedos, and not the more cumbersome Aρ^−2/3₁, we assume ρ₁ = 1, unless stated otherwise. Neglecting the possibility of high porosity, we also refer to A = ρ₁ = 1 as the "minimum mass" case.

The orbit solution of T08 gives μ_Nix = (4.4 ± 3.9) × 10⁻⁵ and μ_Hyd = (2.5 ± 4.9) × 10⁻⁵. The large uncertainties reflect the difficulty of measuring small mass ratios astrometrically. Since the uncertainties extend toward or beyond zero,¹ the T08 solution places upper limits ∼10⁻⁴ on the mass ratios of Nix and Hydra. The equivalent albedo constraint is A ≳ 0.04. Our stability calculations place much tighter constraints of A ≳ 0.3 on Nix and Hydra.

The orbit of P4 is only loosely constrained. The discovery (S11) announces P4 as "consistent with" a circular, coplanar orbit. Its mean motion from HST images spanning ∼20 days implies an orbital period of 32.1 ± 0.3 days, corresponding to a period ratio of 5.03 ± 0.05 with Charon.

For a Keplerian orbit about the PC barycenter (a good approximation at the 1% level), the measured period corresponds to a mean separation of 57, 400 ± 400 km. Figure 1 plots this orbit location as a dark green dot with error bars. The lighter green error bars show the larger range of orbits considered in our numerical simulations.

Zoom In Zoom Out Reset image size

The mean motion is consistent with the measured projected radial separation of 59, 000 ± 2000 km from Pluto (which is 2040 km from the barycenter). The pre-discovery HST images that date back to 2006 (S11), as well as future observations from HST, ALMA, and New Horizons will help constrain the orbit of P4.

As summarized in Table 1, we adopt the orbit solution of T08 for the orbits of Pluto, Charon, Nix, and Hydra and the masses of Pluto and Charon. Ephemerides from JPL's HORIZONS system² are similar, but not identical to the T08 solution. Orbit integrations started with the JPL (specifically the PLU021) ephemerides gave statistically similar results to the T08 ephemerides. Thus, we only present results based the T08 orbit solution.

Table 1. Keplerian Orbital Elements^a and Cartesian State Vectors^b from Tholen et al. (2008)

	Charon	Nix	Hydra		Charon	Nix	Hydra
Mass (kg)c	1.52e21	(5.8e17)	(3.2e17)	μc	0.1166	(4.4e−5)	(2.5e−5)
a (km)	19570.3	49240	65210	x	−0.1677	−0.1744	0.5202
e	0.0035	0.0119	0.0078	y	−0.2551	0.7148	−0.8822
i	96.168	96.190	96.362	z	0.0	5.6384e−5	3.4397e−3
ω	157.9	244.3	45.4	$\dot{x}$	1.5995	−1.0282	1.0625
M	237.0	129.80	129.12	$\dot{y}$	−1.0452	−0.3580	0.4458
Ω	223.054	223.202	223.077	$\dot{z}$	0.0	−3.1683e−3	−1.8596e−4
P (days)	6.3872	25.49	38.85	$P \over P_{\rm PC}$	1	3.99	6.06

Notes. ^aPlutocentric for Charon and Barycentric for Nix and Hydra. The Nix and Hydra masses are the T08 values, not those used in this work. ^bCartesian state vectors are given in code units where G = M_♇ = P_pc = 1. The Pluto–Charon orbit has been rotated into the x–y plane. ^cM_♇ = 1.304e22 kg, μ is mass relative to Pluto.

Download table as:  ASCIITypeset image

Though the T08 (and also JPL) solutions assume specific masses for Nix and Hydra, we vary the masses of Nix and Hydra without varying the orbit solution. The error introduced is hopefully modest compared to the astrometric uncertainty. Future work could establish how orbit solutions vary with uncertain masses.

Salient features of the system's orbital parameters, including circularity, coplanarity, and the proximity of moons to resonances, are discussed in the introduction and evident in Table 1. Though low, the finite eccentricity of the PC orbit, e_C = 3.5 × 10⁻³, has long been considered suspicious. Tidal interaction between Pluto and Charon damps e_C on short, ∼10 Myr, timescales (Dobrovolskis et al. 1997; Lithwick & Wu 2008a). Nix and Hydra cannot excite e_C above ∼10⁻⁵ (LP06); P4 cannot contribute significantly either. External forcing—from solar and planetary tides, Kuiper Belt Objects (KBO) flybys and collisions—appears unable to explain the observed eccentricity (Stern et al. 2003).

After the submission of this manuscript, a lower eccentricity for Charon was announced (Buie et al. 2012). This new solution gives a 1σ limit of "3 km out of round" (Buie et al. 2012), implying e_C ≲ 1.5 × 10⁻⁴ (at 1σ). A preliminary investigation shows that setting e_C = 0 in the T08 orbit solution has a minor effect on orbital stability (Section 4.4).

The stability of a test particle orbiting an inner binary is well studied (Szebehely 1980; Dvorak 1986), including the systematic numerical investigation of Holman & Wiegert (1999, hereafter HW). The HW results show that Nix, Hydra, and P4 are all individually stable about the PC binary. Based on the results of HW (see their Equation (3) and/or their Table 7), period ratios below 2.8 are unstable in the PC system. Nix's period ratio of ≈4 is 41% larger than the stability boundary, and also avoids a possible instability strip near the 3:1 resonance.

We now temporarily ignore the important fact that Pluto and Charon are a binary and consider the three-body stability of two satellites about a central mass that places the combined mass of Pluto and Charon at their barycenter. This approximate problem allows us to simply demonstrate how strong interactions with Nix or Hydra restrict the allowed orbits of P4, as in Figure 1, and also aids the interpretation of numerical results.

The well-known stability criterion for initially circular orbits with semimajor axes a₁ and a₂ = a₁ + Δ requires Δ > 3.5R_H in the restricted three-body problem when one of the satellites is massless (Gladman 1993; which also analyzes two massive satellites). That is, two satellites on circular orbits are stable if their orbits are separated by more than 3.5 Hill radii,

$$\begin{equation} R_{\rm H} = \left(\mu \over 3\right)^{1/3} a , \end{equation} \tag{ 3 }$$

of the massive satellite.

The separation of P4 from Nix or Hydra in terms of Hill radii depends on the assumed masses for Nix and Hydra. For the mass scalings in Equation (2), P4 is separated from Nix by $29 \sqrt{A_{\rm Nix}}$ and from Hydra by $17 \sqrt{A_{\rm Hyd}}$ Hill radii (of Nix and Hydra, respectively). Separation by fewer than 3.5R_H would require either A_Nix < 0.014 or A_Hyd < 0.042. Such low albedos and thus high masses are inconsistent with the T08 (1σ) upper mass limits. The more stringent mass limits that we find ensure that P4 is stable by the Hill stability criterion for circular orbits.

The stability of an eccentric test mass follows from the (approximately) conserved Tisserand parameter

$$\begin{equation} C_{\rm T} = {a^{\prime } \over 2a} + \sqrt{{a \over a^{\prime }}(1-e^2)} \cos i \,, \end{equation} \tag{ 4 }$$

where a, e,  and i (inclination) refer to the test body (P4) and a' refers to the perturber (Nix or Hydra). The stability boundary of circular, coplanar (e = i = 0) orbits at |a − a'| = 3.5R_H defines a critical C_T. (Restricted three-body) stability for arbitrary e and i holds for C_T above the critical value.

Figure 1 plots the critical C_T curves, sometimes called "Tisserand tails." The outer and inner tails of Nix and Hydra, respectively, are relevant for interactions with P4. For low Nix and Hydra masses (the A = 1 curves), the tails would only cross an eccentric P4, e_P4 ≳ 0.1. With higher Nix and Hydra masses (the A = 0.05 curves), Hydra's Tisserand tail intersects plausible P4 orbits at lower e_C, especially if P4 lies on the outer edge of allowed orbits. While a considerable simplification, these considerations from the restricted three-body problem demonstrate how higher masses for Nix and Hydra severely restrict the allowed orbits of P4.

Still ignoring perturbations from the central binary, we now consider three interacting satellites about a central mass. Sharp stability boundaries no longer exist, but the timescale to orbit crossing can be studied numerically.

Most investigations of multi-planet stability consider roughly equal mass companions with evenly spaced orbits, in terms of Hill radii. Chambers et al. (1996, hereafter C96) defined the mutual Hill radius, as

$$\begin{equation} R_{\rm H}^{\prime } = \left(\mu _i + \mu _{i+1} \over 3\right)^{1/3} {a_i + a_{i+1} \over 2} \, \end{equation} \tag{ 5 }$$

for neighboring planets (i and i + 1). This definition does not reduce to the standard Hill radius when μ_i → 0. This deficiency is readily corrected by a mass-weighted averaging of the semimajor axes. This distinction is insignificant when the satellite masses are similar, but is relevant here due to the low mass of P4.

C96 measured the orbit crossing timescale, t_c, for systems of three equal mass planets, with planet–star mass ratios ranging from μ = 10⁻⁹ to μ = 10⁻⁵. Relative to the period of the inner planet, P₁, we fit the data in Figure 4 of C96 as

$$\begin{equation} \log _{10}(t_{\rm c}/P_1) = -9.11 + 4.39 \Delta ^{\prime } \mu ^{1/12} - 1.07 \log (\mu)\,, \end{equation} \tag{ 6 }$$

where Δ' is the orbit separation in mutual Hill radii and the functional form is motivated by the Faber & Quillen (2007) study of systems of >3 planets. The mass scaling Δ'μ^1/12∝μ^−1/4 differs from the μ^−1/3 Hill radius scaling due to three-body resonances (Quillen 2011).

Unequal masses (of P4 in particular) make direct application of these results difficult. To allow simple estimates, we make a range of possible assumptions: μ as the average of all three satellites or of just Nix and Hydra; the mutual Hill radius defined as in Equation (5) or with density weighted a_i. In all cases, the minimum mass A = 1 case is stable, with crossing times >10³⁰ yr. For the higher mass A = 0.05 case, t_c ≲ 3 × 10⁶ yr, implying rapid instability.

Even neglecting binary perturbations, the stability of Pluto's minor moons depends strongly on their assumed masses. Binary perturbations accelerate the destabilization. We show in Section 4 that for the high-mass case with A = 0.05, the simulated crossing time for P4 is ≲ 10³ yr, over a thousand times faster than the above estimate neglecting the central binary.

In this paper, we approximate the Pluto system as a set of isolated point masses. Since Pluto's Hill sphere extends ∼120 times further than Hydra's orbit, solar tides are a minor effect. The protective 3:2 resonance between Neptune and Pluto helps to weaken the dominant planetary perturbation. Nevertheless, follow-up work should test the effect of weak external perturbations on long-term stability.

Collisions with interloping KBOs could significantly perturb the weakly bound outer moons. If the Kuiper Belt was massive in its youth, collisions would unbind ∼100 km class binary KBOs (Nesvorný et al. 2011; Parker & Kavelaars 2012). Pluto's moons, P4 in particular, could be destabilized by collisions that only modestly perturb its eccentricity. Combining the dynamical stability of Pluto's moons with collisional perturbations could be a powerful constraint on both the Pluto system and the collisional environment of the Kuiper Belt.

Tides in the Pluto system are dominated by interactions between Pluto and Charon, which are locked in a dual synchronous state, with both spin periods equaling the orbital period (Dobrovolskis et al. 1997). In principle, the eccentricities of minor moons should be damped by exciting the eccentricity of Charon, which then suffers tidal dissipation. However, Lithwick & Wu (2008b) conclude that this damping mechanism is weak.

We study the stability of P4 with a suite of 4 + N-body integrations. In these simulations, the four massive bodies are Pluto, Charon, Nix, and Hydra. Treating P4 as a massless test particle allows simultaneous investigation of many trial orbits. We follow each test particle until it crosses the orbit of either Nix or Hydra. (No significant perturbations to the orbits of the massive bodies occurs in the simulations.) Integrations were performed using the 15th-order Radau integrator (Everhart 1985), as implemented by the Swifter software package.³ Details concerning code performance are deferred to Appendix B.

Between different sets of simulations, we vary the uncertain masses of Nix and Hydra, keeping their mass ratio fixed at M_Nix/M_Hyd = 0.575, the value implied by their median optical brightness. This mass ratio assumes that Nix and Hydra share the same density and albedo, which should be expected for standard formation scenarios.⁴ The range of Nix and Hydra masses considered corresponds (for a density of 1 g cm⁻³) to albedos from 0.03 to 0.24. Integrations were run with albedos up to 0.4 (μ_Nix  ≈  10⁻⁶), however only 10%–30% of test particles suffered orbit crossing after ∼10⁸ Charon orbits, making systematic investigations too expensive.

The initial conditions for the massive bodies are given in Table 1. The initial orbits for the test particles were chosen by two different methods. The first suite of simulations, described in Section 4.2, populated P4 orbits by randomly sampling Keplerian elements. The second suite of simulations, summarized in Section 4.3 initializes P4 on the "most circular" orbits about the PC binary.

For each adopted mass of Nix and Hydra, we integrate 5000 P4 orbits with initial conditions chosen randomly from a uniform distribution of Keplerian osculating elements.⁵ The semimajor axes are restricted to the range

$$\begin{equation} 56{,}632\,{\rm km} < a < 58{,}832\,{\rm km} , \end{equation} \tag{ 7 }$$

or a/a_PC = 2.89–3.01. Eccentricity and inclination are restricted to e < 0.05 and i < 05. The modest range in i was insignificant for our results and will not be discussed further. Other Keplerian angles (argument of pericenter, longitude of ascending node, and mean longitude) are sampled over the full (2π) range. Figure 1 compares the sampled orbits to the nominal location of mean motion resonances as well as Nix and Hydra crossing trajectories.

4.2.1. Mapping a–e Space

Figure 2 maps the median stability timescale versus initial a and e for several Nix and Hydra masses. As shown by the color bars on each map, crossing times increase significantly as the mass of Nix and Hydra drops from high to low (upper left to lower right plots). Crossing times are always short at higher a and e (upper right of each plot). Interactions with Hydra are the likely culprit, consistent with the Tisserand parameter curves in Figure 1.

Zoom In Zoom Out Reset image size

At low eccentricity, crossing times are longest, and display a complex dependence on a. Resonances play a key role. Indeed resonance overlap is a generic cause of orbital chaos (Mudryk & Wu 2006). From the approximate resonance locations in Figure 1, specific resonances can be implicated.

The 5:1 resonance with Charon helps destabilize the region near a ∼ 5.8 × 10⁴ km (Δa ∼ 8.0[ × 10³ km] as labeled in Figure 2) at lower masses. The 4:5 resonance with Hydra shortens crossing times between 5.65 ≲ a/(10⁴ km) ≲ 5.70 (i.e., 6.5 ≲ Δa/(10³ km) ≲ 7.0), for the highest mass (upper left). This change in the prominence of different resonances is expected as the mass of Nix and Hydra varies relative to PC.

4.2.2. Median Crossing Times: Measured and Extrapolated

Figure 3 shows how crossing times scale with the mass of Nix and Hydra. Median timescales are plotted for three subsets of the initial orbital parameters: (1) all initial orbits, (2) only e < 0.015, and (3) the stable "core" of parameters with both e < 0.015 and 5.70 × 10⁵ km < a < 5.75 × 10⁵ km. For the "core" sample, we also plot the 90th percentile of crossing time (beyond which only 10% of particles survive). At any mass, crossing times increase as cuts become more selective. The "missing" data points for low-mass cases occur where P4 orbits were so stable that the median (or 90th percentile) timescale was not reached after 10⁸ PC periods.

Zoom In Zoom Out Reset image size

The dependence of crossing time on the mass, m, of Nix and Hydra is well described by a power law

$$\begin{equation} t_{\rm c} \propto m^{-\gamma } . \end{equation} \tag{ 8 }$$

The best-fit indices are γ ≈ −3.6, −4.3, −4.4,  and  − 4.6 for the bottom to top curves in Figure 3. The larger scatter about the 90th percentile "core" power law is partly due to Poisson noise, as fewer orbits contribute to this restrictive sample. Physical scatter caused by the shifting locations of the most stable regions could also contribute.

No known theoretical explanation exists for such a power-law scaling. Indeed, if crossing times scale with Hill radius as in Equation (6), the mass dependence is not a simple power law—instead the local power-law index would steepen toward lower masses. However, Duncan & Lissauer (1997) found that orbit crossing timescales for the Uranian moons also exhibit a power-law dependence on satellite mass.

Extrapolation along these power laws allows us to speculate about the stability of P4 for very low masses of Nix and Hydra. While such extrapolation is inherently uncertain, low masses are much more costly to simulate directly due to longer crossing timescales. The goal of extrapolation is to estimate the "allowed" masses of Nix and Hydra, where allowed means that P4 crossing timescales exceed the age of the solar system.

Including all sampled P4 orbits (the blue curve in Figure 3), the allowed masses of Nix and Hydra are implausibly low—below the minimum set by A = 1. However for low eccentricity P4 orbits (all other curves), a range of plausible Nix and Hydra masses are allowed. From extrapolation of the 90th percentile "core" sample (red curve), the lowest allowed albedo is A ∼ 0.5. Lower albedos are possible for special orbits, including the "most circular" orbits that we present next.

The greater stability of low eccentricity P4 orbits shown above motivates the use of the "most circular" orbits about PC as an initial condition. These special orbits are not simply found by setting the Keplerian e = 0. Appendix A describes techniques to find these orbits.

We integrated trial orbits for 48 distinct P4 periods: 41 are uniformly spaced with period ratios (relative to Charon) from 4.79 to 5.31 and the remaining seven provide more dense sampling near the 5:1 resonance. At each orbital period 25 trial orbits were integrated. These 25 orbits are shifted in orbital phase along the most circular orbit by one Charon period each, as described in Appendix A.2.

Figure 4 shows the crossing times of the most circular orbits, for all orbital periods and for different Nix and Hydra masses. The median crossing times at each period are given by square symbols, with error bars spanning the 10th and 90th percentile. Lower limits are given where simulations did not run long enough to determine a timescale.

Zoom In Zoom Out Reset image size

Crossing times become longer as mass decreases, and the dependence on orbital period (or a) is quite complex. This general behavior agrees with the Keplerian initial conditions (Section 4.2). Figure 4 shows further that the period dependence is more varied and irregular for lower masses of Nix and Hydra. While it is difficult to understand all these fluctuations in detail, the trend is intuitively consistent with weaker resonant forcing being more narrowly focused at specific orbital periods.

Figure 4 shows a dramatic drop in crossing times near Charon's 5:1 resonance (C5:1). As also seen in Figure 2, this instability strip appears for lower masses of Nix and Hydra. At higher masses, Nix and Hydra perturbations wash out the influence of C5:1. Longer-lived P4 orbits exist both outside and inside the narrow instability strip at C5:1 Within the observational constraints, shown by the horizontal error bar, our analysis favors locations outside C5:1.

Figure 5 plots crossing times versus the mass of Nix and Hydra, with power-law fits overplotted. We restrict the range of periods to the observational constraint (S11), but with double the uncertainty for inclusiveness. The statistical measures of t_c include: (1) "median," which takes the median of the values given by the symbols in Figure 4 (themselves median values over phases); (2) "max median," the longest phase-median t_c, i.e., highest symbol; (3) "max 90%," the longest 90th percentile t_c, i.e., highest upper error bar; and finally (4) "longest," simply the longest t_c at any phase or period considered.

Zoom In Zoom Out Reset image size

For reference, the median t_c for the relatively stable "core" sample of initial Keplerian parameters is overplotted. Compared to this reference case, the crossing timescales for the most circular orbits are significantly longer, especially toward the (more realistic) lower masses of Nix and Hydra. Longer crossing times equate to higher allowed masses (and lower albedos) for Nix and Hydra.

Extrapolating along the power-law fits in Figure 5 shows that A ≳ 0.3 is needed to achieve t_c > 4 Gyr. This limit is more inclusive than the A ≳ 0.5 found in Section 4.2.2. We cannot definitely rule out even lower albedos. As already discussed, extrapolation could prove misleading.

Characterizing the most stable orbit is difficult. We do not base our estimate on the absolute longest lived orbits, which loosen our constraints. As shown in Figure 5, the longest t_c's do not follow a simple power law and are therefore unreliable for extrapolation. Even without extrapolation, the longest t_c's are highly subject to sampling, especially considering the pronounced period and phase dependence shown in Figure 4. It is also unclear whether P4 is likely to inhabit the most stable orbits, especially if those orbits occupy a tiny volume of phase space. Small neglected effects (such as collisions, see Section 3.4) could easily remove P4 from narrow pockets of parameter space. Thus, we base our constraints not on the absolutely most stable orbit, but on an average of orbits among the most stable.

We include a preliminary investigation of the effect of Nix and Hydra's mass ratio on the orbital stability of P4. Our main set of simulations fixed M_Nix/M_Hyd = 0.575. Here, we also investigate M_Nix/M_Hyd = 0.3 and 1.2, i.e., roughly half and double the fiducial value. For each mass ratio, we consider four values of the fixed total mass equal to the four highest mass cases in our main set of simulations. As in Section 4.3, we initialize the P4 orbits on the most circular orbits, with the same sampling of orbital periods and phases.

Figure 6 plots P4 lifetimes for these runs with different mass ratios. The dependence of orbital stability on the Nix/Hydra mass ratio is relatively minor compared to the stronger dependencies on Nix and Hydra's total mass and P4's orbital period. The decrease in orbital stability near the 5:1 commensurability—probably the most relevant feature—is evident in the lower mass runs for all mass ratios investigated, as are several other trends with period. The most prominent effect of the mass ratio is an increase in P4 crossing times at longer periods (≳ 5.1P_Charon) for larger mass ratios, i.e., smaller Hydra masses. This behavior is unsurprising, and is qualitatively explained by three-body interactions. As shown in Figure 1, a more massive Hydra strongly perturbs orbits exterior to P4's measured orbit.

Zoom In Zoom Out Reset image size

For the two lowest total mass runs, the variation in median crossing timescales was ≲ 50% between the different mass ratios, compared to the fiducial case. Our preliminary investigation suggests that it is difficult for stability studies to precisely constrain Nix and Hydra's mass ratio at present. However, ever-improving orbit determinations should help, both on their own and in combination with stability studies.

We also investigated the effect of setting Charon's eccentricity to zero. We conducted this test with the methods of Section 4.2, i.e., uniform sampling of Keplerian osculating elements. As above, we investigated the four highest mass satellite cases. We find that reducing e_C to zero has a modest stabilizing effect. The median crossing time increased by <50% (40%, 30%, 17%, and 42% for the highest to lowest satellite masses considered). Future work should aim to develop a more systematic understanding of how binary eccentricity, and other orbital parameters, affects stability.

As described in Section 3, the dynamics of the PC system are broadly relevant to our understanding of circumbinary dynamics. We conclude with a brief comparison to the circumbinary planets discovered by Kepler, and emphasize that Pluto is a guide to the circumbinary multi-planet systems that have yet to be discovered.

To make an analogy with planetary systems, we may scale the mass of Pluto to a solar mass. Charon then equates to a 0.12 M_☉ red dwarf. Nix, Hydra, and P4 would then have scaled minimum masses of ∼2M_Mars, ∼3.5M_Mars, and ∼3M_♇, respectively. Kepler 16, 34, and 35 contain planets with masses of 1.1, 0.7, and 0.4 times Saturn (respectively) around binaries of secondary-to-primary mass ratios of 0.3, 0.97, and 0.91, respectively (Welsh et al. 2012). The period ratios of the planets to the inner binaries are 5.6, 10.4, and 6.3 (respectively). Unlike the Pluto system, there is no preference for being near n:1 resonances. However, the Kepler multiples around single stars do show a preference (as yet unexplained) for being near low-order mean motion resonances (Fabrycky et al. 2012).

Kepler's circumbinary planets are relatively close to the circumbinary stability boundary (HW). The period ratio of Kepler 16b is only 14% beyond this limit (Welsh et al. 2012) compared with 41% for Nix. Of course, the observational biases of transit surveys strongly favor shorter period planets (Youdin 2011b).

Kepler has revealed that compact multi-planet systems are common around single stars (Lissauer et al. 2011). However for circumbinary systems, a planet too close to the stability boundary is unlikely to have other planets nearby. Thus, circumbinary planets detected at a safer distance from the stability boundary are more likely to be in circumbinary multi-planet systems—the true Pluto analogs.

Our study shows that multi-planet scattering is more violent around a binary than a single star. Thus, orbital stability will require larger spacing between circumbinary planets, making detection of multi-planet systems more challenging. Nevertheless, Pluto's rich set of companions bodes well for a fruitful search.

In time-periodic potentials—such as an eccentric binary—special orbits form a closed loop when plotted stroboscopically, i.e., when a snapshot is taken once per orbit of the potential. For concreteness, we take this snapshot when the binary orbit is at periapse. These special, stroboscopically closed orbits are the most circular orbits. In the case of resonant orbits, the stroboscopically closed loop breaks into a chain of islands.

Invariant loops in the plane of the binary can be found by iterative methods. Integrations start when the binary is at periapse and the test particle is aligned with the binary. By symmetry, the radial velocity is zero at this location for an invariant loop. Thus at a given radial distance, the only initial condition to vary is the azimuthal velocity, v_o. For the correct choice of v_o, the stroboscopic orbit forms a one-dimensional curve. To iterate toward this orbit, we vary v_o to minimize the radial dispersion of the stroboscopic orbit.

Lithwick & Wu (2008a) measured the radial dispersion in narrow azimuthal bins so that the invariant loop has a nearly constant radial distance in the bin. We found faster and more reliable convergence by fitting the stroboscopic orbits to a (fourth order) cosine series, and measuring the radial dispersion about that best-fit orbit.

In Section 4.3, we present integrations where Nix and Hydra perturb P4 from these most circular orbits. The initial orbital phase along the most circular orbit must be consistent with the initial orbital phase of the binary. This initial phase is not unique as P4 can be advanced along the orbit by an integer number of Charon periods. For different orbital phases, the position of P4 changes relative to Nix and Hydra. Due to the sensitive dependence on initial conditions, stability timescales also change. Thus, our simulations sample many (typically 25) initial phases of P4 for every trial orbital period.

Figure 7 plots the time evolution of the most circular test particle orbits in the vicinity of P4. These orbits only include the perturbations from the PC binary, and not those due to Nix or Hydra. Each curve represents a different orbital period, ranging from 4.8 to 5.2 PC periods. None of these orbits are in a 5:1 resonance with Charon, and the effect of the 5:1 is negligible due to the low e_C, unlike the study of the 4:1 for a high e_C by Lithwick & Wu (2008a).

Zoom In Zoom Out Reset image size

The left panel plots radial distance from the PC barycenter. The non-circularity, given by the magnitude of the radial excursions relative to the mean r, is Δr/(2r) ∼ 2 × 10⁻³. The periodic radial oscillations have well-understood timescales (Lee & Peale 2006). The faster oscillations correspond to the synodic period with Charon (∼5/4 Charon periods) and the first harmonic at half that period. The slower oscillations are on the epicyclic period (∼5 Charon periods) of the test particle itself. The most circular orbits about a circular binary do not have these longer period epicyclic oscillations (Lee & Peale 2006). The epicyclic oscillations are only modestly larger than the synodic timescale variations. When Nix and Hydra are introduced to numerical integrations they rapidly excite a larger epicyclic eccentricity. Thus, Charon's adopted eccentricity is not expected to strongly affect orbital stability. The integrations in Section 4.4 affirm this expectation.

The middle and right panels of Figure 7 plot the osculating (instantaneous) Keplerian a and e, respectively, about the barycenter to demonstrate the non-Keplerian nature of these orbits. The Keplerian a is both systematically larger than the cylindrical radius and varies more significantly with time. The Keplerian e oscillates and reaches values ∼0.015, much larger than the actual radial excursions. This plot demonstrates clearly that setting the osculating e = 0 does not give the most circular orbit, and why e had to be increased above ∼0.01 in Figure 2 to noticeably affect orbital stability.