Inconsistency in solution across inclination

[landenconway]

I am trying to reproduce Figure 3 from https://arxiv.org/abs/1107.2414 and I'm having trouble producing consistent results when varying the individual inclinations. If I use i=64.7 and i3=0.3 I am able to produce the same the plot. If I change the inclination say to i=40 and i3=25, I should get the same results since the relative inclination, i+i3, is unchanged. However, I am getting different results when I do this. I attached a pdf of my ipynb. Any help in understanding where I'm going wrong is much appreciated, thanks!

Rebound_Triple_Evolution.pdf

[hannorein]

I don't see anything wrong with your setup. You are plotting the sum of the inclinations. Is that what you've intended? I could be mistaken, but isn't Figure 3 just showing one of the inclinations?

[landenconway]

You're correct, I've plotted the relative inclination where they just plot i. However, both plots that I'm producing should still be the same since the relative inclination is unchanged. The issue is still the same if I only plot i, the plot for i=40 i3=25 doesn't match the plot for i=64.7 i3=0.3

[hannorein]

The results are different because you're integrating two different systems. One is not just a rotated version of the other. (orbital elements can be counter-intuitive, I suggest you make a few plots of how the two systems looks like)

[landenconway]

I believe these two systems are just rotated versions of one another, they would be related by a -24.7 degree rotation about the x-axis. I show this in the first cell of my linked screenshot. The dynamical equations that govern the evolution of each orbital element should depend only explicitly on the mutual inclination, i+i3, so changing the individual inclinations but maintaining the mutual inclination should return the same dynamics. After evolving these two systems the final state vectors are no longer rotationally related, they have different magnitudes in separation and velocity. This is shown in the second cell. Additionally, plotting the orbital elements over this integration reveals that 4 of the orbital elements are evolving differently, the inclinations and both pericenter angles. The outer binary's pericenter and inclination appears to be not be evolving at all which is even more problematic.

I apologize for all the questions, but I just can't seem to figure out where I'm going wrong in thinking about this.
triple_evolution_issue.pdf

[hannorein]

Fair enough. Then I'm not sure. Is the system just very sensitive to initial conditions? Try running two almost identical simulations that only differ by a small amount initially, say ~1e-8.

Minor: your call to sim2.integrate(...) uses sim1.particles[1].P which will have changed from its initially value after you call sim1.integrate.

[landenconway]

It doesn't seem like its sensitive enough to make any substantial difference with variations at that order. I also found something else peculiar in the inclinations. After defining the system, if I look at the inclination using .inc vs calculating it from the r's and the v's with angular momentum I get two different results. I linked a screenshot of the code showing that.
inclination_issue.pdf
Could you point to where in rebound's code the orbital elements are defined and calculated?

[hannorein]

By default REBOUND is using Jacobi coordinates both when you add particles and when you query for orbital parameters. Maybe you want to use heliocentric coordinates instead? You might want to look at this tutorial: https://rebound.readthedocs.io/en/latest/ipython_examples/OrbitalElements/

PS: It would easier to look at the code if you were to just include it in the markdown instead of attaching a pdf file.

[landenconway]

I apologize for the difficulty of the pdfs, I'm not the most familiar with git or markdown so I have little experience in this area. However, I understand that REBOUND is using Jacobi coordinates, that too is what I use in my work (secular dynamics of triple systems) which makes REBOUND attractive for me to use. I really want to continue to incorporate REBOUND into my work but there just seems to be these inconsistencies that I'm unsure of how to make sense of.

[hannorein]

Note that you're working in heliocentric coordinates when you calculate the inclination using the angular moment. REBOUND will give you the inclination in Jacobi coordinates (unless you tell it otherwise). The routine which calculates orbital elements is here: https://github.com/hannorein/rebound/blob/main/src/tools.c#L1040

It's complicated because it deals with all kind of edge cases. Also, when you call the convenience functions from python, then the primary object will be set according to Jacobi coordinates (again, unless you tell it otherwise). See https://github.com/hannorein/rebound/blob/main/rebound/particle.py#L422

[acpetit]

Hi, jumping on this issue. I agree that both systems are just similar up to a rotation, however, the quantity that you plot (i1+i3) is not the planets mutual inclinations because the xy plane is not the invariant plane for any of these systems.

[hannorein]

Thanks, @acpetit!

[landenconway]

Hi @acpetit , I appreciate the feedback! As I understand it, the inclination of a specific binary orbit is defined as the angle between that binary's angular momentum and the total angular momentum of the system. I agree that the coordinate system that is being used doesn't coincide with the invariant plane, but that shouldn't matter if the inclinations and the corresponding mutual inclination is determined from the angular momentum. The mutual inclination should be invariant of whatever coordinate system you use to calculate it. Additionally, that doesn't remedy the fact that two initial state vectors that are rotationally similar are getting mapped to non-rotationally similar final state vectors when the transformation itself is rotationally invariant.

[acpetit]

The mutual inclination should be invariant of whatever coordinate system you use to calculate it.

Yes but if you are not in the invariant reference frame, the mutual inclination is not given by $i_1+i_3$ instead $$\cos i_{mut}= \cos(i_1)\cos(i_2)+\sin(i_1)\sin(i_2)\cos(\Omega_1-\Omega_2)$$

[landenconway]

Understood. However, I've been using Omega_1 = 0 and Omega_2 = Pi in the setup for my simulations, so that the mutual inclination is just the sum of the two inclinations. Which I think is actually aligning the x-y plane with the invariant plane.

[acpetit]

No. Your invariant plane is slightly inclined in both cases with respect to the xy plane. You should also account for the change of Omega 1 and 2 during the simulation to compute the final mutual inclination. Then as Hanno told you, the slight discrepancies between the values computed using the positions and velocities and from the builtin orbital element functions come from choices of coordinates.

[hannorein]

I'm closing this for now but feel free to reopen if there are additional questions.