orbital elements change with only one planet in the system

[nchoksi4456]

Hi all

I've been running REBOUND simulations of densely-packed systems, starting with N ~ 10 planets and letting them scatter / merge / eject each other. In a small fraction of simulations I've noticed some odd behavior. The system will whittle itself down to just one remaining planet (and the star). The final planet behaves weirdly. See the plot below for an example.

It seems impossible to me that the semi-major axis and eccentricity could change when there are no other planets in the system (it should just be the Kepler two-body problem?). And since the energy E ~ 1/a, energy is apparently not conserved. Checking sim.calculate_energy() verifies this. However, I have verified that angular momentum is conserved in the data plotted above.

I am using MERCURIUS, with a timestep set to 5% of the shortest orbital period (the timestep is updated every 1 yr).

Any ideas would be appreciated. Here I've uploaded the simulation archive for the simulation in the plots above. And here is my little analysis script that generated the plots shown above.

I am continuing to poke around, but I've been unable to find the source of the error so I figured I would post here in case anyone has any good ideas. Thanks!

[hannorein]

Hi Nick,

Good question. Thanks for providing the SimulationArchive! I had quick look at it. I can't find anything wrong.

You're right in your assumption that with just two bodies, it should be easy to solve. I suspect there is something going on with MERCURIUS. There are two possible reasons: 1) It won't be able to solve the two-body problem exactly because of the coordinates it's using. That can be an issue for very eccentric planets. 2) It might switch to and from the close encounter mode every orbit (because it's a very eccentric planet). That will introduce an error every at every switch

I can't see an obvious solution. You can always make the timestep smaller. 5% of the innermost planet's period is a reasonable for circular planets. But for eccentric planets, a better number would be 5% of the shortest pericenter passing time (something like pericenter distance divided by velocity at pericenter).

Sorry for not having a better solution at the moment. I'll keep thinking about it.

Hanno

[nchoksi4456]

Thanks for looking into it Hanno. The systems I'm simulating are super chaotic so it is hard for a single test case to determine whether a shorter timestep is really fixing the issue. I'll have to run an ensemble of systems and see whether these sorts of spurious results are still present in the final distribution.

2) It might switch to and from the close encounter mode every orbit (because it's a very eccentric planet). That will introduce an error every at every switch

Do you mean it is switching into close encounter mode because of close encounters with the star? Before the sudden jump near the end of the integration, the planet has an eccentricity of just ~0.4 and a semi-major axis ~ 1.6 au, so the pericenter is only (1-0.4)*1.6 au ~ 1 au, so it is still not getting super close to the star. Maybe I'm misunderstanding you though, and there's another reason MERCURIUS would jump into close encounter mode.

[hannorein]

Yes, it could be switching to the close encounter mode with the star. If it does that many times, it could slowly increase the eccentricity. I'd need to run a few tests to see what's going on.

Just to state the obvious: this could also be a bug, not related to anything physical.

In any case, I agree, these results are not very satisfactory! I'll look into. I'm not sure how fast I'll be, maybe a couple of days. If you wanted to give me a head start, you could try to integrate a simple system with just 1 planet. Set the eccentricity, semi-major axis, switching radius, timestep to the same values as in your chaotic simulations and see what happens. Can MERCURIUS integrate that?

[nchoksi4456]

Thanks Hanno. I'll try this test and post the results here soon.

[dtamayo]

Another useful test might be to run WHFast with your single planet system parameters. Like Hanno said, it might just be that your global timestep is too long to resolve the pericenter passage

On Sun, May 9, 2021, 4:27 PM nchoksi4456 @.***> wrote:

Thanks Hanno. I'll try this test and post the results here soon.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/hannorein/rebound/issues/533#issuecomment-835881043, or unsubscribe https://github.com/notifications/unsubscribe-auth/AA2ABFY6S2XSFF7W3OSTFPTTM3V47ANCNFSM44OSXYBA .

[nchoksi4456]

I tried the test Hanno suggested: I took the first output from the chaotic simulations when there was one planet remaining (not the same as the first timestep when there was one planet remaining, since I only output to the SimulationArchive every 3000 yr), and extracted the positions and velocities of the planet and the star. I then used these data to set up a new single-planet simulation, also using MERCURIUS.

Here is the result:

And here is the energy error (appreciable) and angular momentum error (~zero):

Qualitatively, the result is similar to that seen in the original simulations (see caveat/question below though): a and e behave ~reasonably until suddenly everything goes haywire (in the plot above, things go haywire at ~0.8 Myr), with a and e both suddenly increasing.

If I cut the timestep down to 2% of the orbital period, I get results that are stable -- at least over the 1 Myr I integrated for.

I'm a little confused why the time for things to go "haywire" in the dt = 0.05 case is different from that seen in my original simulations. In the original simulation (plot at top of this thread), it took (17.5 - 17.3) Myr = 0.2 Myr, whereas in this test case it took ~0.8 Myr. Since I set up everything identically (as much as I could), I would have expected nearly identical results.

[hannorein]

Thanks! I don't think I would expect the result to be exactly the same (unless you really get everything set up exactly the same down to the last bit). The case where something goes wrong might simply be chaotic (for non-physical reasons).

But one other reason why it might not be exactly the same is because each particle's critical radii is calculated at the beginning. If you planet was on a different orbital before all the scattering happened, it will sill have it's origin critical radius.

You could inspect the critical radii with sim.ri_mercurius._dcrit[I].

Can you share how you setup the simulation with 1 planet?

[dmhernan]

Hello,

I wanted to chime in here that, at least in the original formulation of Mercurius, close encounters between planet and star don't introduce any changes in the Hamiltonian splitting (e.g., they don't induce close encounter mode). So in this case, you're essentially just using the Wisdom--Holman method in democratic heliocentric coordinates, and switching errors are not an explanation. Or my memory could be rusty...

However, one thing I wanted to mention is that, whether it's Mercurius or Wisdom--Holman, adapting the timestep can also be problematic and can cause error jumps in energy, but not angular momentum.

David

[nchoksi4456]

However, one thing I wanted to mention is that, whether it's Mercurius or Wisdom--Holman, adapting the timestep can also be problematic and can cause error jumps in energy, but not angular momentum.

Thanks David. I think you are right, updating the timestep is the origin of the issue. I was updating the timestep every ~yr of my simulation by checking the periods of all the remaining planets at the current output and setting the timestep to 5% of the shortest orbital period. The relevant Python code looked like this:

    for time in times: 
        sim.integrate(time, exact_finish_time = 0); 
        particles = sim.particles
                Plist = []
        for j in range(1, len(particles)): # Loop over all planets
            particle = particles[j]
            dist = np.sqrt(particle.x**2 + particle.y**2)
            particle.r = calc_Rplanet(particle.m, dist) # Update particle sized based on planet mass and heliocentric distance
            Plist.append(particle.P); 
                 sim.dt = 0.05*np.amin(Plist) # Update timestep based on shortest orbital period

Now, for the case of a single remaining planet, the orbital period should of course formally stay constant, so dt = 0.05Porb would also be constant. But because the computer is not perfect, the orbital period DOES change by an epsilon between simulation steps, and so my code above also leads to an epsilon change in timestep. Apparently MERCURIUS does not like this epsilon change in timestep. The plot below uses the same single-planet simulation from my previous post, but sets dt = 5% of the INITIAL orbital period and then leaves the timestep FIXED (i.e., just comment out the line sim.dt = 0.005*np.amin(Plist)in the code above). Everything looks great now (compare to dt = 0.05 plots in my last post)!

Does anyone have a better way to update the timestep, one which won't break MERCURIUS? My simulations initially lay down planets between ~10 days and ~few thousand days. I could just set the timestep to a constant 0.05*10 days, but this seems wasteful...The planets at the inner edge often scatter / merge outward. So for example if the system whittles its way down to 2 planets, both at > 100 days, it would be good to update the timestep to something like a few percent of 100 days, rather than few percent of 10 days. Sorry if this is a naive question...

[dmhernan]

Glad that worked out. Just to note that MERCURIUS cannot solve the 1-planet problem though, it gets the period wrong. Otherwise your energy error would be far smaller at machine precision.

Your question is not at all naive. The short answer is that yes, it's possible in many cases, but it would be a mini research project in of itself. You need to make sure you would get an identical set of timesteps if you integrated forwards as if you integrated backwards in time. E.g., make your timestep a symmetric function of the past state before the timestep and future state after the timestep. For more, https://ui.adsabs.harvard.edu/abs/2018MNRAS.475.5570H/abstract

[dmhernan]

P.S.: Do you know why there's a big jump in angular momentum error initially?

[hannorein]

Thanks Nick and David. I think the frequent timestep changes explain everything. The main issue is not that you change the timestep, but that you do it so often (once per year). I don't think there is a perfect solution for this type of problem. Different people have come up with different ways approaches. The easy thing to do would be to continue changing the timestep, but rarely. Then the error that you accumulate each time doesn't add up as much. You could for example check if you need a smaller timestep, like you do now, but only ever decrease the timestep, and only by factors of 2. That way your timestep won't change very often and also doesn't oscillate back and forth. You could sporadically allow the timestep to get larger again depending on what you're trying to simulate. Again, as long as it happens rarely, I wouldn't worry too much about it. You just need to avoid the kind of exponential instability that you initially saw in your simulation.

[nchoksi4456]

I think the initial jump is just because I was plotting (L - L[0])/L[0] on a log-scale. So initially the error is zero by definition, which is off-scale for a log plot?