Coupling to PROTEUS

This page is the theory of how Morrigan plugs into a PROTEUS coupled run: where it sits in the loop, what the wrapper does with its output, how one planet is chosen out of a whole system, and what each impact changes when it lands. For the practical TOML recipe, see the how-to page.

Morrigan runs once, not every iteration

Every other PROTEUS module is called inside the time loop and updates its slice of the planet's state on each pass. Morrigan is the exception. It runs once, before the loop starts, and what it returns is a schedule: a list of impacts with the times they occur and the physical parameters of each. The loop then consumes that schedule as it advances.

   ┌──────────────────────────────────────────────────────────────────────────┐
   │  init_accretion(handler)                    <-- Morrigan runs HERE, once │
   │      build_parameters(config)                   config + star.mass       │
   │      morrigan.run_system(**params)              evolve the whole system  │
   │      select_planet(survivors, config)           pick ONE survivor        │
   │      -> [ImpactEvent, ImpactEvent, ...]         its impact history       │
   │      validate_timeline(events)                  reject an impossible one │
   │      _drop_events_before_start(events, t0)      discard the pre-run tail │
   │                                                                          │
   │  while not done:                                                         │
   │                                                                          │
   │      pending = next_event(events, t)          # the loop reads the       │
   │      interior_o.t_next_impact = pending.time  # schedule, never the model│
   │                                                                          │
   │      run_interior(...)                        # cooling step; the        │
   │                                               # time-stepper shortens dt │
   │                                               # to land on t_next_impact │
   │                                                                          │
   │      hf_row['Time'] += dt                     # advance FIRST            │
   │                                                                          │
   │      for event in due_events(events, t_prev, t_now):   # impact falls due│
   │          apply_impact(handler, event)         # <-- consequences here    │
   │                                                                          │
   │      run_orbit(...) / run_escape(...) / run_outgassing(...) / ...        │
   └──────────────────────────────────────────────────────────────────────────┘

The time advance comes before the impact block, so the orbit, structure and escape steps of the step an impact lands on already see the grown planet. The clamp that shortens dt onto the impact time lives in the interior time-stepper, not in the main loop; the loop only publishes when the next impact is due.

The reason for the split is that the dynamical model and the coupled framework work on different clocks and at different cost. Morrigan integrates a whole system of embryos over hundreds of millions of years in a few seconds; PROTEUS integrates one planet's interior and atmosphere over the same span in hours. Running the dynamics once and replaying its schedule keeps the expensive loop in charge of the timestep.

A consequence worth stating plainly: the coupling is one-way. Nothing the coupled planet does feeds back into the dynamics. If the planet loses half its atmosphere to escape, the impact schedule does not change, because the schedule was fixed before the loop began.

Selecting one planet out of a system

Morrigan evolves a system of embryos and typically leaves several survivors. PROTEUS models one planet. Something has to choose, and that choice is the coupling's main free decision.

select_planet implements four rules, set by accretion.morrigan.selector:

Which end of the history a rule compares is what decides which body you get. match_config matches on the state the embryo started in, so it answers "which of these bodies began as the planet I configured"; mass and semimajoraxis match on the state it ended in, so they answer "which of these bodies became what I am looking for". A body can easily win one and lose the other.

The selected survivor's impact history becomes the schedule. Every other survivor is discarded. This is why num_planets and the initial mass distribution matter even though only one planet is followed: they set the dynamical environment that produced the one you keep.

From Morrigan's record to a PROTEUS event

build_parameters translates the TOML block into the arguments run_system expects:

The stellar mass is deliberately not a Morrigan setting. It is read from config.star.mass so the dynamical model and the rest of PROTEUS cannot disagree about the host star.

Each returned record becomes an ImpactEvent. The one transformation applied on the way is the time axis: Morrigan measures time from disk dispersal, PROTEUS from the start of its own evolution, and accretion.time_offset maps between them. time_offset is added to each impact time, and the result is compared against the run's start time.

Impacts that still land at or before that start time are discarded, with a warning naming the mass they would have added. Their mass is not folded into the initial condition and arrives nowhere: the configured planet.mass_tot and orbit define the starting state on their own. A schedule whose early impacts fall outside the simulated interval therefore grows the planet less than the dynamics described, so choose time_offset to bring the history you care about inside the run.

validate_timeline then rejects a schedule that cannot describe one body: times must increase strictly, each impact's target mass must follow from the previous merged mass, and a body may not gain mass between impacts. A drop of up to 10 % between impacts is allowed, since a consumer may strip an atmosphere in between.

The schedule becomes a property of the run

Once resolved, the schedule is written to impact_timeline.csv in the run's output directory, on the PROTEUS time axis with time_offset already applied. A resumed run reads that file instead of running the dynamics again, provided the file is still there; without it the run derives a fresh schedule.

This is not an optimisation. Re-deriving the schedule at every restart would reproduce the original history only if the dynamical model is bit-reproducible at a fixed seed, which the framework cannot check for a module it merely calls. Reading the file makes the impact history a property of the run rather than of the model's determinism, so a run that is stopped and resumed applies the same impacts it started with.

The same file is a portable artefact. accretion.module = "timeline" replays a timeline from a path, so a Morrigan-derived history can be re-run, shared, or handed to a collaborator who does not have Morrigan installed, and the impact consequences are computed identically either way. That is also the route for reproducing a published history, or for writing one by hand to isolate a single impact.

What an impact does

When the loop reaches a scheduled time, apply_impact applies seven consequences in a fixed order. The order matters: each step reads state the previous one wrote.

1. Size the volatile consequences from the pre-impact state, before anything changes. The erosion fraction, what the target loses, and what the impactor carries are all computed first, so no step sees a half-updated planet.
2. Grow the planet. planet.mass_tot increases by the impactor's rock alone, the impactor's mass less its volatile content. The volatiles arrive through their own channel below, so counting them here would double-count them.
3. Re-solve the structure. The interior module recomputes radius, core size, and pressures at the new mass.
4. Apply the volatile changes. The target loses its stripped fraction, the impactor delivers what survives, and the whole-planet element budget is refreshed.

5. Re-melt the mantle. The interior is reset to its molten initial condition, recomputed for the grown planet, and the next interior solve is told not to clip the resulting temperature jump. On Aragog the heat this injects is measured across the cooled-to-molten entropy jump and recorded as step_dE_impact_J, whenever a pre-impact entropy profile exists to measure it against. That column is added to both sides of the energy budget, so it quantifies the injection under a defined convention: E_residual_cons_J is unchanged across an impact for any booked value, and so cannot validate the magnitude itself. The relative residual is not invariant, because the injection enters its denominator; expect E_residual_cons_frac to spike on an impact row when the two are comparable, and read the absolute residual there instead. The impact's kinetic energy is logged separately; it is a different quantity, not the booked heat.

   Still on Aragog, how molten the reset state is depends on planet.temperature_mode: only liquidus_super solves for a profile guaranteed fully molten, raising if the requested superheat cannot be reached. adiabatic_from_cmb reaches one if tcmb_init is pinned high enough, but neither it nor planet.temperature_mode = "accretion" is checked against the liquidus. The remaining modes are only as molten as the configured initial condition, and an Aragog run warns about those at start-up. The scalar backends ignore temperature_mode entirely: they reset a surface temperature taken from planet.tsurf_init, and warn at the impact if that leaves the mantle short of fully molten.

6. Clear the solidification latch. A mantle that had crystallised is molten again, so the one-way latch is lifted. Without this step outgassing would stay frozen and the volatiles would be treated as locked in a solid mantle for good.

7. Move the orbit. Both orbital elements move by the change this impact made, not by the followed body's absolute values.

The orbit step is worth expanding. Morrigan's absolute orbits belong to its own system, which need not sit where the PROTEUS planet sits. So the coupling applies the ratio a_after / a_before to the semi-major axis and the difference e_after - e_before to the eccentricity: the planet keeps its configured orbit and inherits the dynamical model's kick. A ratio suits the semi-major axis because it is a scale; a difference suits the eccentricity because it is dimensionless and routinely zero, which a ratio cannot express. The result is clamped to a bound orbit. This is why the record reports both eccentricities and not only the post-impact one.

Both are bounded in the histories run_system returns, because a body carrying an unbound eccentricity is removed on the step it acquires one and never reaches the returned history. See the limitations for the case where the two branches of the model disagree about that, which is visible only to a caller reading the raw merger list.

What PROTEUS adds that Morrigan does not do

Morrigan reports bare bodies. Three pieces of physics live entirely on the PROTEUS side:

• Impact atmospheric erosion. Morrigan tracks no atmosphere and its merged masses are exact sums. PROTEUS computes what an impact strips, either as a fixed fraction or through the Kegerreis et al. (2020)¹ scaling law evaluated by zephyrus.collision.mass_loss from the record's own collision parameters. One fraction governs both bodies: the target loses that fraction of its atmosphere, and a volatile-bearing impactor loses the same fraction of its atmospheric part.
• Volatile delivery. Whether an impactor carries volatiles at all is a PROTEUS choice: dry, matching the planet's formation composition, or per-element budgets.
• The thermal consequence. Re-melting the mantle, recording the heat that the reset injects, and lifting the solidification latch is interior physics, not dynamics.

The impact record schema

out['impacts'] is a dictionary keyed by surviving-body id; each value is that body's impacts in time order, and every surviving body is present, with an empty list if it never merged. Each record carries:

time [yr], M_target_before [kg], M_impactor [kg], M_merged_after [kg], v_impact [m/s], v_esc [m/s], impact_parameter, R_target_before [m], R_impactor [m], rho_target [kg/m3], rho_impactor [kg/m3], a_before [m], a_after [m], e_before, e_after, id_target, id_impactor.

Conventions the schema guarantees:

• v_impact is the speed at first contact, \(\sqrt{v_\infty^2 + v_\mathrm{esc}^2}\), the convention the Kegerreis et al. (2020)¹ erosion law expects, and v_esc is the mutual escape speed recomputed from the recorded masses and radii, so the two are exactly consistent.
• M_merged_after is the perfect-merger sum of the two masses. The model sheds nothing, so a consumer applying its own erosion law starts from the unstripped mass.
• Bodies carry no atmosphere at all, in the record or in the model; masses and radii are bulk values, and the densities are recovered from mass and radius, so each record is exactly self-consistent.
• Successive records for one body chain: each M_target_before equals the previous record's M_merged_after to machine precision, since nothing removes mass between impacts.

out['survivors'] lists each surviving body's id together with its initial and final mass and semi-major axis, so the consumer can select a planet and replay its growth.

Any change to the record fields, units, or conventions is a breaking interface change for PROTEUS and must be flagged in the pull request that makes it.

Random-state hygiene

run_system seeds numpy's global random state for its own draws and restores the caller's state afterwards, so an embedding program's random sequence is unaffected by running a system in-process. This matters in a coupled run, where PROTEUS and other modules draw from the same global state.

See also

• Coupling to PROTEUS (how-to) for the TOML recipe, the selector choice, and the pitfalls.
• Model overview for the dynamics that produce the schedule.
• Limitations for what the model does not represent.
• On the PROTEUS side: the configuration listing for every parameter of the block, the model description for how the framework presents its modules, and the coupling loop for where the impact step sits among the others.
• The PROTEUS-side code: src/proteus/accretion/wrapper.py (dispatch and impact consequences), src/proteus/accretion/morrigan.py (this module's adapter), src/proteus/accretion/common.py (the ImpactEvent schema and timeline validation), src/proteus/config/_accretion.py (the config block).

References