[Question]:Transition guide for HNBody users: How to implement Earth-Moon quadrupole attraction?

[astrotuna201]

Hi, first: many thanks for making such nice open-source tools (both Rebound + ReboundX),

I am not an astronomer, but user of a different package (HNBody) so far, mainly doing long-term solar system integrations.

I would like to migrate to an open-source package for reproducibility and ease of hacking, and managed to follow the examples etc., to approximate HNBody's setup (e.g. including postNewtonian GR corrections and sun J2 corrections) ...

HNBody has a some drivers that implements either an Earth-Moon sub integration (ems), or a simpler "lunar" driver, but I am not sure how to translate that into one of the Rebound(X) J2/J4 examples, which seem to be related. Do you have any suggestions?

A separate question about HNBody vs Rebound: I think HNBody implements Kahan Summation in the integration steps, claiming reduced round-off by several orders of magnitude. I looked at the Rebound code and it seems everything is coded in terms of doubles --- would this make a big difference or is this addressed by other means?

Many thanks!

[hannorein]

Hello 🔭🐟,

I'll let Dan answer the questions about the additional forces. Regarding Kahan Summation: REBOUND does include a module for compensated summation. You need to set the gravity method to compensated. But I suspect this will not be useful for your Solar System integrations. Typically other errors dominate (accuracy of ephemeris, etc).

Hanno

[astrotuna201]

Many thanks for the quick answer re. Q2! Amazing work!

[dtamayo]

Hi,

I'm not familiar with HNBody's drivers, so I can't give very concrete advice on that front. One thing to decide is whether you need the moon and Earth in the integration separately. If you don't, then you can to leading order lump their mass at their mutual center-of-mass. This has the advantage of allowing you to take much bigger time steps and use much faster algorithms (WHFast). If you want to include the moon in the integration, then you need to use IAS15 or Mercurius.

As for the Sun's J2, this should be negligible compared to GR, so you should be fine completely ignoring it. The bigger correction if you're putting in the Earth and Moon as a single mass, is that the next order correction looks like a J2 term. Perhaps this is what HNBody is doing?

If you wanted to put this in, you'd have to get several things right. If you think this effect matters I would suggest putting Earth and Moon in separately and just using IAS15 and paying the computational price.

[astrotuna201]

Thanks!

I am actually trying to match the orbital calculations of e.g. Laskar et al. (2004) (paper, data) as best as I can ... Recent papers by Zeebe (e.g. 2017) managed to match for example Earth's eccentricity over the past 50Myr very well, using HNBody. To get a good match, the solar J2 term is important (beyond GR), but the Earth-Moon system can be treated as a a point source if the quadrupole moment of the EM system is included ... (e.g. solution ZB17e in above '17 paper) –– I think this is what you mean by "looks like a J2 term"?

If I use "gr", identical initial positions, masses + velocities, and solar J2, I get below results (black), compared to the near identical La2004 and ZB2018 solutions over 30 Myr (yellow/blue). ...

So, I think your are right, but my questions is then how to adapt the J2 examples in ReboundX to match the lunar.c terms ...

comparison.pdf

My current problem.c is also below ... after doing a coordinate rotation to the solar rotation axis (see Zeebe 2017 and Beck & Giles)

problem.c.txt

[dtamayo]

Interesting! I didn't appreciate that the Sun's J2 term matters. And yes, that is what I meant by the J2 term. If both those J2 terms matter, one problem is that you can align your z axis with the sun's rotation axis, but then your EM J2 will be misaligned. The right way to do this would be to add two parameters to the gravitational_harmonics.c implementation specifying the spin axis direction, and to implement rotation matrices or quaternions to do the rotation. Then every body can have its own orientation.

Looking at that lunar.c it doesn't look like HNBody is accounting for the obliquity between the sun's spin axis and the ecliptic, so it's probably negligible since the perturbation is so small. So then it's just a matter of adding a J2 and radius to the Earth-Moon particle. The relevant equation for the effective J2 is

J2 R_earth^2 = 1/2 a_moon^2 M_moon / Mearth

So try setting R_earth, and calculating the effective J2 as above. Let us know how it goes!

[astrotuna201]

Many thanks! I will see how I get on. The solar J2 limits the solutions to about 40Myr ... but the Earth-Moon system kicks in earlier. The solar J2 is included in HNBody as a separate parameter for the dominant mass; the one in the lunar.c driver is for the EM system only ... Anyway, I'll close this for now, and report back; again many thanks!

[dtamayo]

In REBOUNDx, once you add the gravitational_harmonics effect, you can add a J2 (and radius) to any of the particles in the integration