possible inconsistency in orbital elements

[parkus]

This a fantastic package! Playing around with it, I noticed the orientation of an eccentric orbit is not what I would expect based on the definition of the orbital elements in the documentation. An example:

import exoplanet as xo
import numpy as np
from matplotlib import pyplot as plt
plt.ion()

# make an edge-on highly eccentric orbit with
# 45 deg argument of periapsis
# orbital plane is the X-Z plane
orbit = xo.orbits.KeplerianOrbit(
    period=200,
    t_0=0.0,
    incl=np.pi/2,
    ecc=0.8,
    omega=np.pi/4,
    Omega=0.0,
    m_planet=0.0,
    m_star=1.0,
    r_star=1.0,
)

t = np.linspace(0, 200, 200, endpoint=False)
x, y, z = orbit.get_planet_position(t)

plt.figure()
plt.scatter(x.eval(), z.eval(), c=t)
plt.xlabel('X')
plt.ylabel('Z (+Z is to observer)')
ax = plt.gca()
ax.set_aspect(1)
plt.grid()

Yielding this figure

The docs state that the ascending node is defined to be the point at which the planet passes through the XY plane moving away from the observer, who is looking from the +Z axis. The colors on the plot progress from purple to yellow as time progresses, so the ascending node appears to be at X = 120ish. The argument of periastron is defined as the angle from the ascending note to the periastron in the orbital plane. From the figure it appears has been referenced against the -X axis, opposite the ascending node.

Is this an inconsistency or have I misinterpreted something?

[dfm]

Pinging @iancze

[iancze]

I think this might be an issue of what t_0 corresponds to. If you set t_periastron=0 in xo.orbits.KeplerianOrbit, do things make sense? I think the t_0 is for a reference transit, which may not make sense for a non-transiting orbit (if that is what you're working with).

[parkus]

I wondered the same thing initially, but no. That changes where the planet "starts" in its orbit, but the argument of periastron still seems to be referenced against the node where the planet is "approaching" the observer (moving from -Z to +Z). Maybe I've got my brain flipped?

Edit: added screenshot

[iancze]

I think I might know what be going on... the orbital convention figure should probably be updated. I am assuming that you are talking about this one defined in the data and models tutorial?

In that figure Omega is something like 210 degrees, so an Omega = 0 orbit would indeed be referenced omega from the +X axis (defined as the ascending node). The dashed angle arcs are confusing... it was my best attempt to show the three dimensional nature of the orbit using my limited Inkscape skills at the time, and I think it could easily be misinterpreted as Omega ~ 50 degrees. Hopefully this sorts this as an issue with the documentation and not the underlying conventions of the codebase? Definitely interested to hear if this sorts it for you.

We're mostly following Murray and Correia here, and some more information on these orbital conventions is contained in this document. I've tested the equations as written against some other RV and astrometric timeseries in refereed papers (cited in the document) and the orbital elements checked out (though see the discussion about some papers leaving ambiguity between ascending/descending node and omega_1 vs omega_2). That's not to say that there couldn't be a bug somewhere, though!

[parkus]

That's correct, that definition of the elements on the data and models page (which, by the way, I really appreciated as I was struggling to find anything instructive in the literature) is the foundation of my question.

Yes, I think the issue could just come down to updating the docs. The figure looks to be consistent with the text on the data and models page, which states

In the stellar and exoplanetary fields, there is less agreement on which segment of the line of nodes is defined as the ascending node. We choose to define the ascending node as the point where the orbiting body crosses the plane of the sky moving away from the observer (i.e., crossing from 𝑍 > 0 to 𝑍 < 0)

However, the orbit I plotted seems to be inconsistent with both that statement and the figure.

Actually, the origin of the issue might be revealed in the draft you shared. I see there you state

our main break from Murray & Correia (2010) is that we designate the ascending node as the position on the line of nodes where the secondary star is moving away from the observer.

If you coded things up to be consistent with Murray & Correia, then indeed the ascending node would be where the planet crosses the reference plane approaching the observer.

[iancze]

Thanks for your help digging into these inconsistencies.

We didn't completely follow the conventions in Murray and Correia. In particular, we explicitly chose the ascending node convention (orbiting object crosses the sky plane moving away from observer) to be consistent with the visual binary field, which makes up the largest marketshare of astrometric orbital datasets and provided some of the well-sampled benchmark datasets we were working with at the time. One could argue that we should be using the alternative ascending node convention since the package is called "exoplanet," but I'm not entirely sure whether the definition chosen by Murray and Correia is established as the standard in the exoplanet field---I've seen both node conventions used in direct imaging papers. More information leading to the convention choices is contained in this issue.

As far as I know, the text in the data and models tutorial is correct in describing the rotations (as a function of Omega, inclination, and omega) that are required to orient an orbit in 3D space as shown in the figure (and as implemented in the code). The equations and figure apply equally to whether you are orienting the orbit of the star (primary star) or the orbit of the planet (secondary star). For a two-body orbit, most of the orbital parameters are the same for the star and planet (e.g., e, Omega), but omega_1 = omega_2 + pi (in radians).

However, I think there is some ambiguity around which omega is required by the KeplerianOrbit routine. My understanding is that it is omega_1, i.e., the argument of periapse of the stellar (primary star) orbit.

It would be helpful if others could confirm, but I've seen it common to list omega = omega_1 for radial velocity orbits (i.e., omega corresponds to the star orbit, just like semi-amplitude corresponds to the stellar motion) and omega = omega_2 for astrometric orbits (i.e., omega corresponds to the planet orbit measured relative to the star). For joint RV + astrometric orbits I have seen a (confusing) blend of conventions that seem like omega defaults to whichever dataset the code was written to fit first (RV or astrometric). I am not familiar enough with transit orbits to know what the convention there might be.

Following your example, I made a few plots of X, Y, and Z (our convention is that capital letters refer to the sky-frame and lowercase letters refer to the perifocal frame).

orbit = xo.orbits.KeplerianOrbit(
    period=200,
    t_periastron=0.0,
    incl=0,
    ecc=0.8,
    omega=0.0,
    Omega=0.0,
    m_planet=0.0,
    m_star=1.0,
    r_star=1.0,
)

t = np.linspace(0, 200, 200, endpoint=False)
X, Y, Z = orbit.get_planet_position(t)

Xt, Yt, Zt = X.eval(), Y.eval(), Z.eval()

fig, ax = plt.subplots(ncols=3, figsize=(9, 4))
ax[0].scatter(Yt, Xt, c=t)
ax[0].set_xlabel("Y")
ax[0].set_ylabel("X")
ax[0].invert_xaxis()
ax[0].set_title("Sky")
ax[0].set_aspect(1)

ax[1].scatter(Xt, Zt, c=t)
ax[1].set_xlabel("X")
ax[1].set_ylabel("Z")


ax[2].scatter(Yt, Zt, c=t)
ax[2].set_xlabel("Y")
ax[2].set_ylabel("Z")
for a in ax:
    a.grid()

fig.subplots_adjust(left=0.05, right=0.95)

Here is a repeat of your example with zeroed out geometric parameters, Omega = 0, omega = 0, and incl = 0, which means X=x, Y=y, and Z=z. We see that periapse of the planet is actually at the -X line, because exoplanet is interpreting this as omega_1 = 0, omega_2 = np.pi.

But if we specify omega=np.pi, we get

And then for an inclined orbit (still with Omega=0), with omega_1=200 degrees and omega_2 = 20 degrees we have the following plot of the planet orbit, which I believe correctly follows the conventions in the documentation, including the ascending node convention.

Please let me know if this clarifies the issue that you are having.

[parkus]

Likewise, thank you for putting in so much effort to trying to get to the bottom of this!

We didn't completely follow the conventions in Murray and Correia. In particular, we explicitly chose the ascending node convention (orbiting object crosses the sky plane moving away from observer) to be consistent with the visual binary field,

Right, I think the "data and models" page makes this perfectly clear, and I think the choice is well justified, so no issues there.

However, I think there is some ambiguity around which omega is required by the KeplerianOrbit routine. My understanding is that it is omega_1, i.e., the argument of periapse of the stellar (primary star) orbit.

AHA! That is the answer to my confusion. I was definitely under the impression that the omega specified when instantiating a KeplerianOrbit was that of the secondary object (i.e., the planet). This statement from the data and models page

We choose to define the ascending node as the point where the orbiting body crosses the plane of the sky moving away from the observer (i.e., crossing from 𝑍 > 0 to 𝑍 < 0)

is the source of that.

It appears then that some of the inputs to KeplerianOrbit are for the primary and others for the secondary. The last example above is consistent with omega being omega_1. However, Omega seems to be Omega_2. I.e., the secondary's ascending node (crossing from Z > 0 to Z < 0) happens on the +X axis at Omega=0. Would you agree, @iancze?

[iancze]

Great! I'm glad we identified the source.

Re: Omega, I believe this variable should be the same for both the primary and secondary, since both objects will orbit in the same sense (clockwise vs. counterclockwise) and therefore their ascending nodes will be on the "same side" of the orbit w.r.t. center of mass, and thus have the same position angle. I think omega and a are the main parameters to KeplerianOrbit that are not the same between star and planet (and are not otherwise marked as such).

[parkus]

Re: Omega, I believe this variable should be the same for both the primary and secondary, since both objects will orbit in the same sense (clockwise vs. counterclockwise)

Right you are!

In that case, it seems the solution is either updating the docs to be explicit that omega refers to the primary or updating the code such that the omega users supply to KeplerianOrbit is omega_2. Actually, either way it would be helpful to newbies to exoplanet orbits like me if the docs were explicit as to which omega KeplerianOrbit takes. I would argue for changing KeplerianOrbit to take omega_2, as that is more intuitive to me as a newcomer, but I know backwards compatibility is also desirable.

[iancze]

https://arxiv.org/abs/2212.06966 relevant to this conversation and perhaps explains some of the variation in the literature previously alluded to.