Numerical investigation of Pluto’s system in two tidal evolution scenarios

Pluto’s small moons, Styx, Nix, Kerberos, and Hydra, present a unique probe of the system’s past dynamical evolution. All four are located on circular orbits aligned with the orbital plane of Pluto’s largest satellite Charon, and all reside close to integer mean-motion resonances with Charon, with orbital periods approximately in the ratios of 3:4:5:6 [1]. Styx, Nix, and Hydra are also locked in a Laplace-like three-body resonance [2]. Furthermore, the small satellites have a high albedo, indicating their ice-rich composition and a likely origin within the Pluto system [1].

The system’s present-day orbital configuration was likely affected by past tidal evolution of the central binary. According to the leading formation hypothesis, Charon originated in a collision of its progenitor with an ancient Pluto and the two bodies then evolved by tidal recession onto their present circular mutual orbit, fully synchronising their rotation frequencies with the mean motion [3]. Alternatively, Charon might have been captured as an originally retrograde satellite of Pluto, which would result in an inward orbital evolution and change in Pluto’s rotation direction [4]. While the former scenario provides a clear mechanism for the small satellites’ formation from the post-impact debris [5], the latter scenario may be reconciled with the idea of a subsequent collision of another Kuiper Belt object with Charon [6]. Nevertheless, whether the inward tidal evolution can lead to the observed configuration of the circumbinary satellite system is currently unknown and has not been tested by [4].

In this preliminary study, we combine calculations of the Pluto-Charon binary’s tidally-induced rotational and orbital evolution in both aforementioned scenarios with an N-body simulation of circumbinary satellites, similarly to the model of [7]. For the tidal interaction, we assume that Pluto and Charon are differentiated bodies composed of liquid and viscoelastic solid layers. This assumption leads to different evolution timescales than the constant time lag or constant phase lag tidal models that are usually used. The tidal evolution is treated within the Darwin-Kaula formalism [8,9], where both the tidal potential and the evolution equations (in particular, Lagrange planetary equations and Euler’s second law) are expanded in series of multiple tidal modes at different frequencies. Since the tidal response of viscoelastic bodies is frequency-dependent, we precalculate tables of complex tidal potential Love numbers of both Pluto and Charon at a wide range of frequencies and use the appropriate values for each mode.

The tidal model is then coupled to the N-body integrator Rebound [10]. For testing purposes, we assume that the only massive particles are Pluto and Charon, and we evaluate the stability of possible circumbinary orbits by including a cloud of ~10,000 massless particles extending up to several present-day radii of Hydra’s orbit. The initial conditions for these particles are chosen from a model distribution of Keplerian elements based on a prior test run over 100 orbital periods of Charon, which is initially emplaced on an eccentric orbit below or above its present-day location (at a larger or smaller semi-major axis, depending on the tidal evolution scenario).

During the simulation runs with a tidally-evolving Pluto and Charon, we monitor collisions of the test particles with the two massive bodies as well as ejections from the system. We then assess both the initial and the final orbital elements of the surviving particles and we discuss the configuration of the resulting system.

 

Acknowledgements

This work has been supported by the Czech Science Foundation grant number 23-06513I.

 

References:

[1] Weaver, H. A. et al. (2016), Science, 351(6279), doi: 10.1126/science.aae0030.

[2] Showalter, M. R. and Hamilton, D. P. (2015), Nature, 522(7554):45-49, doi: 10.1038/nature14469.

[3] Canup, R. M. (2005), Science, 307(5709):546-550, doi: 10.1126/science.1106818.

[4] Efroimsky, M. et al. (2026), Planet. Sci. J., 7(3):68, doi: 10.3847/PSJ/ae47f7.

[5] Canup, R. M. (2011), Astron. J., 141(2):35, doi: 10.1088/0004-6256/141/2/35.

[6] Bromley, B. C. and Kenyon, S. J. (2020), Astron. J., 160(2):85, doi: 10.3847/1538-3881/ab9e6c.

[7] Woo, J. M. Y. and Lee, M. H. (2018), Astron. J., 155(4):175, doi: 10.3847/1538-3881/aab367.

[8] Kaula, W. M. (1964), Rev. Geophys. Space Phys., 2(4):661-685, doi: 10.1029/RG002i004p00661.

[9] Boué, G. and Efroimsky, M. (2019), Cel. Mech. Dyn. Astron., 131(7), id. 30, doi: 10.1007/s10569-019-9908-2.

[10] Rein, H. and Spiegel, D. S. (2015), Mon. Not. R. Astron. Soc., 446(2):1424-1437, doi: 10.1093/mnras/stu2164.