Adding 'j2' to Phoebe

[miroslavbroz]

In Phoebe (https://www.phoebe-project.org/), it is necessary to model stellar systems, which are massive, compact, and with rotating components. Generally, this problem is too complex, because each star is distorted by a Roche-like potential (Horvath et al. 2020). Nevertheless, a useful approximation might be that each star is described by oblateness (J2), which acts preferentially on the closest components. Each pair of stars is thus assumed to have the rotatonal axes perpendicular to the orbital plane (e.g., Fabrycky 2010). To this point, one can use an alternative implementation of the J2 force shown here.

[miroslavbroz]

Note: The 'gravitational_harmonics' implementation is not suitable, because the z-axis is fixed. On the contrary, results for i = 0 and 90 deg binaries should be the same. Closely related to https://github.com/miroslavbroz/phoebe2/commit/b5832f03550d91965c5913b29c794763ee944ab9 .

[dtamayo]

Thanks! I've been wanting to generalize the gravitational_harmonics implementation for a long time. I think it would make a big difference for REBOUNDx users if we can find a way to generalize a single implementation of j2/j4 etc (that you would be welcome to be listed as lead author on) rather than having several different effects that do slightly different things.

It would certainly be nice to be able to input a spin vector and have the j2/j4 effects aligned with that spin axis. But in the setup you describe with two stars with rotation axes aligned with the orbit normal, why couldn't one align the orbit normal with the z axis of the simulation and then use the existing implementation? Is the issue that other effects cause the orbit to evolve and you would like the spin axis to follow the orbit normal?

[miroslavbroz]

Dear Dan, it is true that a more general implementation would be ideal! As an intermediate step, I tried to implement this simplified J2, because we definitely need it for Phoebe. In our Keplerian model, we have domega/dt as a free parameter. In our N-body model, we have perturbations from N-1 bodies, relativistic corrections, but we also need oblateness due to internal structure of stars. Otherwise, all simulations would be systematically offset. :-(

The reason why current J2 cannot be used is that we are aiming at misaligned multiple systems (hierarchical, 2+2, ...). I cannot rotate all binaries (one by one) and integrate them separately, everything else (perturbations, relativity) would be lost and the code would be 'horrible'. This is (by far) the simplest solution. One way how to include it in Phoebe is to copy the rebound/reboundx and include it as a modified version. Hmm, I do not like it, because of the maintenance.

Whether a full J2, J4 implementation (w. all axes of all *) is overly complicated or not? I do not know, but...

Clear skies, Mira

[miroslavbroz]

... it is not complicated:

https://github.com/miroslavbroz/reboundx/blob/j2_/src/gravitational_harmonics.c

If you will agree, I shall rewrite the rest (J4, ...). And then submit a separate PR.

Clear skies, Mira

[hannorein]

I think it would be nice to make it backwards compatible: if no vector is given, then the default is the z-axis.

[miroslavbroz]

Done (see #139 ). Hmm, default Omega (z-axis) is at two places.

[miroslavbroz]

Ops! The back-transformation was erroneous. Corrected, so xy, xz, yz planes (z, y, x axes) work the same.