Skip to content

PROTEUS Model Framework

Overview

PROTEUS is a modular framework for simulating the time evolution of small (exo)planets. It is designed to be flexible, reflecting the broad diversity of planetary conditions already discovered, with the view of being updated to incorporate additional physics as the need arises. This approach stands in contrast to common monolithic models in the literature. PROTEUS is free and open-source, which permits external scrutiny of its workings. It is directly based upon the model of Lichtenberg et al. (2021) 1 although the code has evolved substantially from that state.

Design philosophy

George Box famously put that "all models are wrong, but some are useful". PROTEUS therefore leverages a modular and hierarchical modelling approach. Multiple independent models can fill the role of a given module, and each model can be used stand-alone. Hierarchical modelling allows inter-comparison of simple and complex models, taking advantage of the easy comprehension of the simple in order to diagnose and validate the qualitative behaviour of the complex. Simplified dummy modules are not designed for quantitatively meaningful calculations, but only to qualitatively capture end-member behaviours.

System architecture

Although PROTEUS aims to treat the problem of planetary evolution, it must necessarily also handle external processes which act upon the planet (e.g. tidal heating). The framework therefore models the combined system of a planet, its orbital mechanics, and the evolution of its host star. The planet itself is conceptually sub-divided into a vaporised atmosphere component above an interior component containing a silicate mantle and metallic core. PROTEUS facilitates communication between individual software modules which each implement a model for a specific part of the overall system. Conceptually, PROTEUS modules (e.g. the interior) are 'slots' which are filled by specific implementations: the 'models' (e.g. Aragog).

PROTEUS module schematic (light mode) PROTEUS module schematic (dark mode)
Schematic of the PROTEUS modular architecture. Radial domains (labelled boxes) are each handled by interchangeable modules; the star and orbit provide external forcing; arrows show the quantities exchanged at each coupling timestep.

The planet is divided into radial domains: a metallic core, a silicate mantle (solid and partially molten), a surface boundary layer, and a volatile atmosphere. The star and orbit force the system through instellation, XUV flux, and tidal heating. At each coupling timestep, the interior module passes the surface temperature and heat flux upward to the atmosphere, which returns the outgoing longwave radiation; the outgassing module partitions volatiles between the melt and the atmosphere based on the current temperature and melt fraction; and the escape module removes mass from the top of the atmosphere driven by the stellar XUV flux. Multiple implementations are available for each slot (e.g. Aragog or SPIDER for the interior, AGNI or JANUS for the atmosphere), enabling hierarchical model intercomparison.

Module overview

Module Implementations Role
Structure Zalmoxis, dummy Interior structure (radius, density, gravity profiles)
Interior SPIDER, Aragog, boundary, dummy Mantle/core thermal evolution
Atmosphere (climate) AGNI, JANUS, dummy Radiative-convective profile
Atmosphere (chemistry) VULCAN, dummy Chemical kinetics
Star MORS, dummy Stellar evolution and spectrum
Escape ZEPHYRUS, dummy Atmospheric escape
Outgassing CALLIOPE, atmodeller, dummy Volatile exchange between interior and atmosphere
Orbit Obliqua, dummy Orbital evolution and tidal heating
Observations petitRADTRANS, 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.


Interior structure: Zalmoxis

Zalmoxis computes the hydrostatic equilibrium structure of a differentiated planet (metallic core + silicate mantle + volatile envelope). Given a total planet mass, bulk composition, and surface temperature, Zalmoxis integrates the equations of hydrostatic equilibrium inward from the surface using a tabulated equation of state (EOS), returning radial profiles of pressure, density, temperature, and gravitational acceleration. It also computes the core radius, mantle mass, and surface gravity.

Zalmoxis supports several EOS backends, including the PALEOS MgSiO\(_3\) tables and the Wolf & Bower (2018) 2 parameterisation. The PALEOS mantle EOS comes in two forms: PALEOS:MgSiO3 solves the structure on the unified table and derives the solidus from the liquidus with a fixed mushy-zone factor (the Stixrude 2014 solidus-to-liquidus ratio), while PALEOS-2phase:MgSiO3 uses separate solid and liquid tables that set the phase boundaries directly. Both are detailed in the interior configuration reference. The structure solution is used by PROTEUS to initialise the planet's radius and to dynamically update the structure during the simulation when the interior thermal state changes (see Structure-interior coupling below).

Config section: [interior_struct]. Reference: Interior configuration.

Interior energetics: Aragog, SPIDER, boundary

The interior energetics module evolves the mantle temperature and melt fraction in time by solving the energy equation in the planetary interior.

Aragog (Python/JAX) solves the interior energy equation in the temperature-pressure (T-P) formulation. It discretises the mantle on a radial grid and advances the temperature profile using the SUNDIALS CVODE integrator with an analytic Jacobian computed via JAX automatic differentiation. Aragog handles the full mushy-zone (partial melt) regime, including phase-dependent material properties, radiogenic heating, tidal heating, and core cooling. It is the recommended interior module for most configurations.

SPIDER (C) uses the temperature-entropy (T-S) formulation 3. It discretises the mantle on a staggered radial grid and solves the entropy equation using PETSc's implicit time integrator. SPIDER provides an independent cross-check on Aragog for the same physical problem formulated in a different thermodynamic variable. It requires PETSc and is an optional installation component.

Boundary solves a simplified single-node energy balance ODE for the mantle surface temperature, treating the mantle as a single thermal reservoir with Arrhenius or aggregate viscosity. It is useful for rapid exploration of parameter space and for configurations where the full radial resolution of Aragog or SPIDER is not needed.

Config section: [interior_energetics]. Reference: Interior configuration.

Atmosphere climate: AGNI, JANUS

The atmosphere climate module solves the radiative-convective equilibrium structure of the atmosphere given the surface temperature, surface pressure, atmospheric composition, and incoming stellar spectrum.

AGNI (Julia) is a radiative-convective atmosphere model that uses the SOCRATES spectral radiative transfer library (Fortran) as its radiation core. AGNI solves for the self-consistent atmospheric temperature profile by iterating a Newton solver on the radiative-convective flux balance. It supports grey-gas and correlated-k radiative transfer, Rayleigh scattering, clouds, and condensation of volatile species. AGNI returns the outgoing longwave radiation (OLR), Bond albedo, and atmospheric temperature-pressure profile.

JANUS (Python) is a 1D convective atmosphere module that computes a radiative-convective profile using a band-averaged two-stream approach. JANUS is faster than AGNI but makes stronger simplifying assumptions (e.g. no iterative Newton solve; a prescribed adiabatic lower atmosphere). It is useful for parameter sweeps and initial exploration.

Config section: [atmos_clim]. Reference: Atmosphere configuration.

Atmospheric chemistry: VULCAN

VULCAN (Python) is a photochemical kinetics code that computes steady-state atmospheric mixing ratios given a temperature-pressure profile, eddy diffusion profile, and UV stellar spectrum. In PROTEUS, VULCAN can run in two modes: online (at each snapshot during the simulation) or offline (once after the simulation completes). Most configurations use the offline mode because the chemistry calculation is computationally expensive relative to the coupling timestep.

Config section: [atmos_chem]. Reference: Atmosphere configuration.

Stellar evolution: MORS

MORS (Python) provides stellar evolutionary tracks and spectral energy distributions. It supports two track families: Spada 7 (rotation-dependent tracks with activity-calibrated XUV luminosities, suitable for solar-type stars) and Baraffe 8 (mass-luminosity tracks for low-mass stars). MORS interpolates the stellar mass, radius, effective temperature, bolometric luminosity, and XUV luminosity at any stellar age, and synthesises a wavelength-resolved spectrum by scaling a modern reference spectrum (observed MUSCLES data, solar NREL data, or a synthetic PHOENIX spectrum) to the historical luminosity.

Config section: [star]. Reference: Star and orbit configuration.

Atmospheric escape: ZEPHYRUS

The escape module computes atmospheric mass loss rates driven by the stellar XUV flux.

ZEPHYRUS (Python) implements energy-limited escape, computing a bulk mass loss rate from the XUV flux, planet mass, and XUV absorption radius. The bulk rate is then distributed across elements proportionally to their atmospheric abundance (unfractionated escape). ZEPHYRUS also provides a tidal correction factor for planets in close-in orbits.

Config section: [escape]. Reference: Escape and outgassing configuration.

Volatile outgassing: CALLIOPE, atmodeller

The outgassing module computes the thermodynamic equilibrium partitioning of volatiles between the atmosphere, silicate melt, and solid mantle.

CALLIOPE (Python) solves for the equilibrium partial pressures and dissolved volatile concentrations given the mantle temperature, melt fraction, and total element inventories for H, C, N, S, and O 4. It uses parameterised solubility laws and an fO2 buffer (configurable as an IW offset) to compute the redox state 6. CALLIOPE handles the full set of major volcanic gases (H\(_2\)O, CO\(_2\), H\(_2\), CO, N\(_2\), SO\(_2\), S\(_2\), CH\(_4\)).

atmodeller (Python/JAX) is an alternative outgassing backend that uses a real-gas equation of state and a more detailed thermochemical treatment 5. atmodeller provides an independent cross-check on CALLIOPE for the same volatile partitioning problem.

Config section: [outgas]. Reference: Escape and outgassing configuration.

Orbital evolution: Obliqua

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).

Config section: [orbit]. Reference: Star and orbit configuration.

Synthetic observations: petitRADTRANS

petitRADTRANS (Python) is a radiative transfer code for computing exoplanet transmission and emission spectra. PROTEUS uses petitRADTRANS as a forward model to synthesise what an observer would measure given the simulated atmospheric state.

For each active composition source PROTEUS runs two calculations:

  • Transit depth: wavelength-dependent \((R_\mathrm{transit}/R_\star)^2\) in ppm, computed from the atmospheric scale height via Radtrans.calculate_transit_radii.
  • Eclipse depth: wavelength-dependent thermal emission contrast \((F_\mathrm{planet}/F_\star)\) in ppm, computed via Radtrans.calculate_flux.

After computing the full-atmosphere baseline spectrum, PROTEUS can optionally repeat the calculation once per included line species with that species removed. This leave-one-out procedure quantifies the spectral contribution of each gas and produces additional columns in the output CSV files.

The active source set is configured by observe.source (all, outgas, profile, offchem). Leave-one-out spectra are controlled by observe.remove_one_gas. Whether transit, eclipse, or both products are generated is controlled by observe.spectrum_type.

The atmosphere profile is interpolated onto a 100-point log-spaced pressure grid before being passed to petitRADTRANS, and temperature is clipped to the valid table range [100.5 K, 3999.5 K].

Before constructing the Radtrans object, PROTEUS moves the species with the broadest opacity wavelength coverage to the front of the line-species list.

Config section: [observe]. Reference: Observations configuration.

Dummy modules

Every module slot has a dummy implementation for testing, debugging, and hierarchical modelling. Dummy modules replace the full physics with minimal parameterisations that capture qualitative behaviour at negligible computational cost. They require no external solvers, no compiled code, and no reference data, making them suitable for verifying the coupling architecture and for quick parameter exploration.

Module Dummy behaviour
Structure Noack & Lasbleis (2020) 9 analytical scaling laws for interior radius, density, and gravity as a function of planet mass and iron content
Interior energetics Heat-capacity integrator with prescribed solidus/liquidus; cooling driven by atmospheric flux
Atmosphere climate Grey-body model: \(F_\mathrm{OLR} = \sigma [T_\mathrm{surf}(1-\gamma)]^4\). Optionally a fixed-flux mode
Atmosphere chemistry Parameterised vertical profiles with cold trap and approximate photolysis products
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

The Quick start tutorial runs PROTEUS with all modules set to dummy.


Structure-interior coupling

The structure and interior modules serve complementary roles. The interior module (SPIDER or Aragog) evolves the mantle temperature and melt fraction in time, while the structure module (Zalmoxis) solves for the hydrostatic equilibrium of the planet given a temperature profile and equation of state.

When interior_struct.module = 'zalmoxis' and interior_struct.zalmoxis.update_interval > 0, PROTEUS can dynamically recompute the planetary structure during a simulation. Structure updates are governed by a hybrid trigger that combines physics-based criteria with timing constraints:

  • Floor (update_min_interval): minimum time between updates, preventing excessive recomputation during rapid cooling.
  • Ceiling (update_interval): maximum time between updates, guaranteeing periodic recomputation.
  • Temperature trigger (update_dtmagma_frac): fires when the relative change in mantle surface temperature since the last update exceeds the threshold.
  • Melt fraction trigger (update_dphi_abs): fires when the absolute change in global melt fraction since the last update exceeds the threshold.

Setting update_interval = 0 disables dynamic updates entirely; the structure is computed only at initialisation.

Time evolution vs equilibrium

Only the interior and star modules have an explicit notion of time-evolution. All other modules are applied at equilibrium, such that the quantities calculated by these modules are effectively updated instantaneously at each time-step. This assumes that the physical processes handled by these equilibrium modules reach steady-state on time-scales shorter than the physics considered by interior and stellar evolution modules.

Further reading


  1. Lichtenberg, T., Bower, D.J., Hammond, M., et al., Vertically resolved magma ocean-protoatmosphere evolution: H2, H2O, CO2, CH4, CO, O2, and N2 as primary absorbers, Journal of Geophysical Research: Planets, 126, e2020JE006711, 2021. SciX

  2. Wolf, A.S. & Bower, D.J., An equation of state for high pressure-temperature liquids (RTpress) with application to MgSiO3 melt, Physics of the Earth and Planetary Interiors, 278, 59-74, 2018. SciX

  3. Bower, D.J., Sanan, P. & Wolf, A.S., Numerical solution of a non-linear conservation law applicable to the interior dynamics of partially molten planets, Physics of the Earth and Planetary Interiors, 274, 49-62, 2018. SciX

  4. Bower, D.J., Kitzmann, D., Wolf, A.S., et al., Linking the evolution of terrestrial interiors and an early outgassed atmosphere to astrophysical observations, Astronomy & Astrophysics, 631, A103, 2019. SciX

  5. Bower, D.J., Thompson, M.A., Hakim, K., et al., Diversity of low-mass planet atmospheres in the C-H-O-N-S-Cl system, The Astrophysical Journal, 995, 59, 2025. SciX

  6. Nicholls, H., Lichtenberg, T., Bower, D.J. & Pierrehumbert, R., Magma ocean evolution at arbitrary redox state, Journal of Geophysical Research: Planets, 129, e2024JE008576, 2024. SciX

  7. Spada, F., Demarque, P., Kim, Y.C. & Sills, A., The radius discrepancy in low-mass stars: single versus binaries, The Astrophysical Journal, 776, 87, 2013. SciX

  8. Baraffe, I., Homeier, D., Allard, F. & Chabrier, G., New evolutionary models for pre-main sequence and main sequence low-mass stars down to the hydrogen-burning limit, Astronomy & Astrophysics, 577, A42, 2015. SciX

  9. Noack, L. & Lasbleis, M., Parameterisations of interior properties of rocky planets, Astronomy & Astrophysics, 638, A129, 2020. SciX