Kepler orbital elements

A binary’s relative orbit is fully specified by seven orbital elements: the semi-major axis $a$, eccentricity $e$, inclination $i$, longitude of the ascending node $\Omega$, argument of periastron $\omega$, mean anomaly $M$ at epoch, and the orbital period $P$ (or equivalently, the time of periastron passage $t_p$). progenax stores these in the KeplerElements PyTree and provides closed-form mappings between elements and phase-space coordinates $(\mathbf{r}, \mathbf{v})$ for each component star.

This chapter derives the elements-to-phase-space conversion, lists the conventions progenax uses, and documents the differentiability properties that make KeplerElements HMC-compatible.

The seven elements¶

The first six are positional (they specify where the orbit is in space); the seventh specifies the temporal phase. Together they fully determine $(\mathbf{r}_{\mathrm{rel}}(t), \mathbf{v}_{\mathrm{rel}}(t))$ for the relative orbit.

Kepler’s third law¶

The orbital period and semi-major axis are linked by the total mass:

↗ model card

progenax exposes both directions:

• compute_period(a, m_total, *, G) — given $a$, returns $P$.

• period_to_semimajor_axis(P, m_total, *, G) — given $P$, returns $a$.

The choice of $a$ vs $P$ as the parameterised quantity matters in practice: most observational period distributions (Binary period distributions) parameterise $P$ directly, so progenax’s default KeplerElements constructor accepts $P$ and computes $a$ internally via (1).

From elements to phase-space coordinates¶

Converting orbital elements to relative position $\mathbf{r}$ and velocity $\mathbf{v}$ at time $t$ requires four steps:

Step 1 — solve Kepler’s equation for the eccentric anomaly $E$:

↗ model card

This is a transcendental equation in $E$. progenax uses a fixed-iteration Newton solver built on jax.lax.scan (not a Python loop), so the iteration count is static and the whole solve is one JIT-/grad-compatible trace:

import jax
import jax.numpy as jnp

def solve_kepler(M, e, max_iter=50):
    M_wrapped = jnp.mod(M, 2.0 * jnp.pi)
    # Initial guess E0 = M + 0.85 sign(sin M) e (good for moderate e)
    E0 = M_wrapped + 0.85 * jnp.sign(jnp.sin(M_wrapped)) * e

    def newton_step(E, _):
        f = E - e * jnp.sin(E) - M_wrapped
        f_prime = 1.0 - e * jnp.cos(E)
        E_new = E - f / jnp.maximum(jnp.abs(f_prime), 1e-12)
        return E_new, None

    E, _ = jax.lax.scan(newton_step, E0, None, length=max_iter)
    return E

The Newton iteration converges quadratically — typically 3–4 steps for $e \lesssim 0.8$. progenax’s default max_iter=50 (KeplerElements._solve_kepler_equation) is deliberately conservative so the fixed-iteration scheme still reaches machine precision in the slow-converging $e \to 1$ regime; the extra cheap iterations guarantee accuracy across all $e < 1$.

Step 2 — eccentric anomaly to true anomaly $\nu$:

Step 3 — orbit-plane coordinates. Position and velocity in the orbit plane (with periastron along the $x$-axis):

Step 4 — rotate into the reference frame via three Euler-angle rotations (about $z$ by $\Omega$, then $x$ by $i$, then $z$ by $\omega$):

The full pipeline is one JIT trace and fully vmap-able over a population of binaries.

From relative to component coordinates¶

Given $\mathbf{r}_{\mathrm{rel}} = \mathbf{r}_2 - \mathbf{r}_1$ and component masses $m_1, m_2$, the individual particle coordinates are

↗ model card

with the corresponding velocity decomposition for $\mathbf{v}_1, \mathbf{v}_2$. These are then added to the binary’s centre-of-mass position and velocity (which come from the spatial profile and velocity DF — see Binary populations) to produce the absolute particle coordinates.

Sampling: isotropic angles + population distributions¶

For a Monte Carlo binary population, progenax samples the seven elements as follows:

1. $P$ from the population period distribution (Binary period distributions) — typically a log-normal or Sana et al. (2012) empirical fit.

2. $a$ from $P$ and $m_{\mathrm{tot}}$ via (1).

3. $e$ from the population eccentricity distribution (Eccentricity distributions) — thermal, uniform, or Moe & Di Stefano (2017) period-dependent.

4. $\cos i \sim \mathcal{U}(-1, 1)$ — isotropic inclination.

5. $\Omega \sim \mathcal{U}(0, 2\pi)$ — uniform.

6. $\omega \sim \mathcal{U}(0, 2\pi)$ — uniform.

7. $M_0 \sim \mathcal{U}(0, 2\pi)$ — uniform orbital phase.

The mass ratio $q$ comes from Mass-ratio distributions; the binary fraction from Multiplicity statistics (Moe & Di Stefano 2017). The seven orbital elements plus $q$ specify the orbit completely.

Differentiability properties¶

KeplerElements is differentiable in all element values. Specifically:

• $\partial\mathbf{r}/\partial a$, $\partial\mathbf{r}/\partial e$, $\partial\mathbf{r}/\partial(i, \Omega, \omega, M_0)$, and $\partial\mathbf{r}/\partial P$ all flow analytically through (4)–(5).

• The Newton iterations on Kepler’s equation (2) are a fixed-count jax.lax.scan, gradient-compatible.

This matters when fitting binary orbits — e.g. inferring $(a, e, i, \Omega, \omega, M_0)$ from radial-velocity or astrometric observations of a known binary. The full inference loop is one JAX gradient call.

Domain of validity¶

1. Bound orbits only ($e < 1$). Hyperbolic orbits ($e \ge 1$) represent unbound encounters, not binaries; they are not handled by KeplerElements.

2. Two-body approximation. Higher-order multiples (triples, quadruples) are not represented by a single KeplerElements instance. progenax models them as nested binaries (an outer KeplerElements whose “secondary” is itself an inner binary).

3. Newtonian gravity. Relativistic effects (periastron precession, orbital decay via gravitational waves) are not included. For binaries with $a \ll 1$ AU and component masses $\gtrsim 10\,\Msun$, relativistic corrections become observable on Gyr timescales.

Implementation, validation & references¶

• In code: src/progenax/binaries/kepler.py (KeplerElements, the fixed-iteration Kepler-equation solver, and the elements→phase-space conversion), with compute_period / period_to_semimajor_axis in src/progenax/binaries/kepler_period.py. See the binaries API.

• Validated in: binary-aware recovery — the binary regression suite that exercises the orbit machinery (Kepler’s third law to $\sim 10^{-10}$).

• Primary sources: Kepler-element machinery is standard textbook material; Murray & Dermott Solar System Dynamics §2 gives a clean derivation, and modern N-body codes (NBODY6, COSMIC) use the same conventions. progenax’s KeplerElements is JAX-native via Kidger & Garcia (2021)’s PyTree pattern; the Newton solver follows the standard fixed-iteration approach.