A Numerical Method for Determining the Elements of Circumbinary Orbits and Its Application to Circumbinary Planets and the Satellites of Pluto-Charon

Planets and satellites orbiting a binary system exist in the solar system and extrasolar planetary systems. Their orbits can be significantly different from Keplerian orbits, if they are close to the binary and the secondary-to-primary mass ratio is high. A proper description of a circumbinary orbit is in terms of the free eccentricity $e_{\rm free}$ at the epicyclic frequency $\kappa_0$, forced eccentricity $e_{\rm forced}$ at the mean motion $n_0$, and oscillations at higher frequencies forced by the non-axisymmetric components of the binary's potential. We show that accurate numerical values for the amplitudes and frequencies of these terms can be extracted from numerical orbit integrations by applying fast Fourier transformation (FFT) to the cylindrical distance between the circumbinary object and the center of mass of the binary as a function of time. We apply this method to three Kepler circumbinary planets and the satellites of Pluto-Charon. For the satellite Styx of Pluto-Charon, the FFT results for $\kappa_0$ and $e_{\rm free}$ differ significantly from the first-order analytic value and the value reported by Showalter & Hamilton (2015), respectively. We show that the deviation in $\kappa_0$ is likely due to the effect of the 3:1 mean-motion resonance and discuss the implications of the lower value for $e_{\rm free}$.