© ESO 2022

1 Introduction

The Pluto-Charon (PC) binary has a mass ratio of ~0.122 and is currently found in a dual synchronous state, which is the typical endpoint of tidal evolution, over 1–10 Myr for this particular system (e.g. Farinella et al. 1979; Correia 2020). The orbit of Charon is almost circular, as confirmed by a 1σ upper limit of 7.5 × 10⁻⁵ (Buie et al. 2012). In the last 15 yr, the system has gained attention due to the discovery of its four small moons (Styx, Nix, Kerberos, and Hydra) and their complex circumbinary configurations. Several formation mechanisms for this satellite system have been proposed (Kenyon & Bromley 2021, and references therein). However, none of the proposed scenarios of the past evolution of the whole system has provided strong conclusions on whether the four small moons could survive the tidal expansion period of the PC binary.

The orbits of the small moons can be described with respect to the barycentre of the PC binary in the nearly circular and coplanar orbital geometry. Several works have studied the dynamical stability of the small moons, considering the possibility of the existence of putative satellites (e.g. Weaver et al. 2006; Kenyon & Bromley 2019c,a). Some have also considered the moons’ masses to be one order of magnitude larger than those determined by observational data (Tholen et al. 2008; Pires Dos Santos et al. 2011; Youdin et al. 2012). The orbital periods of the small moons place them very close to or even inside (depending on the uncertainties) the N/1 mean motion resonances (MMRs) with Charon, namely, 3/1, 4/1, 5/1, and 6/1 for Styx, Nix, Kerberos, and Hydra, respectively. For example, Brozović et al. (2015) obtained the period ratios of a satellite and Charon of 3.1565, 3.8913, 5.0363, and 5.9810 for Styx, Nix, Kerberos, and Hydra, respectively. While the double synchronous state of the PC binary is an indicator of the tidal evolution of the system, the positions of the moons with respect to the MMRs could be an indicator of a smooth migration process.

It is well accepted that Charon was formed as a result of a giant collision that most likely happened when the population of the Kuiper Belt was much denser than today (Canup 2005, 2011; Ward & Canup 2006; Asphaug et al. 2006; McKinnon et al. 2017; Walsh & Levison 2015). This giant impact could also have originated the very small circumbinary satellites Styx, Nix, Kerberos, and Hydra. There is evidence that the satellites are the same age as Pluto; indeed, the crater-counting data from New Horizons imply that the surface ages of Nix and Hydra are at least 4 billion years (Weaver et al. 2016).

Formation theories for the small moons in the PC system include an ‘intact capture scenario’ and a ‘planetesimal capture scenario’. In the scenario of intact capture (Canup 2005, 2011), a proto-Charon grows rapidly, in about 30 h after the collision event, within a massive debris swarm produced during the collision and extending from 4 R_P to 25 R_P (in units of Pluto’s radius). The commonly assumed initial position of Charon in the disc is around 4 R_P (e.g. Cheng et al. 2014a,b; Woo & Lee 2018). Depending on the characteristics of the impact, the initial orbit of Charon varies its form, from circular to highly eccentric (e_C ~ 0.50). The small moons belong to the debris swarm, and their initial positions are generally considered to be closer to Pluto than today.

To place the moons at their current locations, Ward & Canup (2006) proposed the mechanism of resonant transport. The main idea of the method is that, during the tidal expansion of the PC orbit, the small satellites can be captured into resonances with Charon and migrate outwards together with Charon. Several authors have tested this hypothesis and analysed the probability of capture in the MMRs of the kind N/1 (e.g. Lithwick & Wu 2008a; Cheng et al. 2014b; Woo & Lee 2018; Kenyon & Bromley 2021), which are considered to be strongest around binaries (e.g. Cuello & Giuppone 2019; Gallardo et al. 2021). However, Lithwick & Wu (2008b) found some difficulties in adjusting the values of Charon’s eccentricity, e_C: in order to safely transport Nix to the 4/1 MMR, it should be e_C < 0.024, while in order to transport Hydra to the 6/1 MMR, it should be e_C >0.04. Cheng et al. (2014b) have found stable solutions (that is, without ejections from the system) for the test particles at the 5∕1, 6∕1, and 7∕1 MMRs but none at the 3∕1 and 4∕1 MMRs. Moreover, the orbits of the surviving particles were highly eccentric, in contrast with the currently nearly circular orbits of the small moons. In addition, the authors have found that, when the hydrostatic value of Pluto’s zonal harmonic, J₂, was included in the model, there was no stable transport at the regions near the N/1 MMRs.

To overcome the problems of the resonant transport model, the scenario of planetesimal capture gained more attention. This scenario considers a ring of ejected material that ranges up to 60 R_P, or even 200 R_P; this value is even higher in the intact capture scenario (e.g. Pires dos Santos et al. 2012; Desch 2015; Walsh & Levison 2015; Kenyon & Bromley 2019b). Smullen & Kratter (2017) studied the evolution of a debris disc resulting from the Charon-forming impact, established regions of stability according to the tidal evolution of the PC binary, and characterised the collisions onto Charon’s surface that might leave visible craters. Woo & Lee (2018) studied several possibilities for survival of the test particles in the regions of the known moons during the tidal expansion of the PC binary (in situ formation scenario). Applying different tidal models, the authors found the most promising results when considered a constant-Δt tidal model with a large dissipation coefficient (A ~ 40) and an initially circular orbit for Charon (e_C = 0.0); however, in their simulations they did not consider the impact of the zonal harmonic of Pluto’s oblateness.

The results obtained by either the tidal transport scenario or the in situ formation scenario are still inconclusive due to the poor knowledge of the precise orbital elements and masses of the small moons. The most precise orbital parameters of the PC system, reported in Brozović et al. (2015) and Showalter & Hamilton (2015) (see Table 1), present a dispersion of several orders of magnitude in the eccentricities and inclinations of the small moons. The solution from Showalter & Hamilton (2015) locked Pluto and Charon to the ephemerides of the Jet Propulsion Laboratory (JPL), PLU043. Our data are taken from JPL Horizons and correspond to the pre-computed solution PLU058/DE440, a fit to ground-based Hubble Space Telescope and New Horizons spacecraft encounter astrometry in the interval 1965–2018. To illustrate the uncertainties in orbital fits, we show in Table 1 the in-orbit errors (along the track) reported by Brozović et al. (2015) and also include the amplitude observed in the variation of orbital elements (see also Fig. A.1). The precise orbital dynamics of the small moons depends on the initial osculating orbital elements and the masses of the moons.

With this in mind, we introduce an additional mechanism that could provide a robust explanation for the existence of the four small circumbinary satellites at their current positions. For this, we first give a global view of the dynamics of the current Pluto system of small satellites in Sect. 2. In Sect. 3 we present the model, which describes the tidal interactions between Pluto and Charon, and discuss the choice of the parameter values adopted in this paper. In the next section we analyse how the tidal expansion of the PC binary, starting at different initial configurations, could affect the behaviour of the small moons located at their current positions (Sect. 4). In Sect. 5we analyse,in the frame of the in situ formation scenario, the effects of the tidal evolution of Charon’s orbit on a large grid of parameter values and initial conditions. In that section we also study the impact of the zonal harmonic of Pluto’s oblateness on the moons’ dynamics. To overcome the problem of survival of the moons during the passages through the low-order MMRs with Charon, we investigate the effects of a mass-accreting Charon on the moons’ behaviour in Sect. 6. Finally, we present our conclusions in Sect. 7.

2 Dynamical portrait of the current PC system

It is generally accepted that the current orbital configuration of the PC binary and the four small moons is a product of the evolutionaryhistory of the whole system during its lifetime, since the collision event that gave rise to Pluto’s satellites. In this context, a detailed analysis of the current relative positions of the Pluto system members may provide important constraints on the dynamical past of the whole system.

A representative picture of the PC system is shown in Fig. 1; it is calculated with data from the JPL Horizons site (see Table 1). To construct it, we used a modified Jacobi reference frame, in which orbital elements of the small moon are referred to the centre of mass of the PC system and the mutual interactions between the moons are neglected¹. In this reference frame, the motion of a particle of Styx’s mass is simulated through the numerical integration over 200 yr (~3500 Styx orbital periods) and subsequently analysed using Δe or Δi indicators to better address the structure of the resonances. We calculate Δe as the amplitude of maximum variation in the orbital eccentricity of the satellites during the integrations, $\Delta e=\left(e_{\text{max}}-e_{\text{min}}\right)$. Similarly, $\Delta i=\left(i_{\text{max}}-i_{\text{min}}\right)$ and $\Delta a=\left(a_{\text{max}}-a_{\text{min}}\right)$. Regions with higher Δe lead to chaotic motion, and Δe is a powerful indicator of secular and resonant dynamics that also serves to identify separatrices of MMRs (e.g. Dvorak et al. 2004; Ramos et al. 2015).

The instabilities in the motion of the particles are represented using a colour scale on the maps in Fig. 1, on the a–e plane (left panel) and a–I plane (right panel) of the initial conditions, where a, e, and I are the semi-major axis, the eccentricity, and the inclination of the particle, respectively. In the stability maps, the initial conditions for Charon are given in Table 1, whereas the initial values of the angular variables of the circumbinary moon are those given to Styx in Table 1. We note that all moons lie beyond the black continuous curve in the left panel, which delimits the domain of stable motion, according to the Holman–Wiegert criterion (HWC) derived in Holman & Wiegert (1999) for circumbinary configurations and also noticed by Smullen & Kratter (2017).

Figure 1 shows that the orbits of the moons are nearly circular, nearly planar, and nearly resonant. Although these features of the Pluto system have already been pointed out in previous studies, it is still worth analysing them in more detail. Indeed, as seen on the stability maps on the a–e and a–I planes, the small satellites are confined to the very small eccentricity and inclination regions of regular motion (blue colour). According to secular perturbation theories, very low eccentricities and inclinations are characteristic of the dynamical systems that have an angular momentum deficit (AMD) close to zero (Michtchenko & Rodríguez 2011). Regarding initially excited dynamical systems (in this case, due to the collisional origin of Charon), the low-AMD configurations are generally attained by the systems evolving under weak and slow dissipation.

The interpretation of the near resonant distribution of the small satellites requires the analysis of an additional dynamical mechanism, that is, a capture into MMRs. The stability maps in Fig. 1 show the MMRs produced by Charon as vertical strips of chaotic motion (red colour). The strong low-order MMRs, such as 3/2, 5/3, 2/1, 5/2, and 3/1, are located inside the domain of instability defined by the HWC (black curve in the left panel). In the low-eccentricity region, the 3/1 MMR can be regarded as an inferior limit of stability of the motion, and Styx is stuck to it. This moon lies so close to the 3/1 MMR that uncertainties in the determination of the angular variables of its orbit, such as the mean longitude and the longitude of the pericentre, might put it closer to the chaotic layer of the 3/1 MMR (e.g. in-orbit and radial errors reported in Brozović et al. 2015). The same is also true for Nix with respect to the 4/1 MMR (see the Styx and Nix panels in Fig. A.2) but is most noticeable for Kerberos and Hydra with respect to the 5/1 and 6/1 MMRs, respectively (see the bottom panels of Fig. A.2). This fact suggests that the small moons may have been involved in one of the MMRs at some period of their common history. Therefore, in this paper we analyse the interaction of the small moons with the main MMRs during the tidal expansion of Charon’s orbit.

It is worth emphasising that the uncertainties in the determination of the moons’ masses and the actual positions can provide a general view of the phase space of the PC system, but without specifying whether the objects are evolving near or inside the MMRs.

Fig. 1 Dynamical maps on the a–e plane (left panel) and a–I plane (right panel) of the moon’s initial conditions calculated for the PC binary and a small moon of Styx’s mass with initial Styx orbital elements from JPL in Table 1. The current positions of the moons are shown by red crosses. The locations of the main MMRs are indicated, and the grey line in the right panel indicates the current inclination of Charon. The colour scale varies with Δe and Δi, increasing from blue (regular motion) to red (chaotic motion), as shown in the colour bars on the top of the figure.

3 Tidal model

In this section we describe the tidal model, which we apply to investigate the orbital expansion of Charon and the implications for the orbital motion of the small moons. As described in the previous section, the orbital configuration of Pluto’s satellites is better represented in terms of the modified Jacobi elements. The initial Jacobi orbital elements are transformed to the Jacobi rectangular coordinates and velocities in a three-dimensional reference frame centred on Pluto. The equations of motion are solved using an N-body integrator with adaptable step size (for details, see Rodríguez et al. 2013). We considered the system consisting of Pluto, Charon, and an external small moon, where two large bodies undergo mutual tidal interaction but the external moon is not (directly) tidally affected, due to its small size and mass. We adopted a classical linear tidal model, in which the deformations of the bodies are delayed by a constant Δt, usually referred to as the tidal time lag (Mignard 1979). In this scenario, the PC tidal interaction is modelled following Cheng et al. (2014a)², who assumed the collisional origin of Charon, when Charon emerges with a non-zero orbital eccentricity, e_C, and an initial semi-major axis of a_C = 4 R_P (Ward & Canup 2006).

Both Pluto and Charon would be spinning after the collision, with the spin axes oriented in some directions with respect to the orbital plane of the binary. Due to the relatively small size and mass of Charon, in comparison with Pluto, the contribution of its spin to the total angular momentum is small, and it evolves quickly to the asymptotic tidal spin rate (in less than 10 yr). For the sake of simplicity, in this work we assume zero obliquity for both Pluto and Charon; the model of the PC binary, accounting for the arbitrary obliquities and the Sun’s perturbations, is studied in Correia (2020).

The total angular moment of the PC binary can be written, in the chosen reference frame and with the assumed approximations, as $$L={\mathcal{C}}_{\text{P}}{\Omega }_{\text{P}}+{\mathcal{C}}_{\text{C}}{\Omega }_{\text{C}}+M_{\text{PC}}{n}_{\text{C}}{a}_{\text{C}}^2\sqrt{1-e_{\text{C}}^2},$$(1)

where M_PC = m_Pm_C∕(m_P + m_C) and n_C are the reduced mass of the binary and the mean motion of Charon, respectively. According to classical theories, L is conserved during the tidal expansion of the PC binary. Hence, its amount for the PC system can be obtained knowing the masses of Pluto and Charon, m_P and m_C, and the current values of a_C and n_C. Indeed, at the dual synchronous state of the current PC binary, Charon’s orbit is circularised, that is, we can assume e_C = 0 and Ω_P = Ω_C = n_C.

Then, by substitution of these conditions into Eq. (1), the L value is obtained; in the following, this value is used to obtain the angular velocity of Pluto, Ω_P, resolving Eq. (1) for given starting values of a_C, e_C, and Ω_C, the angular velocity of Charon’s rotation.

Following Cheng et al. (2014a), we define the dissipation parameter, A, as $$A=\frac{k_{2C}}{k_{2P}}\frac{\Delta t_{\text{C}}}{\Delta t_{\text{P}}}\text{Big}{(\frac{m_{\text{P}}}{m_{\text{C}}}\text{Big})}^2\text{Big}{(\frac{R_{\text{C}}}{R_{\text{P}}}\text{Big})}^5,$$(2)

where (k_2C, k_2P) are the second-order Love numbers and (R_P, R_C) are the radii for Pluto and Charon, respectively. We note that, for the masses and the radii of Pluto and Charon (R_C = 606 km), $$\text{Bigg}{(\frac{m_{\text{P}}}{m_{\text{C}}}\text{Bigg})}^2\text{Bigg}{(\frac{R_{\text{C}}}{R_{\text{P}}}\text{Bigg})}^5\simeq 2;$$(3)

thus, the constant A can be closely understood as the ratio of time lags between Charon and Pluto. Cheng et al. (2014a) have also shown that, for an appropriate value of the ratio (A ~ 10), the eccentricity of Charon’s orbit, e_C, remains roughly constant during most of the tidal evolution.

Our current knowledge of the physical properties of the PC system does not allow us to constrain the range of A values (see Eq. (2) and Lainey 2016; Quillen et al. 2017). However, it is easy to show through classical tidal theories that for A << 1 tidal forces on Pluto dominate the tidal evolution of the binary, whereas for A >> 1 the orbital expansion and circularisation are dictated by the tides on Charon. Previous works adopted different ranges for the possible values of A, namely, 1 < A < 18 (Cheng et al. 2014a,b); A = 10 or A = 40 (Woo & Lee 2018); and 1 < A < 16 (Correia 2020). In this work, we initially explored the wide range of 0.01 < A < 100; however, in order to place test moons at the current Styx position, we adopted the range 8 < A < 40.

Cheng et al. (2014a) show the tidal evolution of the PC system for different initial eccentricities of Charon. Here, we summarise the common features observed in the numerical integrations, which we refer to in our discussions below: (i) Independently of the initial values of Charon’s eccentricity, e_C, it presents a peak at some moment during the PC expansion; (ii) the time evolution of Ω_P ∕n is correlated with the peak of e_C and can reach values of ~ 10; (iii) the semi-major axis of Charon attains its current value (16.5 R_P) and the orbit is circularised, whereas both rotations synchronise with the orbital period after around 3 × 10⁶ yr. It is worth emphasising that, according to tidal evolution theories, the double synchronous rotation is only attained when the orbit is completely circularised (see Cheng et al. 2014a). Also, at the end of the evolution, the angular momentum exchange between Pluto and Charon ceases and both orbital and rotational components attain their equilibrium configurations in the system.

4 Dynamics of the moons in a tidally evolving PC system

In this section we investigate the dynamics of the small moons located at their current positions during the tidal expansion of Charon’s orbit. For this, we use the same techniques described in Sect. 2. Considering that Charon’s orbital elements vary during the migration, we mapped the neighbourhood of each moon, applying the conservative model for different values of the semi-major axis, a_C, and the eccentricity, e_C, of Charon’s orbit and the fixed mass and initial conditions of the satellite (Table 1). In this way, we can visualise the dynamical behaviour of each moon at different configurations of the PC binary acquired in the past. Figure 2 shows the dynamical maps on the a_C–e_C planes, obtained for the current locations of Styx, Nix, Kerberos, and Hydra retrieved from JPL Horizons. The black curves show the positions of Charon on these planes as obtained through simulations of its tidal expansion in the past (described in the previous section). Three values were used for Charon’s eccentricity, namely 0.01, 0.1, and 0.2, while the other parameters were fixed at a_C ∕R_P = 4, Ω_C ∕n_C = 2, A =10, and Δt_P = 600 s. By construction, the three paths converge to the current position of Charon, at a_C ≅ 16.5 R_P and e_C ≅ 0.

Now we can investigate the dynamical behaviour of each moon according to the position of Charon at any moment of its history. Indeed, for a chosen position of Charon in one of its paths on the a_C–e_C planes in Fig. 2, we can identify dynamical features that affect the behaviour of the moon at that instant. These features are mainly MMRs, which are identified by plotting the moon’s dynamics with the colour scale. In this way, we can distinguish between the stable (blue-white) and unstable (red) behaviour of satellites. The instabilities of the motion are clearly related to the MMRs, whose locations are indicated by the corresponding ratio in Fig. 2.

The strongest MMRs are of low order, 3/1 and 4/1, and we expect that the dynamical stability of the moons would be significantly affected by these resonances. Figure 2 shows that Charon’s paths (black curves) never cross the 3/1 MMR, and, consequently, none of the four small moons suffer the destabilising effects of this resonance. Currently, Styx stays close to this resonance (left-top panel), but its motion is still stable. It is worth noting here that the 3/1 MMR can be considered as an inferior limit of the stability of the circumbinary motion, as shown in Fig. 1.

According to Fig. 2, the actual position of Styx has been crossed by the 4/1 MMR at some moment in the past, as has that of Nix (right-top panel) at another period. Our simulations of the moons’ dynamics (described in the next section) show that the perturbations produced by this resonance increase with the increasing eccentricity of Charon. For Charon’s paths with e_C = 0.1 and 0.2, Styx and Nix spend some time evolving in the 4/1 MMR, which might be sufficient to strongly excite their eccentricities and even eject them from their current positions. On the contrary, for Charon on a quasi-circular orbit (e_C = 0.01), the passages of Styx and Nix through the 4/1 MMR may occur quickly and without significant excitation of their motions.

The width and, consequently, the dynamical effects of the higher-order MMRs, such as 5/1, 6/1, and, 7/1, rapidly decrease with decreasing order, as seen in Fig. 2. We expect that their passages through the current positions of the small moons would not produce perturbations to destabilise their motion, particularly for Charon evolving on a quasi-circular orbit during these passages. Thus, we conclude that, in order to ensure the simultaneous orbital stability of the four small moons during the orbital expansion of the PC binary, Charon’s eccentricity should be constrained to lower values. This fact points towards a dynamical scenario in which the adequate choice of Charon’s eccentricity at the beginning of the tidal expansion is essential in order to guarantee the orbital stability of the system of moons. In the next section we numerically simulate the impact of a migrating Charon on the orbital motions of Styx, Nix, Kerberos, and Hydra by taking simultaneous effects of the tidal evolution and gravitational perturbations into account.

5 Simulations of the moons’ behaviour under the tidal expansion of Charon’s orbit

In this section we analyse the behaviour of Styx, Nix, Kerberos, and Hydra during the orbital expansion of Charon’s orbit and their passages through the MMRs. Here, we assumed the in situ formation of the small moons, that is, our numerical simulations start with the satellites at their current positions. We worked in the modified Jacobi reference frame, where small moons are orbiting the barycentre of the PC binary. In this frame, we neglected the mutual perturbation between the moons, due to their small masses, which allowed us to solve their orbital motions in a separate way. We solved the exact equations of motion for each satellite according to the tidal model described in Sect. 3 and following Rodríguez et al. (2013). We considered coplanar configurations and the initial mean anomalies, and longitudes of pericentre were set to zero.

We performed a set of numerical simulations, varying the main parameters of the problem, namely, the initial semi-major axis and eccentricity of Charon’s orbit (a_C, e_C), the initial rotation velocities of Pluto and Charon (Ω_C and Ω_P), the time lag of Pluto, Δt_P, and the constant of dissipation, A (see Table 2). The initial orbital angles of Charon were set to zero. We fixed a fiducial set of parameters such that a_C ∕R_P = 4, e_C = 0.001, Ω_C ∕n_C = 2, A =10, and Δt_P = 600s. Then, the simulations were done varying the parameters one by one, choosing the values listed in Table 2. It is worth noting that the initial value of Ω_P is not a free parameter but is determined by the conservation of the angular momentum of the binary through Eq. (1).

The masses and current semi-major axes of the external moons were taken from Brozović et al. (2015) (see Table 1). The initial orbital eccentricities of the moons were assumed to be almost circular (e = 10⁻⁵). We note that in this section we work in the frame of the planar problem, when the small moon evolves on the orbital plane of the PC binary. Physical parameters, such as the masses and radii of Pluto and Charon, are the same as described in Sect. 3. Since the tidal evolution of the PC binary is not affected by the external small moon, we do not show the variations in their orbital and rotational components in this section (see Fig. 2 in Cheng et al. 2014a). In addition, we performed simulations for both Δt_P = 0.06 s and Δt_P = 0.6 s, but the results obtained were discarded since Charon did not reach its current semi-major axis after 100 Myr in these cases. We integrated the equations of motions over 100 Myr for Δt_P = 6 s and Δt_P = 60 s; for Δt_P = 600 s, we considered 10 Myr as the integration timescale.

In the following we summarise the results obtained for each of the four moons, emphasising the cases in which the final orbital configuration of the small moon is similar to its current one. We note that, despite the fact that the moons are currently in nearly resonant configurations, the uncertainties in their masses and orbital parameters could still place them inside the MMRs (see our discussion in Sect. 1 and Fig. A.2). In Table 3 we include the cases for which each of the moons reaches the end of the simulation in stable motion. Additionally, we identify the initial conditions for which the final eccentricity of the moon is smaller than 0.05, which is compatible with the current eccentricities of the moons (see our Fig. A.1).

Fig. 2 Dynamical maps of Styx, Nix, Kerberos, and Hydra on the aC–eC plane of Charon’s initial conditions. On each map, black curves are Charon’s tidal trajectories, obtained for initial eC values of 0.01, 0.1, and 0.2; the initial conditions of the corresponding moon are fixed at those given in JPL from Table 1. The blue tones represent regular motions, while the brick tones correspond to increasing instabilities and chaotic motion; the colour bar on the top shows Δe. The locations of the main MMRs are indicated.

5.1 Styx

Styx is located at a mean distance of 42 650 km (≃35.9 R_P) from the barycentre of the PC binary and is close to the 3/1 MMR with Charon (see Figs. 1 and A.2). Other initial conditions different from those shown in Table 3 tested forthis satellite lead to either ejection of the small moon from the system (Fig. 3, left panel) or a final orbit that is different from the current one (Fig. 3, right panel). In most cases, instabilities occur when Styx approaches the 4/1 or 5/1 MMRs, in this last case, only for a high initial e_C.

Figure 4 illustrates Styx’s orbital evolution obtained for A = 100. We note that, despite an apparently slow increase, the semi-major axis of Styx remains close to the current value. The eccentricity acquires a mean value of 0.02, varying roughly between 0.007 and 0.04, which is also consistent with the oscillation of Styx’s current eccentricity. The final value of the mean-motion ratio is larger than the nominal location of the 3/1 MMR (see Fig. 1); the difference between the pericentre longitudes of Styx and Charon, Δϖ, oscillates around 180°, which is an indicator of the dissipative migration (see Sect. 4). We have verified that, for this particular simulation, all critical angles associated with the 3/1 MMR circulate.

The results obtained for Styx point out that the scenario with an initially quasi-circular orbit of Charon and a high dissipation ratio (large A) favours the stability of Styx’s motion during the orbital expansion of the PC binary. However, we note in Table 3 that only very specific combinations of the parameters allow Styx to remain on a stable orbit as the 4/1 MMR with Charon is crossed. For most of the explored initial conditions, Styx’s motion is strongly excited by the passage through this resonance, resulting in the escape from the system. In Sect. 6 we test the orbital stability of Styx (and the other small moons) assuming anon-time accretion of Charon’s mass.

Fig. 3 Two typical examples of the Styx orbital evolution resulting in either instability (left panel) or the final configuration, which is significantly different from its current orbit (right panel). Left panel: initial conditions used are aC ∕RP= 5, eC = 0.001, ΩC ∕nC = 2, A = 10, and ΔtP = 600 s. Styx is ultimately ejected from the system when the system crosses the 4/1 MMR between Charon and Styx (not shown), resulting in a strong increase in Styx’s eccentricity. Right panel: initial conditions used are aC ∕RP= 4, eC = 0.001, ΩC ∕nC = 2, A = 10, and ΔtP = 6 s. The orbital eccentricity is excited, due to the crossing of the 4/1 and 7/2 MMRs (not shown), to attain the mean value ≃ 0.1, which is substantially larger than the current Styx eccentricity.

5.2 Nix

Nix lies at a mean distance of 48 690 km (≃41 R_P) from the PCbarycentre and is currently close to the 4/1 MMR with Charon. The values of the initial conditions resulting in final configurations similar to the current orbit of Nix are shown in Table 3. Figure 5 shows examples of the orbital evolution obtained for a_C∕R_P = 3 and a_C∕R_P = 7. For a_C ∕R_P = 3 (orange curves), the resulting orbit clearly deviates from the current orbital configuration of Nix. In this case, the trapping into the 4/1 MMR excites Nix’s eccentricity, up to e = 0.16. The final semi-major axis is larger than the current value, with a difference of ≃ 2 R_P. Moreover, the critical angle associated with the 4/1MMR, ϕ_4∕1 = 4λ − λ_C − 3ϖ_C, librates around 180°, as do all other angles associated with the same MMR, while Δϖ oscillates around 0° (Fig. 5, bottom right). In this case, e_C reaches ≃4 × 10⁻⁵ at 10⁷ yr in one of the simulations due to perturbations from Nix (not shown here).

On the other hand, for a_C∕R_P = 7 (blue curves), the orbital evolution of Nix occurs with the eccentricity and semi-major axis close to their current values. All critical arguments associated with the 4/1 MMR circulate in the case a_C ∕R_P = 7, indicating that the moon is out of the resonance, although it is still very close to it. In this case, the mean-motion ratio oscillates around a mean value slightly smaller than the nominal position of the 4/1 MMR (Fig. 5, bottom left), in accordance with the current position of Nix given in Brozović et al. (2015) (see Figs. 1 and 2).

The results shown in Table 3 indicate that, in comparison with Styx, almost 40% of the tested initial conditions result in the final orbital configurations, which are similar to the current orbit of Nix.

Fig. 4 Examples of initial conditions for Styx. Top row: time evolution of Styx’s semi-major axis (left) and eccentricity (right) during the tidal expansion of the PC binary. The initial conditions used in the simulations are aC ∕RP= 4, es C = 0.001, ΩC ∕nC = 2, A = 100, and ΔtP = 600 s. The final values of the orbital elements remain close to the current ones. Bottom row: time evolution of the mean-motion ratio (left) and Δϖ (right). We note that the mean value of nC∕n is larger than 3. The secular angle, Δϖ, oscillates around 180°.

Fig. 5 Examples of initial conditions for Nix. Top row: time variations in Nix’s semi-major axis (left) and eccentricity (right) for two different initial values of aC. The initial conditions are aC∕RP = 3 and aC ∕RP = 7, eC = 0.001, ΩC ∕nC = 2, A = 10, and ΔtP = 600 s. For aC ∕RP = 3, the eccentricity of Nix is excited to a mean value close to 0.16. Bottom row: time evolution of the mean-motion ratio (left) and two angles for Nix (right), namely, the secular angle, Δϖ, and the 4/1 MMR critical angle ϕ4∕1 = 4λ − λC − 3ϖC. The mean value of nC∕n is slightly smaller than 4 for aC∕RP = 7, whereas ϕ4∕1 and Δϖ librate around 180° and 0°, respectively, for aC∕RP = 3.

5.3 Kerberos

Kerberos is the third small satellite orbiting the PC binary, at a mean distance of 48.6 R_P, which places it close to the 5/1 MMR with Charon. Table 3 shows that all numerical simulations result in good agreement with the current orbit of this small moon.

Figure 6 shows the orbital evolution of Kerberos during the tidal expansion of the PC binary for two values of the initial Charon eccentricity, e_C = 0.001 (grey curves) and e_C = 0.1 (red curves). For both initial conditions, the output of the simulations satisfactory fits the current orbit of the moon. All critical angular combinations of the 5/1 MMR circulate, and the ratio n_C ∕n oscillates around a mean value slightly larger than 5, which places the moon outside, but still close to, the 5/1 MMR (see Fig. 2). We show the angle ϕ_5∕1 = 5λ − λ_C − ϖ_C − 3ϖ, obtained for e_C = 0.001, in the bottom-right panel of Fig. 6. At around 5 × 10⁶ yr, this angle changes its behaviour from circulation to libration around 0°.

Fig. 6 Examples of initial conditions for Kerberos. Top row: orbital evolution of Kerberos for two initial conditions of eC, namely, 0.001 and 0.1. The other parameters are aC∕R = 4, ΩC ∕nC = 2, A = 10, and ΔtP = 600 s. These simulations are successful according to the adopted criteria (see Sect. 5). Bottom row: time variations in the mean-motion ratio (left) and ϕ5∕1 = 5λ − λC − ϖC − 3ϖ (right). We show the angle for the specific case of eC = 0.001. The mean value of nC∕n is larger than 5.

5.4 Hydra

Hydra is the most distant known satellite of the PC binary, with a current barycentric semi-major axis of a ~ 54.5 R_P. This small satellite is close to the 6/1 MMR with Charon. Most of the simulations reproduce the current orbit of Hydra, as shown in Table 3. Figure 7 shows the evolution for Δt_P = 6 s (brown curves) and Δt_P = 60 s (grey curves). We can note that, for Δt_P = 60 s, Hydra is trapped in the 6/1 MMR, with the critical angle ϕ_6∕1 = 6λ − λ_C − ϖ_C − 4ϖ librating around 0°. Moreover, the ratio n_C ∕n is slightly smaller than 6 for this specific trapping. In the case with Δt_P = 6 s, all critical angles of the 6∕1 MMR are circulating.

As in the case of Kerberos, all initial conditions result in a stable motion of the test moons, with e < 0.05, indicating good agreement with the current orbit of Hydra.

5.5 Discussion of the moons’ behaviour and proximity to MMRs

Analysing the parameter values in Table 3, we realise that it is unlikely that we can collect a specific set of parameters and initial conditions that would result in a final configuration of all four satellites that is consistent with the current one. This fact is mainly due to the difficulty of finding Styx in a stable configuration.

Since the determination of the satellite’s current orbits is affected by uncertainties, it is not clear whether they are captured or not in the MMRs. Thus, for the sake of completeness, we show in Table 4 the values of the parameters that provide the final orbits of the moons trapped in the nearest MMR with Charon with at least one of the associated critical angles librating. It is remarkable that the motion of Styx in the 3/1 MMR is always unstable, which confirms our suggestion that the 3/1 MMR delimits the domain of stable motion in the PC binary.

For the sake of completeness, we compare our results with those previously reported by Woo & Lee (2018). Their model considered A = 10 and A = 40 and an initial e_C = 0 and e_C = 0.2, in the case of tidal model of constant Δt. For a large A and an initial e_C = 0, all test particles survive the entire simulation (see their Table 3, last row). Moreover, the particles near the MMRs with Charon have final eccentricity values larger than the current ones. On the other hand, for small and large A and initial e_C = 0.2, the final orbits do not match the current ones. These results are in good agreement with ours; however, Woo & Lee (2018) only considered fixed values of Δt = 600 s and Ω_P ∕n = 5.65 and an initial distance for Charon of a_C = 4 R_P, while our study presents a wider range of initial parameters and tidal dissipation.

Fig. 7 Examples of initial conditions for Hydra. Top row: orbital variation in the semi-major axis (left) and eccentricity (right) of Hydra according to the adopted criteria. The initial conditions are aC ∕RP= 4, eC = 0.001, ΩC ∕nC = 2, A = 10, and ΔtP = 6 s. Bottom row:time variation in the mean-motion ratio (left) and ϕ6∕1 = 6λ − λC − ϖC − 4ϖ (right) for two values of ΔtP, namely, 6 s and 60 s. We note that, for ΔtP = 60 s, ϕ6∕1 librates witha large amplitude around 0°. For ΔtP = 6 s, all critical angles associated with the 6/1 MMR circulate.

5.6 Survival time of the small moons under J₂ evolution

In order to refine our model, we introduced an additional force due to the oblateness of Pluto. Pluto’s polar oblateness, parameterised by the zonal harmonic, J₂, is modelled as (Cheng et al. 2014b) $$J_2=\frac{k_{\text{fP}}{R}_{\text{P}}^3{\Omega }_{\text{P}}^2}{3G{m}_{\text{P}}},$$(4)

where k_{f P} is the second-order fluid Love number of Pluto and Ω_P is the angular rotation velocity of Pluto. Recently, Correia (2020) derived from the shape of Pluto a static value of J₂ = 6 × 10⁻⁵. In this work, we use both formulations, adding them to the model described in Sect. 3.

We applied the model to study the behaviour of a ring composed of 300 particles, each with a very small mass of 1.49× 10¹⁵ kg (thus, we can neglect theirmutual perturbations and contributions in the total angular momentum of the system), at different positions chosen from the interval of 10 R_P < a < 55 R_P of the barycentric semi-major axis. We used the constant Δt model, A = 10, and the different initial values of Charon’s eccentricity, e_C. It is worth noting that the value of e_C constrains Ω_P (and, equivalently, the initial rotation period of Charon). For example, starting at a_C = 4R_P, with e_C = 0.001 and Ω_C ∕n_C=2 (i.e. P_C = 0.38136 d) and using Eq. (1), we calculate the corresponding initial rotation period of Pluto as being P₀ = 0.13858 d (P₀ = 2π∕Ω_P).

We worked with two kinds of disc particles: one is coplanar with the orbital plane of the PC binary, and the other is inclined, with the angle randomly chosen from the interval − 0.01° < i < 0.01°. In the coplanar disc, the particles start on nearly circular orbits, with eccentricities uniformly distributed in the interval 0.0 < e < 10⁻³. For Charon’s initial configuration, we considered a_C = 4 R_P and three values of the eccentricity, e_C, namely, 0.001, 0.1, and 0.2. In the case of the inclined disc, the particles remain on nearly circular orbits, while a_C = 4 R_P and e_C = 0.001. In both cases, we set random initial conditions for the angular variables of Charon’s and the particles’ orbits. The system composed of the PC binary and the ring of particles was integrated over 3 × 10⁶ yr (approximately 150 × 10⁶ periods of the binary)³.

The case of the initially inclined disc of particles was considered due to recently published results (Ćuk et al. 2020) that demonstrate that the tidal evolution of particles at low inclinations can induce resonant trapping and the excitation of the particles’ inclination, which is incompatible with the observed low inclinations of the moons in the PC system. The tests with an initially eccentric orbit for Charon were inspired by the results of the smooth particle hydrodynamic simulations reported in Canup (2005, 2011).

Figure 8 shows the survival times of the particles from the different discs. The results were obtained considering the zonal harmonic, J₂, which evolves as the PC binary tidally expands, according to Eq. (4). We also investigated the behaviour of the disc considering the cases of J₂ = 0 and J₂ = 6 × 10⁻⁵, the latter value derived from Correia (2020); the results obtained are discussed in the appendix. We note in Fig. 8 that the particles can survive the tidal expansion of the PC binary (that is, remain in the PC system after 3 × 10⁶ yr) only when Charon starts on a nearly circular orbit (e_C = 0.001) and only in the region beyond Nix’s orbit. Also, there are no significant differences between the coplanar and low inclined discs (blue crosses and orange dots, respectively). Contrarily, in the case of an initially eccentric Charon orbit, all particles (green crosses and red dots) are ejected from the system in less than 10³ yr due to capture in the MMRs and the subsequent excitation of the eccentricities, as described in Cheng et al. (2014a) and Kenyon & Bromley (2019b).

For the quasi-circular Charon orbit, the low-order 3/1 and 4/1 MMRs, at the positions of Styx andNix, respectively, are responsible for the rapid ejection of the particles (less than 10⁴ yr). In this case, the survivors are observed only beyond Nix’s orbit. Moreover, even farther, in the regions around the current positions of Kerberos and Hydra (5/1 and 6/1 MMRs, respectively), most of the fictitious moons do not survive more than 10⁵ yr. It is clear that these times are incompatible with the hypothesis of primordial moons after the system settled into its present configuration through an impact on Charon (Bromley & Kenyon 2020).

It is worth observing that, applying the simplified model, which considers the constant value of J₂ (zero or J₂ = 6 × 10⁻⁵), we obtain that the particles starting on nearly circular and nearly planar orbits at the current positions of Styx and Nix can survive the tidal expansion of the PC binary (for details, see the appendix). However, in the more realistic model, which accounts for the zonal harmonic, J₂, evolving during the expansion, the additional perturbations destroy the stability of the particles in the range 10 R_P < a < 55 R_P. Moreover, no particles remain at the current positions of Styx and Nix. We repeated these experiments with a higher dissipation ratio, A = 20, and qualitatively obtain the same result.

Fig. 8 Survival times of the small particles, for the different discs extending from 10 RP to 55 RP, during the tidal expansion of the PC binary, with the dissipation parameter A = 10, under the effect of the evolving zonal harmonic, J2. The colour symbols are used to identify the particles from the different discs (see text).

6 Accreting Charon’s mass

In this section we introduce a scenario that considers the mass accretion of proto-Charon from the disc of debris left soon after the giant impact. The growth of Charon occurs simultaneously with the tidal expansion of its orbit around Pluto and could affect the behaviour of the small circumbinary moons described in the previous sections.

First, we analyse the evolution of the expanding binary with a growing Charon in Sect. 6.1. Then we introduce an additional small moon (either Styx or Nix) to assess the parameters of the tidal evolution in Sect. 6.2. Finally, we simulate the behaviour of the disc of particles in the PC system in Sect. 6.4. In this section we choose the plane of Charon’s orbit as a reference plane.

6.1 Tidal evolution of the PC binary with the accreting Charon mass

We started by considering the evolution of the PC pair described in Sect. 3, now including the accretion of Charon from the disc of debris. The duration of the accretion process was arbitrarily fixed at 10⁴ yr after a giant impact, starting with Charon’s mass at 60% of its current value, which is in accordance with the results in Canup (2005) and Kokubo et al. (2000). Following Kokubo et al. (2000), we adopted a logarithmic law for the mass accretion, but we also considered an oligarchic growth defined as m ∝ t^1∕3. The evolution of the increasing mass of Charon is shown in the left panel of Fig. 9 for both models. Assuming Charon to be a homogeneous sphere, we calculated its radius for a constant density of ρ = 1.70 g cm⁻³ and show its time evolution in the middle panel of Fig. 9. Finally, the right panel shows the time evolution of the zonal harmonic, J₂, of Pluto’s oblateness during the tidal migration of a growing Charon. According to Eq. (4), J₂ is a function of the orbital spin of Pluto Ω_P, whose starting value was chosen as described below. For the chosen Ω_P, J₂ starts at ~ 0.17 and reaches 8 × 10⁻⁵ at the end of the migration, for both logarithmic and oligarchic laws of the mass accretion. We note that both models represent almost the same evolution of J₂. We also show in the right panel of Fig. 9 the constant J₂ of 6 × 10⁻⁵ derived in Correia (2020).

It is worth emphasising that, in the growing mass scenario, the total angular momentum of the PC binary given in Eq. (1) is no longer conserved. Thus, we need to obtain the initial value of Pluto’s spin, Ω_P, and the dissipation ratio, A, in such a way that both can guarantee the current double synchronous configuration of the binary with orbital parameters compatible with those from Table 1.

Using a logarithmic law for Charon’s accretion and fixing A = 10 and e_C = 0.001, we varied the initial rotation period of Pluto, P₀ = 2π∕Ω_P, in the range from 0.12 d to 0.21 d in order to determine the value that assures, at the end of the migration, the current configuration of the PC binary. We chose a low initial eccentricity for Charon due to constraints already described in Sect. 5.6. We simulated the tidal evolution of the PC binary and show the results in Fig. 10, where the time evolution of Charon’s semi-major axis, a_C (in units of R_P), is presented in panel (a). The horizontal dashed line, which shows the final state of the current system, corresponds to the initial value P₀ = 0.149 d. Smaller initial values of P₀ lead to wider PC binaries, while larger values of P₀ lead to morecompact systems, more similar to the current PC configuration. The eccentricity of the binary is more efficiently damped for higher initial values of P₀, as shown in panel b in Fig. 10. Panels c and d show the evolution of the angular velocities of Pluto and Charon, Ω_P and Ω_C (in units of the mean motion n_C), respectively;it is clear that the system reaches the double synchronous rotation more rapidly for higher values of P₀.

The time evolution of the accreting system for the different values of the dissipation parameter, A, is presented in Fig. 11. The values of A are chosen from the range [8 − 21], while the initial rotation period of Pluto and the initial Charon eccentricity are fixed at P₀ = 0.149 d and e_C = 0.001, respectively.Charon’s eccentricity, shown in the top panel of Fig. 11, is damped more efficiently as A increases. The evolution of Ω_C (in units of the mean motion, n_C) is shown in the bottom panel; we note that the synchronous state of Charon is more quickly attained for higher values of the dissipation parameter, A. The variations in the semi-major axis of Charon and the rotation of Pluto, Ω_P, with the increasing A values are insignificant and are not shown in Fig. 11.

The tests described above allowed us to define the parameter values of P₀ = 0.149 d and A for the interval [8 − 21], which are compatible with the current configuration of the PC binary. The tests, done with the oligarchic mass accretion, indicate an initial value of P₀ = 0.157 d; however, the corresponding simulations have shown results that are qualitatively similar to those obtained with the logarithmic mass growth model.

Fig. 9 Time evolution of Charon’s physical properties during the mass accretion: the mass, the radius, and the zonal harmonic, J2, from left to right. Colours represent the results obtained from the different accretion models. The dashed horizontal line (right panel) presents the constant J2 obtained in Correia (2020).

Fig. 10 Tidal evolution of the PC binary, for A = 10 and a logarithmic mass accretion of Charon during the first 104 yr, starting with the initial eccentricity eC = 0.001. Panels show the orbital semi-major axis, aC, in units of RP (panel a), the orbital eccentricity, eC (panel b), and the spin angular velocities of Pluto, ΩP (panel c) and Charon, ΩC (panel d) in units of the mean motion, nC. Each colour represents a different initial P0, in units of days.

6.2 PC binary and a small moon

In this section the model is expanded with the additional of a small moon with the mass and orbital elements of either Styx or Nix. As described in the previous section, the behaviour of these two satellites is significantly affected by the passages through the strong 4/1 and 5/1 MMRs. We explored the parameter space of the problem defined by A and m_i, where m_i is the mass of the small moon. For the mass, we chose the extreme values of the masses of Styx and Nix, given by the uncertainties reported in Brozović et al. (2015) and Showalter & Hamilton (2015). We set the initial eccentricity of the small moon at e = 0.001, the semi-major axis at a = 44413 km (a ~ 37.4 R_P) and a = 50690 km (a ~ 42.7 R_P) for the Styx and Nyx clones, respectively, and the angular elements equal to zero. As shown below, these initial values of a were chosen to obtain the current positions of the small moons. All simulations in this section were done with the accretion time of 10⁴ yr.

We started analysing the behaviour of the Styx clones during the tidal evolution of a growing Charon. We considered first a Styx mass of 1.49× 10¹⁵ kg, which corresponds to the 1σ estimate given in Brozović et al. (2015). Figure 12 shows the time evolution of the Styx orbital elements obtained for the values of the dissipation coefficient, A, ranging from 11 to 20; the elements are the semi-major axis (top panel) and the eccentricity (middle panel). We note that the test moon smoothly migrates inwards during the first ~ 10⁴ yr, a period during which Charon’s mass is increasing. After the satellite reaches its current position (continuous horizontal line), it evolves with a nearly constant semi-major axis up to an encounter with some important MMRs with Charon.

The interaction of the moon with the MMRs is generally accompanied by the excitation of the moon’s eccentricity. We can observe this phenomenon in the middle panel of Fig. 12, which shows the time evolution of the moon’s eccentricity during the migration of a growing Charon. When Charon’s mass grows to 90% of its current value (in about 3 × 10³ yr), the moon’s eccentricity is excited due to the approximation to the resonance domain, in this case, of the 9/1 MMR. From that event, the eccentricity can suffer consecutive excitations during the passages of the moon through the resonances, even being ejected from the system, depending on the value of A. The currentconfiguration of Styx is acquired for around 30% of the tested values of A, namely 16, 17, and 18. As observed in the bottom panel of Fig. 12, there are two ejections due to the 5/1 MMR and four due to the 4/1 MMR, and one is associated with the 7/2 MMR.

We repeated the simulations described above but this time analysing the behaviour of the Nix clones during the tidal expansion of an accreting Charon. We considered a value for Nix’s mass of 1.5 × 10¹⁶ kg, which is smaller than the mass range obtained in the orbital fits reported in Brozović et al. (2015) and Showalter & Hamilton (2015). The results obtained for Nix are shown in Fig. 13, which is similar to Fig. 12, but in this case we observe fewer ejections: one from the 6/1 MMR and another from the 4/1 MMR (see the bottom panel). We can observe in the middle panel of Fig. 13 that the particles experience consecutive eccentricity excitations during the passages through the MMRs, resulting in the higher dispersion of the final values when compared to the case of the Styx clones.

We redid the simulations described above for the Styx and Nix clones, considering now larger masses of 1.5 × 10¹⁶ kg and 4.5 × 10¹⁶ kg, respectively.We could not observe significant differences when comparing the new results to those obtained for the smaller masses of the test moons. Thus, we conclude that, in the range of values considered in this paper, the masses of the moons are not an important factor for achieving favourable configurations of the satellite system. On the other hand, the simulations conducted with the values of the dissipation coefficient, A, in the range 11− 20 did not allow definitive conclusions to be drawn regarding the influence of this parameter on the moon’s behaviour. In the next section we return to explore this issue.

Fig. 11 Time evolution of the PC binary, eccentricity (top panel) and ΩC∕n (bottom panel), for the fixed initial P0 = 0.149 d and differentvalues of A, which are represented with different colours.

Fig. 12 Time evolution of the orbital elements of Styx, with the mass mS = 1.49 × 1015 kg, for the different dissipation values of A. The top panel shows Styx’s semi-major axis, the middle panel shows the eccentricity, and the bottom panel shows nC ∕n. The continuous horizontal line in the top panel indicates Styx’s current position, at a∕Rp = 35.7, while the dashed lines in the bottom panel identify the 3/1, 4/1, and 5/1 MMRs.

6.3 Dependence on the accretion time

The accretion time of Charon, t_acc, is an important parameter since it affects the time evolution of Charon’s orbital elements, particularly Charon’s eccentricity. As shown in Sect. 4, the ejection of a small circumbinary moon during the passage through one of the main MMRs is most likely for higher Charon eccentricity values. On the other hand, the rate of Charon’s mass accretion determines the damping rate of its eccentricity. Indeed, Fig. 14 shows the evolution of e_C during the tidal expansion of Charon’s orbit. We use different colours to identify the tested values of the accretion time, namely 10², 10³, 10⁴, and 10⁵ yr; the curve corresponding to the simulation with the ‘no-growing’ Charon, starting with its current mass, is also shown in Fig. 14. All simulations were performed with the initial value of e_C = 0.001 and the A value fixed at 16. It is worth cautioning that the accretion rate determines the initial value of Pluto’s rotation period, P₀; for the tested times of 10² yr, 10³ yr, 10⁴ yr, and 10⁵ yr, we obtain optimal values of P₀ = 0.1400 d, P₀ = 0.1437 d, P₀ = 0.149 d, and P₀ = 0.1548 d, respectively (see Sect. 6.1 for details).

Figure 14 shows that the damping of Charon’s eccentricity is more efficient for an increasing t_acc. Indeed, we can clearly observe that Charon’s eccentricity systematically increases with decreasing accretion time and reaches its maximum value in the case of the non-accreting process. The vertical dashed lines mark the main MMRs and allow us to evaluate e_C when one of these MMRs crosses the current position of Styx. Hence, during the passages of the small moons through the MMRs, their motions will be more excited by Charon accreting the mass more rapidly. This fact would explain why we could not obtain successful results by simulating the evolution of the PC system without considering the accretion of Charon’s mass (see Sect. 5.6).

As shown in the previous sections, another important parameter of the problem is the dissipation coefficient, A. In order to link information on the dissipation parameter, A, to information on the accretion times, we prepared the parametric A–t_acc plane, which is shown in Fig. 15. For a given point on the plane, we show the value of e_C (in logarithmic scale) at the instant at which the 5/1 MMR, originated by a migrating Charon, crosses the current position of Styx. For this, we use the colour scale palette shown at the top of the figure. We note that the zero value on the y axis corresponds to the non-accreting Charon mass simulations. The red domains are associated with the highest e_C values, while the dark blue domains are related to the lowest values. Observing the plane, we can conclude that the efficient damping of e_C during the tidal expansion occurs for higher diffusion parameter values and the slower rate of accretion that favours the survival of the small satellites.

Fig. 13 Same as Fig. 12 except for Nix, with the mass mN = 1.5 × 1016 kg at the position a∕Rp = 40.98 (continuous horizontal line).

Fig. 14 Eccentricity of a growing Charon during the tidal expansion of Charon’s orbit. Different colours are used to present the different values of the accretion time. The vertical dashed lines show the nominal locations of the main MMRs originated by a migrating Charon at Styx’s current position; each MMR is identified by the corresponding label.

Fig. 15 Parametric A–tacc plane showing, in colour scale, the values of Charon’s eccentricity, eC (in logarithmic scale) at the instant when the 5/1 MMR crosses Styx’s current location.

6.4 Ring of particles in the accreting PC binary

In this section we consider the behaviour of the circumbinary thin disc of 500 particles during the tidal expansion of the accreting PC binary. Following Kenyon & Bromley (2019b), we first chose the disc, which is initially extended from 34 R_p to 60 R_p, with the eccentricities of the particles uniformly distributed between 0.5 × 10⁻³ and 1.5 × 10⁻³, no inclinations with respect to the orbital plane of the PC binary, and randomly chosen angular elements. We worked with particles with masses of either 1.49 × 10¹⁵ kg or 4.5 × 10¹⁶ kg; these values are still smaller than the mass of the PC binary; thus, the particles’ interactions can be neglected. Also, we selected random values of the dissipation parameter, A, in the interval [14−20] because these values favour the survival of the Styx and Nix clones, as described in the previous section. Finally, the time of the accretion of Charon’s mass is 10⁴ yr.

Figure 16 shows the results obtained for the disc of particles with masses of 4.5 × 10¹⁶ kg as a function of the initial positions of the particles. The top panel shows the survival times of the test particles. In the regions beyond the Nix position, almost all particles survive the period of the tidal expansion of the accreting PC binary, while the region between Styx and Nix presents more instabilities (black dots in the top panel). In the middle panel of Fig. 16, we plot the final positions of the particles (orange) and Charon (blue). The final positions of the particles are generally closer to the binary than the initial conditions (the points lie slightly below the line a_i = a_f, plotted as a continue blue line); this is a consequence of the increasing mass of the binary.

The bottom panel of Fig. 16 shows the final eccentricities of the PC binary (blue dots) and the particles (orange dots). For the particles with an initial a chosen in the range 34 R_p−45 R_p, the eccentricities are slightly higher than those from the range 45 R_p−60 R_p. The shaded region represents the extreme values of the eccentricities of the known moons, according to Table 1; the output of several simulations ends in this shaded region. The PC binary is almost circular, with e_C ~ 10⁻⁵, which is compatible with its orbital fits (Brozović et al. 2015; Showalter & Hamilton 2015). We do not find any preference for A values in the range [14−20] in the region, which is interior to the Nix position; around 35% of particles survive in this region. We repeated the experiment choosing A values in the range [9−13] and found that 99% of the clones located inside the Nix’s orbit were ejected.

Several works have shown that the particles that start with low initial inclinations could be captured in a secular resonance, with a consequent excitation of their inclinations (e.g. Ćuk et al. 2020). We ran some experiments that considered a ring of particles similar to that described above but now inclined with an angle whose value was uniformly distributed in the range − 0.01° < i₂ < 0.01°. We find no qualitative differences with the results presented in Fig. 16. Figure 17 shows the result obtained for the inclinations of the particles after 3 × 10⁶ yr of the tidal evolution of the system. We observe some excitation of the inclination in the domain surrounding Styx and Nix, which is expected during the passages through the strong 4/1 MMR. These results are in agreement with the orbital fits reported in Showalter &Hamilton (2015).

The results obtained for the disc of the particles with a mass of 1.49× 10¹⁵ kg are similar to those described above. Thus, we can conclude that the simple model, which considers the accreting of Charon’s mass during the tidal expansion of the PC binary, is effective to explain the current locations of the small moons of Pluto without having to necessarily introduce the resonant transport mechanism.

Fig. 16 Evolution of the ring of particles, each with a mass of 4.5 × 1016 kg, for an A in the range [14–21]. Top panel: survival times for particles. Middle panel: final semi-major axis of the small moons and Charon. Bottom panel: final eccentricity of Charon and the particles. The mass growth of Charon is set to 104 yr. Blue is used to indicate orbital parameters of Charon and orange for the particles. Orange dots represent survivors (the particles that remain close to initial configurations during the integration time of 3 × 106 yr). The black points in the top panel identify those particles that escape while the PC binary evolves. The shaded grey region indicates the range of values according to diverse authors (see Table 1).

Fig. 17 Distribution of inclination for the disc of particles after the dual synchronisation of the PC system with mass growth. Theorbits that do not survive the entire integration time are identified with black dots. The shaded grey region indicates the range of values according to diverse authors (see Table 1).

7 Conclusions

We have analysed the stability and dynamical evolution of the small moons in the PC system, namely, Styx, Nix, Kerberos, and Hydra, during the tidal orbital expansion of Charon. We study the current architecture of the system and investigate how the tidal expansion of the PC binary system affects the orbital evolution of the small moons, extending previous studies.

The resonant structure around the PC binary exhibits N∕1 and N∕2 resonances, although most of them are weak in the low-eccentricity and low-inclination regimes. The 3/1 MMR is the strongest circumplanetary resonance, even for a low eccentricity and inclination of the moons (see Fig. 1). The current orbital configuration of the system does not show dynamical preferences of stable motion within the N/1 resonances despite the four small moons currently being located close to the 3/1, 4/1, 5/1, and 6/1 MMRs.

Using dynamical maps in the a_C–e_C plane, we reconstructed the tidal evolution of the PC system with a small moon initially placed at its current position (see Fig. 2). We adopted a tidal model parameterised by the constant time lag, Δt. When the PC system evolves assuming aninitial value of e_C ~ 0.001, the resonances are weaker than when the PC system evolves from an initial value of e_C ~ 0.1 or 0.2. Thus, the evolution from initially quasi-circular orbits avoids chaotic regions associated with N∕1 resonances, constraining the initial eccentricity of Charon. In addition, for small moons located near their current position, their semi-major axes do not secularly change unless they are captured in an MMR due to the tidal expansion of the PC binary.

The results of our numerical simulations indicate that it is possible to reproduce the current orbits of the four small moons of Pluto during the early tidal orbital expansion of Charon if we ignore the contribution of zonal harmonic, J₂, and for a specific combination of initial conditions, such as Charon’s eccentricity, semi-major axis, and dissipation parameter, A. Hence, under the assumption of some specific values for the physical parameters and for the initial semi-major axis and eccentricity of Charon, these results point towards some specific scenarios (see Table 3). For example, for Styx and Nix, an initial almost circular orbit of Charon is needed in order to obtain a final stable orbit of the small satellites. In addition, a high dissipation in Pluto (Δt_P) allows for final orbital configurations of Styx that are in good agreement with their current orbits.

In the case of accreting Charon mass and considering an evolving J₂ with dissipation coefficient A ≳ 14, we find that the small moons can survive the rapid evolution of the PC system and that the unlikely resonant transport of particles is no longer necessary. The ad hoc mass evolution of Charon might allow the moons to cross some resonances and not be trapped in specific positions of resonant motion. We also find that the eccentricity of Charon, e_C, when Charon crosses the 5/1 MMR with the actual position of Styx might be used to estimate the accretion time or dissipation coefficient in order to obtain survival moons as the system tidally evolves. Our numerical simulations that consider a ring of particles between 34 R_C and 60 R_C show that a significant amount of particles remain in stable motion, with final orbits similar to the present small moons. Thus, the detection of additional small moons in future missions to the Pluto system should not be ruled out.

Further work is necessary to constrain the mass accretion in a more self-consistent scenario (e.g. Kokubo et al. 2000; Kokubo & Ida 1998). In addition, our model could be improved by considering mutual perturbations between the small moons(Showalter & Hamilton 2015; Lee & Peale 2006), by using a more realistic model for tidal dissipation within icy bodies (see Ferraz-Mello 2013; Folonier et al. 2018), or by considering the obliquity of Charon and perturbations from the Sun and Neptune (see e.g. Correia 2020).

Acknowledgements

N-body computations were performed at Mulatona Cluster from CCAD - UNC, which is part of SNCAD - MinCyT, Argentina. The constructionof the dynamical maps has made use of the facilities of the Laboratory of Astroinformatics (IAG/USP, NAT/Unicsul), whose purchase was made possible by FAPESP (grant 2009/54006-4) and the INCT-A. T.A.M. was supported by the Brazilian agencies Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), and São Paulo Research Foundation (FAPESP, grant 2016/13750-6). We acknowledge the program Professor Visitante do Exterior (PVE) from Pró-Reitoria de Pós-Graduação - Universidade de São Paulo. Edital PRPG 02/2019. The research done in this project made use of the Astroquery (Ginsburg et al. 2019), a community-developed core Python package for Astronomy (Astropy Collaboration 2013, 2018). The data presented in this paper are original from the authors, and are available upon reasonable request.

Appendix A Dynamical evolution of small moons

Here we provide complementary information to construct Table 1. In Fig. A.1 we show N-body integration of a three-body problem for a time span of 40 yr (grey lines), together with ephemerides taken from JPL Horizons. Amplitude oscillations of the orbital elements are resumed in Table 1 and agree with the colour scale presented in Fig 1. It is interesting to look first at the JPL ephemerides(blue lines in the figure). The left panels indicate the evolution of the semi-major axis, which presents amplitude of oscillations, from top to bottom, of 1761 km, 1086 km, 637 km, and 443 km (roughly Δa = 1.5R_P, 0.9R_P, 0.5R_P, and 0.4R_P) for Styx, Nix, Kerberos, and Hydra, respectively. We note that orbital elements are given with respect to the barycentre of the PC system. Osculating eccentricities present higher variations for closer moons and lower amplitudes for exterior moons (0.034, 0.025, 0.018, and 0.016, respectively). These variations are of the same order as the variations presented in Fig. 2 in Woo & Lee (2018) and Fig. 1 in Woo & Lee (2020), although the authors used the eccentricity for Charon from the JPL solution PLU043, which is older and an order of magnitude lower than the current determinations (solution PLU058). Finally, inclination evolution (right panels of Fig. A.1) shows oscillations around the midplane of the PC binary (denoted by an orange line, which is nearly coincident between N-body integrations and JPL query data). It is evident that the small moons can be either above or below the plane with different phases, and this inspired us in Sects. 5.6 and 6 to choose the PC orbit as reference plane, thus, for example, randomly defining inclinations around zero degrees.

Our N-body integrations (grey lines in Fig.A.1) show almost the same amplitudes as the JPL query (see also Table 1). However, there is an evident shift in the evolution of the semi-major axis of Nix and Hydra, although Δa is almost the same. We believe that the differences in the evolution of the semi-major axis in our simulations and JPL data could be due to the fact that the more sophisticated JPL model considers a significant number of perturbations, including the Sun and Neptune.

Figure A.2 shows four dynamical maps in the plane semi-major axis and mean anomaly, constructed using the mean exponential growth of nearby orbits (MEGNO) colour scale (e.g. MEGNO=⟨Y ⟩; Cincotta & Simó 2000, Y^* = log₁₀(⟨Y ⟩− 2)). MEGNO allows us to identify the separatrices of each resonance as chaotic (reddish colour of each panel). The MMR 3/1 is evidently more asymmetric than the 4/1, 5/1, and 6/1 MMRs because of its proximity to the PC binary. For example, for the Styx plane the separatrices show that the initial mean anomaly modifies the location of the resonance in ~1000 kms. Also, the MMRs 3/1 and 4/1 are wider than 5/1 and 6/1 MMRs because the moons experience greater excursions in eccentricity. These maps also help to understand that stability maps in planes a − e strongly depend on the chosen value of mean anomaly (in Fig. 1 these correspond to those from Styx, M_S = 18.79°).

Fig. A.1 Orbital evolution of Styx, Nix, Hydra, and Kerberos around the PC binary, from top to bottom, respectively. We show the evolution of the semi-major axis (left panels), eccentricity (middle panels), and inclination (right panels) for 40 yr. Blue lines correspond to data from JPL and grey lines to our N-body integrations. Orange lines correspond to the integrated inclination of Charon.

Fig. A.2 Stability maps on the a–M plane using the MEGNO colour scale for the four small moons. All other orbital parameters are set to the JPL values from Table 1. The current positions of the moons are shown by red crosses. The locations of the main MMRs are indicated. Small deviations from initial conditions within the estimated errors from PLU043 can place the moons inside the MMR regime.

Appendix B Survival times of the disc described by simplified models

Here we present complementary information to that discussed in Sect. 5.6. Now we investigate the survival times of the disc of small moons, applying the model that uses fixed values of the zonal harmonic of Pluto’s oblateness, either J₂ = 0 or J₂ = 6 × 10⁻⁵ (the latter value was determined from the shape of Pluto in Correia 2020). We remember that our model consists of the tidally evolving PC binary and a disc of small moons extending from 10R_P to 55R_P (for more details, see Sect. 5.6).

Figure B.1 shows the case when the effect of the zonal harmonic is not considered (J₂ = 0). For Charon starting on a nearly circular orbit, the survival times of the test moons on the planar orbits (blue crosses) monotonously increase (in log scale) with the increasing semi-major axis, reaching 10³ yr for particles in the internal disc (10R_P < a < 20R_P). The particles at Styx’s position andbeyond it (a > 35R_P) show survival times greater than 3 × 10⁶ yr. Moreover, the numerical integrations of the individual orbits show that the moon’s motions are stable for more than 1 × 10⁸ orbital periods of the PC binary.

When we consider a thin disc of small moons with inclinations uniformly distributed between − 0.01° and 0.01° (Fig. B.1, orange dots), we observe that survival times are similar to those obtained for a coplanar disc. A difference is observed essentially for the moons initially located at Charon’s position. We analyse the final inclinations of the particles in the thin disc in Fig. B.2 and observe that the orbits of the surviving moons around 20 R_P reach inclinations of up to 25°. The inclinations of the particles located inside the limit given by the HWC are dispersed around the mean value of 5°. Also, some portion of the test particles near the Nix position is excited in inclination up to 5° but survive the time span of the simulations.

Figure B.1 also presents the survival times obtained for the coplanar disc, but with higher initial values of Charon eccentricities, namely, e_C = 0.1 and e_C = 0.2 (green crosses and red dots, respectively). These values put constraints on the initial rotational periods of Pluto, now P₀ = 0.13791 d and P₀ = 0.13593 d, respectively. For e_C = 0.1, only particles beyond the Kerberos orbit survive for 10⁶ yr. Setting the initial values of e_C = 0.2 causes a faster depletion of the initial disc, < 10⁵ yr. We do not find differences when we initially choose aligned or anti-aligned longitudes for the particle’s pericentre and the pericentre of Charon’s orbit.

Figure B.1 can be directly compared with Fig. 10 in Woo & Lee (2018), where almost all the particles in the range a > 30R_P survive. However, their model was constructed for initial e_C = 0, massless particles, and A = 40. We can observe in our simulations that moderate values (A = 10) allow the survival of particles. However, Woo & Lee (2018) explore neither cases with an almost initially circular Charon (e_C ~ 0) nor cases with inclined particles.

The results obtained when setting a fixed value of J₂ are presented in Fig. B.3. We remember that the static value of J₂ = 6 × 10⁻⁵ was derived from the shape of Pluto in Correia (2020). We can observe that, for all initial eccentricities of Charon, the particles evolve almost in the same way as shown in Fig. B.1.

Fig. B.1 Survival times of the test particles from different discs, all extending from 10RP < a < 55RP, for the dissipation parameter A = 10. We use different symbols and colours for initially coplanar and inclined discs as well as for different initial eccentricities of Charon’s orbit. We do not consider J2 effects, and the Charon mass is conserved during the integrations (see the text for more details).

Fig. B.2 Distribution of the final inclinations of the thin disc of particles shown in Fig. B.1. The particles ejected fromthe PC system are identified by black dots.

Fig. B.3 Same as in Fig. B.1, except considering the effect of a fixed zonal harmonic, J2 = 6 × 10−5.

References

All Tables

All Figures

Fig. 1 Dynamical maps on the a–e plane (left panel) and a–I plane (right panel) of the moon’s initial conditions calculated for the PC binary and a small moon of Styx’s mass with initial Styx orbital elements from JPL in Table 1. The current positions of the moons are shown by red crosses. The locations of the main MMRs are indicated, and the grey line in the right panel indicates the current inclination of Charon. The colour scale varies with Δe and Δi, increasing from blue (regular motion) to red (chaotic motion), as shown in the colour bars on the top of the figure.

In the text

Fig. 2 Dynamical maps of Styx, Nix, Kerberos, and Hydra on the aC–eC plane of Charon’s initial conditions. On each map, black curves are Charon’s tidal trajectories, obtained for initial eC values of 0.01, 0.1, and 0.2; the initial conditions of the corresponding moon are fixed at those given in JPL from Table 1. The blue tones represent regular motions, while the brick tones correspond to increasing instabilities and chaotic motion; the colour bar on the top shows Δe. The locations of the main MMRs are indicated.

In the text

Fig. 3 Two typical examples of the Styx orbital evolution resulting in either instability (left panel) or the final configuration, which is significantly different from its current orbit (right panel). Left panel: initial conditions used are aC ∕RP= 5, eC = 0.001, ΩC ∕nC = 2, A = 10, and ΔtP = 600 s. Styx is ultimately ejected from the system when the system crosses the 4/1 MMR between Charon and Styx (not shown), resulting in a strong increase in Styx’s eccentricity. Right panel: initial conditions used are aC ∕RP= 4, eC = 0.001, ΩC ∕nC = 2, A = 10, and ΔtP = 6 s. The orbital eccentricity is excited, due to the crossing of the 4/1 and 7/2 MMRs (not shown), to attain the mean value ≃ 0.1, which is substantially larger than the current Styx eccentricity.

In the text

Fig. 4 Examples of initial conditions for Styx. Top row: time evolution of Styx’s semi-major axis (left) and eccentricity (right) during the tidal expansion of the PC binary. The initial conditions used in the simulations are aC ∕RP= 4, es C = 0.001, ΩC ∕nC = 2, A = 100, and ΔtP = 600 s. The final values of the orbital elements remain close to the current ones. Bottom row: time evolution of the mean-motion ratio (left) and Δϖ (right). We note that the mean value of nC∕n is larger than 3. The secular angle, Δϖ, oscillates around 180°.

In the text

Fig. 5 Examples of initial conditions for Nix. Top row: time variations in Nix’s semi-major axis (left) and eccentricity (right) for two different initial values of aC. The initial conditions are aC∕RP = 3 and aC ∕RP = 7, eC = 0.001, ΩC ∕nC = 2, A = 10, and ΔtP = 600 s. For aC ∕RP = 3, the eccentricity of Nix is excited to a mean value close to 0.16. Bottom row: time evolution of the mean-motion ratio (left) and two angles for Nix (right), namely, the secular angle, Δϖ, and the 4/1 MMR critical angle ϕ4∕1 = 4λ − λC − 3ϖC. The mean value of nC∕n is slightly smaller than 4 for aC∕RP = 7, whereas ϕ4∕1 and Δϖ librate around 180° and 0°, respectively, for aC∕RP = 3.

In the text

Fig. 6 Examples of initial conditions for Kerberos. Top row: orbital evolution of Kerberos for two initial conditions of eC, namely, 0.001 and 0.1. The other parameters are aC∕R = 4, ΩC ∕nC = 2, A = 10, and ΔtP = 600 s. These simulations are successful according to the adopted criteria (see Sect. 5). Bottom row: time variations in the mean-motion ratio (left) and ϕ5∕1 = 5λ − λC − ϖC − 3ϖ (right). We show the angle for the specific case of eC = 0.001. The mean value of nC∕n is larger than 5.

In the text

Fig. 7 Examples of initial conditions for Hydra. Top row: orbital variation in the semi-major axis (left) and eccentricity (right) of Hydra according to the adopted criteria. The initial conditions are aC ∕RP= 4, eC = 0.001, ΩC ∕nC = 2, A = 10, and ΔtP = 6 s. Bottom row:time variation in the mean-motion ratio (left) and ϕ6∕1 = 6λ − λC − ϖC − 4ϖ (right) for two values of ΔtP, namely, 6 s and 60 s. We note that, for ΔtP = 60 s, ϕ6∕1 librates witha large amplitude around 0°. For ΔtP = 6 s, all critical angles associated with the 6/1 MMR circulate.

In the text

	Fig. 8 Survival times of the small particles, for the different discs extending from 10 RP to 55 RP, during the tidal expansion of the PC binary, with the dissipation parameter A = 10, under the effect of the evolving zonal harmonic, J2. The colour symbols are used to identify the particles from the different discs (see text).
In the text

	Fig. 9 Time evolution of Charon’s physical properties during the mass accretion: the mass, the radius, and the zonal harmonic, J2, from left to right. Colours represent the results obtained from the different accretion models. The dashed horizontal line (right panel) presents the constant J2 obtained in Correia (2020).
In the text

Fig. 10 Tidal evolution of the PC binary, for A = 10 and a logarithmic mass accretion of Charon during the first 104 yr, starting with the initial eccentricity eC = 0.001. Panels show the orbital semi-major axis, aC, in units of RP (panel a), the orbital eccentricity, eC (panel b), and the spin angular velocities of Pluto, ΩP (panel c) and Charon, ΩC (panel d) in units of the mean motion, nC. Each colour represents a different initial P0, in units of days.

In the text

	Fig. 11 Time evolution of the PC binary, eccentricity (top panel) and ΩC∕n (bottom panel), for the fixed initial P0 = 0.149 d and differentvalues of A, which are represented with different colours.
In the text

Fig. 12 Time evolution of the orbital elements of Styx, with the mass mS = 1.49 × 1015 kg, for the different dissipation values of A. The top panel shows Styx’s semi-major axis, the middle panel shows the eccentricity, and the bottom panel shows nC ∕n. The continuous horizontal line in the top panel indicates Styx’s current position, at a∕Rp = 35.7, while the dashed lines in the bottom panel identify the 3/1, 4/1, and 5/1 MMRs.

In the text

	Fig. 13 Same as Fig. 12 except for Nix, with the mass mN = 1.5 × 1016 kg at the position a∕Rp = 40.98 (continuous horizontal line).
In the text

	Fig. 14 Eccentricity of a growing Charon during the tidal expansion of Charon’s orbit. Different colours are used to present the different values of the accretion time. The vertical dashed lines show the nominal locations of the main MMRs originated by a migrating Charon at Styx’s current position; each MMR is identified by the corresponding label.
In the text

	Fig. 15 Parametric A–tacc plane showing, in colour scale, the values of Charon’s eccentricity, eC (in logarithmic scale) at the instant when the 5/1 MMR crosses Styx’s current location.
In the text

Fig. 16 Evolution of the ring of particles, each with a mass of 4.5 × 1016 kg, for an A in the range [14–21]. Top panel: survival times for particles. Middle panel: final semi-major axis of the small moons and Charon. Bottom panel: final eccentricity of Charon and the particles. The mass growth of Charon is set to 104 yr. Blue is used to indicate orbital parameters of Charon and orange for the particles. Orange dots represent survivors (the particles that remain close to initial configurations during the integration time of 3 × 106 yr). The black points in the top panel identify those particles that escape while the PC binary evolves. The shaded grey region indicates the range of values according to diverse authors (see Table 1).

In the text

	Fig. 17 Distribution of inclination for the disc of particles after the dual synchronisation of the PC system with mass growth. Theorbits that do not survive the entire integration time are identified with black dots. The shaded grey region indicates the range of values according to diverse authors (see Table 1).
In the text

	Fig. A.1 Orbital evolution of Styx, Nix, Hydra, and Kerberos around the PC binary, from top to bottom, respectively. We show the evolution of the semi-major axis (left panels), eccentricity (middle panels), and inclination (right panels) for 40 yr. Blue lines correspond to data from JPL and grey lines to our N-body integrations. Orange lines correspond to the integrated inclination of Charon.
In the text

	Fig. A.2 Stability maps on the a–M plane using the MEGNO colour scale for the four small moons. All other orbital parameters are set to the JPL values from Table 1. The current positions of the moons are shown by red crosses. The locations of the main MMRs are indicated. Small deviations from initial conditions within the estimated errors from PLU043 can place the moons inside the MMR regime.
In the text

Fig. B.1 Survival times of the test particles from different discs, all extending from 10RP < a < 55RP, for the dissipation parameter A = 10. We use different symbols and colours for initially coplanar and inclined discs as well as for different initial eccentricities of Charon’s orbit. We do not consider J2 effects, and the Charon mass is conserved during the integrations (see the text for more details).

In the text

	Fig. B.2 Distribution of the final inclinations of the thin disc of particles shown in Fig. B.1. The particles ejected fromthe PC system are identified by black dots.
In the text

	Fig. B.3 Same as in Fig. B.1, except considering the effect of a fixed zonal harmonic, J2 = 6 × 10−5.
In the text