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

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¶

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) ² 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 ³. 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 ⁹ (rotation-dependent tracks with activity-calibrated XUV luminosities, suitable for solar-type stars) and Baraffe ¹⁰ (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 ⁴. It uses parameterised solubility laws and an fO2 buffer (configurable as an IW offset) to compute the redox state ⁶. 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 ⁵. atmodeller provides an independent cross-check on CALLIOPE for the same volatile partitioning problem.

LavAtmos (Python) is an optional addition, enabled with outgas.vapourise = true, that computes the thermodynamic vapourisation equilibrium of the surface melt ⁷ and adds the resulting rock-vapour species to the volatile outgassing ⁸. Vapourisation is a distinct process from outgassing. It uses the FastChem equilibrium chemistry solver and requires ThermoEngineLite for melt thermodynamics.

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

Volatile outgassing and rock vapourisation are combined¶

With outgas.vapourise = true the following steps run sequentially, at each PROTEUS iteration:

• The volatile outgassing backend runs first on its own
• The vapourisation step then re-equilibrates the volatile atmosphere together with the rock vapour
• Returns one combined composition

 reset rock-vapour species and rock element reservoirs; M_vaps = 0
                          |
                          v
 [1] run_outgassing  (CALLIOPE / atmodeller / dummy)
     volatile equilibrium at fO2 = IW + outgas.fO2_shift_IW
     -> <gas>_bar, <gas>_vmr, <gas>_kg_*, atm_kg_per_mol
     -> P_surf = volatile total ;  P_vol = P_surf ;  P_vap = 0
     -> M_atm  = sum over gases of <gas>_kg_atm
     -> fO2_shift_IW_derived                       (volatile redox column)
                          |
          Phi_global < params.stop.solid.phi_crit ? ---- yes ---> stop here
                          | no                                   (no vapour)
                          v
 [2] run_vapourisation  (LavAtmos + FastChem)
     input  = melt at T = max(T_magma, outgas.lavatmos.T_min),
              P_volatile = P_vol from [1],
              element number fractions of: 
              (H, C, N, S, O + noble gases, from <e>_kg_atm),
              fO2 initial guess = 10 ** fO2_vapourise_derived (previous step)
     output = FastChem boa_chem: total pressure Pbar, mean mass mu,
              mixing ratios of ALL species
                          |
                          v
 [3] combine, in this order:
     P_vap  = max(Pbar - P_vol_from_[1], 0)        rock-vapour part
     P_surf = Pbar                                 combined total
     P_vol  = P_surf - P_vap                       volatile part (identity holds)
     M_atm  = Pbar * 4*pi*R_int^2 / gravity        hydrostatic column
     atm_kg_per_mol = mu                           combined MMW
     per gas:      <gas>_bar    = vmr * Pbar
                   <gas>_kg_atm = vmr * M_atm * W_gas / mu
                   <gas>_kg_total = 0 for rock-vapour species (kept in M_vaps),
                                    else atm + liquid + solid
     per element:  <e>_kg_atm = frac_e * M_atm * W_e / mmw_elements
                   O_kg_atm stays at the volatile-step
                   value and the vapour-derived excess goes to O_vapourised_kg
                   rock-forming elements accumulate into M_vaps
     fO2_vapourise_derived          = log10( vmr_O2 * (P_vol + P_vap) )
                                      absolute, log10 bar
     fO2_vapourise_shift_IW_derived = fO2_vapourise_derived - IW(T_magma)
                                      relative to the buffer, dex

Some notable consequences of step 3:

• Pressures add, compositions do not. P_surf = P_vol + P_vap, with P_vap defined as the excess of the combined FastChem total. The per-gas and per-element amounts are taken from the combined equilibrium, not summed across the two steps since adding them would count the volatile inventory twice.

• The vapourisation step reports its redox state twice, absolute and buffer-relative. The two must not be used interchangeably:

  Column	Meaning	Units
  fO2_vapourise_derived	Absolute fO2 of the combined atmosphere: the FastChem O2 mixing ratio times the total pressure P_vol + P_vap. Zero means 1 bar of O2.	log\(_{10}\) bar
  fO2_vapourise_shift_IW_derived	Relative to the iron-wustite buffer at T_magma: fO2_vapourise_derived - IW(T_magma). Zero means "exactly on the buffer", so it moves as T_magma evolves, even at fixed absolute fO2.	dex

  The absolute value is fed back to LavAtmos as the fO2 initial guess for the next iteration's LavAtmos call.

• The vapourisation fO2 is separate from the volatile-chemistry fO2. fO2_shift_IW_derived is the redox buffer offset obtained from the volatile step (e.g. CALLIOPE or atmodeller rather than by LavAtmos). It is not expected to equal fO2_vapourise_shift_IW_derived since the two come from different solvers.

Whole-planet mass is not conserved when vapourisation is enabled¶

• Rock vapour is added to the atmosphere without being removed from the interior. The vapourised rock mass is accumulated into M_vaps and enters M_atm, but no matching mass is subtracted from the interior: M_int, M_mantle.
• Rock-vapour elements dilute the escape outflow but are not depleted by it. Rock-forming elements take part in the unfractionated escape partitioning and reduce the escape rate available to H/C/N/O/S. This applies when escape.reservoir = "outgas", where the partitioning weights are atmospheric masses. With escape.reservoir = "bulk" the rock-forming whole-planet totals are reset to zero on every outgassing step, so they carry no weight and do not dilute the outflow.

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¶

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.

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¶

• Coupling loop: how the modules exchange data at each timestep
• Code architecture: source code layout and module patterns
• Configuration reference: all configuration parameters
• Observations configuration: synthetic spectrum options
• Tutorials: worked examples from dummy to production
• Bibliography: published references for PROTEUS and its modules