diff --git a/benches/LADDER.md b/benches/LADDER.md
new file mode 100644
index 00000000..36701a09
--- /dev/null
+++ b/benches/LADDER.md
@@ -0,0 +1,231 @@
+# Comprehensive Benchmarking
+
+This benchmark separates kernel, shared-memory, process, network, and scaling effects. Define:
+
+- `N`: nodes
+- `C`: usable cores per node
+- `R`: MPI ranks per node
+- `P`: in-process partitions per rank
+- `P_total = N × R × P`: total partitions
+
+Choose `R` for the platform, normally one rank per NUMA domain, then set `P = C / R`. Keep `R`
+and `P` fixed wherever the ladder requires a direct comparison.
+
+> **Example platform.** All flags, term counts, timings, and memory figures below were measured on
+> code based on `main` commit `97f95f762dfc7174243ccd59dd4ecbbb9775b610`, on Deucalion x86
+> nodes with `C=128`, using `R=8` and `P=16` for MPI runs. They are calibration and job-sizing
+> examples, not portable targets but memory should be relatively similar.
+
+## Ladder
+
+| rung | shape | purpose |
+| --- | --- | --- |
+| L1 | `N=1, R=1, P=1` | isolate single-thread kernel performance |
+| L2a | `N=1, R=1, P=C` | measure full-node shared-memory performance |
+| L2b | `N=1, R=R_platform, P=C/R` | isolate process and MPI overhead without a network hop |
+| L3 | `N>1`, same `R` and `P` as L2b | isolate the network |
+| L4 | vary `N`, hold `R` and `P` fixed | measure strong and weak scaling |
+
+L2a, L2b, and L3 use the same model flags. Only the execution shape changes. Term counts are
+topology-independent, so a mismatch means the problem definition changed.
+
+## Running a row
+
+```bash
+uv sync --all-groups --all-extras          # once
+just bench L1-hubbard --hubbard-cutoff=10 --hubbard-lower-atol=4.2e-05 \
+    -k "test_model_propagate and hubbard"
+```
+
+Name `LABEL` by rung, row, and arm, for example `L1-hubbard-branch`. A run writes
+`results/time-<label>.json`, `results/<label>.json`, and rebuilds `benches/results/REPORT.md`.
+Use separate labels for repetitions and leave `--bench-rounds=1`; later rounds can overlap the
+previous round's live propagator and inflate memory.
+
+Run `build_graph`, `energy`, and `gradient` together when all three are needed. They share one
+graph. Run `propagate` separately because it owns another operator.
+
+## Declaring the shape
+
+The launcher sets `N` and `R`; export `P` explicitly:
+
+```bash
+export monoprop_PARTITIONS=<partitions-per-rank>
+export monoprop_NUM_THREADS=<cores-per-rank>
+```
+
+For multi-rank runs this is mandatory. The engine otherwise defaults to one partition per rank.
+Before comparing results, check `ranks`, `nodes`, `ranks_per_node`, `partitions_env`, and
+`monoprop_threads` in the report's Configuration table.
+
+## Calibrating model size
+
+| model | preferred size control | constraint |
+| --- | --- | --- |
+| Hubbard | `--hubbard-lower-atol` | `2 × num_sites <= MAX_NUM_MODES` |
+| Pauli | `--pauli-lower-atol` | the supplied heavy-hex model requires 127 qubits |
+| random, Heisenberg | `--obs-terms` | duplicates make this an upper bound; fix `--seed` |
+| random, Schrödinger | `--num-generators` | also sets gradient length and work |
+
+
+Calibrate each row on the target platform before running the ladder:
+
+1. Fix the model structure and seed.
+2. Sweep the preferred size control until the full row fits the desired memory budget.
+3. Size graph rows against the largest peak of `build_graph`, `energy`, and `gradient`, not the
+   graph build alone.
+4. Reuse the chosen flags unchanged in L2a, L2b, and L3.
+
+In the Schrödinger picture, sweep number of qubits or `--num-generators` instead.
+
+## Example problem set
+
+Examples 128-core platform.
+
+### L1 example
+
+Set `monoprop_PARTITIONS=1` and `monoprop_NUM_THREADS=1`.
+
+| row | flags | selector | terms | time | peak RSS |
+| --- | --- | --- | ---: | ---: | ---: |
+| Hubbard `propagate` | `--hubbard-cutoff=10 --hubbard-lower-atol=4.2e-05` | `test_model_propagate and hubbard` | 9.95M | 30 s | 1 GiB |
+| Pauli `propagate` | `--pauli-cutoff=12 --pauli-lower-atol=1.22e-04` | `test_model_propagate and pauli` | 10.07M | 20 s | 1 GiB |
+| random `gradient` | `--num-generators=1000 --num-modes=142 --cutoff=6 --obs-terms=295000` | `test_random_gradient and heisenberg` | 19.90M | 11 s | 2.5 GiB |
+| random `gradient`, pared | same, plus `--pare-threshold=1e-10` | `test_random_gradient and heisenberg` | 19.90M | 0.8 s | 2.8 GiB |
+
+### L2a example
+
+Set `monoprop_PARTITIONS=C` and `monoprop_NUM_THREADS=C`. On the example platform, `C=128`. These are medians of three repetitions.
+
+
+| row | flags | selector | terms | time | peak RSS |
+| --- | --- | --- | ---: | ---: | ---: |
+| Hubbard `propagate` | `--hubbard-cutoff=10 --hubbard-lower-atol=3.38e-06` | `test_model_propagate and hubbard` | 1.002B | 140 s | 60 GiB |
+| Pauli `propagate` | `--pauli-cutoff=14 --pauli-lower-atol=8.9e-06` | `test_model_propagate and pauli` | 0.986B | 90 s | 80 GiB |
+| random Heisenberg `propagate` | `--num-generators=1000 --num-modes=142 --cutoff=6 --obs-terms=14750000` | `test_random_propagate and heisenberg` | 0.949B | 50 s | 95 GiB |
+| random Heisenberg `build_graph` | same | `test_random_build_graph and heisenberg` | 0.949B | 50 s | 109 GiB |
+| random Heisenberg `energy` | same | `test_random_energy and heisenberg` | 0.949B | 9.9 s | 102 GiB |
+| random Heisenberg `gradient` | same | `test_random_gradient and heisenberg` | 0.949B | 34 s | 110 GiB |
+| random Schrödinger `build_graph` | `--num-generators=5800 --num-modes=142 --cutoff=6 --obs-terms=200000` | `test_random_build_graph and schrodinger` | 0.597B | 66 s | 73 GiB |
+| random Schrödinger `energy` | same | `test_random_energy and schrodinger` | 0.597B | 11 s | 79 GiB |
+| random Schrödinger `gradient` | same | `test_random_gradient and schrodinger` | 0.597B | 39 s | 100 GiB |
+
+
+## L2b: process overhead
+
+Run the L2a problems on one node with `R × P = C`. Use the same `R` and `P` in L3. One rank per
+NUMA domain is a portable starting point, not a guaranteed optimum.
+
+Example results at `N=1, R=8, P=16`:
+
+| row | time | RSS sum per node | worst-rank RSS |
+| --- | ---: | ---: | ---: |
+| Hubbard `propagate` | 169 s | 59.5 GiB | 7.5 GiB |
+| Pauli `propagate` | 102 s | 85.6 GiB | 10.7 GiB |
+| random Heisenberg `propagate` | 48.4 s | 160 GiB | 20.2 GiB |
+| random Heisenberg `build_graph` | 53.2 s | 174 GiB | 21.8 GiB |
+| random Heisenberg `energy` | 10.4 s | 131 GiB | 16.4 GiB |
+| random Heisenberg `gradient` | 30.9 s | 139 GiB | 17.4 GiB |
+| random Schrödinger `build_graph` | 67.3 s | 77.5 GiB | 9.8 GiB |
+| random Schrödinger `energy` | 11.7 s | 81.3 GiB | 10.2 GiB |
+| random Schrödinger `gradient` | 39.5 s | 102 GiB | 12.8 GiB |
+
+The node sum is an upper bound because shared pages are charged to every rank. Use it for
+provisioning, but do not compare it directly with a single process's peak. Per-rank overhead is
+model-dependent, so size L2b from its own measurements.
+
+## L3: network overhead
+
+Run the L2b shape on multiple nodes without changing `R`, `P`, model flags, or selectors. The
+L2b-to-L3 difference then isolates the network.
+
+Example results at `N=4, R=8, P=16`:
+
+| row | time | RSS sum per node | worst-rank RSS |
+| --- | ---: | ---: | ---: |
+| Hubbard `propagate` | 35 s | 18 GiB | — |
+| Pauli `propagate` | 23 s | 25 GiB | — |
+| random Heisenberg `propagate` | 13 s | 100 GiB | — |
+| random Heisenberg `build_graph` | 14 s | 105 GiB | 13.1 GiB |
+| random Heisenberg `energy` | 2.3 s | 61 GiB | 7.7 GiB |
+| random Heisenberg `gradient` | 7.6 s | 63 GiB | 7.9 GiB |
+| random Schrödinger `build_graph` | 20 s | 26 GiB | 3.3 GiB |
+| random Schrödinger `energy` | 1.9 s | 26 GiB | 3.3 GiB |
+| random Schrödinger `gradient` | 8.9 s | 36 GiB | 4.5 GiB |
+
+## L4: strong and weak scaling
+
+Use one representative operation and hold `R` and `P` fixed.
+
+- **Strong scaling:** hold the problem fixed and vary `N`.
+- **Weak scaling:** hold work per node approximately fixed while varying `N`.
+
+For strong scaling relative to the smallest run `(N0, t0)`, efficiency is
+`N0 × t0 / (N × t(N))`. For weak scaling relative to one node, it is `t(1) / t(N)`.
+
+The example uses Hubbard `propagate` with `--hubbard-cutoff=10`. This platform-specific sweep
+approximately doubles the term count at each step:
+
+| k | `--hubbard-lower-atol` | terms |
+| ---: | --- | ---: |
+| 0 | `1.25e-05` | 96,981,051 |
+| 1 | `8.8e-06` | 184,124,520 |
+| 2 | `5.9e-06` | 377,482,074 |
+| 3 | `3.9e-06` | 781,669,404 |
+| 4 | `2.6e-06` | 1,569,152,761 |
+| 5 | `1.73e-06` | 3,104,527,573 |
+| 6 | `1.14e-06` | 6,125,805,627 |
+| 7 | `7.35e-07` | 12,255,330,837 |
+| 8 | `4.68e-07` | 24,419,998,198 |
+| 9 | `2.94e-07` | 48,317,129,677 |
+| 10 | `1.82e-07` | 94,684,031,363 |
+
+For a weak-scaling series that doubles `N`, advance one `k` per point. For strong scaling, choose
+one `k` that fits the smallest node count and keep it fixed. On another platform, build a fresh
+size sweep rather than treating this table as a target.
+
+## Launch templates
+
+For L1 and L2a, request one node and one rank. Set the core and partition counts from the granted
+CPU affinity:
+
+```bash
+export monoprop_PARTITIONS="$C"
+export monoprop_NUM_THREADS="$C"
+just bench L2a-hubbard-branch --hubbard-cutoff=10 --hubbard-lower-atol=3.38e-06 \
+    -k "test_model_propagate and hubbard"
+```
+
+For L2b, L3, and L4, request `R` ranks per node and `P` cores per rank. This Slurm example assumes
+the scheduler variables have already been set consistently:
+
+```bash
+export monoprop_PARTITIONS="$P"
+export monoprop_NUM_THREADS="$P"
+export monoprop_BENCH_LABEL=L3-hubbard-branch
+export monoprop_BENCH_RESULTS=benches/results
+
+srun --ntasks-per-node="$R" --cpus-per-task="$P" --cpu-bind=cores \
+    .venv/bin/python -m pytest benches -o filterwarnings=default \
+    --benchmark-json="benches/results/time-$monoprop_BENCH_LABEL.json" \
+    -k "test_model_propagate and hubbard" \
+    --hubbard-cutoff=10 --hubbard-lower-atol=3.38e-06
+
+uv run --no-sync monoprop-bench-report benches/results
+```
+
+Run the report outside the MPI launcher. Invoke the virtual-environment interpreter directly on
+each rank to avoid repeated dependency resolution on shared filesystems.
+
+## Validation and common errors
+
+- Repeat the rank count, cores per rank, and binding options on the launcher command; scheduler
+  allocation directives may not propagate them.
+- Keep `P` no larger than the CPUs visible to each rank. Otherwise thread placement is disabled.
+- Use the scheduler affinity mask, such as `os.sched_getaffinity(0)`, to verify visible CPUs;
+  `nproc` may honor unrelated thread limits.
+- Keep CPU binding enabled unless the target platform has been measured otherwise.
+- Interpret peak RSS as a process high-water mark. Allocator capacity and temporary reallocation
+  can dominate live object size.
+- For grouped graph operations, compare the maximum peak across operations, not their sum.
+- Verify term counts and the recorded shape before interpreting timing differences.
diff --git a/benches/README.md b/benches/README.md
index 32eb8dec..852d1ac2 100644
--- a/benches/README.md
+++ b/benches/README.md
@@ -11,7 +11,9 @@ documentation for detailed instructions.
 ## What lives where
 
 This directory holds only monoprop's own benchmarks — `conftest.py` (the fixtures
-and the results schema), `bench_random.py`, `bench_models.py`, and `results/`.
+and the results schema), `bench_random.py`, `bench_models.py`,
+[`LADDER.md`](LADDER.md) (the benchmark ladder for sensitive pull requests: the
+groups, their flags and their shapes) and `results/`.
 
 Both bench modules measure the same four operations — `build_graph`, `propagate`,
 `energy` and `gradient` — so a number means the same thing whichever problem
diff --git a/benches/bench_models.py b/benches/bench_models.py
index a5df09c5..8f21c4c9 100644
--- a/benches/bench_models.py
+++ b/benches/bench_models.py
@@ -130,12 +130,19 @@ def run(built, n_steps):
 @pytest.mark.slow
 @pytest.mark.parametrize("model", list(MODELS))
 def test_model_energy(
-    benchmark, model_graph, model, model_configs, bench_comm, bench_rounds, op_memory
+    benchmark,
+    model_graph,
+    model,
+    model_configs,
+    bench_comm,
+    bench_rounds,
+    op_memory,
+    pare_threshold,
 ):
     """Benchmark evaluating a fixed model's expectation-value functional."""
     skip_if_graph_will_not_fit(model, MODELS[model][2](model_configs[model]))
     propagator, parameters = model_graph(model)
-    functional = propagator.expectation_value_functional()
+    functional = propagator.expectation_value_functional(pare_threshold)
     op_memory.open()
     result = benchmark.pedantic(
         barriered(functional, bench_comm),
@@ -151,12 +158,19 @@ def test_model_energy(
 @pytest.mark.slow
 @pytest.mark.parametrize("model", list(MODELS))
 def test_model_gradient(
-    benchmark, model_graph, model, model_configs, bench_comm, bench_rounds, op_memory
+    benchmark,
+    model_graph,
+    model,
+    model_configs,
+    bench_comm,
+    bench_rounds,
+    op_memory,
+    pare_threshold,
 ):
     """Benchmark evaluating a fixed model's expectation-value-and-gradient functional."""
     skip_if_graph_will_not_fit(model, MODELS[model][2](model_configs[model]))
     propagator, parameters = model_graph(model)
-    functional = propagator.expectation_value_and_gradient_functional()
+    functional = propagator.expectation_value_and_gradient_functional(pare_threshold)
     op_memory.open()
     _value, gradient = benchmark.pedantic(
         barriered(functional, bench_comm),
diff --git a/benches/bench_random.py b/benches/bench_random.py
index 5c00b820..bdb4e4fd 100644
--- a/benches/bench_random.py
+++ b/benches/bench_random.py
@@ -26,6 +26,7 @@ def test_random_build_graph(
     bench_rounds,
     op_memory,
     record_opsize,
+    publish_graph,
 ):
     """Benchmark building the propagation graph from a fresh propagator."""
     last = []
@@ -47,8 +48,11 @@ def build(built):
         iterations=1,
     )
     op_memory.close(last[0])
-    assert record_opsize(last[0]) > 0
+    terms = record_opsize(last[0])
+    assert terms > 0
     assert last[0].graph_layers > 0
+    # Hands the timed graph to energy/gradient, so one process builds once.
+    publish_graph(last[0], terms)
 
 
 def test_random_propagate(
@@ -83,10 +87,16 @@ def run(built):
 
 
 def test_random_energy(
-    benchmark, built_graph, random_problem, bench_comm, bench_rounds, op_memory
+    benchmark,
+    built_graph,
+    random_problem,
+    bench_comm,
+    bench_rounds,
+    op_memory,
+    pare_threshold,
 ):
     """Benchmark evaluating the expectation-value functional."""
-    functional = built_graph.expectation_value_functional()
+    functional = built_graph.expectation_value_functional(pare_threshold)
     op_memory.open()
 
     # The entry barrier must stay untimed, so it lives in a setup, and a setup that returns args forbids args=.
@@ -104,10 +114,16 @@ def setup():
 
 
 def test_random_gradient(
-    benchmark, built_graph, random_problem, bench_comm, bench_rounds, op_memory
+    benchmark,
+    built_graph,
+    random_problem,
+    bench_comm,
+    bench_rounds,
+    op_memory,
+    pare_threshold,
 ):
     """Benchmark evaluating the expectation-value-and-gradient functional."""
-    functional = built_graph.expectation_value_and_gradient_functional()
+    functional = built_graph.expectation_value_and_gradient_functional(pare_threshold)
     op_memory.open()
 
     def setup():
diff --git a/benches/conftest.py b/benches/conftest.py
index b2af8acd..40ab612c 100644
--- a/benches/conftest.py
+++ b/benches/conftest.py
@@ -123,21 +123,27 @@ def _spread(comm: Any, value: int) -> dict[str, int]:
     return {"sum": _reduce_sum(comm, value), "max": _reduce_max(comm, value)}
 
 
+# (name, type, default, help).
 _RANDOM_OPTIONS = (
-    ("gen-length", 4, "Majorana operators per generator."),
-    ("obs-terms", 10000, "Observable terms."),
-    ("num-generators", 100, "Random generators (circuit gates)."),
-    ("num-modes", 128, "Fermionic modes."),
-    ("cutoff", 6, "Truncation cutoff."),
-    ("seed", 0, "Random seed."),
-    ("bench-rounds", 1, "Fixed timing rounds (MPI-safe)."),
+    ("gen-length", int, 4, "Majorana operators per generator."),
+    ("obs-terms", int, 10000, "Observable terms."),
+    ("num-generators", int, 100, "Random generators (circuit gates)."),
+    ("num-modes", int, 128, "Fermionic modes."),
+    ("cutoff", int, 6, "Truncation cutoff."),
+    ("seed", int, 0, "Random seed."),
+    ("bench-rounds", int, 1, "Fixed timing rounds (MPI-safe)."),
+    ("pare-threshold", float, None, "Edge-retention cutoff for the graph functionals."),
 )
 
+# Picture -> built graph, so one process can time build_graph, then energy, then gradient
+# against the graph it just built rather than building a second one.
+_GRAPH_CACHE: dict[str, Any] = {}
+
 _RESULTS: dict[str, Any] = {
     "meta": {},  # run configuration (ranks, threads, host, ...)
     "params": {},  # resolved random-problem hyperparameters
-    "memhwm": {},  # node id -> summed peak RSS, whole test, setup() included
-    "memhwm_max": {},  # node id -> worst-rank peak RSS, whole test, setup() included
+    "memhwm": {},  # pytest node id -> summed-over-ranks peak RSS, whole test, setup() included
+    "memhwm_max": {},  # pytest node id -> worst-rank peak RSS, whole test, setup() included
     "opsize": {},  # picture / model / node id -> {"terms": n}
     "configs": {},  # fixed model -> config dataclass fields
     "opmem": {},  # fixed model -> per-field operator memory split (bytes)
@@ -168,8 +174,8 @@ def _results_path() -> Path | None:
 def pytest_addoption(parser: pytest.Parser) -> None:
     """Register benchmark configuration options for the random benchmarks."""
     group = parser.getgroup("monoprop-bench", "monoprop random benchmark sizing")
-    for name, default, help_text in _RANDOM_OPTIONS:
-        group.addoption(f"--{name}", type=int, default=default, help=help_text)
+    for name, kind, default, help_text in _RANDOM_OPTIONS:
+        group.addoption(f"--{name}", type=kind, default=default, help=help_text)
 
     models = parser.getgroup("monoprop-models", "monoprop fixed-model overrides")
     for model, (config_cls, _builder, _steps) in MODELS.items():
@@ -226,9 +232,8 @@ def _meta(nodes: int, ranks_per_node: int) -> dict[str, Any]:
     }
     # No nanobind property exposes the engine's resolved partition count (see bindings/binder.h),
     # so record the requested env var under its own name rather than claim it is the effective value.
-    partitions_env = os.environ.get("monoprop_PARTITIONS")  # noqa: SIM112
-    if partitions_env is not None:
-        meta["partitions_env"] = partitions_env
+    # Recorded even when absent: "unset" and a missing field render identically in the report.
+    meta["partitions_env"] = os.environ.get("monoprop_PARTITIONS", "unset")  # noqa: SIM112
     return meta
 
 
@@ -236,10 +241,28 @@ def _params(config: pytest.Config) -> dict[str, Any]:
     """Return the resolved random-problem hyperparameters (defaults included)."""
     return {
         name.replace("-", "_"): config.getoption(f"--{name}")
-        for name, _default, _help in _RANDOM_OPTIONS
+        for name, _kind, _default, _help in _RANDOM_OPTIONS
     }
 
 
+def _require_shape() -> None:
+    """Fail a multi-rank run that has not declared its partition count.
+
+    ``resolve_partition_count_`` defaults to ``ranks == 1 ? cores : 1``, so an unset knob
+    under MPI measures one partition per rank -- a single-threaded run at a plausible wall
+    time, with nothing in the timing to say so. Every rank raises, so no rank is left in a
+    collective. Serial runs take the engine default and are unaffected.
+    """
+    if _size() > 1 and not os.environ.get("monoprop_PARTITIONS"):  # noqa: SIM112
+        msg = (
+            "monoprop_PARTITIONS is unset on a run of "
+            f"{_size()} ranks. The engine would use one partition per rank. Export the "
+            "shape you intend, e.g. `export monoprop_PARTITIONS=16 monoprop_NUM_THREADS=16` "
+            "for 8 ranks per 128-core node."
+        )
+        raise pytest.UsageError(msg)
+
+
 @pytest.hookimpl(trylast=True)
 def pytest_configure(config: pytest.Config) -> None:
     """Record run metadata on rank 0; silence the other ranks under MPI.
@@ -248,6 +271,7 @@ def pytest_configure(config: pytest.Config) -> None:
     ``trylast`` so the terminal reporter exists before non-root ranks unregister it.
     """
     nodes, ranks_per_node = _nodes()  # collective; every rank must call this
+    _require_shape()
     if _rank() == 0:
         _RESULTS["meta"] = _meta(nodes, ranks_per_node)
         _RESULTS["params"] = _params(config)
@@ -289,6 +313,12 @@ def bench_rounds(request: pytest.FixtureRequest) -> int:
     return int(request.config.getoption("--bench-rounds"))
 
 
+@pytest.fixture(scope="session")
+def pare_threshold(request: pytest.FixtureRequest) -> float | None:
+    """Return the edge-retention cutoff for the graph functionals (``None`` keeps the exact graph)."""
+    return request.config.getoption("--pare-threshold")
+
+
 @pytest.fixture(scope="session")
 def model_configs(request: pytest.FixtureRequest) -> dict[str, Any]:
     """Return each fixed model's config, every field resolved from the CLI.
@@ -520,6 +550,9 @@ def built_graph(
     read-only graph benchmarks (``energy``, ``gradient``), and records the operator
     size for this picture while the graph is resident.
     """
+    if picture in _GRAPH_CACHE:
+        return _GRAPH_CACHE[picture]
+
     mp, circuit = build_random_propagator(
         random_problem, comm=bench_comm, schrodinger=picture == "schrodinger"
     )
@@ -528,9 +561,25 @@ def built_graph(
     # Under MPI the operator is partitioned, so sum the partitions.
     _record("opsize", picture, {"terms": _reduce_sum(bench_comm, mp.size())})
 
+    _GRAPH_CACHE[picture] = mp
     return mp
 
 
+@pytest.fixture(scope="session")
+def publish_graph(picture: str) -> Callable[[Any, int], None]:
+    """Return ``publish(propagator, terms)``, handing a timed build to ``built_graph``.
+
+    Takes the term count rather than recomputing it: ``record_opsize`` has already
+    reduced it, and a second collective reduce over a 1B-term operator is not free.
+    """
+
+    def _publish(mp: Any, terms: int) -> None:
+        _GRAPH_CACHE[picture] = mp
+        _record("opsize", picture, {"terms": terms})
+
+    return _publish
+
+
 @pytest.fixture(scope="session")
 def model_graph(
     model_configs: dict[str, Any], bench_comm: Any
diff --git a/docs/content/docs/benchmarks.mdx b/docs/content/docs/benchmarks.mdx
index 36326af1..d72f679b 100644
--- a/docs/content/docs/benchmarks.mdx
+++ b/docs/content/docs/benchmarks.mdx
@@ -294,11 +294,24 @@ across ranks. Memory is the sum of the per-rank peak RSS values. This is an uppe
 job-wide peak: shared pages are counted once per rank, and rank-local peaks are summed even when
 they occur at different times.
 
+Above one rank, `monoprop_PARTITIONS` is mandatory: the engine's own default is
+`ranks == 1 ? cores : 1`, so an unset knob would silently measure one partition per rank, and the
+benchmark session refuses to start without it. The shape it declares is recorded in each run's
+`meta` and shown in the report's configuration table.
+
 ```bash
 just bench-build-mpi                           # build once (MPI on)
-monoprop_NUM_THREADS=2 just bench-mpi r5t2 5 --map-by slot:PE=2 --bind-to core
+monoprop_PARTITIONS=2 monoprop_NUM_THREADS=2 just bench-mpi r5t2 5 --map-by slot:PE=2 --bind-to core
 ```
 
+### Benchmarking ladder
+
+`benches/LADDER.md` holds the ladder used for large-scale MPI benchmarking during development: five
+groups of configurations, the flags and `-k` selector that produce each row, and the shapes to run
+them at. Each rung differs from the one below it in exactly one thing -- partitions (L1 to L2a),
+processes (L2a to L2b), the network (L2b to L3), then node count (L3 to L4) -- so a regression that
+appears at one rung and not the one below it names its own cause.
+
 ### Continuous benchmarking
 
 The `bench_main.yml` workflow benchmarks every commit on `main`. [Bencher](https://bencher.dev/)
diff --git a/justfile b/justfile
index da08ee19..b0f2055a 100644
--- a/justfile
+++ b/justfile
@@ -349,9 +349,10 @@ bench LABEL *ARGS:
     uv run --no-sync monoprop-bench-report "{{ bench_results }}"
 
 # Needs an MPI build (`just bench-build-mpi`) -- a non-MPI build is rejected by the
-# preflight. Extra args are passed to mpiexec for pinning (and, as root, add
+# preflight. `monoprop_PARTITIONS` is mandatory above one rank (benches/conftest.py refuses
+# the run without it). Extra args are passed to mpiexec for pinning (and, as root, add
 # `--allow-run-as-root`), e.g.
-#   monoprop_NUM_THREADS=2 just bench-mpi r5t2 5 --map-by slot:PE=2 --bind-to core
+#   monoprop_PARTITIONS=2 monoprop_NUM_THREADS=2 just bench-mpi r5t2 5 --map-by slot:PE=2 --bind-to core
 
 # Run under MPI: RANKS ranks recorded as one LABEL column.
 bench-mpi LABEL RANKS *MPIARGS:
@@ -360,7 +361,9 @@ bench-mpi LABEL RANKS *MPIARGS:
     label="$1"; ranks="$2"; shift 2; \
     monoprop_BENCH_LABEL="$label" monoprop_BENCH_RESULTS="{{ bench_results }}" \
         uv run --no-sync mpiexec -n "$ranks" \
-        -x monoprop_BENCH_LABEL -x monoprop_BENCH_RESULTS "$@" \
+        -x monoprop_BENCH_LABEL -x monoprop_BENCH_RESULTS \
+        ${monoprop_PARTITIONS+-x monoprop_PARTITIONS} \
+        ${monoprop_NUM_THREADS+-x monoprop_NUM_THREADS} "$@" \
         python -m pytest benches -o filterwarnings=default \
         --benchmark-json="{{ bench_results }}/time-$label.json"
     uv run --no-sync monoprop-bench-report "{{ bench_results }}"
diff --git a/packages/monoprop-bench-tools/src/monoprop_bench_tools/models.py b/packages/monoprop-bench-tools/src/monoprop_bench_tools/models.py
index 1c9ae567..a95c1c36 100644
--- a/packages/monoprop-bench-tools/src/monoprop_bench_tools/models.py
+++ b/packages/monoprop-bench-tools/src/monoprop_bench_tools/models.py
@@ -28,6 +28,7 @@
 import numpy as np
 
 from monoprop import (
+    MAX_NUM_MODES,
     Circuit,
     ExpGate,
     MajoranaPropagator,
@@ -120,6 +121,28 @@ def parameters(self) -> np.ndarray:
         return np.asarray(self.circuit.parameters, dtype=float)
 
 
+def _check_gen_length(gen_length: int, option: str) -> None:
+    """Reject a generator length whose real-coefficient monomials are not Hermitian."""
+    if gen_length % 4 not in (0, 1):
+        msg = (
+            f"{option} gives length {gen_length}. A product of k Majoranas obeys "
+            f"(g_1...g_k)^dag = (-1)^(k(k-1)/2) g_1...g_k, so it is Hermitian only for "
+            f"k % 4 in (0, 1); at k={gen_length} the real coefficients drawn here are "
+            "anti-Hermitian and the propagator's constructor rejects them."
+        )
+        raise ValueError(msg)
+
+
+def _check_modes(num_modes: int, option: str) -> None:
+    """Reject a size above the extension's compile-time mode limit."""
+    if num_modes > MAX_NUM_MODES:
+        msg = (
+            f"{option} gives {num_modes} modes, above this build's MAX_NUM_MODES "
+            f"({MAX_NUM_MODES}). Rebuild with -Dmonoprop_MAX_NUM_MODES to go higher."
+        )
+        raise ValueError(msg)
+
+
 def _random_terms(
     rng: np.random.Generator,
     num_terms: int,
@@ -166,7 +189,9 @@ def make_random_problem(
     """Build a random observable and generator circuit for benchmarking.
 
     Args:
-        gen_length: Number of Majorana operators in each random generator.
+        gen_length: Number of Majorana operators in each random generator. Must satisfy
+            ``gen_length % 4 in (0, 1)``, the lengths at which a real-coefficient
+            Majorana product is Hermitian.
         obs_terms: Number of terms in the random observable.
         num_generators: Number of random generators in the circuit.
         num_modes: Number of fermionic modes (Majorana indices = ``2 * num_modes``).
@@ -177,6 +202,8 @@ def make_random_problem(
     Returns:
         A :class:`RandomProblem` bundling the observable, circuit, and cutoff.
     """
+    _check_modes(num_modes, "--num-modes")
+    _check_gen_length(gen_length, "--gen-length")
     rng = np.random.default_rng(seed)
     num_majorana_indices = 2 * num_modes
 
@@ -321,6 +348,15 @@ def build_hubbard_problem(
         ``(propagator, circuit)`` for one Trotter step, re-applied by the driver.
     """
     config = config or HubbardConfig()
+    _check_modes(config.num_qubits, "--hubbard-num-sites")
+    # A lattice position, not a constant: shrinking the lattice without scaling it shifts the light cone.
+    if not 0 <= config.observable_site < config.num_sites:
+        msg = (
+            f"--hubbard-observable-site {config.observable_site} is outside the "
+            f"{config.num_sites}-site lattice. Scale it with --hubbard-num-sites "
+            "(the default 46 sits 46/60 along the lattice)."
+        )
+        raise ValueError(msg)
     fermi_gates = [ExpGate(term) for term in _hubbard_fermion_terms(config)]
     parameters = [config.trotter_dt] * len(fermi_gates)
     occupied = _neel_occupied_modes(config.num_sites, config.neel_start_spin)
@@ -501,12 +537,15 @@ def build_hubbard_problem(
     (125, 126),
 ]
 
+# The map is fixed, so it also fixes the register: qubit counts are checked against it.
+HEAVY_HEX_NUM_QUBITS: int = 1 + max(j for _, j in HEAVY_HEX_TOPOLOGY)
+
 
 @dataclass(frozen=True, slots=True)
 class KickedIsingConfig:
     """Configuration for the static Pauli-basis (kicked-Ising) benchmark."""
 
-    num_qubits: int = 127
+    num_qubits: int = HEAVY_HEX_NUM_QUBITS
     num_layers: int = 20
     observable_qubit: int = 62
     theta: float = np.pi / 4
@@ -551,6 +590,19 @@ def build_kicked_ising_problem(
         ``(propagator, circuit)`` ready for in-place propagation.
     """
     config = config or KickedIsingConfig()
+    if config.num_qubits != HEAVY_HEX_NUM_QUBITS:
+        msg = (
+            f"--pauli-num-qubits must be {HEAVY_HEX_NUM_QUBITS}: HEAVY_HEX_TOPOLOGY is the "
+            "fixed IBM Eagle coupling map, so any other count either indexes past the "
+            "operator or leaves qubits uncoupled. Size this model with --pauli-lower-atol."
+        )
+        raise ValueError(msg)
+    if not 0 <= config.observable_qubit < config.num_qubits:
+        msg = (
+            f"--pauli-observable-qubit {config.observable_qubit} is outside the "
+            f"{config.num_qubits}-qubit register."
+        )
+        raise ValueError(msg)
     gate_angles: list[tuple[ExpGate, float]] = []
     for _ in range(config.num_layers):
         gate_angles.extend(_xlayer(config.num_qubits, config.theta / 2))
diff --git a/packages/monoprop-bench-tools/tests/test_models.py b/packages/monoprop-bench-tools/tests/test_models.py
index 9c54865e..02a6c6be 100644
--- a/packages/monoprop-bench-tools/tests/test_models.py
+++ b/packages/monoprop-bench-tools/tests/test_models.py
@@ -18,6 +18,7 @@
 
 import inspect
 
+import pytest
 from monoprop_bench_tools.models import build_random_propagator, make_random_problem
 
 
@@ -47,3 +48,31 @@ def test_built_graph_is_populated(serial_comm) -> None:
     propagator.build_graph(circuit)
     _n_cos_indices, n_cycles = propagator.graph_size()
     assert n_cycles > 0
+
+
+@pytest.mark.parametrize("gen_length", [4, 5, 8, 9])
+def test_hermitian_gen_lengths_are_accepted(gen_length: int) -> None:
+    problem = make_random_problem(
+        gen_length=gen_length,
+        obs_terms=3,
+        num_generators=2,
+        num_modes=8,
+        cutoff=6,
+        seed=0,
+    )
+    assert problem.circuit.system_size == 8
+
+
+@pytest.mark.parametrize("gen_length", [2, 3, 6, 7])
+def test_anti_hermitian_gen_lengths_raise_naming_the_option(gen_length: int) -> None:
+    # Unvalidated, the length reaches the extension and surfaces as a RuntimeError from the
+    # propagator's constructor, naming neither --gen-length nor the parity rule.
+    with pytest.raises(ValueError, match="--gen-length"):
+        make_random_problem(
+            gen_length=gen_length,
+            obs_terms=3,
+            num_generators=2,
+            num_modes=8,
+            cutoff=6,
+            seed=0,
+        )