diff --git a/.github/actions/setup-proteus/action.yml b/.github/actions/setup-proteus/action.yml
index 44766bc80..c797877b9 100644
--- a/.github/actions/setup-proteus/action.yml
+++ b/.github/actions/setup-proteus/action.yml
@@ -16,7 +16,7 @@ inputs:
julia-version:
description: Julia version
required: false
- default: '1.12.6'
+ default: '1.13.0'
install-editable-submodules:
description: >
Whether to clone MORS / aragog / JANUS / CALLIOPE / ZEPHYRUS / Zalmoxis
diff --git a/.github/workflows/ci-nightly.yml b/.github/workflows/ci-nightly.yml
index 27fa94458..d9abd8a97 100644
--- a/.github/workflows/ci-nightly.yml
+++ b/.github/workflows/ci-nightly.yml
@@ -262,6 +262,13 @@ jobs:
os: ubuntu-latest
files: >-
tests/integration/test_slow_agni_aragog.py
+ # ---- shard: orbit-evection (pure orbit.satellite physics,
+ # no real submodule needed; Linux-only, ~25-40 min) ----
+ - shard: orbit-evection
+ tier: extended
+ os: ubuntu-latest
+ files: >-
+ tests/integration/test_slow_orbit_evection_ctl.py
# ---- shard: grid (dummy-backend parameter grid; a single
# module-scoped fixture runs the whole grid, ~27 min. Linux-only
# extended because it needs no macOS parity) ----
diff --git a/.gitignore b/.gitignore
index 2aedede6b..d4907f7cf 100644
--- a/.gitignore
+++ b/.gitignore
@@ -86,6 +86,13 @@ lovepy
Lovepy
LOVEPY
LovePy
+obliqua
+Obliqua
+OBLIQUA
+Obliqua_data
+platon
+Platon
+PLATON
/prt
petitradtrans
petitRADTRANS
diff --git a/docs/Explanations/code_architecture.md b/docs/Explanations/code_architecture.md
index 1dded235f..a1b468d93 100644
--- a/docs/Explanations/code_architecture.md
+++ b/docs/Explanations/code_architecture.md
@@ -12,7 +12,7 @@ coupled planetary evolution simulation:
- [`atmos_chem/`](https://github.com/FormingWorlds/PROTEUS/tree/main/src/proteus/atmos_chem): atmospheric photochemistry (VULCAN, dummy)
- [`escape/`](https://github.com/FormingWorlds/PROTEUS/tree/main/src/proteus/escape): atmospheric mass loss (ZEPHYRUS, dummy)
- [`outgas/`](https://github.com/FormingWorlds/PROTEUS/tree/main/src/proteus/outgas): volatile partitioning (CALLIOPE, atmodeller, dummy)
-- [`orbit/`](https://github.com/FormingWorlds/PROTEUS/tree/main/src/proteus/orbit): orbital evolution and tides (Obliqua/LovePy, dummy)
+- [`orbit/`](https://github.com/FormingWorlds/PROTEUS/tree/main/src/proteus/orbit): tidal response (Obliqua/LovePy, dummy) and orbital evolution (native)
- [`star/`](https://github.com/FormingWorlds/PROTEUS/tree/main/src/proteus/star): stellar evolution and spectra (MORS, dummy)
Most modules follow a common pattern: a `wrapper.py` defining the dispatch
diff --git a/docs/Explanations/dummy_modules.md b/docs/Explanations/dummy_modules.md
index 5cd90dbb9..3c2ef03ca 100644
--- a/docs/Explanations/dummy_modules.md
+++ b/docs/Explanations/dummy_modules.md
@@ -139,7 +139,7 @@ from Kepler's third law. A configurable tidal heating amplitude
(`H_tide`) is applied to mantle layers where the melt fraction exceeds
a threshold (`Phi_tide`), providing a simple parameterised heat source
for testing the interior module's response to tidal power without
-running the full LovePy viscoelastic solver.
+running the full LovePy or Obliqua viscoelastic solvers.
---
diff --git a/docs/Explanations/model.md b/docs/Explanations/model.md
index 7b711dcdc..9652bdbe1 100644
--- a/docs/Explanations/model.md
+++ b/docs/Explanations/model.md
@@ -45,7 +45,8 @@ atmosphere), enabling hierarchical model intercomparison.
| Star | [MORS](https://proteus-framework.org/MORS/), dummy | Stellar evolution and spectrum |
| Escape | [ZEPHYRUS](https://github.com/FormingWorlds/ZEPHYRUS), dummy | Atmospheric escape |
| Outgassing | [CALLIOPE](https://proteus-framework.org/CALLIOPE/), [atmodeller](https://github.com/djbower/atmodeller), dummy | Volatile exchange between interior and atmosphere |
-| Orbit | [Obliqua](https://github.com/FormingWorlds/Obliqua), dummy | Orbital evolution and tidal heating |
+| Tides | [Obliqua](https://proteus-framework.org/Obliqua), [Lovepy](https://github.com/nichollsh/LovePy), dummy | Tidal response of the planet (Love numbers, heating) |
+| Orbit | PROTEUS (internal) | Orbital evolution (semi-major axis, eccentricity, spin) |
| Observations | [petitRADTRANS](https://petitradtrans.readthedocs.io/), none | Synthetic transit and eclipse spectra |
Each module is maintained in its own repository and can be used as a standalone package outside of PROTEUS. The following sections describe each module's physical role and how PROTEUS couples to it.
@@ -189,11 +190,17 @@ Some notable consequences of step 3:
`utils.coupler.assert_mass_conservation` therefore checks two things separately. `M_atm <= M_planet` is enforced under `outgas.vapourise = false`. `M_vol_atm` equals the sum of the per-species atmospheric masses; rock vapour is excluded from `M_vol_atm` by definition.
-## Orbital evolution: Obliqua
+## Tidal evolution: Obliqua, Lovepy
-**[Obliqua](https://github.com/FormingWorlds/Obliqua)** (Julia) evolves the orbital semi-major axis and eccentricity under the influence of tidal dissipation. The tidal response of the planet is computed from its interior structure and rheology using a viscoelastic love-number solver (LovePy). Tidal heating power is distributed radially across the mantle and fed back into the interior energy equation. Obliqua also computes the spin-orbit evolution and checks for dynamical stability (Roche limit, Hill sphere).
+**[Obliqua](https://proteus-framework.org/Obliqua)** (Julia) computes the multi-phase tidal response of the planet. The tidal love-numbers of the planet is computed from its interior structure and rheology using various viscoelastic rheological models. Tidal heating power is distributed radially across the mantle and fed back into the interior energy equation. The model is valid for arbitrary eccentricity and spin-orbit misalignment, as it models arbitrary tidal degrees and modes.
-Config section: `[orbit]`. Reference: [Star and orbit configuration](../Reference/config/star_orbit.md).
+**[Lovepy](https://github.com/nichollsh/LovePy)** (Julia) computes the solid-only tidal response of the planet using a Maxwell rheology. Tidal heating power is distributed radially across the mantle and fed back into the interior energy equation. It assumes spin-orbit synchronisation and a small eccentricity, as only the dominant degree-2 tidal modes are considered.
+
+## Orbital evolution: PROTEUS (internal)
+
+**[Orbital evolution](https://github.com/FormingWorlds/PROTEUS/tree/main/src/proteus/orbit)** (Python) computes the time evolution of the independent orbital parameters (semi-major axis, eccentricity, spin, apsidal precession) under the influence of tidal dissipation in both the primary and the perturbing body. The model combines spin-orbit dynamics with eccentricity evolution, based on a vectorial approach expressed in Hansen coefficients. The model allows for angular momentum draining through the evection resonance, but conserves it in all other cases.
+
+Config section: `[orbit]`. Reference: [Star and orbit configuration](../Reference/config/star_orbit.md). For the full set of star-planet and planet-satellite models (sp0d/sp1d/ps0d/ps1d/ps1d_evec), the evection resonance, and the satellite Love-number lookup workflow, see [Orbital dynamics and tides](orbit.md).
## Synthetic observations: petitRADTRANS
@@ -240,7 +247,7 @@ architecture and for quick parameter exploration.
| Star | Fixed effective temperature and luminosity; Planck-function spectrum at a user-specified $T_\mathrm{eff}$ |
| Escape | Constant bulk mass loss rate (user-specified kg/s), distributed proportionally across elements |
| Outgassing | Melt-fraction-dependent volatile partitioning with fixed stoichiometry, no equilibrium chemistry |
-| Orbit | Fixed semi-major axis and eccentricity; configurable parameterised tidal heating |
+| Orbit & Tides | Fixed semi-major axis and eccentricity; configurable parameterised tidal heating |
The [Quick start tutorial](../Tutorials/quick_start_dummy.md) runs PROTEUS
with all modules set to dummy.
diff --git a/docs/Explanations/orbit.md b/docs/Explanations/orbit.md
new file mode 100644
index 000000000..fa576e1f4
--- /dev/null
+++ b/docs/Explanations/orbit.md
@@ -0,0 +1,403 @@
+
+
+
+
+
+
+
+# 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](../Reference/config/star_orbit.md). For the
+one-paragraph summaries of each tidal module, see
+[Tidal evolution: Obliqua, Lovepy](model.md#tidal-evolution-obliqua-lovepy)
+and [Orbital evolution: PROTEUS (internal)](model.md#orbital-evolution-proteus-internal).
+
+## Two independent evolution families
+
+PROTEUS evolves an orbit around exactly one body at a time:
+
+
+
+- 
+
+ **Star-planet**
+
+ Evolve the planet's own orbit around its host star.
+
+ [Star-planet models](#star-planet-models-orbitstar_planet_model){ .md-button .md-button--primary }
+
+- 
+
+ **Planet-satellite**
+
+ Evolve a satellite's orbit around the planet.
+
+ [Planet-satellite models](#planet-satellite-models-orbitplanet_satellite_model){ .md-button .md-button--primary }
+
+
+
+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`).
+
+| Module | What it computes | Assumptions |
+|---|---|---|
+| `dummy` | A fixed heating rate `H_tide` applied where the local melt fraction satisfies an inequality (`Phi_tide`, e.g. `"<0.3"`); returns a fixed `Imk2`. | No physical rheology; a parameterised heat source for testing the interior's response. |
+| `lovepy` | Solid-only viscoelastic (Maxwell) Love number from the topmost region above a viscosity threshold. | Degree-2 only, small eccentricity, spin-orbit synchronisation. |
+| `obliqua` | Multi-phase (solid/mushy/fluid) Love-number spectrum from the full interior profile, for arbitrary tidal degree/mode (`orbit.obliqua.n`/`m`) and eccentricity. | The only module that can compute a satellite-side response (see [Satellite Love-number lookup](#satellite-love-number-lookup-obliqua-only) below); requires `orbit.perturber` set explicitly. |
+
+??? note "dummy in a nutshell - van Dijk et al. (2026)[^cite-vandijk2026]"
+ Grew out of a study of the Hadean Earth-Moon system, asking how long
+ tidal heating could keep a magma ocean from fully solidifying without
+ modelling the rheology in detail: heating is simply switched on below
+ a melt-fraction threshold and scaled linearly with the remaining
+ solid fraction. Sweeping that heating rate reveals quasi-steady
+ "global radiative equilibrium" epochs, where interior heating and
+ atmospheric cooling balance.
+
+??? note "lovepy in a nutshell - Nicholls et al. (2025)[^cite-nicholls2025lovepy]"
+ Solves for the planet's actual viscoelastic (Maxwell) response by
+ propagating the tidal deformation through radial layers, rather than
+ prescribing a heating rate. Applied to the L 98-59 system, it revealed
+ a self-limiting "radiation-tide-rheology" feedback: as tidal heating
+ softens the mantle, dissipation efficiency drops too, capping heating
+ at levels up to two orders of magnitude below earlier estimates -
+ while still being enough to sustain magma oceans for billions of
+ years.
+
+??? note "obliqua in a nutshell"
+ Obliqua generalises the same viscoelastic idea beyond `lovepy`'s
+ single solid layer and low-eccentricity limit: it resolves solid,
+ mushy, and fluid regions together, at arbitrary tidal degree, mode,
+ and eccentricity. The dummy-module study above hinted at how much
+ tidal heating can matter for early evolution, but only for one
+ fixed, simplified regime; Obliqua exists to track the tidal response
+ self-consistently across the much wider range of thermal and
+ orbital states real exoplanets occupy.
+
+!!! warning "Dynamic-tide resonances (Obliqua)"
+ Setting `orbit.obliqua.solid.inertial_terms` to `true` solves the full
+ finite-frequency problem instead of the quasi-static ($\omega \to 0$) approximation.
+ When tidal forcing matches a normal-mode frequency of the body, it produces a
+ physically real, bounded peak in the Love number—not a numerical grid artifact.
+
+ * **When they occur:** Only when parts of the mantle are molten or mushy.
+ * **Control knob:** Use `params.dt.mushy_maximum` to tighten timesteps during solidification.
+ This forces PROTEUS to re-sample Obliqua frequently enough to resolve resonance crossings.
+ * **Safety cap:** `orbit.obliqua.cap_LN` clamps each mode's Love number to a fixed multiple
+ of the classical fluid limit for its degree n. This prevents extreme heating spikes while
+ macro-steps are too large to fully resolve the resonance timescale.
+
+## Star-planet models (`orbit.star_planet_model`)
+
+| Model | Evolves | Reference | Notes |
+|---|---|---|---|
+| `sp0d` | `semimajorax`, `eccentricity` | Driscoll & Barnes (2015)[^cite-driscoll2015], Eq. 15-16 | Closed-form two-ODE system in `(a, e)` only; no spin dynamics, so it is **not** angular-momentum-conserving by construction. |
+| `sp1d` | `axial_period`, `semimajorax`, `eccentricity`, `plan_star_am` | Correia & Valente (2022)[^cite-correia2022] | Vectorial, Hansen-coefficient formulation restricted to planetary tides (star assumed non-dissipative). Genuinely angular-momentum-conserving; verified by dedicated tests. |
+
+??? note "sp0d in a nutshell - Driscoll & Barnes (2015)"
+ Written for rocky planets around M dwarfs, where the habitable zone
+ sits close enough in that tides matter. Treats the planet as a
+ passive, non-rotating "equilibrium tide" bulge dragged slightly
+ behind (or ahead of) the star: that lag drains eccentricity and
+ trades orbital energy for heat inside the planet. No spin, no
+ resonances -- just a slow circularisation clock coupled to whatever
+ the interior does with the heat.
+
+??? note "sp1d in a nutshell - Correia & Valente (2022)"
+ Instead of one lumped tidal bulge, the tidal potential is decomposed
+ into its individual Fourier harmonics (Hansen coefficients), each
+ oscillating at its own forcing frequency and dissipating
+ independently. This removes the low-eccentricity assumption baked
+ into classical tidal theory, and it means spin and orbit are evolved
+ together as one system, exchanging angular momentum internally.
+
+Both integrate with `scipy.solve_ivp` (`orbit.solver.*` controls method
+and tolerances).
+
+## Planet-satellite models (`orbit.planet_satellite_model`)
+
+| Model | Evolves | Reference | Notes |
+|---|---|---|---|
+| `ps0d` | `semimajorax_sat`, `axial_period` | Korenaga (2023)[^cite-korenaga2023], Eq. 58-60 | No eccentricity evolution, no satellite-side tide. Uses the `M_sat << M_planet` limit of the orbital angular-momentum term (~1.2% error for Earth-Moon). |
+| `ps1d` | `axial_period`, `axial_period_sat`, `semimajorax_sat`, `eccentricity_sat`, `plan_sat_am` | Correia & Valente (2022)[^cite-correia2022] | Same vectorial approach as `sp1d`, extended to track both planet-raised and satellite-raised tidal contributions separately. Requires satellite-side Love-numbers (see below). |
+| `ps1d_evec` | Everything `ps1d` evolves, plus `evection_angle` | `ps1d` physics plus Rufu & Canup (2020)[^cite-rufu2020] evection-resonance terms | Adds a J2-driven apsidal-precession term and a resonant forcing term. See [Evection resonance](#evection-resonance-ps1d_evec) below. |
+
+??? note "ps0d in a nutshell - Korenaga (2023)"
+ Built to explain why the Moon's magma ocean stayed molten for so
+ long: rather than solving the tidal potential in detail, it tracks
+ one number, the system's total (spin + orbital) angular momentum,
+ and lets the planet's tidal dissipation rate spend it. As the
+ planet's spin winds down, the satellite's orbit must expand to keep
+ the ledger balanced - a bookkeeping model, not a torque model, so
+ it is cheap and exactly momentum-conserving, at the cost of no
+ eccentricity evolution.
+
+??? note "ps1d in a nutshell - Correia & Valente (2022)"
+ The same Hansen-coefficient decomposition as `sp1d`, but with two
+ dissipating bodies instead of one: both the planet's and the
+ satellite's tidal responses pull on the shared orbit, so each of
+ their spins, the semi-major axis, and the eccentricity all evolve
+ together, coupled through one exchange of angular momentum.
+
+??? note "ps1d_evec in a nutshell - Rufu & Canup (2020)"
+ As a tidally-receding moon's orbit expands, its slow apsidal
+ precession can fall into step with the star's apparent yearly
+ motion - a secular resonance. Falling into that resonance is like
+ pushing a swing at just the right moment: it pumps up the moon's
+ eccentricity long after ordinary tides alone would have damped it
+ flat, which is the paper's proposed route to the Moon's present-day
+ orbital tilt.
+
+!!! warning "Satellite Love-number lookup"
+ Both `ps1d` and `ps1d_evec` need the satellite's own Love-number spectrum as a
+ function of forcing frequency, which only `orbit.module='obliqua'` can
+ supply (via [`LN_from_lookup`](#satellite-love-number-lookup-obliqua-only)).
+ Using `ps1d`/`ps1d_evec` unconditionally populates the satellite's tidal
+ parameters in `tides_o` through Obliqua's `lookup_from_interior` at the
+ start of the run.
+
+## 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:**
+
+| `orbit.module` | `Imk2` | `tides_o` | `hf_row['F_tidal']` |
+|---|---|---|---|
+| `dummy` | Yes | No | Yes |
+| `lovepy` | Yes | Yes | yes |
+| `obliqua` | Yes, only when `orbit.obliqua.n == [2]` (`0.0` otherwise) | Yes, planet always, satellite too when `orbit.perturber='satellite'` | yes |
+
+**What each orbit model reads:**
+
+| Model | Reads | Compatible `orbit.module` |
+|---|---|---|
+| `sp0d` | `hf_row['Imk2']` | `dummy`, `lovepy`, `obliqua` (requires `orbit.obliqua.n == [2]`) |
+| `sp1d` | `tides_o`, (`primary='planet', perturber='star'`) | `lovepy`, `obliqua` |
+| `ps0d` | `hf_row['F_tidal']` | `dummy`, `lovepy`, `obliqua` |
+| `ps1d` | `tides_o`, (both `primary='planet', perturber='satellite'` and `primary='satellite', perturber='planet'`) | `lovepy`, `obliqua` |
+| `ps1d_evec` | Same as `ps1d`, plus `evection_angle` | `lovepy`, `obliqua` (Note that `lovepy` breaks down at high eccentricities, so it is not recommended for this case) |
+
+!!! warning "Note on `*1d` models"
+ `sp1d`, `ps1d`, and `ps1d_evec` are rejected at config load when
+ `orbit.module` is not `'obliqua'` or `'lovepy'`. Prefer
+ `orbit.module='obliqua'` for any `*1d` orbit model.
+
+---
+
+### 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.
+
+
+
+
+
+
+Example evection-resonance episode. The satellite starts outside
+the resonance band, evolving freely; capture into the band locks the
+evection angle to the resonant condition and pumps up the eccentricity;
+escape from the band later returns the system to free, non-resonant
+precession.
+