Issue 
A&A
Volume 670, February 2023



Article Number  A152  
Number of page(s)  9  
Section  Celestial mechanics and astrometry  
DOI  https://doi.org/10.1051/00046361/202244717  
Published online  21 February 2023 
Orbital analysis of the PlutoCharon moon system’s mutual interactions and forced frequencies
Department of Physics, University of Patras,
Patras,
26504 Rio, Greece
email: dgakis@upnet.gr; kngourg@upatras.gr
Received:
9
August
2022
Accepted:
4
January
2023
Context. The orbits of the four small moons in the PlutoCharon system, Styx, Nix, Kerberos, and Hydra, are circumbinary, as Pluto and Charon form a binary dwarf planet. Consequently, the orbit of each moon is characterized by a number of frequencies, arising from the central binary and the mutual gravitational interactions.
Aims. In this work, we identify the most prominent of these forced frequencies using fast Fourier transforms.
Methods. Two methods were implemented, a semianalytic and a numerical one, and comparisons are made.
Results. The results indicate that as a first approximation, moon orbits may well be modeled as the superposition of a series of inevitable oscillations induced by Pluto and Charon, deviating from circular orbits, even if the eccentricity is set to zero. Moreover, the mutual gravitational effects are significant in their longterm evolution, especially for the lighter moons Styx and Kerberos, activating modes that dominate the lowfrequency region of the power spectrum. This becomes evident through the comparison of simulations where only one moon is included along with the binary dwarf planet and simulations of the entire sixbody system. These modes become noticeable over long integration times and may affect the orbits of the lighter moons of the system.
Key words: planets and satellites: dynamical evolution and stability / Kuiper belt objects: individual: PlutoCharon / celestial mechanics
© The Authors 2023
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
1 Introduction
Pluto’s moon system is a dynamical treasure. As the mass ratio of the dwarf planet Pluto to its largest moon Charon is 8:1 (Stern et al. 2015), they are in fact a binary dwarf planet. Along with the central binary, with at present four known moons orbiting the system’s center of mass, namely Styx, Nix, Kerberos, and Hydra, this structure is valuable for studying circumbinary orbits in depth. Thus, studying the motions of these small moons is of particular interest, as the potential arising from Pluto and Charon forces them into orbits that deviate significantly from the standard elliptical ones.
Circumbinary orbits differ greatly from those described by Keplerian orbital elements. Lee & Peale (2006) developed a theoretical solution to modeling orbits around a zeroeccentricity binary system. Their theory, which holds for point masses on circumbinary coplanar orbits, yields that a circumbinary orbit is the superposition of a circular orbit around the center of mass, an epicyclic motion caused by the binary and a vertical component. Leung & Lee (2013) generalized this theory to include eccentric orbits of the central binary as well.
Another study, by Bromley & Kenyon (2021), revisited the above theory and provided quantitative tools to apply it in practice. A “most circular” circumbinary orbit is defined, corresponding to a circular orbit around a single mass. Deviations from the most circular orbit are quantified using the free eccentricity (e_{free}). They tested outcomes for the eccentricity damping of tracer particles in the PlutoCharon system, along with other extrasolar planetary systems, achieving their objective satisfactorily. Nevertheless, it is acknowledged that more precise techniques are required to analyze the actual moon orbits instead of test particles. In a followup study, Kenyon & Bromley (2022) further examined and set improved constraints on the dynamical behavior and masses of the smallest moons by performing an array of simulations.
Woo & Lee (2020) used fast Fourier transforms (FFTs) on numerical simulations of the dynamical system to estimate the exact values of the amplitudes and frequencies that outline the peculiar orbits of the moons of Pluto and Charon. Although they confirmed the accuracy of this method, they found significant deviations from the expected orbit of Styx. By adopting the Hamiltonian approach of Lithwick & Wu (2008), they propose that at least a part of these deviations can be explained by the 3:1 mean motion resonance.
Some other theoretical solutions for circumbinary orbits exist as well. For instance, Georgakarakos & Eggl (2015) used perturbations of the Runge–Lenz vector to study the shortterm of the evolution of loweccentricity orbits in a hierarchical triple system. Another example appears in the work of Sutherland & Kratter (2019), who proposed using empirical geometric orbital elements to search for active resonances in orbits around a binary system. These studies do result in compatible solutions to those based on the Lee & Peale (2006) theory, so we do not discuss them further.
In a previous paper (Gakis & Gourgouliatos 2022), we examined the moon motions within the dynamical system of Pluto and Charon. Despite attempting to define orbits that are as close as possible to circular, moons, nonetheless appeared to deviate from such orbits, and the barycentric distance varied significantly. Our conclusion was that the timedependent, nonaxisymmetric potential of Pluto and Charon induces irregular patterns in the orbits. We also inspected that major discrepancies concerning Keplerian orbital elements between different studies so far do not primarily reflect limitations and inaccuracies in measurements, observations and calculations, but are in fact a result of the underlying specifications of the actual system.
This work is the second part of our analysis concerning the circumbinary orbits within the gravitational system of Pluto and Charon. Here, our goal is to give a more quantitative view on the dynamical specifications of the system, focusing on the effect of mutual interactions between the moons, in addition to the impact of the central binary that we have already studied. To this end, we utilize both semianalytic and numerical approaches, and infer the most prominent amplitudes and frequencies. Specifically, we apply FFTs to identify the exact frequencies of the many oscillations that moons perform. This way we quantify the impact of the mutual interactions between the moons that force additional frequencies in the orbits.
The structure of this paper is the following. In Sect. 2, we describe our calculations, which are both semianalytic (Sect. 2.1) and numerical (Sect. 2.2). Section 3 contains our main results. Our final conclusions are summarized in the last section (Sect. 4).
Orbital parameters for Pluto’s moon system.
2 Calculations
The orbital elements of the dynamical system, used in this work, are presented in Table 1. There are notable differences between the data set of Showalter & Hamilton (2015), and another prominent study, by Brozović et al. (2015). These deviations are most probably explained by the intrinsic behavior of the system, as illustrated by Gakis & Gourgouliatos (2022). In particular, observations at different time periods inevitably produced distinct outcomes because of the variations in the relative positions of the bodies.
2.1 Semianalytic approach
A short description of the model for circumbinary orbits, introduced by Lee & Peale (2006) and reconsidered by Bromley & Kenyon (2021), on which our semianalytic approach is principally based, is given below. Adopting cylindrical coordinates and assuming that the barycenter lies on the origin O (0, 0, 0), the potential caused by the central binary system at a point P (R, ϕ, z = 0) may be approximated by a cosine series: (1)
The coefficients Φ_{k} are related to the binary properties (masses M_{P} and M_{C} and separation a_{PC}) and the orbital radius of the moon, R_{S}. The orbital frequency of Pluto and Charon is given by (2)
Solving the equations of motion, the solution yields the following expressions for the position of a pointmass body initially located at P, as a function of time: (3) (4)
Here n_{syn} stands for the synodic frequency, i.e. n_{syn} = n_{PC} − n_{S} and i is the inclination of the moon orbit with respect to the PlutoCharon orbital plane. The moon’s mean motion n_{S}, the epicyclic frequency υ_{e} and the vertical frequency υ_{i} are defined as follows (see Appendix in Bromley & Kenyon 2021): (5) (6) (7)
where M = M_{P} + M_{C} and μ = (M_{P} · M_{C})/(M_{P} + M_{C}) are the total and reduced mass of the binary, respectively, and (the subscripts denote the respective objects).
The factor C_{k} represents the amplitudes of the oscillations: (8)
Assuming that the whole system’s barycenter coincides with the PlutoCharon’s barycenter, the orbital distance of a small moon would be: (9)
We assign the nominal eccentricities (Table 1) of each moon as e_{free}. Table 2 presents the calculated values of the major frequencies at which the moons oscillate. The central binary frequency is 0.1566 days^{−1} (according to Table 1).
The first major oscillatory frequencies, in (2π days)^{−1}, for all small moons.
2.2 Numerical approach
We approach the orbits numerically by implementing an nbody symplectic integrator in a Python 3.9.6 IDLE environment. The nbody simulation code utilizes the kickdrift technique to solve the differential equations representing gravitational interactions. The desired accuracy of the code is validated, since the total calculated energy of the system is kept constant.
Simulations of the 6body system in question, with different initial data sets each time, have been performed using the same nbody code by Gakis & Gourgouliatos (2022). We reexamine the basic situation of these simulations, in order to give a more thorough perspective on the concepts discussed in this paper. Specifically, initial conditions are taken from Brozović et al. (2015), where a table is provided (Table 8 therein) of measured 3D vectors of positions and velocities for every object. We analyzed the evolution of the system forward in time using this data set.
The numerical integration timestep is fixed to ∆t = 5000 s, which maintains computational times under manageable limits, and at the same time keeps uncertainties below 0.1%, as determined in Gakis & Gourgouliatos (2022). Besides, Kenyon & Bromley (2019a,b) propose at least ∆t ⪅ 13 500 s for reliable integrations. Timesteps like these are used by various studies on the orbits within the PlutoCharon system; for example the numerical calculations performed by Woo & Lee (2020) have a ∆t = 3000 s, whereas Lee & Peale (2006) use a larger timestep, ∆t = 10 000 s. Nevertheless, we also ran numerical tests with smaller timesteps, not identifying however any noticeable changes. The gravitational effect of the Sun is significant at distances of about ten times larger than Hydra’s average orbital distance (Michaely et al. 2017); therefore, we focus on a restricted 6body problem in our simulations, neglecting the Sun and other Solar System bodies.
3 Results
Several algorithms are implemented in order to study the orbits in several dynamical systems. In our analysis, we chose to adopt FFT to decompose the orbits in the PlutoCharon system. Although FFT may be less accurate than other methods, such as the Frequency Map Analysis (FMA; Laskar 1999), it still produces reliable outcomes significantly quickly, and hence is often applied to analyze circumbinary orbits (e.g. Woo & Lee 2020; Gakis & Gourgouliatos 2023). Our results indicate that the resolution provided by FFT is suitable to detect and separate the forced frequencies caused by the central binary and some of their harmonics, as well as the main trends of the reciprocal effects by each of the other moons. Besides, since our goal is to identify the frequencies of the moon’s oscillations rather than studying chaos in the system, which has been explored thoroughly in past studies (e.g. Kenyon & Bromley 2019a,b, 2022), the FFT method will suffice.
The general formula used to convert a sequence x[n] of length N into a new one y[k] using a Fourier transform is: (10)
More precisely, we convert a distance domain into a domain of frequencies. Fast Fourier transforms are performed using the Python scipy routine fft.
3.1 Central binary effects
At first, we apply FFT of r(t) for the outcomes of the semianalytic model (Fig. 1). The timestep adopted for the semianalytic calculations was the same with the nbody integration timestep, and the total duration was set to 10^{6} days. Unlike Woo & Lee (2020), here we examine the vertical motion, since our calculations include the proper inclinations. The most outstanding frequencies arising are those defined when computing Eqs. (3) and (4), as expected. The red vertical dotted line in each periodogram corresponds to the value of υ_{e}, whereas the gray lines correspond to the harmonics k(n_{PC} − n_{S}). The spikes vary in height, as anticipated by the factor C_{k}. In other words, the relative size of each peak would give us a comparison of the different amplitudes of each frequency.
As far as the vertical frequency υ_{i} is concerned, there also appears to be a minuscule peak, though not visible in the frequency spectra of Fig. 1. Having zoomed into the lowfrequency area of each spectrum and identifying a corresponding (barely visible) formation, we advocate that its apparent absence is not a problem of the frequency resolution nor with the wide frequency range. Instead, this is a result of the factor i R_{S} in Eq. (4) outlining the vertical frequency, which much smaller than the values of R_{S} in Eq. (3), which dominates in the configuration of the strength of each peak in the frequency spectrum (i < 1° for all moons).
There are also some other secondary spikes, not corresponding to the values of Table 2. They originate from the vertical component of the motion and the sinusoidal products deriving from Eq. (9). Namely, R^{2}(t) and z^{2}(t), along with the square root, give a number of harmonic cross terms that eventually result in frequencies of forms like υ_{e} ± k n_{syn} (and multiples), 2υ_{e}, 2υ_{i} and so on. In general, many combinations of υ_{e}, k n_{syn} and υ_{i} arise from Eq. (9), which are present in the periodograms of Fig. 1. Most of the secondary peaks have a quite small amplitude, which makes them only clearly visible when using a highly magnified image.
There are some distinctive differences once the nbody simulation is employed. The resulting power spectra are shown in Fig. 2. The system is let to evolve for 10^{4} and 10^{6} days. Again, in this figure, vertical lines show the main expected frequencies, as they have been computed in Table 2. The primary peaks are observed at these frequencies in this case as well. Nevertheless, a large number of other minor peaks are also visible. In fact, when we increase the simulated time, additional frequencies appear or, alternatively, the alreadypresent frequencies are enhanced. Furthermore, the increase of the total timespan reduces unwanted noise; the values of the power spectrum are diminished.
We notice that the vertical lines (i.e. the findings of the semianalytic model) do not exactly match with the peaks of the periodogram in Fig. 2. This is not largely visible on large scale (deviations scale to ~0.5%), but may be observed when zooming in on each individual spike. This distinction is evidence of the unavoidable inconsistency between the two approaches that we follow. Any differences could safely be attributed to approximations made in the semianalytic model, as justified in Sect. 2. For example, we neglect the nonlinear terms that definitely rule the motions of the moons. Styx is the most striking example and reveals the limitations of the Lee & Peale (2006) theory in distances close to the binary.
Yet, apart from the peaks dominating the region around zero, which we discuss later, the lowamplitude spikes appearing in between the vertical lines can clearly be understood within the semianalytic model. As we discussed earlier for Fig. 1, these peaks are the result of the multiple sinusoidal products of Eq. (9). Thus, they are not in principal caused by numerical errors, but are undeniably anticipated by the theoretical model.
Bromley & Kenyon (2021) and Woo & Lee (2020) found further formations in their power spectra, unable to be explained by the epicyclic theory. In the first paper, the authors found a residual signal at the epicyclic frequency, when simulating a mostcircular orbit for Nix (purple curve in their Fig. 2). We believe that attributing this residual solely to numerical challenges seems unlikely, as it occurs in the exact frequency of υ_{e} and appears prominently when increasing e_{free}. We argue, instead, that this behavior most probably reflects the practical impossibility of defining a zeroeccentric circumbinary orbit, as validated by Gakis & Gourgouliatos (2022). In other words, although the linearized theory allows an orbit exempt of free eccentricity, even adopting such initial conditions inevitably yields that some frequencies will couple to υ_{e}. Accordingly, forced frequencies of the form n_{syn} − υ_{e} are expected, as we explained earlier. In the latter study, we inspect that peaks near the values of n_{S} and υ_{e} in Figs. 5 and 6 of Woo & Lee (2020) might truly appear at some extent from higherorder terms, but perhaps are generated by the merging of the frequencies υ_{e}, k (n_{PC} − n_{S}), and υ_{i}, as quantified by Eq. (9).
Fig. 1 FFT power spectrum for all small moons using the semianalytic approach. In each panel the red vertical line represents the epicyclic frequency of the moon (υ_{e}) and the gray lines the synodic frequency and its harmonics (k n_{syn}), as shown in Table 2. 
Fig. 2 FFT power spectrum for all small moons by varying the total simulated time of 6body integrations (10^{4} and 10^{6} days). In each panel the red vertical dotted line represents the epicyclic frequency of the moon (υ_{e}) and the gray lines the synodic frequency and its harmonics (k n_{syn}), as shown in Table 2. 
3.2 Mutual interactions
The remaining peaks in Fig. 2 at low frequencies are caused by the mutual interactions, corresponding to resonances between the moons. Considering that there are 6 bodies constituting the system, it is understandable to expect that several synodic periods (and their harmonics) can be found. In this section, we determine and identify these mutual frequencies caused by one moon to another. Some frequencies of the perturbations, for the case of Kerberos, were identified by Showalter & Hamilton (2015, extended data Fig. 3 therein). In order to separate the effects by the binary system from those by the other moons, the authors chose to merge Pluto and Charon into a single central body and compare the harmonics of resonances with the peaks of their power spectrum. In this work, we identify the mutual gravitational effects between moons by collating the simulation of the system, accounting for all objects, with a set of runs of a fictitious system where we consider the motion of each moon, accounting only for the gravitational attraction by Pluto and Charon.
Especially in frequencies near zero (i.e., long periods), an immense number of secondary peaks are visible, implying that the orbits also include longperiod frequencies. Actually, in this region, the spectrum is occupied by a forest of very long frequencies. As noted above, this behavior is evident in the 6body, longtimescale simulations. The situation is more prominent for the two lightest bodies, Styx and Kerberos. Their low masses explain their comparatively broader susceptibility to perturbations caused by the other bodies.
To quantify this effect, we provide a comparison of the power spectrum coming from the 6body integrations and the one deriving from 3body integrations (when simulating only the binary and one of the moons). This is shown in Fig. 3. Discrepancies are indeed more significant for Styx and Kerberos. Additionally, Fig. 3 confidently establishes that the large number of peaks that appear in Fig. 2 do not indicate numerical noise, but appear because of the mutual gravitational effects. That is, the longterm periodicity terms of one moon to each other cover almost entirely the frequency region near zero. In 3body simulations, especially this area of the spectrum lacks peaks because only moonbinary interactions are considered.
In order to study more precisely the mutual gravitational interactions, we zoom in the area of low frequencies (<0.20 days^{−1}) of the frequency spectra in Fig. 3. These magnified spectra are presented successively in Figs. 4–7 for Styx, Nix, Kerberos and Hydra, respectively. For a more efficient visualization, we divide the lowfrequency area of each moon into four panels, corresponding to frequencies ≤0.0025 days^{−1} (panels a), between 0.0025 and 0.050 days^{−1} (panels b), between 0.050 days^{−1} and 0.10 days^{−1} (panels c), and between 0.10 and 0.20 days^{−1} (panels d). That way we managed to study in detail the reciprocal effects that gradually become weaker compared to the forced oscillations as the frequency rises. Obviously, when improving the resolution of the periodogram, separate oscillations arise from the seemingly almost continuous spectrum of frequencies in Fig. 3. Some of these frequencies, though immense in number, drop in favor of a few more dominant peaks.
At first, all four small moons are placed close to mean motion resonances (MMRs) with Charon. Specifically, the ratios of their Keplerian orbital periods are about 1:3:4:5:6, similar to the Laplace configuration in the Galilean moons of Jupiter. Although they do not definitely belong in the resonance according to the currently accepted orbital elements for the four moons (Brozović et al. 2015; Showalter & Hamilton 2015), adopting the range of their uncertainties could certainly place them well inside the MMRs (Giuppone et al. 2022). Additionally, a resonant term could affect the tidal damping of a moon even if it is not located in the actual position of the resonance, as shown, for example, for the 3:1 resonance of Nix in Lithwick & Wu (2008). For that reason, we examine this type of resonance in our analysis. In the same study it is shown that even other MMRs of the form N:1 with Charon could be significant for some moons (e.g., the 2:1 MMR of Nix with Charon). Nevertheless, in order to maintain our analysis feasible and focused on the most strongest mutual interactions, we examine only the most prominent 1:3:4:5:6 resonance with Charon. For each moon, the anticipated positions of the harmonics of this resonances is shown in purple vertical dotted lines in the frequency spectra of Figs. 4–7.
More complex resonances can also be defined. After an extensive search for potential resonances, Showalter & Hamilton (2015) found two such major resonances implicating three moons. The strongest resonance identified was Φ = 3λ_{S} − 5λ_{N} + 2λ_{H} ≈ 180°, which implies that the synodic period of NixHydra divided by that of StyxNix is 3:2, i.e., 3S_{NH} = 2S_{SN}, where the subscripts note the respective moons. A second resonance, now involving Styx, Nix and Hydra, was found as 42S_{NK} ≈ 43S_{SN}. We calculate the frequencies induced by the above resonances and indicate their harmonics with yellow and cyan dotted lines in the frequency spectra.
Of course, apart from any resonances we note the strong interaction of one with another over a time period equal to synodic period of them. Therefore, we calculate the respective frequencies by the synodic frequencies S_{SN}, S_{SK}, S_{SH}, S_{NK}, S_{NH}, and S_{KH}. The positions of their harmonics are shown in blue, green, orange, lime, olive and brown vertical dotted lines in Figs. 4–7, respectively.
To avoid further contamination of the images with many vertical lines representing expected oscillatory modes, in these figures we present only the frequencies by the mutual interactions (the forced frequencies by the binary system are shown clearly in Figs. 1 and 2). Peaks where the 6body spectra match with the 3body spectra mark the modes by Pluto and Charon, as explained in Sect. 3.1. As far as the longperiod resonance 42S_{NK} ≈ 43S_{SN} is considered, we show only its first 20 harmonics because additional lines could heavily fill the figures.
Even though each mode would ideally appear as delta function (negligible width), in practise many of them have a significant width. As is evident from the vertical lines in Figs. 4–7, this is not primarily an effect of the limited accuracy of the FFT method. Instead, it is a result of a number of oscillations with very close frequencies.
As noted by Gakis & Gourgouliatos (2022), Keplerian osculating elements are not sufficient to describe circumbinary orbits. Hence, obtaining the synodic periods of one moon with another or the locations of the resonances is uncertain, as it is based on their deduced orbital periods. Future observations may increase the accuracy of the measurements providing more robust estimates for the masses and the 3D position and velocity vectors of the objects. This is yet another reason for possible minor discrepancies between a peak in the periodogram and the anticipated position of its respective frequency (no more than ~0.2%). Apart from that, these orbital elements themselves have not been decisively determined yet, as we find significant differences in previous data tables. For our calculations, we adopted the values in Table 1, deduced by Showalter & Hamilton (2015).
Styx is the moon most heavily affected by mutual perturbations, along with Kerberos. The strongest mode in Styx’s spectrum (Fig. 4) is the one induced by Nix, which is the closest moon. The interaction with Hydra has a similar strength, since Hydra is the most massive moon (though most distant to Styx). It is then evident that the resonance 3S_{NH} = 2S_{SN} produces the most important perturbations in the orbital pattern of Styx. On the other hand, due to its low mass, Kerberos forces much weaker modes to Styx, about 3 times weaker compared to the perturbations by Nix and Hydra. The purple lines in Fig. 4, showing the resonance 1:3:4:5:6, correspond to small peaks either, as Styx does not belong well within the resonance. On the contrary, Styx experiences the effects of the resonance 42S_{NK} ≈ 43S_{SN} strongly at low frequencies, but as the frequency of their harmonics rises, their strength gradually drops. This is, however, the reason for the large number of peaks in between the most prominent ones, by 3S_{NH} = 2S_{SN}, though they are not shown in detail in the diagrams.
Nix appears to have a similar behavior, though the number and strength of peaks is not as large as Styx’s. Again, the resonance 3S_{NH} = 2S_{SN} is the dominant one, while there is a slight preference to the modes by Hydra rather than those by Styx. Despite their proximity, the gravitational interaction between Nix and Kerberos is not powerful; Hydra has a much clearer imprint in the frequency spectrum of Nix than Kerberos. The resonance 42S_{NK} ≈ 43S_{SN} mainly affects the lowfrequency region, as noted by the cyan vertical lines in Fig. 5, whereas the 4:1 MMR with Charon does not seem to produce an intense peak.
The strongest peaks in the spectrum of Kerberos (see Fig. 6) are undeniably the harmonics of the synodic frequency S_{KH}. They are followed by the interactions of Kerberos with Styx, while the mutual effects between Kerberos and Nix seem to be minuscule compared to the above two. However, the 5:1 MMR with Charon has a more substantial impact on Kerberos than the effect N:1 MMR has on Styx or Nix. The resonance 42S_{NK} ≈ 43S_{SN} remains quite strong even for higher frequencies as well, forming evident peaks in between the main components of the synodic frequencies.
Hydra is the most massive moon and has the largest distance from the system’s barycenter. Evidently the effects by either Pluto–Charon or the other three moons are the lowest, and its motion is the closest approximation to a typical Keplerianelliptic orbit in the system. As Fig. 7 suggests, Hydra is mainly disturbed by the resonance 3S_{NH} = 2S_{SN}. The motion of Kerberos also induces some perturbations, though of a smaller scale. Apart from this, Hydra’s resonance with Charon is evident through some intermediate peaks.
For each one of the moons there are a number of loweramplitude peaks that are not explained by the resonances examined above. This effect is more obvious for the lightest moons, Styx and Kerberos. We argue that their origin lies in the N:1 MMRs with Charon. As mentioned, a specific type of resonance could have a gravitational effect on a moon even if it is not located in the exact position of the resonance. Consequently, our conclusion is that these smallerscale modes are created by nearresonance kicks other than 1:3:4:5:6.
Lastly, a significant remark on the power spectra of the four moons is that the peaks by mutual interactions have a comparable size to the binaryinduced peaks for Styx and Kerberos, at least in the lowfrequency region. This is not the case for the heavier Nix and Hydra, where the binary effects are in general stronger in the entire frequency region. This implies that the reciprocal effects for the less massive moons may influence the orbits just as much as the binary system, or even dominate over them for long periods. Consequently, we speculate that modeling the orbits of these moons would not be realistic when ignoring the dynamics of the rest of the objects in the system.
Fig. 3 FFT power spectrum for all small moons for 6body (violet) and 3body integrations (dark blue). The simulated time is 10^{6} days (~2700 yr). 
Fig. 4 FFT successive power spectra for Styx, magnified in low frequencies. Panel a includes frequencies ≤0.0025 days^{−1}, panel b between 0.0025 days^{−1} and 0.050 days^{−1}, panel c between 0.050 and 0.10 days^{−1} and panel d ≥0.20 days^{−1}. Violet plots represent the 6body integrations and dark gray blue plots the 3body ones. Vertical dotted lines mark the expected positions of the main mutual gravitational interactions. 
Fig. 5 FFT successive power spectra for Nix, magnified in low frequencies. Panel a includes frequencies ≤0.0025 days^{−1}, panel b between 0.0025 days^{−1} and 0.050 days^{−1}, panel c between 0.050 days^{−1} and 0.10 days^{−1} and panel d ≥0.20 days^{−1}. Violet plots represent the 6body integrations and dark gray blue plots the 3body ones. Vertical dotted lines mark the expected positions of the main mutual gravitational interactions. 
Fig. 6 FFT successive power spectra for Kerberos, magnified in low frequencies. Panel a includes frequencies ≤0.0025 days^{−1}, panel b between 0.0025 and 0.050 days^{−1}, panel c between 0.050 days^{−1} and 0.10 days^{−1} and panel d ≥0.20 days^{−1}. Violet plots represent the 6body integrations and dark gray blue plots the 3body ones. Vertical dotted lines mark the expected positions of the main mutual gravitational interactions. 
Fig. 7 FFT successive power spectra for Hydra, magnified in low frequencies. Panel a includes frequencies ≤0.0025 days^{−1}, panel b between 0.0025 days^{−1} and 0.050 days^{−1}, panel c between 0.050 and 0.10 days^{−1} and panel d ≥0.20 days^{−1}. Violet plots represent the 6body integrations and dark gray blue plots the 3body ones. Vertical dotted lines mark the expected positions of the main mutual gravitational interactions. 
4 Conclusions
In this paper we identified the most prominent frequencies at which the moons of Pluto and Charon swing. To achieve this, we employed FFTs to computed orbital elements by a semianalytic and an arithmetic process. It is confirmed that forced oscillations caused by the rotating nonaxisymmetric components of the central binary extend to values set by the C_{k} term and occur at frequencies k(n_{S} − n_{PC}). We also notice that although our adopted linearized theory allows it, minimizing such frequencies is impossible even when demanding e_{free} = 0 for the circumbinary orbits, because of the ubiquitous nonlinear terms of gravity.
By collating outcomes from restricted 3body simulations along with the case in which every moon is present, we also managed to assess the mutual gravitational effect. We deduce that the lowmass moons Styx and Kerberos are more intensely affected by the other objects, and over long integration times the mutual interactions have a significant effect on their orbital frequencies. In fact, the mutual effects concerning these two moons have comparable (in strength) perturbations to the orbits with the central binary, at least in the lowfrequency region.
Specifically, we found that the strongest mutual gravitational interactions are caused by the resonances 3S_{NH} = 2S_{SN} and 42S_{NK} ≈ 43S_{SN}. In the first case, strong kicks are sparsely produced, while the second, longerperiod, type of resonance fills the region in between. Nonetheless, the N:1 MMR with Charon is not as powerful, which is attributed to the proximity of the moons to the resonance.
Acknowledgements
The authors thank Prof. Alain Vienne for useful suggestions that enhanced this work. The numerical code used in this work was branched from an nbody code (https://github.com/pmocz/nbodypython) created by Philip Mocz. K.N.G. acknowledges support by grant University of Patras, ELKE81641.
References
 Bromley, B. C., & Kenyon, S. J. 2021, AJ, 161, 25 [Google Scholar]
 Brozović, M., Showalter, M. R., Jacobson, R. A., & Buie, M. W. 2015, Icarus, 246, 317 [CrossRef] [Google Scholar]
 Buie, M. W., Tholen, D. J., & Grundy, W. M. 2012, AJ, 144, 15 [CrossRef] [Google Scholar]
 Gakis, D., & Gourgouliatos, K. N. 2022, Celest. Mech. Dyn. Astron., 134, 14 [NASA ADS] [CrossRef] [Google Scholar]
 Gakis, D., & Gourgouliatos, K. N. 2023, MNRAS, 519, 3832 [NASA ADS] [CrossRef] [Google Scholar]
 Georgakarakos, N., & Eggl, S. 2015, ApJ, 802, 94 [NASA ADS] [CrossRef] [Google Scholar]
 Giuppone, C. A., Rodríguez, A., Michtchenko, T. A., & de Almeida, A. A. 2022, A & A, 658, A99 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kenyon, S. J., & Bromley, B. C. 2019a, AJ, 158, 69 [NASA ADS] [CrossRef] [Google Scholar]
 Kenyon, S. J., & Bromley, B. C. 2019b, AJ, 157, 79 [NASA ADS] [CrossRef] [Google Scholar]
 Kenyon, S. J., & Bromley, B. C. 2022, AJ, 163, 238 [NASA ADS] [CrossRef] [Google Scholar]
 Laskar, J. 1999, in Hamiltonian Systems with Three or More Degrees of Freedom (Berlin: Springer), 134 [CrossRef] [Google Scholar]
 Lee, M. H., & Peale, S. J. 2006, Icarus, 184, 573 [CrossRef] [Google Scholar]
 Leung, G. C. K., & Lee, M. H. 2013, ApJ, 763, 107 [Google Scholar]
 Lithwick, Y., & Wu, Y. 2008, ArXiv eprints [arXiv:0802.2939] [Google Scholar]
 Michaely, E., Perets, H. B., & Grishin, E. 2017, ApJ, 836, 27 [NASA ADS] [CrossRef] [Google Scholar]
 Showalter, M. R., & Hamilton, D. P. 2015, Nature, 522, 45 [NASA ADS] [CrossRef] [Google Scholar]
 Stern, S. A., Bagenal, F., Ennico, K., et al. 2015, Science, 350, aad1815 [Google Scholar]
 Sutherland, A. P., & Kratter, K. M. 2019, MNRAS, 487, 3288 [NASA ADS] [CrossRef] [Google Scholar]
 Woo, J. M. Y., & Lee, M. H. 2020, AJ, 159, 277 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
All Figures
Fig. 1 FFT power spectrum for all small moons using the semianalytic approach. In each panel the red vertical line represents the epicyclic frequency of the moon (υ_{e}) and the gray lines the synodic frequency and its harmonics (k n_{syn}), as shown in Table 2. 

In the text 
Fig. 2 FFT power spectrum for all small moons by varying the total simulated time of 6body integrations (10^{4} and 10^{6} days). In each panel the red vertical dotted line represents the epicyclic frequency of the moon (υ_{e}) and the gray lines the synodic frequency and its harmonics (k n_{syn}), as shown in Table 2. 

In the text 
Fig. 3 FFT power spectrum for all small moons for 6body (violet) and 3body integrations (dark blue). The simulated time is 10^{6} days (~2700 yr). 

In the text 
Fig. 4 FFT successive power spectra for Styx, magnified in low frequencies. Panel a includes frequencies ≤0.0025 days^{−1}, panel b between 0.0025 days^{−1} and 0.050 days^{−1}, panel c between 0.050 and 0.10 days^{−1} and panel d ≥0.20 days^{−1}. Violet plots represent the 6body integrations and dark gray blue plots the 3body ones. Vertical dotted lines mark the expected positions of the main mutual gravitational interactions. 

In the text 
Fig. 5 FFT successive power spectra for Nix, magnified in low frequencies. Panel a includes frequencies ≤0.0025 days^{−1}, panel b between 0.0025 days^{−1} and 0.050 days^{−1}, panel c between 0.050 days^{−1} and 0.10 days^{−1} and panel d ≥0.20 days^{−1}. Violet plots represent the 6body integrations and dark gray blue plots the 3body ones. Vertical dotted lines mark the expected positions of the main mutual gravitational interactions. 

In the text 
Fig. 6 FFT successive power spectra for Kerberos, magnified in low frequencies. Panel a includes frequencies ≤0.0025 days^{−1}, panel b between 0.0025 and 0.050 days^{−1}, panel c between 0.050 days^{−1} and 0.10 days^{−1} and panel d ≥0.20 days^{−1}. Violet plots represent the 6body integrations and dark gray blue plots the 3body ones. Vertical dotted lines mark the expected positions of the main mutual gravitational interactions. 

In the text 
Fig. 7 FFT successive power spectra for Hydra, magnified in low frequencies. Panel a includes frequencies ≤0.0025 days^{−1}, panel b between 0.0025 days^{−1} and 0.050 days^{−1}, panel c between 0.050 and 0.10 days^{−1} and panel d ≥0.20 days^{−1}. Violet plots represent the 6body integrations and dark gray blue plots the 3body ones. Vertical dotted lines mark the expected positions of the main mutual gravitational interactions. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.