N-body Integrator

The N-body integrator I have developed for exomoon problems uses Mercury (Chambers 1999) as the basis. It treats planets realistically as oblate, and the planet spin is permitted to evolve using a non-secular equation of motion, as required by the orbital evolution on short time scales in the planet-planet scattering phase.

Treating planets as oblate ensures realistic orbital behavior of the moons close to their parent planets. Without planet oblateness, unrealistic instability effect could happen to the moons (Hong et. al 2015).

The evolution of the planet spin can affect the stability of moons (Tremain 2009), and in systems with planetary close encounters, planets could acquire very high obliquities (the angle between the orbit and spin angular momentum of the planet), considering this, the evolution of planet spin is incorporated in the code. In addition, in such systems, the orbits of planets evolve on short timescales, therefore requiring a non-secular equation of motion for the evolution of planet spin.

The integrator uses the Bulirsch-Stoer algorithm, which is a symplectic N-body integrator that can retain high accuracy over long time scales with close encounters between bodies. Since stable moons always orbit their parent planets closely, this algorithm suits moons problems best, and of course as well as planetary close encounters.