diff --git a/docs/adapters.rst b/docs/adapters.rst
index 0592f3e..21cfdb4 100644
--- a/docs/adapters.rst
+++ b/docs/adapters.rst
@@ -129,6 +129,10 @@ wrote, but always what the run used:
   omitted (four edges on a 2D lattice, six on a 3D one), and an edge whose
   condition takes parameters gets them (``u0 = 0.05`` for a
   ``velocity-inlet``).
+* a ``[multiphase]`` section whose EOS has a coexistence solver
+  (Carnahan–Starling) gets any omitted ``rho_l``/``rho_v`` from the Maxwell
+  construction of its own parameters, rounded to six significant figures; a
+  configured density is kept.
 * ``equilibrium`` defaults by phase: ``standard_equilibrium`` single-phase,
   ``wb`` multiphase; the source term is the one it registered as its pair.
 * the ``[hysteresis]``, ``[chemical_step]`` and ``[*_force]`` sections get
@@ -144,8 +148,9 @@ collision scheme (with ``k_diag`` required for MRT), the equilibrium,
 ``save_fields`` membership, every boundary-condition name and parameter key,
 obstacle geometry, each force against its schema, the wetting, hysteresis
 and chemical-step sections, and — for multiphase — the common keys, that
-``rho_l > rho_v``, and the further parameters of the selected EOS (the four
-Carnahan–Starling ones).  Only one gravity force may be configured.  A 3D
+``rho_l > rho_v`` (both required for a double-well EOS), and the further
+parameters of the selected EOS (the four Carnahan–Starling ones).  Only one
+gravity force may be configured.  A 3D
 lattice refuses what is still two-dimensional: wetting, hysteresis, chemical
 step, the electric force, MRT and an obstacle.
 
diff --git a/docs/notes/REFACTOR_PLAN.md b/docs/notes/REFACTOR_PLAN.md
deleted file mode 100644
index 3466b8e..0000000
--- a/docs/notes/REFACTOR_PLAN.md
+++ /dev/null
@@ -1,113 +0,0 @@
-# Refactor plan: feature-slot architecture
-
-Structure-only refactor. Internal physics (stencils, contact-angle maths, hysteresis optimiser) is untouched.
-
-## 1. Decisions
-
-| # | Decision | Outcome |
-|---|---|---|
-| 1 | Static vs dynamic wetting | Wall parameters always live in State, even for static runs. One kernel path serves static wetting, hysteresis and trial steps. |
-| 2 | Wetting physics | A modification of the differential operators' padded region, implemented as `ghost_fill`. |
-| 3 | Coupling mechanism | Fixed, typed slots. Features register like operators (kind `feature`); `build_setup` merges every active feature's contribution into the slots. No resolver, no shared context object. |
-| 4 | Step | One skeleton: `kernel` plus `step`. `sim_type` no longer selects step functions; the features present in the config decide. |
-| 5 | Coupling cost | Plain operators and features: a new file with its registry decorator (auto-discovered). A feature or force with its own TOML section also needs one field on `SimulationConfig`, as forces do today (removing that cost is a follow-up, see open items). `build_setup` is edited only to add a new kind of slot. |
-| 6 | State | Core fields plus `ext: dict[str, pytree]` keyed by feature name. Each feature owns its substate. |
-| 7 | Reference pressure | Rule, documented but not implemented in this refactor: pressure comes from rho inside equilibrium and source, never from a force. The zhang2022 files (split off in 4c36d1d to `.claude/split/zhang2022_files`, planned `feat/zhang2022`, including `ViscosityParams`) are out of scope and get rebuilt on the new slots later. |
-| 8 | Chemical step | Part of the wetting feature (hysteresis parameters vary with contact-line location). |
-| 9 | Config | Structure unchanged. `SimulationConfig` parses and validates; all guarding stays there. Section fields (`wetting_config`, `hysteresis_config`, `chemical_step_config`, `*_force`) and the derived properties read by analysis (`wetting_defaults`, `chemical_step_wall`, `restored_wetting`) stay as they are; features read them from the config. |
-| 10 | Scope | Structure only. Internal physics untouched. |
-| 11 | Registry | Stays. Only the kinds the feature slots replace are removed: `update_timestep`, `wetting` (hysteresis variants), `extra_state`. Everything else (plotting/analysis discovery, dimensionless numbers, EOS capability checks, force config schemas, CLI listings) keeps using it. |
-| 12 | `sim_type` | Kept, with a reduced meaning. All five values still load (saved configs, analysis, CLI and adapters read it). `single_phase` vs `multiphase*` selects the macroscopic operator (and the equilibrium/source pair where relevant). Wetting, hysteresis and chemical step are switched on only by their config sections. |
-
-## 2. Target architecture
-
-```
-config.toml -> SimulationConfig (parse and validate) -> build_setup (recipe) -> Setup -> run (scan over step)
-```
-
-**Vocabulary**
-- **Feature:** a module that adds behaviour touching several parts of the simulation (wetting, electric potential). Returns a `Contribution`, or `None` when its config section is absent.
-- **Slot:** one of the four fixed places a feature can contribute to.
-- **Plain operator:** a homogeneous, selectable component (BC, EOS, collision, force, initialiser). Needs no coupling.
-
-**Contribution slots**
-
-| Slot | Shape | Wetting example |
-|---|---|---|
-| `ghost_fill` | `(padded, ext) -> padded` | Rewrites the wall ghost nodes between padding and the stencil. |
-| `forces` | tuple of force operators | Electric force. |
-| `post_step` | `(setup, prev_state, new_state) -> new_state` | Contact-angle measurement, hysteresis optimiser, electric potential update. |
-| `init_state` | `(setup) -> pytree` | Initial wetting parameters and angles. |
-
-**Feature registration.** A feature registers with a decorator under the registry kind `feature` and is auto-loaded like any operator, so the CLI lists it with no extra work. `build_setup` calls `feature.build(config, lattice)` for every registered feature; it returns a `Contribution`, or `None` when the feature's config section is absent. Merge rules: `ghost_fill` functions are composed, `forces` are concatenated, `init_state` results are keyed into `ext` by feature name, and `post_step` hooks run in name order. Features must be independent of each other. If one needs another's output, that is a real coupling and needs its own slot (an edit to `build_setup`).
-
-**Why `post_step` receives both states.** The electric potential is a second lattice Boltzmann equation (Zhang & Kwok 2005, Eqs. 15-18): `h(t)` is collided with `tau_h = 3*sigma + 0.5`, where sigma comes from the density at the same time level, then streamed to `h(t+1)`. The existing plugin does this with `prev_state.f` and `prev_state.h`. The kernel leaves `h` untouched, so `h(t)` would still be in `new_state`, but `f(t)` is gone; taking sigma from `new_state.f` would pair `h(t)` with `sigma(t+1)` and change the result. Wetting needs the other side: it measures contact angles on `new_state.rho` and reads `prev_state.wetting`. One signature serves both. Hooks read forces from `new_state.force_ext`; there is no `**context`.
-
-No `pre_step` slot for now: hysteresis parameters live in `ext` and the kernel reads them. Add the slot only when a feature needs it.
-
-**Differential.** One builder, `build_differential(lattice, pad_modes, ghost_fill=None)`, returning `op(field, ext)`. The plain stencil ignores `ext`. The package never imports `wetting`.
-
-**Kernel and step**
-- `kernel(setup, f, force_ext, ext) -> (f, rho, u, force_tot, pressure)`: macroscopic, equilibrium, collision + source, streaming, obstacle, BC, macroscopic. Passes `ext` through to the differential operators; knows no feature by name.
-- `step(setup, state)`: compute forces, call kernel, write new State, run the merged `post_step` hooks (features in name order, see feature registration).
-- The hysteresis hook calls `setup.kernel` with `ext | {"wetting": candidate}` for its trial steps. Hooks receive `setup` at call time, so there is no circular dependency.
-- Hysteresis stays post_step to preserve current behaviour (optimiser runs on post-step populations and force, writes parameters for the next step).
-
-**Registry stays.** It is the catalog of plain operators and is load-bearing: the CLI lists and selects from it (`cli/display.py`, `visualise.py`, `analysis.py`), `SimulationConfig` validates names and reads force metadata from it, and capability checks use kind membership (e.g. an EOS has a `pressure` or `surface_tension` entry or it does not). Decorator registration plus `auto_load_operators` is what makes a new BC, EOS or force a one-file addition. Only the entries that stop being selectable operators leave it (see section 3).
-
-**build_setup documentation requirement.** A header comment states the pattern: plain operators and features are auto-discovered through the registry, a feature fills the four slots, features are merged in name order and must be independent, and `build_setup` is edited only to add a new kind of slot. Each step of the recipe (lattice, feature merge, differential, forces, operators, step) gets a one-line comment saying what it does and which slots it consumes. There is no per-feature line, because there is no feature list.
-
-## 3. Deleted
-
-- Differential: `build_diff_ops`, `build_wetting_gradient_fn`, `build_wetting_laplacian_fn`, `_gradient_wetting.py`, `_laplacian_wetting.py`.
-- Step: `step_multiphase_wetting`, `step_multiphase_hysteresis`, `step_multiphase_hysteresis_chemical_step`, `with_default_wetting`, `_make_wetting_ops`.
-- Setup: any equilibrium/source special-casing (`isinstance` checks, by-name lookups) and `partial` binding hacks still present on the branch; verify in step 0.
-- State: typed `wetting` and `h` fields, extra-state plugin registry, `pipeline/state/_extra_state.py`.
-- Registry kinds that stop being selectable operators: `update_timestep` (step variants), `extra_state` (plugins), and the `wetting` hysteresis-variant kind; `simulation_type` pending a check of its consumers. The registry itself, `_loader.py` and the rest of the kinds stay.
-- The ~20 `if x is None: raise TypeError` guards across the step modules.
-
-## 4. Execution order
-
-Each step is completed fully, including removal of replaced code, before the next. No tests are run until step 6.
-
-**Step 0: Preparation**
-- Read files not yet seen: `_hysteresis.py`, `_applicator.py`, `simulation_config.py`, `runner.py`.
-- Capture baseline outputs from the current code, one case per path: single-phase, multiphase, static wetting, hysteresis, chemical step, electric. This is the only run before the refactor.
-
-**Step 1: State and step skeleton**
-- Add `ext` to `State`; remove typed `wetting` / `h` fields.
-- Write `kernel` and `step`.
-- Delete all step variants and `_common.py`.
-
-**Step 2: Features and ghost_fill**
-- Define `Contribution`.
-- Wetting becomes a feature: applicator -> `ghost_fill`; contact-angle measurement and hysteresis optimiser -> `post_step`; chemical step folded in; `init_state` builds initial parameters.
-- Electric force and potential become a feature.
-- Collapse `differential` into one builder.
-
-**Step 3: Dropped.** The zhang2022 files are not on this branch. The rule "pressure comes from rho inside equilibrium/source, never from a force" goes into the docs (step 7). zhang2022 is rebuilt on the new slots later.
-
-**Step 4: build_setup and registration**
-- Recipe in fixed order: lattice, active features (loop over the registry kind `feature`, merge slots), differential, forces (config forces plus feature forces), collision/equilibrium/streaming/BC/obstacle, step.
-- Header and per-step comments as specified above.
-- Remove the obsolete registry kinds listed in section 3; the registry otherwise stays as is.
-
-**Step 5: Config**
-- `SimulationConfig` keeps parsing and validation (BCs, forces, multiphase, MRT coupling, obstacle, and so on).
-- No change to section fields, derived properties, sweeps, overrides or adapters.
-- `sim_type` stops selecting step functions; it keeps only the phase meaning (decision 12). The check "hysteresis config requires a hysteresis `sim_type`" goes away, since config sections switch the features on.
-
-**Step 6: Tests and verification**
-- Delete tests that only cover the old structure.
-- Rewrite the rest for the new design and invariants (slot wiring, `ext` structure, kernel purity).
-- Run everything, compare against step 0 baselines, fix only refactor-introduced issues.
-
-**Step 7: Docs and cleanup**
-- Update `docs/architecture.rst` and `docs/operators.rst`, including the rule that pressure is resolved from rho inside equilibrium/source and never supplied by a force.
-- Check whether `build/lib` is tracked; remove if it is a stale duplicate.
-
-## 5. Open items
-
-- **Hook order:** extra-state plugins currently run alphabetically; the feature merge runs hooks in name order. Verify that wetting measurement and the optimiser (same feature, ordered inside it) keep their present order, and that the current features (wetting, electric) are independent.
-- **Follow-up, out of scope (option B):** a generic `sections` mechanism so features and forces declare their TOML section at registration and `SimulationConfig` needs no new field. It touches sweeps (`array_expansion`), `cli/overrides.py`, the adapters, the saved-config round trip and the analysis readers, so it is a separate project.
-- **Tooling:** a shell on the user's machine is needed for baselines (step 0) and tests (step 6).
diff --git a/docs/notes/REVIEW_FIX_PLAN.md b/docs/notes/REVIEW_FIX_PLAN.md
deleted file mode 100644
index 863f82f..0000000
--- a/docs/notes/REVIEW_FIX_PLAN.md
+++ /dev/null
@@ -1,174 +0,0 @@
-# Plan: code-review fixes (review notes of 23 and 26 October 2026)
-
-Follows the feature-slot refactor (`docs/notes/REFACTOR_PLAN.md`), which is currently uncommitted on `fix/bubble_simulation` (65 files). The project rules still apply: structure first, behaviour preserved unless a change is stated as intended, replaced code removed in the same step, tests adapted and run only once a phase is complete.
-
-## 1. Decisions
-
-| # | Topic | Decision |
-|---|---|---|
-| 1 | Config ownership | Centralised in the `config` package, organised to mirror the operator packages: one module per section (`config/boundary.py`, `config/multiphase.py`, `config/forces.py`, `config/wetting.py`, ...) holding that section's schema, validation and derived objects. Operator packages hold no config knowledge. `SimulationConfig` holds the core fields and one field per section, and runs all validation when it is constructed. |
-| 2 | Old saved configs | Converted on load by one isolated legacy adapter (`src/config/legacy.py`), called from one place. Removing it at the next version bump is deleting that file and its one call. |
-| 3 | Defaults | Every default value appears once, in `src/config/defaults.py`. The section modules read it; no operator or addon package declares a default. |
-| 4 | 3D | Build 3D support (phase 5, scope below). |
-| 5 | Equilibrium / source term | Selected by name from the config, not hard-coded by phase in `build_setup`. Which source term belongs to which equilibrium is physics, so an equilibrium declares its paired source term at registration; the default equilibrium per phase is a value in `config/defaults.py`. |
-| 6 | Aggregators | The obstacle becomes the first entry of the boundary chain (same `(f_stream, f_col) -> f` signature). Forces keep their own aggregator: they are summed from the state, a different operation. |
-| 7 | Protocols | One file, `src/protocols.py` (moved from `operators/protocols.py`), holds every interface in the codebase: the operator protocols, the addon hook types (`GhostFill`, `PostStepHook`, `InitStateHook`, now in `features/contribution.py`), `DimensionlessNumberOperator` (`analysis/physical_parameters/_inputs.py`) and the optax shim `_OptaxLike` (`wetting/_hysteresis.py`). The unused `SimulationRepository` and `ConfigReader`, and the `PlotOperator` Protocol that duplicates the plotting ABC, are deleted. The four ABCs (`PlotOperator`, `AnalysisPlot`, `ConfigAdapter`, `OutputWriter`) move there as whole classes, shared implementation included; the file therefore holds some implementation, accepted. Constraint: `src/protocols.py` imports only the standard library, `jax` and `numpy` at runtime; everything else (matplotlib, `SimulationConfig`, `State`, ...) only under `TYPE_CHECKING` or inside methods, so importing the core never pulls in plotting and no import cycle forms through this file. |
-| 8 | `sim_type` | Removed. A run is multiphase when it has a `[multiphase]` section. Old configs carrying `sim_type` are handled by the legacy adapter. |
-| 9 | Naming (D1) | Feature becomes **addon** (`src/addons/`, `@addon_operator`, registry kind `addon`); `Contribution` becomes **`Hooks`**. |
-| 10 | Differential registry (D2) | Not added. There is one stencil; a registry is added when a second one exists. |
-| 11 | Legacy adapter writes (D4) | Yes: an old-layout `config.toml` is converted on load and the converted file is written back to its results folder. |
-| 12 | General initialiser (D5) | `multiphase_bubbles` (`operators/initialise/_multiphase_bubbles.py`). |
-| 13 | 3D scope (D6) | Single-phase and multiphase bulk runs in 3D. Wetting and electric refuse 3D with a clear error. |
-| 14 | One-line builders (D3) | Kept. Each package keeps its named `build_x` entry point typed with its protocol, including `build_plot_operator` and `build_analysis_plot` in plotting; this is the same pattern as every operator package, so it stays consistent. |
-
-## 2. Review notes mapped to the plan
-
-| Note | Where |
-|---|---|
-| Differential registry | Not added (decision 10) |
-| Plotting builders / protocol location | Builders kept (decision 14); all protocols in `src/protocols.py`, duplicates removed, phase 2 (decision 7) |
-| `front`/`back` missing from BC application order | Phase 5 |
-| `build_boundary_edges` in config; config mirroring operators; chemical step in wetting | Phase 3: `config/boundary.py`, `config/wetting.py` (with hysteresis and chemical step) |
-| Variable input params (multiphase) | Phase 3: `config/multiphase.py` holds a per-EOS parameter table |
-| `@property`s on `SimulationConfig` | Phase 3 |
-| `wetting_defaults` and other constants in one overview | Phase 3 (decision 3) |
-| Initialise legacy examples | Phase 4 |
-| Two functions in `_multiphase.py` | Phase 2 |
-| Central aggregator | Phase 2 (decision 6) |
-| `GhostFill` typed with `Any` | Phase 1 |
-| Naming of feature / Contribution | Phase 1 (decision 9) |
-| `architecture.rst` extension points | Phase 7 |
-| Fixed `sim_type` list | Phase 3 (decision 8) |
-| Forces named in the config | Phase 3: one `forces` field holding every `[*_force]` section, schemas in `config/forces.py` |
-| `lattice` passed to feature builders | Phase 1: removed, both builders ignore it |
-| `forces` slot in `Contribution` | Phase 1: removed, no feature fills it |
-| Eq/source split by phase limits new operators | Phase 2 (decision 5) |
-| Static `rho_l`/`rho_v` in wetting | Phase 6 |
-| `State` storing derived fields | Already done: `State` is `(f, t, ext)` and `observe()` derives the rest. No work. |
-
-Found during the review, not in the notes: `config/__init__.py` `from_toml` has its lazy import inside a comment (a literal `\n`); fixed in phase 2.
-
-## 3. Target config design (phase 3)
-
-The config package is the only place that knows about config: section layout, schemas, defaults, validation and the objects derived from a section. It is organised to mirror the operator packages, so finding the config for an operator means opening the module with the same name.
-
-```
-src/config/
-  simulation_config.py  SimulationConfig: core fields + one field per section; __post_init__ runs every section's validation
-  defaults.py           every default value, once
-  boundary.py           [boundary_conditions]: schema, validation, edges in application order, pad modes, periodic axes
-  obstacle.py           [obstacle]: schema, validation, mask
-  multiphase.py         [multiphase]: schema, per-EOS parameter table, MultiphaseParams, phase references
-  forces.py             [*_force]: one schema per force, active forces
-  wetting.py            [wetting], [hysteresis], [chemical_step]: schemas, cross-section rules, wall parameters, chemical-step wall, restored state
-  legacy.py             old-layout upgrade (see below)
-  adapter_*.py, array_expansion.py, run_config.py, config_overview.py   unchanged roles
-```
-
-**SimulationConfig** holds:
-- Core fields: `simulation_name`, `lattice_type`, `grid_shape`, `tau`, time stepping, collision, initialisation, output.
-- One field per section, named like its TOML section: `boundary_conditions`, `obstacle`, `multiphase`, `wetting`, `hysteresis`, `chemical_step`, and `forces` (a mapping holding every `[*_force]` section). The old `*_config` names and the per-force fields go away.
-- No `@property`s. The derived objects are functions in the section modules, e.g. `config.boundary.edges(cfg)`, `config.wetting.chemical_step_wall(cfg)`, `config.multiphase.params(cfg)`. Setup, analysis (`wetting_overlay`, `setup_figure`) and the CLI call those.
-
-**Schemas and validation.** Each section module defines its schema (required, optional, positive keys) and its cross-section rules (hysteresis needs wetting, the chemical-step edge must be the wetting wall, only one gravity force), reading defaults from `defaults.py`. The force schemas now on the `@force_model` registrations move to `config/forces.py`. The config still asks the registry which names exist (BCs, EOS, forces, initialisers) to validate a configured name; that is the only link from config to the operators.
-
-**What moves where:**
-
-| Today | Becomes |
-|---|---|
-| `config/boundary_edges.py`, the `boundary_edges`/`pad_modes`/`periodic_axes` properties | `config/boundary.py` |
-| `config/obstacle_mask.py`, the `obstacle_mask` property | `config/obstacle.py` |
-| `config/multiphase_params.py`, the multiphase fields, `phase_references`, the Carnahan–Starling special case in `_validate_multiphase` | `config/multiphase.py`; `a_eos`…`t_eos` only in the Carnahan–Starling row of the per-EOS table |
-| `gravity_force`/`electric_force`/`gravity_masked_force` fields, force metadata on `@force_model`, `active_forces` | `forces` field; `config/forces.py` |
-| `config/wetting_defaults.py`, `config/chemical_step.py`, `_HYSTERESIS_DEFAULTS`, `_CHEMICAL_STEP_DEFAULTS`, the wetting/hysteresis/chemical-step validation, `wetting_defaults`/`chemical_step_wall`/`restored_wetting` | `config/wetting.py`; values to `defaults.py` |
-| `sim_type`, `is_multiphase`, `is_single_phase` | multiphase means a `[multiphase]` section is present (`config.multiphase.is_multiphase(cfg)`) |
-
-**Sweeps, adapters, overrides.** Sections keep being nested-sweepable fields; with field names equal to section names, the `cli/overrides.py` alias map goes away. The adapter routes every `[*_force]` TOML section into `forces`, and writes it back the same way.
-
-**Extension cost, stated plainly.** Adding a force or an addon with its own TOML section touches its operator/addon file plus the config package (its schema in the section module, its values in `defaults.py`). That is the price of keeping all config in one place.
-
-**Legacy adapter:** `src/config/legacy.py`, `upgrade(raw) -> raw`. It detects the old layout (`sim_type`, `*_config` and per-force keys, the old multiphase keys) and rewrites it to the new one. It is called once, in the adapter's load path, and nothing else imports it. When it converts a file, the loader writes the converted `config.toml` back to the same path (decision 11), so each old results folder is converted once.
-
-## 4. Execution phases
-
-Each phase is completed fully before the next. Order rationale: renames first so files are touched once; config before 3D because the BC application order moves into `config/boundary.py`; the one intended physics change (phase 6) last, so phases 1–5 can be checked against unchanged baselines.
-
-**Phase 0: Baseline**
-- Commit the current refactor so this plan starts from a clean tree.
-- Capture baselines, one case per path: single-phase, multiphase, static wetting, hysteresis, chemical step, electric, obstacle, and one sweep config.
-- Keep one results folder saved with the old config layout, to test the legacy adapter.
-
-**Phase 1: Feature slots**
-- Rename feature to addon and `Contribution` to `Hooks` (decision 9) across code, registry kind, docs and the `build_setup` header comment.
-- Remove the `forces` slot and the merge code for it in `build_setup`.
-- Remove the `lattice` parameter from feature builders.
-- Typed ghost fill: a feature's fill receives only its own substate (`GhostFill[S] = (padded, S) -> padded`); the merge layer passes `ext[name]`. Electric's substate becomes a NamedTuple (`ElectricState`) instead of `{"h": ...}`. This also enforces the independence rule for the fill.
-
-**Phase 2: Operator cleanups**
-- `macroscopic/_multiphase.py`: one function. The multiphase parameters and EOS are bound at build time, so both macroscopic operators share one signature and the phase branch in `step._moments` goes away.
-- Equilibrium and source term selected by config name with registered pairing (decision 5); `build_setup` names no equilibrium or source.
-- Obstacle joins the boundary chain (decision 6); `operators/obstacle` keeps only its geometry.
-- Protocols (decision 7): move `operators/protocols.py` to `src/protocols.py`; move in the addon hook types, `DimensionlessNumberOperator`, `_OptaxLike` and the four ABCs (`plotting/base.py`, `config/adapter_base.py`'s `ConfigAdapter`, `output_data/base.py`'s `OutputWriter`); delete the emptied modules and the three unused protocols; update every import; check the import constraint by importing `src.pipeline` without matplotlib loaded.
-- Fix the `from_toml` comment line.
-
-**Phase 3: Centralised config** (section 3)
-1. `config/defaults.py`: move every default into it (wetting scalars, hysteresis, chemical step, force defaults, BC defaults, default equilibrium per phase).
-2. Section modules, one at a time, each with its schema, validation and derived functions: `boundary.py`, `obstacle.py`, `multiphase.py`, `forces.py`, `wetting.py`. Delete the old `config/*` modules each one replaces.
-3. `SimulationConfig`: one field per section, `__post_init__` calls the section validations; remove `sim_type` and every `@property`.
-4. Remove the schema metadata from `@force_model` registrations.
-5. Sweeps, adapters, overrides on the new field names; the adapter's `[*_force]` routing.
-6. Update readers in setup, addons, analysis and CLI to the section-module functions.
-7. Legacy adapter, tested on the old-layout results folder from phase 0.
-8. Check that two-level sweeps (`forces.gravity_force.force_g = [...]`) expand; `array_expansion` handles one nesting level today.
-
-**Phase 4: Initialisers**
-- For each legacy initialiser (`wetting`, `wetting_drop_top`, `wetting_chem_step`, `multiphase_bubble_top`), show that the general initialiser reproduces its field (array equality on a few parameter sets). The general one is `multiphase_bubbles` (decision 12).
-- Remove those that are reproduced; migrate the example configs, tests and `cli/wetting_init.py` to the general initialiser.
-- Any that cannot be reproduced stay, with the reason recorded here.
-- Outcome: all four removed; none stay. `tests/operators/initialise/test_legacy_equivalence.py` pins array equality against each removed formula on three parameter sets. The three wall droplets need `interface_width` halved in `[initialisation]`, because their tanh was twice as steep. Other grids and densities differ by up to 2e-15, from two rounding causes: a radius fraction times `min(nx, ny)` cannot always hit the old radius, and `avg - amp` can land an ulp below `rho_v`, which `multiphase_bubbles` clamps. No run under `~/TUD_LBM_data` used any of the four, so the legacy adapter maps nothing.
-
-**Phase 5: 3D support** (scope fixed by decision 13)
-Today 2D-only: BC application order (4 edges), every BC, streaming periodicity (x, y), padding (4 faces), the gradient and Laplacian stencils (`to_2d`, `grad_core_2d`), the well-balanced source term, the electric force, all wetting code, and the wetting initialisers.
-Scope:
-- `front`/`back` in the application order and in every BC.
-- z-periodicity in streaming.
-- 6-face padding and 3D gradient/Laplacian stencils.
-- Single-phase and multiphase bulk runs in 3D, including the well-balanced source.
-- Wetting, electric and the 2D initialisers refuse 3D with a clear error instead of silently running.
-Verification: 3D Poiseuille flow against the analytic profile; a 3D droplet against the Laplace law.
-
-**Phase 6: Phase densities in wetting** (intended behaviour change)
-Today `rho_l`/`rho_v` are the config values in three places: the applicator, the `rho_mean` contact-angle threshold, and the hysteresis clamps (`_hysteresis.py`). Options:
-- (a) keep the config values;
-- (b) the phase densities measured off the initial field (what `phase_references` already does), static for the run;
-- (c) measured every step and carried in the wetting substate.
-Run (b) and (c) against the phase 0 wetting and hysteresis baselines, show the differences, and decide with you before changing the default.
-
-*Outcome (2026-11-04): (b), decided with you.* The densities enter in two places, not three: the `rho_mean` threshold (applicator anchors, addon measurement, hysteresis `_wall_constants`) and the applicator's contrast gate. The hysteresis clamps depend only on `W`. (b) differs from (a) only on an `init_from_file` run, so each case was equilibrated with fixed wetting for 10,000 steps, then restarted under each mode (driver and results in `.claude/review_fix/phase_densities.py` and `.claude/review_fix/phase6/`).
-
-| Case (3000 steps; real: 10,000) | max \|Δθ\| (b / c) vs (a) | max \|Δcll\| (b / c) | inclusion mass, relative (b / c) |
-|---|---|---|---|
-| wetting | 0.03° / 0.03° | 0.02 / 0.02 cells | 0 / 0 |
-| hysteresis | 0.50° / 0.56° | 0.00 / 0.00 | 4.5e-4 / 4.5e-4 |
-| hysteresis bubble | 0.29° / 0.34° | 0.04 / 0.05 | 1.2e-3 / 1.2e-3 |
-| chemical step | 0.25° / 0.19° | 0.02 / 0.03 | 3.2e-4 / 3.2e-4 |
-| CS chem-step bubble, 201x101, t ≥ 50000 (up to 56000) | 6.5° / 5.4° | 1.8 / 0.7 | 5.2e-2 / 2.6e-2 |
-
-Total mass drift stayed at round-off in every run. The real bubble (init densities 0.0083/11.645 against configured 0.015/12.18) dissolved in all three modes on the same schedule (area 1624 → ~100 by t = 57000): the closed-box hydrostatic loss, which no choice of threshold changes. (c) ran 6.5x slower there (26,788 s against ~4,100 s), and its references lost their meaning once the inclusion detached (vapour reference 5.97). Implemented as `config.wetting.wall_rho_mean` (the midpoint of `phase_references`), read by the addon, the hysteresis optimiser and the contact-angle overlay; the applicator takes `phase_references` directly. Nothing of (c) entered `src`. Baselines are unchanged bit for bit, because none of them restarts from a file.
-
-**Phase 7: Tests, docs, verification**
-- Delete tests that only cover replaced structure; adapt the rest to the new design.
-- Run everything; compare phases 1–5 against the phase 0 baselines exactly.
-- `architecture.rst`: extension-points table gains rows for an addon, a force with its config section, an EOS with its parameters, and a config section, each naming the operator file and the config-package edits; describe the config package layout, `defaults.py` and the legacy adapter.
-- `operators.rst`: equilibrium/source pairing, obstacle in the boundary chain.
-
-*Outcome (2026-11-04).* Docs: `architecture.rst` (layout, extension-points rows for an addon, a force, an EOS, a BC and a config section; a section on the config package, `defaults.py` and the legacy adapter), `operators.rst` (kinds table, equilibrium/source pairing, obstacle in the boundary chain, the symmetry fix, the force API, the wall threshold, 3D scope), `adapters.rst` (section mapping, defaults, old layouts, three-level sweeps) and `quickstart.rst` (no simulation type, `observe`, the in-memory trajectory's memory cost). CLAUDE.md rewritten to match, including the drift it carried before this plan: it said there was no `macroscopic` kind, that `MacroscopicOperator` is generic and also the pressure type (pressure is `EosOperator` of `rho`), and it referred to the deleted `step/_common.py`. Found and fixed during the sweep:
-- The `typing.cast` ban CLAUDE.md describes was never configured. `TID251` is now on (`[tool.ruff.lint.flake8-tidy-imports.banned-api]`), and its three users are fixed at the source: `tud-lbm run` takes explicit typed parameters instead of `**cli_kwargs: object` (the narrowing helpers `optional_int`/`str_tuple` are deleted), `regime-map` annotates its `click.Choice` values with their `Literal` types, and `run_parallel` types its parameter list as a covariant `Sequence`.
-- The velocity-inlet `u0` default moved into `defaults.py` (`BOUNDARY_PARAMETERS`) and is filled into the edge's parameter key, so the saved `config.toml` shows it. The MRT free rates stay in `_mrt.py`: the config refuses MRT without `k_diag`, so no run reaches them, and the conformance test calls the unbound operator.
-- `observe` is now exported from `src`; the package docstring described `State` as `(rho, u, f)`.
-- Tests that still named `sim_type`, `bc_config` or `force_enabled` were renamed to what they check; nothing pinned replaced structure. Nothing deselects the `slow` marker, so every gate run already included those tests.
-- Left as found: an explicit `save_fields` still always gains `"f"`, though nothing reads a snapshot's `f` any more (restarts use `rho`/`u`, pressure comes from `rho`). Removing it changes what every run writes, so it is reported, not changed.
-
-## 5. Open decisions
-
-None. All decisions are recorded in section 1 (decisions 9–14 settle D1–D6).
diff --git a/docs/operators.rst b/docs/operators.rst
index 9ec75a6..5594301 100644
--- a/docs/operators.rst
+++ b/docs/operators.rst
@@ -56,6 +56,7 @@ Kind                        Decorator                       Resolved by
 ``eos``                     ``@eos_operator``               ``multiphase.eos``
 ``pressure``                ``@pressure_operator``          ``multiphase.eos`` (if it has one)
 ``surface_tension``         ``@surface_tension_operator``   ``multiphase.eos`` (if it has one)
+``coexistence``             ``@coexistence_operator``       ``multiphase.eos`` (if it has one)
 ``boundary_condition``      ``@boundary_condition``         ``config.boundary_conditions`` per face
 ``force``                   ``@force_model``                the keys of ``config.forces``
 ``stream``                  ``@stream_operator``            fixed (``"standard"``)
@@ -174,6 +175,32 @@ Equation of state — ``eos``
     densities and positive in the unstable region between them — the
     double-well character is in the shape, not the sign.
 
+Coexistence densities — ``coexistence``
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+An EOS whose liquid and vapour densities follow from its own parameters
+registers a solver here under its own name; membership is the capability
+check, and :func:`~src.operators.macroscopic.eos.coexistence_densities`
+returns ``(ρ_v, ρ_l)``, or ``None`` for an EOS outside the kind.
+
+``carnahan-starling``
+    The Maxwell construction: equal bulk pressure and equal chemical
+    potential, i.e. equal areas under ``p₀(1/ρ)`` about the tie line
+    ``p_sat``.  Solved by nested bisection under ``jax.jit`` from the two
+    spinodals of ``dp₀/dρ``.  A ``t_eos`` at or above the critical
+    temperature raises ``ValueError``.  ``rho_l``/``rho_v`` may be omitted
+    from ``[multiphase]``; the config then fills them with this pair, rounded
+    to six significant figures.  A configured value is never overwritten.
+
+``double-well`` registers no solver: ``rho_l`` and ``rho_v`` are its
+parameters, so both are required.
+
+A multiphase run with a solver writes ``coexistence/data/data.json`` (the
+pair, ``p_sat`` and the configured densities) and
+``coexistence/plots/maxwell.png`` into its run directory at start, and
+``physical_parameters.txt`` gets a *Maxwell coexistence* row.
+``tud-lbm analyse config.toml --coexistence`` writes the same without a run.
+
 Boundary conditions — ``boundary_condition``
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
diff --git a/docs/quickstart.rst b/docs/quickstart.rst
index d0b7804..4aee1cc 100644
--- a/docs/quickstart.rst
+++ b/docs/quickstart.rst
@@ -178,7 +178,9 @@ and they combine:
 
 ``[multiphase]`` needs ``eos``, ``kappa``, ``rho_l``, ``rho_v`` and
 ``interface_width``; ``eos = "carnahan-starling"`` additionally needs
-``a_eos``, ``b_eos``, ``r_eos`` and ``t_eos``.  The wetting sections need
+``a_eos``, ``b_eos``, ``r_eos`` and ``t_eos``, and may omit ``rho_l``/``rho_v``:
+they are then the Maxwell coexistence pair of those four (omit them to sweep
+``t_eos``).  The wetting sections need
 ``[multiphase]`` and a ``"wetting"`` face.  A 3D run (``lattice_type =
 "D3Q19"``) supports the bulk physics only: see :doc:`operators`.
 
diff --git a/src/cli/commands/analysis.py b/src/cli/commands/analysis.py
index f2aa4dd..5836530 100644
--- a/src/cli/commands/analysis.py
+++ b/src/cli/commands/analysis.py
@@ -206,6 +206,32 @@ def _analyse_surface_tension(config: SimulationConfig, out_dir: Path) -> None:
     console.print(f"[bold green]Calibration outputs saved to:[/bold green] {surface_tension_dir(out_dir)}")
 
 
+def _analyse_coexistence(config: SimulationConfig, out_dir: Path) -> None:
+    """Write the Maxwell-construction artefacts for *config* and report the coexistence pair."""
+    if not is_multiphase(config):
+        msg = "coexistence requires a multiphase configuration; got no [multiphase] section"
+        raise ValueError(msg)
+
+    from src.config.jax_config import configure_jax
+
+    configure_jax()
+
+    from src.simulation_io.analysis.coexistence import calibrate_coexistence
+    from src.simulation_io.analysis.coexistence import coexistence_dir
+
+    pair = calibrate_coexistence(config, out_dir)
+    if pair is None:
+        eos = config.multiphase["eos"] if config.multiphase else None
+        console.print(
+            f"[yellow]EOS {eos!r} has no coexistence solver: its rho_l/rho_v are its own parameters.[/yellow]"
+        )
+        return
+
+    rho_v, rho_l = pair
+    success(f"Maxwell coexistence: rho_v = {rho_v:.6g}, rho_l = {rho_l:.6g}")
+    console.print(f"[bold green]Coexistence outputs saved to:[/bold green] {coexistence_dir(out_dir)}")
+
+
 def _analyse_length_scale(config: SimulationConfig, config_toml: str, out_dir: Path) -> None:
     """Plot the region behind the Bo/Oh length scale and refresh the parameter overview."""
     from src.simulation_io.analysis.physical_parameters import write_length_scale_figure
@@ -238,6 +264,13 @@ def _analyse_length_scale(config: SimulationConfig, config_toml: str, out_dir: P
     help="Report the measured lattice surface tension of the configured fluid and write its "
     "Young-Laplace calibration figures (a fluid not yet measured prints how to stage the droplet sweep).",
 )
+@click.option(
+    "--coexistence",
+    "coexistence",
+    is_flag=True,
+    help="Report the Maxwell coexistence densities of the configured EOS beside the configured ones and "
+    "write the equal-area figure (an EOS without a coexistence solver prints that its densities are parameters).",
+)
 @click.option(
     "--length-scale",
     "length_scale",
@@ -253,7 +286,13 @@ def _analyse_length_scale(config: SimulationConfig, config_toml: str, out_dir: P
     help="Directory for analysis outputs (default: the config file's directory).",
 )
 @cli_command(title="Analysis", interrupt_message="Analysis interrupted by user.")
-def analyse(config_toml: str, surface_tension: bool, length_scale: bool, out_dir: str | None) -> None:
+def analyse(
+    config_toml: str,
+    surface_tension: bool,
+    coexistence: bool,
+    length_scale: bool,
+    out_dir: str | None,
+) -> None:
     """Run standalone analyses for the configuration in CONFIG_TOML.
 
     --surface-tension reads the calibration cache for any multiphase EOS,
@@ -264,14 +303,17 @@ def analyse(config_toml: str, surface_tension: bool, length_scale: bool, out_dir
         # Report the measured surface tension for the configured fluid
         tud-lbm analyse config.toml --surface-tension
 
+        # Compare the configured densities with the EOS's own Maxwell pair
+        tud-lbm analyse config.toml --coexistence
+
         # Check the droplet region the Bond number is built from
         tud-lbm analyse run_dir/config.toml --length-scale
 
         # Write outputs somewhere other than the config's directory
         tud-lbm analyse config.toml --surface-tension --out-dir results/
     """
-    if not (surface_tension or length_scale):
-        msg = "select at least one analysis: --surface-tension or --length-scale"
+    if not (surface_tension or coexistence or length_scale):
+        msg = "select at least one analysis: --surface-tension, --coexistence or --length-scale"
         raise click.UsageError(msg)
 
     target_dir = Path(out_dir) if out_dir is not None else Path(config_toml).resolve().parent
@@ -280,5 +322,7 @@ def analyse(config_toml: str, surface_tension: bool, length_scale: bool, out_dir
     console.print()
     if surface_tension:
         _analyse_surface_tension(config, target_dir)
+    if coexistence:
+        _analyse_coexistence(config, target_dir)
     if length_scale:
         _analyse_length_scale(config, config_toml, target_dir)
diff --git a/src/cli/execution.py b/src/cli/execution.py
index 9c43868..cd3fe2b 100644
--- a/src/cli/execution.py
+++ b/src/cli/execution.py
@@ -62,9 +62,11 @@ def _run_simulation(config: SimulationConfig) -> str:
         simulation_name=config.simulation_name,
     )
 
+    from src.simulation_io.analysis.coexistence import calibrate_coexistence
     from src.simulation_io.analysis.surface_tension import record_surface_tension
 
     config = record_surface_tension(config, io.run_dir)
+    calibrate_coexistence(config, io.run_dir)
 
     console.print("[bold green]Starting simulation...[/bold green]")
     console.print(f"[dim]Results directory: {io.run_dir}[/dim]")
diff --git a/src/config/multiphase.py b/src/config/multiphase.py
index 3870b68..88950e2 100644
--- a/src/config/multiphase.py
+++ b/src/config/multiphase.py
@@ -3,6 +3,12 @@
 Holds the section's schema, the per-EOS parameter table, the JIT-static
 :class:`MultiphaseParams` and the densities the run's phases actually sit at.
 Adding an EOS with its own parameters is one row of :data:`EOS_PARAMETERS`.
+
+``rho_l``/``rho_v`` may be omitted for an EOS that registers a coexistence
+solver (kind ``"coexistence"``, the Maxwell construction):
+:func:`complete_multiphase` fills them from the EOS's own parameters. A
+configured density always wins. To sweep ``t_eos``, omit them, since a
+configured density does not follow the temperature.
 """
 
 from __future__ import annotations
@@ -16,7 +22,8 @@
 if TYPE_CHECKING:
     from src.config.simulation_config import SimulationConfig
 
-#: Every EOS needs these.
+#: Every EOS needs these. ``rho_l``/``rho_v`` are filled by :func:`complete_multiphase`
+#: when omitted for an EOS with a registered coexistence solver.
 REQUIRED: tuple[str, ...] = ("kappa", "rho_l", "rho_v", "interface_width", "eos")
 #: Keys any EOS may carry.
 OPTIONAL: tuple[str, ...] = ("g",)
@@ -51,8 +58,39 @@ def is_multiphase(config: SimulationConfig) -> bool:
 
 
 def complete_multiphase(section: dict[str, Any] | None) -> dict[str, Any] | None:
-    """*section* without its null keys: a null means the key is not set, and TOML cannot write one."""
-    return None if section is None else {name: value for name, value in section.items() if value is not None}
+    """*section* without its null keys, and with omitted densities filled from the Maxwell construction.
+
+    A null means the key is not set, and TOML cannot write one. When ``rho_l``
+    or ``rho_v`` is missing and the EOS registers a coexistence solver with
+    every parameter of its :data:`EOS_PARAMETERS` row present, the missing
+    density is the Maxwell value. A configured density is never overwritten;
+    anything else is left to :func:`validate_multiphase`.
+
+    Raises:
+        ValueError: If the EOS parameters admit no liquid-vapour coexistence.
+    """
+    if section is None:
+        return None
+    completed = {name: value for name, value in section.items() if value is not None}
+    eos = completed.get("eos")
+    if not isinstance(eos, str) or {"rho_l", "rho_v"} <= completed.keys():
+        return completed
+    row = EOS_PARAMETERS.get(eos)
+    if row is None or not set(row) <= completed.keys():
+        return completed
+
+    # Imported here: the EOS package imports this module.
+    from src.operators.macroscopic.eos import coexistence_densities
+
+    pair = coexistence_densities(eos, completed)
+    if pair is not None:
+        rho_v, rho_l = pair
+        # Six significant figures: the saved config.toml stays readable, and the
+        # surface-tension cache key (which holds rho_l/rho_v) cannot differ in the
+        # last bits between machines.
+        completed.setdefault("rho_v", float(f"{rho_v:.6g}"))
+        completed.setdefault("rho_l", float(f"{rho_l:.6g}"))
+    return completed
 
 
 def validate_multiphase(config: SimulationConfig) -> None:
diff --git a/src/operators/macroscopic/eos/__init__.py b/src/operators/macroscopic/eos/__init__.py
index d64f216..6947f45 100644
--- a/src/operators/macroscopic/eos/__init__.py
+++ b/src/operators/macroscopic/eos/__init__.py
@@ -1,15 +1,17 @@
 """Equation-of-state operators for multiphase macroscopic computation.
 
 Public API: build_eos_fn(), build_pressure_fn(), analytical_surface_tension(),
-has_analytical_surface_tension()
+has_analytical_surface_tension(), coexistence_densities()
 
 Implementation modules (_double_well.py, _carnahan_starling.py) are internal;
 use the factories to access them. Each module registers its EOS under the
-``"eos"`` kind and, if it has them, its bulk pressure under ``"pressure"`` and
-its closed-form surface tension under ``"surface_tension"``, using the same
-name in each. To ask whether an EOS has one of those, query the registry —
+``"eos"`` kind and, if it has them, its bulk pressure under ``"pressure"``,
+its closed-form surface tension under ``"surface_tension"`` and its Maxwell
+coexistence solver under ``"coexistence"``, using the same name in each. To ask
+whether an EOS has one of those, query the registry —
 ``eos in get_operator_names("pressure")``,
-:func:`has_analytical_surface_tension` — rather than any list maintained here.
+:func:`has_analytical_surface_tension`, :func:`coexistence_densities` returning
+``None`` — rather than any list maintained here.
 
 Example:
     from src.operators.macroscopic.eos import build_eos_fn, build_pressure_fn
@@ -19,16 +21,21 @@
 
     pressure_fn = build_pressure_fn(mp)
     p_0 = pressure_fn(rho)
+
+    pair = coexistence_densities(section["eos"], section)  # (rho_v, rho_l), or None for double-well
 """
 
 from __future__ import annotations
 from typing import TYPE_CHECKING
+from typing import Any
+from src.config.multiphase import EOS_PARAMETERS
 from src.config.multiphase import multiphase_params
 from src.operators._loader import auto_load_operators
 from src.registry import get_operator_names
 from src.registry import get_operators
 
 if TYPE_CHECKING:
+    from collections.abc import Mapping
     from src.config.multiphase import MultiphaseParams
     from src.protocols import EosOperator
 
@@ -101,10 +108,30 @@ def analytical_surface_tension(mp: MultiphaseParams) -> float | None:
     return None if sigma is None else float(sigma)
 
 
+def coexistence_densities(eos: str, parameters: Mapping[str, Any]) -> tuple[float, float] | None:
+    """Maxwell coexistence densities ``(rho_v, rho_l)`` of *eos*, or ``None``.
+
+    ``None`` when *eos* registers no coexistence solver: it is then
+    parameterised *by* its densities (the double-well's wells), so there is
+    nothing to compute. *parameters* is any mapping holding the EOS's row of
+    :data:`~src.config.multiphase.EOS_PARAMETERS` — the raw ``[multiphase]``
+    section or ``mp._asdict()`` alike; other keys are ignored.
+
+    Raises:
+        ValueError: If the parameters admit no liquid-vapour coexistence
+            (at or above the critical temperature).
+    """
+    if eos not in get_operator_names("coexistence"):
+        return None
+    solver = get_operators("coexistence")[eos].target
+    return solver(**{name: parameters[name] for name in EOS_PARAMETERS[eos]})
+
+
 __all__ = [
     "analytical_surface_tension",
     "build_eos_fn",
     "build_pressure_fn",
+    "coexistence_densities",
     "has_analytical_surface_tension",
     "multiphase_params",
 ]
diff --git a/src/operators/macroscopic/eos/_carnahan_starling.py b/src/operators/macroscopic/eos/_carnahan_starling.py
index 8c6217b..80179c7 100644
--- a/src/operators/macroscopic/eos/_carnahan_starling.py
+++ b/src/operators/macroscopic/eos/_carnahan_starling.py
@@ -1,15 +1,25 @@
-"""Carnahan-Starling EOS: bulk chemical potential and bulk pressure."""
+"""Carnahan-Starling EOS: bulk chemical potential, bulk pressure and coexistence densities."""
 
 from __future__ import annotations
+import math
 from typing import TYPE_CHECKING
+import jax
 import jax.numpy as jnp
+from src.registry import coexistence_operator
 from src.registry import eos_operator
 from src.registry import pressure_operator
 
 if TYPE_CHECKING:
+    from collections.abc import Callable
     from src.config.multiphase import MultiphaseParams
     from src.protocols import EosOperator
 
+#: Halvings per bisection. Far below float64 resolution on every bracket used
+#: here, and a fixed count keeps each loop a ``fori_loop`` with a static trip count.
+_N_BISECT = 200
+#: Samples of ``dp_0/drho`` across ``(0, 4/b)`` that bracket the two spinodals.
+_N_SCAN = 4096
+
 
 def _cs_params(mp: MultiphaseParams) -> tuple[float, float, float, float]:
     """Return the validated ``(a, b, r, t)`` scalars.
@@ -53,6 +63,99 @@ def _pressure_carnahan_starling(
     return ideal - a * rho**2
 
 
+def _bisect(f: Callable[[jnp.ndarray], jnp.ndarray], lo: jnp.ndarray, hi: jnp.ndarray) -> jnp.ndarray:
+    """Root of *f* on ``[lo, hi]``, across which *f* changes sign, after ``_N_BISECT`` halvings."""
+
+    def halve(
+        _: int,
+        bracket: tuple[jnp.ndarray, jnp.ndarray, jnp.ndarray],
+    ) -> tuple[jnp.ndarray, jnp.ndarray, jnp.ndarray]:
+        low, high, f_low = bracket
+        mid = 0.5 * (low + high)
+        f_mid = f(mid)
+        keep_low = jnp.sign(f_mid) != jnp.sign(f_low)
+        return jnp.where(keep_low, low, mid), jnp.where(keep_low, mid, high), jnp.where(keep_low, f_low, f_mid)
+
+    low, high, _ = jax.lax.fori_loop(0, _N_BISECT, halve, (lo, hi, f(lo)))
+    return 0.5 * (low + high)
+
+
+@jax.jit
+def _maxwell_construction(
+    a: float,
+    b: float,
+    r: float,
+    t: float,
+) -> tuple[jnp.ndarray, jnp.ndarray, jnp.ndarray, jnp.ndarray]:
+    """``(rho_v, rho_l, rho_s1, rho_s2)``: the coexistence pair and the two spinodals.
+
+    Coexistence is equal pressure and equal chemical potential. For a trial
+    pressure ``p*`` between the spinodal pressures, each stable branch has one
+    density at ``p*``; ``p_sat`` is the ``p*`` at which their chemical
+    potentials agree. That residual falls strictly with ``p*`` (Gibbs-Duhem:
+    its derivative is ``1/rho_l - 1/rho_v < 0``), so the root is unique. Equal
+    chemical potential is the Maxwell equal-area rule, since
+    ``∫(p - p*) dv`` from ``v_l`` to ``v_v`` equals ``mu_l - mu_v``.
+
+    The lower pressure bracket is clamped positive: on a low isotherm the
+    liquid spinodal pressure is negative, where the vapour branch has no
+    density. Above the critical temperature ``dp_0/drho`` never changes sign,
+    both spinodal indices fall to 0, and ``rho_s1 == rho_s2`` signals it.
+    """
+
+    def pressure(rho: jnp.ndarray) -> jnp.ndarray:
+        return _pressure_carnahan_starling(rho, a, b, r, t)
+
+    def potential(rho: jnp.ndarray) -> jnp.ndarray:
+        return _eos_carnahan_starling(rho, a, b, r, t)
+
+    slope = jax.grad(pressure)
+    pole = jnp.asarray(4.0 / b)
+
+    # Spinodals: the first fall of dp/drho, then the first rise after it.
+    rho = jnp.linspace(1e-6 * pole, (1.0 - 1e-6) * pole, _N_SCAN)
+    rising = jax.vmap(slope)(rho) > 0.0
+    falls = rising[:-1] & ~rising[1:]
+    i_vapour = jnp.argmax(falls)
+    rises = ~rising[:-1] & rising[1:] & (jnp.arange(_N_SCAN - 1) > i_vapour)
+    i_liquid = jnp.argmax(rises)
+    rho_s1 = _bisect(slope, rho[i_vapour], rho[i_vapour + 1])
+    rho_s2 = _bisect(slope, rho[i_liquid], rho[i_liquid + 1])
+
+    p_hi = pressure(rho_s1)
+    p_lo = jnp.maximum(pressure(rho_s2), 1e-12 * p_hi)
+
+    def vapour(p: jnp.ndarray) -> jnp.ndarray:
+        return _bisect(lambda x: pressure(x) - p, 1e-14 * rho_s1, rho_s1)
+
+    def liquid(p: jnp.ndarray) -> jnp.ndarray:
+        return _bisect(lambda x: pressure(x) - p, rho_s2, (1.0 - 1e-12) * pole)
+
+    p_sat = _bisect(lambda p: potential(liquid(p)) - potential(vapour(p)), p_lo, p_hi)
+    return vapour(p_sat), liquid(p_sat), rho_s1, rho_s2
+
+
+@coexistence_operator(name="carnahan-starling")
+def carnahan_starling_coexistence(*, a_eos: float, b_eos: float, r_eos: float, t_eos: float) -> tuple[float, float]:
+    """Maxwell coexistence densities ``(rho_v, rho_l)`` of the Carnahan-Starling EOS.
+
+    The keywords are the EOS's row of ``EOS_PARAMETERS``, so a caller passes
+    exactly that row. Precision relies on x64, which ``src`` enables at import.
+
+    Raises:
+        ValueError: If ``t_eos`` is at or above the critical temperature, where
+            there is no liquid-vapour coexistence.
+    """
+    rho_v, rho_l, rho_s1, rho_s2 = (float(x) for x in _maxwell_construction(a_eos, b_eos, r_eos, t_eos))
+    if not (math.isfinite(rho_s1) and math.isfinite(rho_s2) and rho_s1 < rho_s2):
+        msg = (
+            f"t_eos = {t_eos} is at or above the critical temperature of the Carnahan-Starling EOS "
+            f"(a_eos = {a_eos}, b_eos = {b_eos}): no liquid-vapour coexistence"
+        )
+        raise ValueError(msg)
+    return rho_v, rho_l
+
+
 @eos_operator(name="carnahan-starling")
 def build_carnahan_starling_eos(mp: MultiphaseParams) -> EosOperator:
     """Return ``eos_fn(rho)`` for the Carnahan-Starling EOS using bound params."""
diff --git a/src/pipeline/parallel_runner.py b/src/pipeline/parallel_runner.py
index b82a62d..542e937 100644
--- a/src/pipeline/parallel_runner.py
+++ b/src/pipeline/parallel_runner.py
@@ -129,9 +129,11 @@ def run_single_simulation(
             output_format=config.output_format,
         )
 
+        from src.simulation_io.analysis.coexistence import calibrate_coexistence
         from src.simulation_io.analysis.surface_tension import record_surface_tension
 
         config = record_surface_tension(config, io.run_dir)
+        calibrate_coexistence(config, io.run_dir)
 
         # Run simulation with streaming I/O
         _, _ = run_fn(
diff --git a/src/registry.py b/src/registry.py
index 0938a7c..667147e 100644
--- a/src/registry.py
+++ b/src/registry.py
@@ -254,6 +254,20 @@ def surface_tension_operator(
     return register_operator("surface_tension", name=name, **meta)
 
 
+def coexistence_operator(
+    *,
+    name: str | None = None,
+    **meta: object,
+) -> Callable[[_OT], _OT]:
+    """Register a liquid-vapour coexistence solver (kind ``"coexistence"``).
+
+    Names match the ``"eos"`` kind: an EOS registers here when its coexistence
+    densities ``(rho_v, rho_l)`` follow from the Maxwell
+    construction. Absence means the EOS depends on the supplied density.
+    """
+    return register_operator("coexistence", name=name, **meta)
+
+
 def initialise_operator(
     *,
     name: str | None = None,
diff --git a/src/simulation_io/analysis/coexistence.py b/src/simulation_io/analysis/coexistence.py
new file mode 100644
index 0000000..8eaf888
--- /dev/null
+++ b/src/simulation_io/analysis/coexistence.py
@@ -0,0 +1,147 @@
+"""Maxwell-construction calibration: the coexistence densities an EOS predicts for itself.
+
+For an EOS that registers a coexistence solver (kind ``"coexistence"``,
+Carnahan-Starling), the liquid and vapour densities at coexistence follow from
+its own parameters. :func:`calibrate_coexistence` computes them and writes them
+into the run directory beside the configured densities, so a gap between the two
+is visible. Unlike the Young-Laplace calibration of
+:mod:`~src.simulation_io.analysis.surface_tension` it needs no sweep and no
+cache: it is a scalar root-find that takes milliseconds.
+
+The artefacts sit under ``<run_dir>/coexistence/``, in the ``data/`` + ``plots/``
+shape a run directory itself has:
+
+``data/data.json``
+    ``eos``, the Maxwell ``rho_v``/``rho_l``, ``p_sat`` and the configured densities.
+``plots/maxwell.png``
+    The bulk pressure against specific volume, the tie line at ``p_sat`` and the
+    two equal areas, with the configured densities marked.
+"""
+
+from __future__ import annotations
+import json
+from pathlib import Path
+from typing import TYPE_CHECKING
+import jax.numpy as jnp
+import numpy as np
+from src.config.multiphase import multiphase_params
+from src.config.run_config import DATA_DIRNAME
+from src.config.run_config import PLOTS_DIRNAME
+from src.operators.macroscopic.eos import build_pressure_fn
+from src.operators.macroscopic.eos import coexistence_densities
+
+if TYPE_CHECKING:
+    from collections.abc import Callable
+    from src.config.simulation_config import SimulationConfig
+
+_OUTPUT_DIRNAME = "coexistence"
+_DATA_FILENAME = "data.json"
+_PLOT_FILENAME = "maxwell.png"
+#: Samples of the isotherm drawn in the figure.
+_N_CURVE = 2000
+
+
+def coexistence_dir(run_dir: str | Path) -> Path:
+    """Return the subdirectory of *run_dir* holding the coexistence artefacts."""
+    return Path(run_dir) / _OUTPUT_DIRNAME
+
+
+def calibrate_coexistence(config: SimulationConfig, run_dir: str | Path) -> tuple[float, float] | None:
+    """Compute the Maxwell pair of *config*'s EOS and write its artefacts into *run_dir*.
+
+    Returns:
+        ``(rho_v, rho_l)``, or ``None`` for a single-phase run or an EOS that
+        registers no coexistence solver (its densities are its own parameters).
+
+    Raises:
+        ValueError: If the EOS parameters admit no liquid-vapour coexistence.
+    """
+    params = multiphase_params(config)
+    if params is None:
+        return None
+    pair = coexistence_densities(params.eos, params._asdict())
+    if pair is None:
+        return None
+    rho_v, rho_l = pair
+    pressure_fn = build_pressure_fn(params)
+
+    def pressure(rho: np.ndarray) -> np.ndarray:
+        return np.asarray(pressure_fn(jnp.asarray(rho)))
+
+    p_sat = float(pressure(np.asarray(rho_v)))
+    out = coexistence_dir(run_dir)
+    _save_data(
+        out / DATA_DIRNAME / _DATA_FILENAME,
+        {
+            "eos": params.eos,
+            "rho_v": rho_v,
+            "rho_l": rho_l,
+            "p_sat": p_sat,
+            "rho_v_configured": float(params.rho_v),
+            "rho_l_configured": float(params.rho_l),
+        },
+    )
+    _save_plot(
+        out / PLOTS_DIRNAME / _PLOT_FILENAME,
+        pressure,
+        maxwell=(rho_v, rho_l),
+        configured=(float(params.rho_v), float(params.rho_l)),
+        p_sat=p_sat,
+    )
+    return rho_v, rho_l
+
+
+def _save_data(path: Path, payload: dict[str, str | float]) -> None:
+    """Write the coexistence values into the run directory."""
+    path.parent.mkdir(parents=True, exist_ok=True)
+    path.write_text(json.dumps(payload, indent=2), encoding="utf-8")
+
+
+def _save_plot(
+    path: Path,
+    pressure: Callable[[np.ndarray], np.ndarray],
+    *,
+    maxwell: tuple[float, float],
+    configured: tuple[float, float],
+    p_sat: float,
+) -> None:
+    """Draw ``p_0(1/v)`` with the tie line at *p_sat* and the two equal areas shaded."""
+    import matplotlib as mpl
+
+    mpl.use("Agg")
+    import matplotlib.pyplot as plt
+
+    rho_v, rho_l = maxwell
+    v_liquid, v_vapour = 1.0 / rho_l, 1.0 / rho_v
+    volume = np.geomspace(0.8 * v_liquid, 2.0 * v_vapour, _N_CURVE)
+    p = pressure(1.0 / volume)
+    inside = (volume >= v_liquid) & (volume <= v_vapour)
+    p_inside = p[inside]
+    span = float(p_inside.max() - p_inside.min())
+
+    fig, ax = plt.subplots(figsize=(7, 5))
+    ax.plot(volume, p, color="black", lw=1.5, label="$p_0(1/\\rho)$")
+    ax.hlines(p_sat, v_liquid, v_vapour, colors="tab:red", lw=1.5, label=f"$p_{{sat}}$ = {p_sat:.6g}")
+    # The isotherm crosses p_sat only at v_l, once between, and at v_v, so each area is one contiguous slice.
+    below, above = inside & (p < p_sat), inside & (p > p_sat)
+    ax.fill_between(volume[below], p[below], p_sat, color="tab:blue", alpha=0.3, label="liquid-side area")
+    ax.fill_between(volume[above], p[above], p_sat, color="tab:orange", alpha=0.3, label="vapour-side area")
+    for rho, label, colour in (
+        (rho_v, f"Maxwell $\\rho_v$ = {rho_v:.6g}", "tab:red"),
+        (rho_l, f"Maxwell $\\rho_l$ = {rho_l:.6g}", "tab:red"),
+        (configured[0], f"configured $\\rho_v$ = {configured[0]:.6g}", "tab:green"),
+        (configured[1], f"configured $\\rho_l$ = {configured[1]:.6g}", "tab:green"),
+    ):
+        ax.axvline(1.0 / rho, color=colour, ls=":", lw=1.0, label=label)
+    ax.set_xscale("log")
+    ax.set_ylim(p_inside.min() - 0.1 * span, p_inside.max() + 0.1 * span)
+    ax.set_xlabel("specific volume $1/\\rho$ [lattice units]")
+    ax.set_ylabel("$p_0$ [lattice units]")
+    ax.set_title("Maxwell construction: areas equal in $v$ (the log axis distorts them)")
+    ax.legend(fontsize="small")
+    ax.grid(visible=True, alpha=0.3)
+
+    fig.tight_layout()
+    path.parent.mkdir(parents=True, exist_ok=True)
+    fig.savefig(path, dpi=150, bbox_inches="tight")
+    plt.close(fig)
diff --git a/src/simulation_io/analysis/physical_parameters/physical_parameters.py b/src/simulation_io/analysis/physical_parameters/physical_parameters.py
index 293de1b..410a854 100644
--- a/src/simulation_io/analysis/physical_parameters/physical_parameters.py
+++ b/src/simulation_io/analysis/physical_parameters/physical_parameters.py
@@ -28,6 +28,7 @@
 from src.config.multiphase import multiphase_params
 from src.operators.macroscopic.eos import analytical_surface_tension
 from src.operators.macroscopic.eos import build_pressure_fn
+from src.operators.macroscopic.eos import coexistence_densities
 from src.operators.macroscopic.eos import has_analytical_surface_tension
 from src.registry import get_operator_names
 from src.registry import get_operators
@@ -670,9 +671,13 @@ def _add_multiphase_section(lines: list[str], config: SimulationConfig) -> None:
             _row("kappa:", params.kappa),
             _row("rho_liquid:", params.rho_l),
             _row("rho_vapour:", params.rho_v),
-            _row("Interface width:", params.interface_width),
         ]
     )
+    # The pair the EOS predicts for itself, beside the configured one, so a gap shows.
+    if (pair := coexistence_densities(params.eos, params._asdict())) is not None:
+        rho_v, rho_l = pair
+        lines.append(_row("Maxwell coexistence:", f"rho_v = {rho_v:.6g}, rho_l = {rho_l:.6g}"))
+    lines.append(_row("Interface width:", params.interface_width))
     _add_measured_density_rows(lines, config)
     if params.g is not None:
         lines.append(_row("g (gravity):", params.g))
diff --git a/tests/cli/test_calibration_commands.py b/tests/cli/test_calibration_commands.py
index d1561a4..59ed4ad 100644
--- a/tests/cli/test_calibration_commands.py
+++ b/tests/cli/test_calibration_commands.py
@@ -29,10 +29,10 @@ def _config(**overrides: object) -> SimulationConfig:
             "rho_l": 0.4,
             "rho_v": 0.02,
             "interface_width": 4,
-            "a_eos": 0.5,
-            "b_eos": 4.0,
+            "a_eos": 0.00031459670905604266,
+            "b_eos": 0.1490857142857143,
             "r_eos": 1.0,
-            "t_eos": 0.05,
+            "t_eos": 0.00039808421247983624,
         },
         "simulation_name": "drop",
         "grid_shape": (48, 32),
diff --git a/tests/cli/test_cli.py b/tests/cli/test_cli.py
index 89e03d9..cfd3b89 100644
--- a/tests/cli/test_cli.py
+++ b/tests/cli/test_cli.py
@@ -2714,6 +2714,42 @@ def test_analyse_reports_calibrated_surface_tension(self, tmp_path):
         assert result.exit_code == 0
         assert "0.0123" in result.output
 
+    def test_coexistence_reports_the_maxwell_pair(self, tmp_path):
+        cfg_toml = self._write_multiphase_toml(tmp_path)
+        content = (tmp_path / "config.toml").read_text(encoding="utf-8")
+        cs_section = (
+            'eos = "carnahan-starling"\n'
+            "a_eos = 0.00031459670905604266\n"
+            "b_eos = 0.1490857142857143\n"
+            "r_eos = 1.0\n"
+            "t_eos = 0.00039808421247983624\n"
+        )
+        content = content.replace("rho_l = 1.0\nrho_v = 0.001\n", "").replace('eos = "double-well"\n', cs_section)
+        (tmp_path / "config.toml").write_text(content, encoding="utf-8")
+
+        result = CliRunner().invoke(cli, ["analyse", cfg_toml, "--coexistence"])
+
+        assert result.exit_code == 0
+        assert "rho_v = 0.0168045, rho_l = 12.1833" in result.output
+        assert (tmp_path / "coexistence" / "data" / "data.json").exists()
+
+    def test_coexistence_of_a_double_well_names_its_parameters(self, tmp_path):
+        cfg_toml = self._write_multiphase_toml(tmp_path)
+
+        result = CliRunner().invoke(cli, ["analyse", cfg_toml, "--coexistence"])
+
+        assert result.exit_code == 0
+        assert "own" in result.output
+        assert "parameters" in result.output
+        assert not (tmp_path / "coexistence").exists()
+
+    def test_coexistence_of_a_single_phase_config_errors(self, tmp_path):
+        cfg_toml = tmp_path / "config.toml"
+        cfg_toml.write_text("[simulation_type]\ntau=0.8\nnt=10\ngrid_shape=[8, 8]\n", encoding="utf-8")
+        result = CliRunner().invoke(cli, ["analyse", str(cfg_toml), "--coexistence"])
+        assert result.exit_code == 1
+        assert "multiphase" in result.output
+
     def test_length_scale_writes_figure_and_refreshes_overview(self, tmp_path):
         cfg_toml = self._write_multiphase_toml(tmp_path)
 
diff --git a/tests/config/test_multiphase_coexistence.py b/tests/config/test_multiphase_coexistence.py
new file mode 100644
index 0000000..61e8b2e
--- /dev/null
+++ b/tests/config/test_multiphase_coexistence.py
@@ -0,0 +1,54 @@
+"""An omitted ``rho_l``/``rho_v`` is filled with the Maxwell pair of an EOS that has a coexistence solver."""
+
+from __future__ import annotations
+import pytest
+from src.config.multiphase import multiphase_params
+from src.simulation_io.readers import DictAdapter
+
+_CS = {
+    "eos": "carnahan-starling",
+    "kappa": 0.015,
+    "interface_width": 4,
+    "a_eos": 0.00031459670905604266,
+    "b_eos": 0.1490857142857143,
+    "r_eos": 1.0,
+    "t_eos": 0.00039808421247983624,
+}
+#: The Maxwell pair of _CS, rounded as the config rounds it.
+_RHO_V, _RHO_L = 0.0168045, 12.1833
+
+
+def _densities(section: dict) -> tuple[float, float]:
+    config = DictAdapter().load({"grid_shape": (16, 16), "nt": 10, "multiphase": section})
+    params = multiphase_params(config)
+    assert params is not None
+    return params.rho_v, params.rho_l
+
+
+def test_omitted_densities_are_the_maxwell_pair():
+    assert _densities(dict(_CS)) == (_RHO_V, _RHO_L)
+
+
+def test_configured_densities_are_kept():
+    assert _densities(_CS | {"rho_v": 0.015, "rho_l": 12.18}) == (0.015, 12.18)
+
+
+def test_only_the_missing_density_is_filled():
+    assert _densities(_CS | {"rho_v": 0.015}) == (0.015, _RHO_L)
+
+
+def test_double_well_still_requires_both_densities():
+    adapter = DictAdapter()
+    raw = {
+        "grid_shape": (16, 16),
+        "multiphase": {"eos": "double-well", "kappa": 0.015, "rho_v": 0.1, "interface_width": 4},
+    }
+    with pytest.raises(ValueError, match="'rho_l' is required"):
+        adapter.load(raw)
+
+
+def test_supercritical_fluid_without_densities_raises():
+    adapter = DictAdapter()
+    raw = {"grid_shape": (16, 16), "multiphase": _CS | {"t_eos": 1.0}}
+    with pytest.raises(ValueError, match="critical temperature"):
+        adapter.load(raw)
diff --git a/tests/io/test_coexistence_calibration.py b/tests/io/test_coexistence_calibration.py
new file mode 100644
index 0000000..2c3dd5c
--- /dev/null
+++ b/tests/io/test_coexistence_calibration.py
@@ -0,0 +1,52 @@
+"""``calibrate_coexistence`` writes the Maxwell pair into a run directory, for a CS fluid only."""
+
+from __future__ import annotations
+import json
+import pytest
+from src.config.run_config import DATA_DIRNAME
+from src.config.run_config import PLOTS_DIRNAME
+from src.simulation_io.analysis.coexistence import calibrate_coexistence
+from src.simulation_io.analysis.coexistence import coexistence_dir
+from src.simulation_io.readers import DictAdapter
+
+_CS = {
+    "eos": "carnahan-starling",
+    "kappa": 0.015,
+    "rho_l": 12.18,
+    "rho_v": 0.015,
+    "interface_width": 4,
+    "a_eos": 0.00031459670905604266,
+    "b_eos": 0.1490857142857143,
+    "r_eos": 1.0,
+    "t_eos": 0.00039808421247983624,
+}
+
+
+def _config(multiphase: dict | None):
+    raw: dict = {"grid_shape": (16, 16), "nt": 10}
+    if multiphase is not None:
+        raw["multiphase"] = multiphase
+    return DictAdapter().load(raw)
+
+
+def test_carnahan_starling_writes_the_pair_and_the_figure(tmp_path):
+    pair = calibrate_coexistence(_config(_CS), tmp_path)
+
+    assert pair == pytest.approx((0.01680445791, 12.18327187), rel=1e-6)
+    out = coexistence_dir(tmp_path)
+    data = json.loads((out / DATA_DIRNAME / "data.json").read_text(encoding="utf-8"))
+    assert data["eos"] == "carnahan-starling"
+    assert (data["rho_v"], data["rho_l"]) == pair
+    assert data["p_sat"] == pytest.approx(6.61753622e-6, rel=1e-6)
+    assert (data["rho_v_configured"], data["rho_l_configured"]) == (0.015, 12.18)
+    assert (out / PLOTS_DIRNAME / "maxwell.png").stat().st_size > 0
+
+
+@pytest.mark.parametrize(
+    "multiphase",
+    [None, {"eos": "double-well", "kappa": 0.015, "rho_l": 1.0, "rho_v": 0.1, "interface_width": 4}],
+    ids=["single-phase", "double-well"],
+)
+def test_no_solver_returns_none_and_writes_nothing(tmp_path, multiphase):
+    assert calibrate_coexistence(_config(multiphase), tmp_path) is None
+    assert not coexistence_dir(tmp_path).exists()
diff --git a/tests/io/test_physical_parameters.py b/tests/io/test_physical_parameters.py
index 694bbaa..8c1964e 100644
--- a/tests/io/test_physical_parameters.py
+++ b/tests/io/test_physical_parameters.py
@@ -279,3 +279,26 @@ def test_inclusion_mask_selects_the_bubble_not_its_ambient():
 def test_inclusion_mask_returns_none_for_a_single_phase_field():
     rho = np.full((40, 20, 1, 1, 1), 0.5)
     assert inclusion_mask_from_rho(rho, rho_mean=0.75) is None
+
+
+def test_overview_reports_the_maxwell_pair_for_carnahan_starling():
+    cfg = _mp_config(
+        multiphase={
+            "eos": "carnahan-starling",
+            "rho_l": 12.18,
+            "rho_v": 0.015,
+            "a_eos": 0.00031459670905604266,
+            "b_eos": 0.1490857142857143,
+            "r_eos": 1.0,
+            "t_eos": 0.00039808421247983624,
+        },
+    )
+
+    text = build_overview(cfg)
+
+    assert "Maxwell coexistence:" in text
+    assert "rho_v = 0.0168045, rho_l = 12.1833" in text
+
+
+def test_overview_has_no_maxwell_row_for_double_well():
+    assert "Maxwell coexistence:" not in build_overview(_mp_config())
diff --git a/tests/io/test_surface_tension_calibration.py b/tests/io/test_surface_tension_calibration.py
index b549cdd..b541da5 100644
--- a/tests/io/test_surface_tension_calibration.py
+++ b/tests/io/test_surface_tension_calibration.py
@@ -37,10 +37,10 @@ def _stub_config(**overrides):
             "rho_l": 0.4,
             "rho_v": 0.02,
             "interface_width": 4,
-            "a_eos": 0.5,
-            "b_eos": 4.0,
+            "a_eos": 0.00031459670905604266,
+            "b_eos": 0.1490857142857143,
             "r_eos": 1.0,
-            "t_eos": 0.05,
+            "t_eos": 0.00039808421247983624,
         },
         "grid_shape": (64, 64, 1),
     }
@@ -643,10 +643,10 @@ def _cs_config(**overrides):
             "rho_l": 0.4,
             "rho_v": 0.02,
             "interface_width": 4,
-            "a_eos": 0.5,
-            "b_eos": 4.0,
+            "a_eos": 0.00031459670905604266,
+            "b_eos": 0.1490857142857143,
             "r_eos": 1.0,
-            "t_eos": 0.05,
+            "t_eos": 0.00039808421247983624,
         },
         "grid_shape": (32, 32),
         "tau": 0.99,
diff --git a/tests/io/test_surface_tension_policy.py b/tests/io/test_surface_tension_policy.py
index 9232bee..023baf4 100644
--- a/tests/io/test_surface_tension_policy.py
+++ b/tests/io/test_surface_tension_policy.py
@@ -99,10 +99,10 @@ def test_uncalibrated_csv_reports_analytical_source(tmp_path: Path):
         "rho_l": 0.4,
         "rho_v": 0.02,
         "interface_width": 4,
-        "a_eos": 0.5,
-        "b_eos": 4.0,
+        "a_eos": 0.00031459670905604266,
+        "b_eos": 0.1490857142857143,
         "r_eos": 1.0,
-        "t_eos": 0.05,
+        "t_eos": 0.00039808421247983624,
     }
 }
 
diff --git a/tests/operators/macroscopic/test_eos.py b/tests/operators/macroscopic/test_eos.py
index b565e12..5b4b6ab 100644
--- a/tests/operators/macroscopic/test_eos.py
+++ b/tests/operators/macroscopic/test_eos.py
@@ -14,6 +14,7 @@
        - ``compute_macroscopic_multiphase`` with CS EOS: finite outputs
        - Uniform density field → zero interaction force
        - ``build_setup`` with full CS config succeeds
+    4. Coexistence (Maxwell construction) for CS; none for double-well
 """
 
 from __future__ import annotations
@@ -25,7 +26,8 @@
 
 # ---------------------------------------------------------------------------
 # Shared CS EOS parameters (from examples/config_cs_simple.toml)
-# Subcritical temperature; coexistence at rho_v≈0.015, rho_l≈12.18
+# Subcritical temperature (T_r = 0.5); Maxwell coexistence at rho_v≈0.0168, rho_l≈12.18
+# (the configured rho_v = 0.015 sits below it)
 # Singularity at bρ = 4  →  ρ_sing ≈ 26.8 — well above rho_l
 # ---------------------------------------------------------------------------
 
@@ -545,3 +547,69 @@ def test_surface_tension_kind_names_match_the_eos_that_implement_one(self):
 
         assert get_operator_names("surface_tension") == {"double-well"}
         assert get_operator_names("surface_tension") <= get_operator_names("eos")
+
+
+# The Julia Maxwell script's preamble: rho_c = 3.5, p_c = 0.001, T_r = 0.65.
+_JULIA_B = 0.5218 / 3.5
+_JULIA_A = _JULIA_B**2 * 0.001 / (0.3773**2 * 0.4963)
+_JULIA_T_CRITICAL = 0.3773 * _JULIA_A / _JULIA_B
+
+
+class TestCoexistence:
+    """The Carnahan-Starling coexistence pair is the Maxwell construction of its own parameters."""
+
+    @staticmethod
+    def _pair(t: float = _T, a: float = _A, b: float = _B) -> tuple[float, float]:
+        from src.operators.macroscopic.eos import coexistence_densities
+
+        pair = coexistence_densities("carnahan-starling", {"a_eos": a, "b_eos": b, "r_eos": _R, "t_eos": t})
+        assert pair is not None
+        return pair
+
+    def test_shipped_fluid_matches_the_reference_pair(self):
+        rho_v, rho_l = self._pair()
+        assert rho_v == pytest.approx(0.01680445791, rel=1e-6)
+        assert rho_l == pytest.approx(12.18327187, rel=1e-6)
+
+    def test_pair_has_equal_pressure_and_chemical_potential(self):
+        from src.operators.macroscopic.eos._carnahan_starling import _eos_carnahan_starling
+        from src.operators.macroscopic.eos._carnahan_starling import _pressure_carnahan_starling
+
+        rho_v, rho_l = self._pair()
+        p_v, p_l = (float(_pressure_carnahan_starling(jnp.asarray(r), _A, _B, _R, _T)) for r in (rho_v, rho_l))
+        mu_v, mu_l = (float(_eos_carnahan_starling(jnp.asarray(r), _A, _B, _R, _T)) for r in (rho_v, rho_l))
+        assert p_l == pytest.approx(p_v, rel=1e-10)
+        assert mu_l == pytest.approx(mu_v, rel=1e-10)
+
+    def test_areas_are_equal_in_the_pressure_volume_plane(self):
+        from scipy.integrate import quad
+        from src.operators.macroscopic.eos._carnahan_starling import _pressure_carnahan_starling
+
+        rho_v, rho_l = self._pair()
+        p_sat = float(_pressure_carnahan_starling(jnp.asarray(rho_v), _A, _B, _R, _T))
+
+        def excess(volume: float) -> float:
+            return float(_pressure_carnahan_starling(jnp.asarray(1.0 / volume), _A, _B, _R, _T)) - p_sat
+
+        area, _ = quad(excess, 1.0 / rho_l, 1.0 / rho_v, limit=200)
+        assert area == pytest.approx(0.0, abs=1e-12)
+
+    def test_julia_defaults_reproduce_the_script(self):
+        rho_v, rho_l = self._pair(t=0.65 * _JULIA_T_CRITICAL, a=_JULIA_A, b=_JULIA_B)
+        assert rho_v == pytest.approx(0.14985544, rel=1e-6)
+        assert rho_l == pytest.approx(10.258465, rel=1e-6)
+
+    def test_supercritical_temperature_raises(self):
+        from src.operators.macroscopic.eos._carnahan_starling import carnahan_starling_coexistence
+
+        t_eos = 1.05 * _JULIA_T_CRITICAL
+        with pytest.raises(ValueError, match="critical temperature"):
+            carnahan_starling_coexistence(a_eos=_JULIA_A, b_eos=_JULIA_B, r_eos=_R, t_eos=t_eos)
+
+    def test_kind_membership_is_the_capability_check(self):
+        from src.operators.macroscopic.eos import coexistence_densities
+        from src.registry import get_operator_names
+
+        assert get_operator_names("coexistence") == {"carnahan-starling"}
+        assert get_operator_names("coexistence") <= get_operator_names("eos")
+        assert coexistence_densities("double-well", {}) is None