Discontinuity in argument of pericenter
Closed this issue · 3 comments
Environment
Which version of REBOUND are you using and on what operating system?
- REBOUND Version: 4.4.1
- API interface: Python
- Operating System (including version): Linux
Hi everybody,
Whilst studying a two-planet system, I have looked at the time evolution of several orbital elements and I want to do a frequency analysis on these time series.
I have noticed a discontinuity in the time evolution of the argument of pericenter of the outer planet when computing the orbits with respect to the system's center of mass, and was wondering why could that be.
Please notice that each point corresponds roughly to one orbit.
I also want to stress that this happens when taking the center of mass as the primary as in sim.particles[j].orbit(primary = com).omega.
This does not happen when using sim.particles[j].orbit(primary = sim.particles[0]).omega, possibly beacuse using the star as reference means imposing on the parameters several frequencies given by the star's orbit around the COM.
The point at t = 0 is an outlier since the two orbits start as circular, so the argument of pericenter cannot be defined at that time. I have already checked that at any t>t0 the eccentricity is greater than zero, so the arg. of pericenter should be defined.
I have also checked the way REBOUND computes the arg. of pericenter and looked at the time evolution of the longitude of ascending node: the variation of the long. of the ascending node is at most of order 10^-4, so it's not connected to the discontinuity, which is of order 1.
Another test I tried is to check for time reversibility and the behaviour is the same wether I integrate forwards or backwards.
To my untrained eye it could seem that the argument of the outer planet precesses with period pi, instead of 2pi (if you fold the upper sections by subtracting 2pi), but I wouldn't know why...
As a comparison this is the time series for the inner planet
Hm. I'm pretty sure this is not a bug in the orbital elements routine, but a question on how to interpret the result. Maybe calculating the planet's orbital elements with the respect to the system's centre of mass (including the planet itself) is not what makes sense physically...
Did you check e is not 0 also relative to the COM? It might be useful to make a plot of (e cos omega, e sin omega) to see whether that's discontinuous. I think those are usually better variables for frequency analysis than the angles that have coordinate singularities
Thank you both for your answers!
The eccentricity for the outer planet is indeed very small at the times of the discontinuity (order 10^-5), so, even if it's never zero, I guess that's why.
Looking at the eccentricity and the (e cos omega, e sin omega) plot also gave me a possible physical answer to what is happening: the orbit becomes progressively less eccentric until it circularizes and starts to get more eccentric in the other direction (i.e. the COM shifts to being the the other focus of the ellipse), that would also explain the difference of pi between above and below the jump; after that the orbit precesses linearly and the arg. of periastron decreases linearly untill it wraps back to 2pi and keeps on decreasing until another jump during the next time e~0.
Let me know if it makes sense to you!
Also, thank you @dtamayo for the suggestion on using e cos(omega) and e sin(omega) for the frequency analysis and @hannorein for the reflection on the COM, I'll need to think about that.