Anomalous MEGNO values

[jrlivesey]

Hello,

I’m simulating a two-planet system to evaluate the MEGNO, using a range of eccentricities for both planets. Some cases exhibit a strange behavior in which the orbits are clearly quasi-periodic, and the MEGNO seems to converge to 2, but then after so long the MEGNO diverges sharply in the negative direction. Below is my simulation set-up.

import rebound
import numpy as np
import astropy.units as u

sim = rebound.Simulation()
sim.integrator = 'whfast'
sim.ri_whfast.safe_mode = False
sim.ri_whfast.corrector = 11

# Primary
sim.add(m=0.12)

# Planet d
d_mass = float(0.26 * u.earthMass/u.solMass) / np.sin(133.)
sim.add(m=d_mass, a=0.029, e=0.61, inc=0., Omega=0., pomega=0.)

# Planet b
b_mass = float(1.07 * u.earthMass/u.solMass) / np.sin(133.)
sim.add(m=b_mass, a=0.049, e=0.36, inc=0., Omega=np.pi, pomega=np.pi)

sim.dt = 0.05 * sim.particles[1].P
sim.move_to_com()
sim.init_megno(seed=0)

Integration for $\sim 10^6$ orbits of the outer planet yields the following plot.

This happens using both the WHFAST and IAS15 integrators. I cannot explain why this might occur, so any help is appreciated!

[hannorein]

Which version of REBOUND do you use?

[hannorein]

Alright. I can reproduce this with the latest version of REBOUND. That means it has nothing to do with any rescaling issue or running out of the floating point range.

You seem to be doing everything right in your setup. I can never rule it out 100%, but I don't think this is a bug in REBOUND either. So I suspect this is a real feature of the system. You're right that it looks well behaved over 30kyr. But it could be that the system is nevertheless chaotic on slightly longer timescales and that MEGNO picks that up. I would explore that idea a bit: Run a longer integrations and see if anything interesting happens. MEGNO is sometimes a bit hard to interpret, so you could also manually add a set of variational equations and see how nearby trajectories diverge (there are several examples that show you how to do that).

[jrlivesey]

Thank you for your reply, @hannorein. Integration for 5 Myr ($\sim 3 \times 10^8$ orbits of the outer planet) yields the following plot: rebound-issue-5Myr.pdf. I have incorporated one of the variational co-ordinates, which does grow exponentially in this simulation. The orbital evolution remains regular (the eccentricities oscillate at only one frequency), but the MEGNO dives to extremely negative values. To me, these two results appear inconsistent. Could you tell me whether I'm missing something here?

[hannorein]

Again, I don't think you're doing anything wrong, nor do I think this is a bug in REBOUND. It's just a matter of how to interpret the MEGNO values. I'm just speculating here, but maybe the orbital phases get out of sync after ~30kyr. That might not affect the stability or further evolution of the system because it evolves secularly. But MEGNO might have a hard time with it because it works on the cartesian coordinates rather than orbital elements. If you want to gain some confidence on how to interpret MEGNO, try running a grid of simulations (see e.g. this example). Choose some reasonable limits on the axes so that you are guaranteed to see some effects (MMR, orbit crossing, etc). Then play with different integration lengths, 10kyr, 30kyr, etc, and see where the map looks cleanest.

[hannorein]

I'm closing this issue for now but feel free to reopen it if you think there still is a problem.