Orbital dynamics¶

This page describes the orbital dynamics included within PROTEUS, how its config options combine, and how the physics is validated. For the config field tables themselves, see Star and orbit configuration. For the one-paragraph summaries of each tidal module, see Tidal evolution: Obliqua, Lovepy and Orbital evolution: PROTEUS (internal).

Two independent evolution families¶

PROTEUS evolves an orbit around exactly one body at a time:

These are mutually exclusive (a config error is raised if both are set). A satellite can still be tracked (orbit.satellite.include_satellite) without its own evolution model, in which case its semi-major axis and eccentricity stay fixed at their configured initial values while the planet's orbit around the star evolves independently. The reverse also works: planet_satellite_model can evolve the satellite while the planet-star orbit stays fixed at its initial semimajoraxis/eccentricity.

Tidal response modules (orbit.module)¶

The tidal module supplies two things every evolution model needs: a heating profile to add to the interior's energy equation, and a Love-number description of the dissipation (hf_row['Imk2'] and/or the full mode spectrum in tides_o).

Star-planet models (orbit.star_planet_model)¶

sp0d and sp1d integrate with scipy.solve_ivp (orbit.solver.* controls method and tolerances). parameterized solves nothing: after time_migration it sets the orbit from a closed form in time, with parameters chosen by the user.

Configured under [orbit.parameterized]:

The track starts from a_0 = orbit.semimajoraxis and orbit.eccentricity, which seed the orbit at the initial condition as for every other star-planet model; the track is evaluated from the first step after it. The initial condition spans Time <= 1 yr, so time_migration must be greater than 1 yr. Setting the orbit from a target flux (orbit.instellation_method = 'inst') is rejected at config load, since it would give the run a second starting orbit.

A law writes only the orbital elements it sets: instant and sigmoid the semi-major axis, high_ecc the semi-major axis and the eccentricity. none, and every law before time_migration, leave the orbit as the previous step left it.

sigmoid holds the orbit until time_migration, carries it to sma_final over the following tau_migration along the cubic 3u^2 - 2u^3, and holds it there afterwards.

high_ecc circularises at constant orbital angular momentum (Postolec et al. 2026)⁷: it excites the eccentricity to sqrt(1 - sma_final / a_0) at the migration epoch and then decays it, until the orbit reaches sma_final. The eccentricity jumps discontinuously at time_migration from orbit.eccentricity to its excited value. That step is physical, since a scattering or Kozai event is fast compared with the orbital evolution that follows.

The migration window must also be resolved by the timestep. sigmoid and high_ecc are sampled wherever the coupled loop happens to step, and nothing aligns a step to time_migration. If tau_migration spans fewer than three timesteps the track is sampled at little more than its endpoints and silently degenerates to instant; the orbit module logs a warning when that happens.

Spin, stellar flux and impacts on a prescribed track¶

No torque acts on the spin, so the planet stays synchronous: axial_period is set to the current orbital period at every step, and AGNI and the breakup check read that value. orbit.axial_period must be unset, for sp0d as well, which evolves no spin either.

On an eccentric orbit the stellar flux is averaged over the orbit as <1/r^2> = 1 / (a^2 sqrt(1 - e^2)), which is the flux at the distance a (1 - e^2)^(1/4). Every module that scales a flux by distance uses that distance: the bolometric and XUV instellation, the stored stellar spectrum, the eclipse depth, VULCAN's star.dat and orbit_radius, and the stellar-surface flux petitRADTRANS recovers from the stored spectrum. F_ins is refreshed every params.dt.starinst and the stored spectrum every params.dt.starspec. While the orbit evolves, petitRADTRANS and VULCAN undo the latest stored spectrum at the current distance, so their stellar flux is off by the square of the ratio of the current distance to the one the file was written at. AGNI is not affected: it takes only the spectral shape from the file and its heating from F_ins. The time-averaged separation a (1 + e^2 / 2) stays in use for geometry only (the Roche-limit checks and the orbit plots). This applies to every eccentric run, not only to parameterized.

Giant impacts are rejected at config load (accretion.module = 'none' is required): after time_migration the track sets the semi-major axis from its own parameters every step, so an impact's new semi-major axis would be overwritten while its change in eccentricity persisted.

Where the orbital energy goes¶

The track changes the orbital energy E = -G M_star M_planet / (2 a) but deposits that energy nowhere: no tidal heating reaches the interior and the energy balance of the planet does not include it. The helpfile column dEdt_orb [W] records the rate the track implies, dE/dt = G M_star M_planet (da/dt) / (2 a^2), negative while the orbit shrinks. It is zero for none, outside the migration window and for instant, whose step releases its energy at a single time.

For high_ecc the rate is largest in magnitude at time_migration, where it equals 2 dE / tau_migration, with dE = E(sma_final) - E(a_0) = -(G M_star M_planet / 2) (1/sma_final - 1/a_0) the whole energy change (negative for inward migration), and it then decays as exp(-2 (t - time_migration) / tau_migration). The eccentricity step at time_migration changes the orbital angular momentum instantly while leaving the energy unchanged, since a is still a_0 there.

Visualizing the four parameterized regimes¶

Each regime was run as a dummy PROTEUS simulation and compared against the closed form in src/proteus/orbit/parameterized.py:

Semi-major axis (top) and eccentricity (bottom) for the four regimes, with orbit.semimajoraxis = 2.0 au, sma_final = 0.8 au, time_migration = 1e3 yr and tau_migration = 1e4 yr. The dashed vertical line marks the migration epoch and the shaded band spans one tau_migration after it. The dotted horizontal lines in the top panel mark the starting and final orbits, a_0 = orbit.semimajoraxis and a_f = sma_final.

Planet-satellite models (orbit.planet_satellite_model)¶

Compatibility between orbit models and tidal modules¶

A tidal module makes up to two things available: the scalar hf_row['Imk2'], and/or the full per-mode spectrum in tides_o. Which one an orbit model reads is exactly what its 0d/1d suffix tracks -- a 0d model reads the scalar path, a 1d model reads tides_o directly.

What each tidal module provides:

What each orbit model reads:

Satellite Love-number lookup (Obliqua only)¶

Unlike the planet, whose interior structure evolves and is re-queried every coupling step, the satellite's interior is treated as static for the lifetime of a run. orbit.obliqua.lookup_from_interior builds a full frequency-spectrum Love-number table once, from a fixed satellite interior description (orbit.satellite.love_number_sat, a JSON initial condition read by a simplified 0-D solid/fluid Obliqua configuration), and writes it to a NetCDF file (sat_tides.nc). Alternatively, the user can provide their own pre-computed table (orbit.satellite.love_number_sat, a NetCDF file), which will be used instead of the one generated by lookup_from_interior. Every subsequent coupling step, LN_from_lookup computes the satellite's own forcing frequencies from its current spin and orbital state and interpolates the satellite's Love numbers from that fixed table (linear in frequency, per tidal degree).

Evection resonance (ps1d_evec)¶

Evection resonance happens when a moon’s elongated orbit rotates at the exact same speed that the central planet orbits its star.; capture into it can pump the satellite's eccentricity well above what tides alone would produce. ps1d_evec detects proximity to the resonance location a'_res (Rufu & Canup 2020, Eq. 12) with a debounced, hysteretic band detector (separate entry/exit margins, orbit.solver.resonance_margin_enter/resonance_margin_exit, avoid chattering at the band edge) and gates only the oscillating resonant forcing term on that detector. The secular apsidal-precession term and the evection angle's own evolution are always active regardless of band status. Setting the gate to zero decouples the resonant forcing term, reducing ps1d_evec to plain ps1d dynamics.

While in or near the band, two additional controls apply:

• Rate cap (params.dt.evection_*): bounds the next PROTEUS coupling step so the fractional change in eccentricity_sat stays near evection_target_rel_de, since Obliqua's own adaptive spectrum window is chosen once per step from the eccentricity at that step's start.
• Growth limiter: caps how fast the step size can grow relative to the previous step while inside the band, or for evection_cooldown_iters iterations after leaving it, so the coupling step does not snap back to its ordinary size the instant the band is exited.

Both are folded into the single exported column evection_dt_cap_yr, which interior_energetics.timestep.next_step applies as one of several caps on the next main-loop step.

Reproducing capture and peak eccentricity against Rufu & Canup (2020) Figure 3 is validated in tests/integration/test_slow_orbit_evection_ctl.py (@pytest.mark.slow): the real ps1d_evec model, driven by a Mignard constant-time-lag tidal spectrum, reproduces the paper's resonance-capture timing window and peak-eccentricity location to better than 0.1% in semi-major axis.

The "three clocks"¶

Fine-grained diagnostics for ps1d_evec (fine_evection_data.csv) distinguish three different notions of "step":

1. PROTEUS main-loop clock: one evolve_orbit_satellite call per coupling iteration, spanning interior_o.dt years.
2. Solver clock: the adaptive-substep controller's own internal dt_yr steps within one main-loop call (see below), and, within a single substep, solve_ivp's own internal integration points.
3. Storage clock: how densely those solver-clock samples are written to disk. Every solver-clock sample is kept while inside the evection band; outside the band, samples are thinned to a target spacing (orbit.solver.fine_csv_target_rel_dt, a fraction of the current main-loop step) so the file does not grow unbounded over a long, quiescent run.

Hansen coefficients¶

The sp1d/ps1d/ps1d_evec vectorial tidal models expand the tide-raising potential in Hansen coefficients X_k^{n,m}(e), evaluated at every ODE substep. A direct FFT evaluation is too slow for that (implicit solvers probe many micro-varying eccentricities per step), and nearest-neighbor caching would introduce discontinuities that break implicit solvers. orbit.hansen instead pre-tabulates, once per run, the eccentricity-dependent mode window [k_min, k_max] (since the number of significant modes grows from ~10 near e=0 to several hundred above e=0.8) and then the coefficient values themselves on that window, both linearly interpolated in e thereafter. The tables are warmed up once by orbit.wrapper.run_orbit at Time<=1; a hot-path call lazily builds them if warm-up was skipped.

Adaptive substep controller¶

sp1d, ps0d, ps1d, and ps1d_evec all integrate through the same accept/reject controller, orbit.common.run_adaptive_orbit_substeps. For each attempted internal step it stages the tentative result, checks it for unphysical values (negative semi-major axis, eccentricity outside [0, 1), non-finite spin) and for excessive relative change in tracked quantities (orbit.solver.max_rel_*), then either merges it in and grows the step, or discards it and shrinks the step (orbit.solver.growth/shrink).

The same call also keeps the planet's moment of inertia (C_int, from interior_energetics.common.get_C_planet) consistent with the live interior structure: rather than jumping to the freshly computed value once per call (which would put a discontinuity in any quantity that depends on the planet's spin rate, such as ps1d_evec's oblateness-driven precession), the controller ramps C_int linearly across the call's accepted substeps, rescaling axial_period at each one to conserve C_int * Omega_p (angular momentum).

ps0d additionally gets a cumulative drift cap: ps0d has no eccentricity or spin feedback of its own, so many small, individually-legal substeps can compound into a large silent migration within a single call.

Termination criteria¶

Orbital and rotational state feed three physical stopping conditions (params.stop.*, checked in utils.terminate):

• Disintegration (params.stop.disint/disint_sat): the planet or satellite orbiting within its partner's Roche limit, or spinning faster than its breakup rate.
• Satellite escape (params.stop.satellite): the satellite's semi-major axis exceeding sma_max.

Model selection guide¶

• Want a cheap, non-physical heat source for testing the interior's response to tides? orbit.module = 'dummy', no evolution model.
• Want a self-consistent solid-body tidal response with minimal setup, spin-orbit synchronised, low eccentricity? orbit.module = 'lovepy'.
• Want the planet's orbit and spin to evolve self-consistently with its own interior structure, at arbitrary eccentricity? orbit.module = 'obliqua', orbit.perturber = 'star', orbit.star_planet_model = 'sp1d'.
• Want a fast, closed-form estimate of star-planet tidal circularisation without resolving spin? orbit.module.dummy supplying Imk2 plus orbit.star_planet_model = 'sp0d'.
• Want a satellite's orbit (e.g. a moon) to evolve, including its own tidal response? orbit.module = 'obliqua', orbit.perturber = 'satellite', orbit.planet_satellite_model = 'ps1d' (or 'ps0d' for a cheaper, eccentricity-frozen estimate).
• Want to also study capture into, and eccentricity pumping by, the evection resonance? orbit.module = 'obliqua', orbit.perturber = 'satellite', orbit.planet_satellite_model = 'ps1d_evec'.

Testing¶

• Unit tests (tests/orbit/*.py, @pytest.mark.unit) mock Julia calls for lovepy/obliqua and pin closed-form results for sp0d, sp1d, ps0d, ps1d against their source papers, each with a discrimination guard against the nearest plausible wrong formula (wrong exponent, wrong sign, or the wrong mass in a prefactor). Angular momentum conservation for sp1d/ps1d is checked directly, both as a raw total and as a targeted spin-vs-orbital exchange test sized to catch a bug the raw total alone would miss.
• Slow/integration tests (tests/integration/test_slow_orbit_evection_ctl.py, @pytest.mark.slow) drive the real ps1d_evec model end-to-end and compare against the published Rufu & Canup (2020) evection trajectory.
• Per-source-file test inventories, references, and re-derivations are tracked under docs/Validation/orbit/: orbit.py, satellite.py, wrapper.py, obliqua.py.

See also: Model description | Star and orbit configuration | Execution and output configuration | Validation: orbit