Stability of Laplace resonance

Question

Someone is making a mod for the video game Kerbal Space Program (KSP), which implements propper N-body with point masses into the game (and possibly in the future also other dynamics), instead of the stock patched conics approximation. The "solar system" in KSP is a simplified and scaled down version of our own solar system. Someone did a N-body numerical simulation of this system, in which it can be seen that not all orbits are actually stable, namely two moons around the gas giant Jool, called Vall and Bop.

The gas giant Jool has, similar to Jupiter, three inner Moons which have a 1:2:4 resonance, however the middle one Vall gets kicked out very early on in the simulation.

This raises a question about when 1:2:4 resonance is stable/self-correcting. The main difference I noticed between the KSP and the real solar system is that in KSP the moons are much heavier compared to their parent body. For Io, Europa and Ganymede the mass ratios vary roughly between $2.5\ 10^{-5}$ and $7.8\ 10^{-5}$, while for Laythe, Vall and Tylo (their KSP equivalents) the mass ratios vary roughly between $7.3\ 10^{-4}$ and $1.0\ 10^{-2}$. The orbital elements and physical parameters of the celestial bodies in KSP can be found here. So this makes me wonder if there is a known limit (rough estimate) to this.

I have done N-body simulations with point masses myself with only two of the moons in 1:2 resonance with their default masses (so Laythe and Vall or Vall and Tylo) in which case their orbits seems to remain stable. Here are a few results, in which I calculated the semi-major axis, $a$, and eccentricity, $e$, as follows, $$ a = \frac{\mu r}{2\mu - v^2 r}, $$

$$ e = \sqrt{1 + \frac{v_\theta^2 r}{\mu} \left(\frac{v^2 r}{\mu} - 2\right)}. $$

For a reference in terms of time scale on the results the orbital period of the most inner moon Laythe is equal to roughly $5.3\ 10^4$ seconds. But I also have done simulations with all three moons in which I decreased the masses of Laythe and Tylo by a factor 10 (Vall already has a mass an order of magnitude smaller), in which case the their orbits also seem to remain stable. These results can be viewed here. This does looks more chaotic than any of the two moons configuration.

This is a related question. It doesn't has an answer yet, but one of its comments does link to a document, but this does not mention mass ratios as a criteria for whether three body resonance would be stable or not.

To be clear, I am only considering gravitational interaction between point masses, so for now no tidal forces or nodal precession. I am looking for some limit in terms of mass ratios between the parent and its satellites, similar as for example the Hill sphere which gives an approximate limit for the distance at which stable orbits are possible. Even though the N-body problem is considered chaotic I still would suspect that there would be some way to determine this ratio. Intuitively I would think that with only gravitational interaction such a system still can stable and self-correcting.

For example if the middle satellite would get perturbed at the closest approach of the inner satellite, such that the opposite side of its orbit would get raised, then that would mean that its period would increase and the points of closest approach with the other two satellites would start to drift backwards in the orbit. This initially would amplify the initial perturbation, since when the inner satellite would be at the initial position of closest approach it would be slightly ahead of the middle satellite, pulling it along slightly, increasing the opposite side of the orbit of the middle satellite even more. But once the points of closest approach would have drifted back far enough that the separation of would be to large to induce any significant pull. Once this drift of the points of closest approach would nearly have looped around 360° the separation distance with the inner satellite would decrease again, but his times each closest approach the inner satellite would pull the middle satellite back, counteracting the initial perturbations.

Answer 1

First off, are you capturing all of the dynamics? In particular, are you modeling

• Jupiter's second dynamic form factor J₂.
  This is fairly simple. It just means you need to take a tiny step beyond point mass gravitational models. Since Jupiter has a very large J₂ (0.014733), you had better be modeling it.

• Jovian tidal dissipation, which transfers angular momentum from Jupiter's rotation to Io's orbit.
  This is a tougher task. You will need to model it if you want to study long-term behavior (e.g., Yoder). That means you'll need values for Jupiter's k₂ Love number and Jupiter's tidal dissipation quality factor Q. For another, there's a huge disagreement in the scientific medium over the value of Jupiter's tidal dissipation quality factor Q, ranging from 36,000 (Lainey et al.) to 10⁹ (Wu). Modeling this is essential. This transfer is often cited as the reason the innermost Galilean moons came into Lagrange resonance. Whether this is a stabilizing influence is debated. Some (e.g. Lainey et al.) think that this resonance is rather transitory.

• Ionian tidal dissipation, which tends to circularize Io's orbit.
  This is an even tougher task. Multiple papers show a hysteresis loop where Io alternates between a very nearly circular orbit with very little tidal heating to a more eccentric (but still close to circular) orbit with significant tidal heating. Suppose Io's orbit is very, very close to circular. The resonances with Europa and Ganymede tend to make this circular orbit eccentric. This makes Io subject to tidal heating by Jupiter. This in turn distorts the shape of Io, which in turn makes Jupiter act to circularize Io's orbit. Rinse and repeat. The thermal time lag means that the eccentricity oscillates rather than stabilizes to some steady value.

The latter two items listed above are non-conservative forces, which draws into question the relevance of using a symplectic integrator. That KSP's target usage is to simulate space voyages also calls that choice into question. There's long been a debate amongst developers of N body simulators for small N (e.g., our solar system) whether its better to use a lower accuracy symplectic integrator versus a non-symplectic integrator that offers greater accuracy and stability over the short haul. Those who develop planetary ephemerides that span thousands of years or more tend to prefer accuracy over symplecticity. It's the other way around with those who study the long-term (many millions of years) stability of the solar system. One thing is certain: You can't have both. It just doesn't work that way.

I don't know how far KSP goes with regard to modeling spacecraft fidelity, but if you're developing a simulation that has flight software in the loop, you can generally kiss the conservation laws and accuracy goodbye.

References:

Hussmann, et al. "Implications of rotation, orbital states, energy sources, and heat transport for internal processes in icy satellites," Space Science Reviews 153.1-4 (2010): 317-348.

Lainey, et al. "Strong tidal dissipation in Io and Jupiter from astrometric observations," Nature 459.7249 (2009): 957-959.

Peale, "Origin and evolution of the natural satellites," Annual Review of Astronomy and Astrophysics 37.1 (1999): 533-602.

Wu, "Origin of tidal dissipation in Jupiter. II. The value of Q," The Astrophysical Journal 635.1 (2005): 688.

Yoder, "How tidal heating in Io drives the Galilean orbital resonance locks," Nature 279 (1979): 767-770.

Answer 2

I'm eggrobin, I happened across this while doing research for my mod.

Regarding the question:

The simulations by Matt Roesle that showed Vall being ejected from the Jool system were computed by interpreting KSP's orbital elements as body-centric, which means that when the system was integrated, the bodies ended up in significantly different orbits around the barycentre of the system, so that they were not in resonance to start with. Scott Manley [1] has performed simulations where the orbital elements were interpreted so as to put the bodies in resonance (this was before the introduction of Pol, but it is unlikely that Pol changes things that much), in that case the system was found to be stable for at least a thousand years (Scott did not mention his methods, but since he worked with Chambers on that I suspect the integration was done using the Mercury package).

Regarding the answer:

I am taking $J_2$ into account when modeling our solar system, not doing so would be quite silly: the influence of $J_2$ on the Galilean moons is two orders of magnitude greater than the influence the other Galilean moons [2, p. 69]. I am not modeling the nonconservative effects mentioned here at the moment.

The concerns about spacecraft are irrelevant. While it is obvious that the spacecraft are subject to nonconservative forces, and that, at least during manoeuvres, they must be integrated differently (and in the current state of the mod, they are, the timestep used for the vessel is shorter by several orders of magnitude during manoeuvres), the whole system still has to be simulated (ephemerides for fictional solar systems are hard to come by), and long stretches of flight without correction could legitimately be symplectically integrated.

The question of symplectic vs. nonsymplectic methods is definitely relevant, and I am investigating it now. It should be noted that going for symplecticity does not necessarily hamper precision, for instance, Blanes, Casas and Ros [3, fig. 1] pit good old Dormand-Prince against a processing symplectic RKN method on the Kepler problem (usual $T + V$ splitting), and the symplectic method is more precise (average position error) for the same computational cost by at least four orders of magnitude.

One can note that most papers focusing on symplectic integration of the solar system do not take moons into account with the occasional exception of ours, which is generally bizarrely treated. One of the reason seems to be that the usual splittings are only efficient without close encounters, with some integrators handling close encounters by special treatment (for instance [4]), and that moons are effectively permanent close encounters. I recently came across the splitting given [5], which, while designed for stellar systems, may apply to this (possibly with individual timesteps in the style of [6], and certainly with higher-order methods from [7]). I may experiment with other splittings.

While the remark on low $N$ is valid in the current state of affairs, if I were to try to model asteroid belts $N$ would cease to be low.

References, not in alphabetical order because I'm lazy, no links because this is my first post around here:

[1] Scott Manley, private communication. Yes, this kind of citation is bad.

[2] Chandler (1973), Determination of the dynamical properties of the jovian system by numerical analysis.

[3] Blanes, Casas and Ros (2001), High-order Runge–Kutta–Nyström geometric methods with processing.

[4] Chambers (1999), A hybrid symplectic integrator that permits close encounters between massive bodies.

[5] Beust (2002), Symplectic integration of hierarchical stellar systems.

[6] Saha and Tremaine (1994), Long-term planetary integration with individual time steps.

[7] Farrés, Laskar, Blanes, Casas, Makazaga and Murua (2013), High precision symplectic integrators for the Solar System.