Adding 'j2' to Phoebe (generalisation)

[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 the J2, which acts on all other components. Each component has its own rotation axis, with respect to which the acceleration must be computed (and transformed back to the inertial frame).

[gitnotebooks]

Found 1 changed notebook. Review the changes at https://app.gitnotebooks.com/dtamayo/reboundx/pull/139

[dtamayo]

Hi Mira,

Thank you so much for all these improvements. One quick question:

This should work if I don't line up the spin vector with the z axis and have inclined orbits, right? I tried taking the python example, changing the orbit to inc=0.5:

sim.add(m = 3.e-6, a = 1., e = 0.01, inc = 0.5)

and changing Omega to

sim.particles[0].params["Omega"] = [0.0, 0.5, np.sqrt(3)/2]

The integration then gives answers that disagree with the predictions from Oplistilova below. Is this an issue with the code, or do the predictions assume that the z axis aligns with Omega?

Dan

[miroslavbroz]

Dear Dan, I think it is all about the reference plane and projection effects. In Oplistilova, it is assumed to be the equator of the central body. Precession (of Omega, omega) then looks like a circulation (from 0 to 360 deg). In the numerical integration, it is xy; all orbital elements are computed w.r.t. xy. However, if the spin is tilted, the precession is projected to the xy plane and only then the elements are computed. Precession then might look like a libration (i.e., an oscillation around an arbitrary value). Two possible solutions: i) add a warning; ii) a transformation of the coordinate system to the equatorial plane prior to using Oplistilova (but this is equivalent to having spin aligned with the z-axis). Of course, i) is much simpler than ii). Clear skies, Mira

[dtamayo]

Hi Mira,

I think that should be the expected behavior, so I don't think we even have to add a warning...I was just wondering, given that the main change is allowing for a tilted spin axis, whether there's a way to test the outputs of the code in the tilted case to verify that it is working as expected.

Could we test it together with the updated potential function by checking that the total energy is conserved to machine precision?

Dan

[miroslavbroz]

Regarding the conservation, it is actually better than the previous implementation: https://sirrah.troja.mff.cuni.cz/~mira/tmp/reboundx/test_conservation/energy.png

[miroslavbroz]

Regarding the spin axis, I have tests for xy, xz, yz planes: https://sirrah.troja.mff.cuni.cz/~mira/tmp/reboundx/, see directories test_xy, test_xz, test_yz, please. All of them exhibit the same precession rate.

[dtamayo]

Hi Mira,

I've gone through the J2.ipynb example and added some calculations to match up the tilted case. Just sent you a pull request. It's all a bit tricky, so hopefully having it all in the example will help people out and lead to fewer questions on how to use the effect and interpret the results.

Thanks so much for this, it's great to finally have this implemented. We put in some functionality in REBOUND to do those kind of rotations using quaternions, but not worth changing something that's already working.

Let me know if the changes look good to you. If so, once you merge it should show up here and we can merge into main and update the version number.

Thanks again,

Dan

[miroslavbroz]

Dear Dan, a very nice analysis, indeed! Thank you too for checking the J2 computations independently. Clear skies, Mira