Running MEGNO questions

[zhexingli]

Hello again,

• When running MEGNO for a planetary system of say three planets, is there a way to run that simultaneously for all planets at once, or I can only do one at a time (let one of the planet's starting e and P vary within certain range while fixing starting positions for the other two planets)?

• Also, when exploring different parameter space of the same planetary system, does changing the time steps affect the consistency of the result? For example, I want to see if a planet will be stable within 1 AU and I'd set a time step of 0.05 year. But when checking if a planet will be stable between 10-15 AU, a time step of 0.5 year will probably enough. So is there a consistency problem I need to worry about when exploreing different orbital distances for the same system?

Thank you!

[hannorein]

I don't understand your first question. MEGNO gives you one value per system.

Regarding your second question. If you have a chaotic system, yes, the results of an individual integration will change when you change the timestep. However, any quantity which you measure over an ensemble of simulations should be the same within statistical error bars. If not, then the simulations are not converged.

[zhexingli]

Thx, is there any way that I can tell if the simulations have converged?

For my first question, the example given here https://rebound.readthedocs.io/en/latest/ipython/Megno.html The example system has two planets. But when running MEGNO here, only the outer planet varies its e and a while the inner planet's parameters are fixed (is it fixed throughout the simulation or just the starting position?). And the result shows the stability for the outer planet. In this case, only one planet's parameters are varied, I'm wondering if I could vary e and a for more than one planet if I want to see the stability map for the inner planet as well? Or I have to run another simulation with inner planet's parameters float while keeping the outer ones fixed.

[hannorein]

You can tell if a set of simulations is converged by changing the timestep. If the results in an statistically large enough ensemble depend on the timestep, then the timestep is not small enough.

You can very whichever parameters you like. MEGNO measures if a system is chaotic or not. As I mentioned above, this is one value per system, not one per planet. It does not make sense to speak about the stability of each planet individually in a two planet system.

[zhexingli]

So in this example, we're only varying the outer planet's e and a. Are the parameters of the inner planet fixed throughout the simulation or it's just the starting point that's fixed? And the final plot tells us if the system will be stable when the outer planet sits between 7 and 10 AU with eccentricity between 0 and 0.5 right?

If I want to check if the system is stable when I vary the inner planet within some ranges, can I do that at the same time as while I'm varying the outer planet's orbit? Or I can only vary one planet at at time? Such as fix inner vary outer, then second time fix outer vary inner?

[hannorein]

It's just an example. You can vary any parameter you like.

[zhexingli]

Ok...thx!

[zhexingli]

Hello, I've got some more questions I need to double check with the expert here. So I'm trying to run MEGNO for 47 UMa, specifically trying to see if I can put an Earth-sized planet in betwwen planet c and d and still make the system dynamically stable. So here's the result of the MEGNO:

As you can see there's a very tiny region of greenish area that reprensents stability. I doubled check that by running stability code just for the new planet at P = 20 and e = 0.06 and the system is stable. But as you can see from the MEGNO plot, at P = 20 and e = 0.00, it shows red which means it's unstable. But when running stability code for that case, it indicates stability, as shown here:

So I'm wondering why does MEGNO indicate the system would be unstable if I add a planet at P=20 years with e = 0.00 but the stability code thinks it's stable?

Another question, I think I asked last time what those white dots in the MEGNO plot mean and I was told that's likely a close encounter right? I checked at those white dots locations with stability code and found out they're indeed unstable. But just curious, why don't they just show up as red points like the rest of the field in the MEGNO plot? Why is it white?

Thank you.

[hannorein]

The most important thing is that MEGNO (or equally the Lyapunov exponent) tells you if the system is chaotic or not. That's not the same as stability.

The parameter space you're exploring seems quite large. For example, I'd be very surprised if any planet with an eccentricity of 0.9 would be stable in such a system. You might see more structure if you zoom in a bit more in your MEGNO plot, i.e. eccentricities which do not lead to orbit crossing).

I'd suggest you run the same grid of simulations but without MEGNO. Integrate for a short time so they finish in a reasonable time. Then plot whether the system is stable over that short timescale or not.

Ideally you'd like to just integrate these simulations for a billion years. Then you would know if they are stable or not. But that takes too long, which is why you might want to use MEGNO. You can use it to predict if the simulation is chaotic or not. Chaotic simulations might be more likely to become unstable over long timescales. But it's not the same criteria.

I'd also suggest running these tests first with a very small mass for your additional planet (or even a test particle). The results will be easier to interpret. If you find a stable island, you can then slowly increase the mass of the additional planet and see if it remains stable.

Regarding the white dots. It's probably a NAN or INF value. You might want to do something special to capture simulations which lead to close encounters or ejections, then mark them with a special colour.

[zhexingli]

That sounds great! This is super helpful! Thank you!

[zhexingli]

I ran the above MEGNO simulation with a smaller eccentricity range, but it doesn't show any more finer features. Is there a reason why at P = 20, e = 0 it says chaotic but for stability simulation it shows stable?

[hannorein]

Maybe it is both chaotic and stable?

[zhexingli]

It sounds very counter-intuitive... but the stability simulation shows the variation is very periodic rather than random changes through time. Doesn't that mean non-chaotic if it doesn't show any chaotic behavior? What are the typical features that belong to chaotic behavior when we run stability simulations?

[hannorein]

It might be very hard to see by eye. I'd encourage you to run a grid of stability simulations so that you can make a plot where you see some structure. For example, you should definitely be able to see orbit crossing orbits, mean motion resonances, etc. Even for relatively short integration times. Try to understand those plots first. Then see if you can get the same results with a stability indicator such as MEGNO, but much faster. Then you can start some "production runs" with higher resolution for your grid, more parameters to explore, etc.

[zhexingli]

Sounds good. thank you.

[zhexingli]

Is there a unit system in REBOUND that uses 'day' as its time unit? Currently I'm using the 'yr', 'au','Msun' unit system and I'd like to switch the 'yr' to 'day', but REBOUND seems to have a problem with that. I changed the orbital period in days and also did that for the time steps and integration time. But REBOUND is saying the time step exceeds the orbital period, even though my time step is set to sim.dt = 0.4 and the shortest period is 8.7 days. I'm not sure why what's going on...

[hannorein]

The by far easiest way to handle units is to use G=1 and adjust all your units accordingly (mass in solarmass, length in AU, time in year/2pi).

[zhexingli]

Alright. Thx.

[zhexingli]

Just double checking: when running in rebound without specifying inclination, is it assuming all planets are orbiting in the same orbital plane?

[hannorein]

Yes.

[zhexingli]

ty!

[zhexingli]

Hello, this has nothing to do with previous postings in this thread, but falls under the same category as the title. We're running some fun simulations about Jupiter's moon system stability by injecting a particle into various locations within the system. We tried using simply N body to see the time evolutin of orbital elements as well as MEGNO. But we're getting different results. For MEGNO, we tested the particle's location from Roche limit to 1/10 Hill's radii with only 1000 year integration time and the result already indicates chaos for all locations tested. But when checking the time evolution of orbital elements of the system bodies for 1 million years when particle is at 0.03 AU from Jupiter with zero inclination, everything looks fine. Nothing got ejected. (There's variation in semi-major axis evolution for all bodies)

So MEGNO says the system with injected particle would be chaotic after 1000 year run, but N body with even 1 million years shows nothing is ejected. Are these results contradicting themselves or is this one of the cases where the test system is chaotic but stable due to different criteria for chaos and stability? Thank you.

[hannorein]

A guess: you might need to run the simulations with megno for longer. At the moment, you're effectively extrapolating from 1e3 to 1e6 years. That's hard!

[zhexingli]

Thanks for the quick reply! I thought if the system is chaotic at the 1000 year mark, it'll continue to be chaotic for the longer term, is that assumption not right? Can a system return from being chaotic to not chaotic? (Meanwhile I'm putting a new megno sim with 1e6 years for one grid to see what megno returns for this (0.03 au and 0 inclination).)

[dtamayo]

To my eye the results are consistent. The eccentricity evolution in purple of the test particle definitely does not look oscillatory/regular, it is wandering around, so definitely chaotic. Unfortunately, I think the problem of calculating how long it takes a chaotic system to go unstable is not at all obvious/straightforward. You might have to run your Nbody integrations for longer to see the instabilities (which quickly becomes painful!)

[hannorein]

It's not always easy to estimate how long you have to integrate for to get a reasonable estimate for MEGNO. If you have an a priori estimate of the Lyapunov timescale, that's a good order of magnitude to keep in mind.

[zhexingli]

Ok, thanks Hanno and Dan. So one should always be extremely careful when interpreting megno results? Always better to check system bodies through N body time evolution for stability because it could still be stable within short term but is chaotic?

[hannorein]

Yes, it can be tricky to interpret what the MEGNO value means. Keep in mind that chaos is not the same as instability. MEGNO can help you to get a rough overview of a large parameter space, to identify resonances, stable islands, etc.

[zhexingli]

Thank you!

[zhexingli]

Sorry, one more question. The result of the MEGNO simulation I put in earlier as Hanno suggested to run MEGNO for a longer time came out. I ran the simulation with the same setup as the N body time evolution sim shown above for 1e6 years when the particle is at 0.03 au and 0 inclination, but MEGNO returned NaN value. Why is that? Is it because even though within 1e6 years the system seems stable, but the irregularities in the orbital evolution is making the chaos indicator impossible to predict and yield a value?

[hannorein]

The distance between nearby orbits can grow exponentially fast. This is probably what's causing your NaNs.

[zhexingli]

Hmm.. but I'm not seeing this growth in between orbits from the semi-major axis, eccentricity, and inclination time evolution plots..

[hannorein]

I'm not sure I can offer any useful insight. From everything you've described, I don't see anything wrong with the way you use REBOUND, nor do I think there is a bug.

As a way forward, I'd suggest to re(read) the MEGNO papers. Try to reproduce some simple examples and get a feeling for when MEGNO is the right tool to answer a dynamics questions and when not. If you want me to help you out with this specific problem, please post some code that I can run to reproduce the inclination versus semi-major axis plot from above (remove everything that is not essential but keep enough so I can run the code).

[zhexingli]

Yep. I'll close it for now. Thanks for the help.