diff --git a/tools/simulators/tsi/tsi-sim-mc/.gitignore b/tools/simulators/tsi/tsi-sim-mc/.gitignore
new file mode 100644
index 0000000..c4ee506
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim-mc/.gitignore
@@ -0,0 +1,18 @@
+# Generated artifacts
+runs/
+results/
+figures/
+!results/.gitkeep
+!figures/.gitkeep
+
+# Python
+__pycache__/
+*.py[cod]
+*.egg-info/
+.venv/
+venv/
+.pytest_cache/
+.mypy_cache/
+.ruff_cache/
+*.parquet
+*.csv
diff --git a/tools/simulators/tsi/tsi-sim-mc/Makefile b/tools/simulators/tsi/tsi-sim-mc/Makefile
new file mode 100644
index 0000000..4521f9c
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim-mc/Makefile
@@ -0,0 +1,57 @@
+.PHONY: install smoke sweep sweep-fullscale figures verify test test-slow test-all lint clean
+
+VENV ?= .venv
+PY   := $(VENV)/bin/python
+PIP  := $(VENV)/bin/pip
+
+# Pin BLAS/OpenMP to a single thread so joblib's N worker processes don't oversubscribe
+# the cores (numpy in each worker would otherwise each spawn a full BLAS thread pool).
+export OMP_NUM_THREADS := 1
+export OPENBLAS_NUM_THREADS := 1
+export MKL_NUM_THREADS := 1
+export NUMEXPR_NUM_THREADS := 1
+
+$(VENV):
+	python3 -m venv $(VENV)
+	$(PIP) install --upgrade pip
+
+install: $(VENV)
+	$(PIP) install -e ".[dev]"
+
+# Every sweep writes results + figures into a fresh dated folder under runs/ so new
+# runs never overwrite old ones (runs/<YYYY-MM-DD_HHMMSS>_<label>/).
+
+# Fast end-to-end check (tiny scaled-k grid); uses all cores across configs.
+smoke: install
+	$(PY) scripts/run_sweep.py --config configs/smoke.yaml --label smoke
+
+# Full scaled-k parameter sweep (multicore across configs; auto-generates figures).
+sweep: install
+	$(PY) scripts/run_sweep.py --config configs/default.yaml --label default
+
+# Full-scale (true k) grid; small + heavy, so load-balance one config at a time.
+sweep-fullscale: install
+	$(PY) scripts/run_sweep.py --config configs/fullscale.yaml --label fullscale --batch-size 1
+
+# Re-render figures from a run's parquet into a fresh dated figures/ folder:
+#   make figures RESULTS=runs/<dir>/results.parquet
+figures: install
+	$(PY) scripts/make_figures.py --results $(RESULTS)
+
+# Analytic sanity checks (simulator vs closed-form theory; replicates run across cores).
+verify: install
+	$(PY) scripts/verify.py
+
+test: install          # fast subset (addopts already excludes slow)
+	$(PY) -m pytest
+test-slow: install
+	$(PY) -m pytest -m slow
+test-all: install
+	$(PY) -m pytest -m ''
+
+lint: install
+	$(VENV)/bin/ruff check src scripts tests
+
+clean:
+	rm -rf runs/* results/*.parquet figures/* .pytest_cache .ruff_cache .mypy_cache
+	find . -name __pycache__ -type d -prune -exec rm -rf {} +
diff --git a/tools/simulators/tsi/tsi-sim-mc/README.md b/tools/simulators/tsi/tsi-sim-mc/README.md
new file mode 100644
index 0000000..45fef8f
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim-mc/README.md
@@ -0,0 +1,102 @@
+# tsi-sim-mc — Cryptarchia TSI simulator (multicore build)
+
+> **This is the multicore-optimised, reviewed copy of `../tsi-sim/`.** The original is left
+> untouched. Versus the original it adds: a **sparse lottery sampler** (~30–100× faster per
+> epoch), a **hardened multicore sweep** (loky + single-thread BLAS), an **opt-in parallel
+> chunked lottery**, an optional **spec fixed-point mode**, corrected figures, config
+> validation, and a larger test suite. See "Performance & reproducibility" below.
+
+Monte-Carlo simulation framework used to choose safe values for the **uncle-reference**
+parameters of Cryptarchia's Total Stake Inference (TSI):
+
+- `U` — max uncles referenced per block (`MAX_UNCLES`); `U=0` is the no-uncle baseline.
+- `W` — uncle reference window in slots (spec default 300).
+- swept against network size `N`, stake distribution `S` (uniform / Pareto), and
+  network **latency** `L` (in slots — deliberately *not* `D`, which denotes the stake estimate).
+
+It measures how well the inferred total active stake `D` tracks the true total stake, and
+whether uncle references recover the active slots that network latency loses to forks.
+
+> This lives under a `raw/` docs path, so it is invisible to the repository's
+> markdown-lint CI. It is a standalone Python package with its own tooling.
+
+## Model
+
+Reduced **canonical-chain-with-orphans** model: we simulate the global winning-slot
+sequence (stake-weighted φ lottery), build a real block tree with latency- and
+multi-winner-induced forks, resolve the canonical chain (honest longest-chain), let
+canonical blocks reference uncles per the spec's selection rules, and count TSI density
+`m = honest-chain blocks + deduplicated referenced uncles` in the measurement window.
+All honest nodes converge to the same deep chain (k-finality), so a single per-epoch
+`D` is faithful. A full per-node model is the planned next phase (`per_node_dest` flag
+scaffolds it).
+
+See the sibling spec `../` and `../../cryptarchia-total-stake-inference.md` for the math.
+
+## Quick start
+
+```bash
+make install         # create .venv and install (editable) with dev deps
+make test            # unit tests + fast theory checks
+make verify          # simulator vs closed-form analytic checks
+make smoke           # tiny scaled-k sweep + figures (end-to-end smoke test)
+make sweep figures   # full scaled-k parameter sweep + academic figures
+```
+
+Outputs: `results/*.parquet` (one row per config×epoch) and `figures/*.{pdf,png}`
+(both git-ignored).
+
+## Scale
+
+True constants (`k=2160`, `f=1/30`) give 648,000-slot epochs. Mean accuracy is provably
+`k`-invariant (only variance scales `~1/T`), so sweeps use a **scaled `k`**
+(`configs/default.yaml`); the final accuracy/variance figures re-run at true `k`
+(`configs/fullscale.yaml`, `make sweep-fullscale`). `configs/smoke.yaml` is a tiny dev grid.
+
+## Performance & reproducibility
+
+- **Sparse lottery (the main win).** The number of slots a node wins is `Binomial(n_slots,
+  p_i)` and the won slots are a uniform distinct subset — distributionally identical to an
+  independent Bernoulli per slot, but without the dense `(n_nodes, n_slots)` array that was
+  ~95% of runtime. Full-scale (`k=2160`) epoch: **~3.3 s → ~0.1 s (~30×)**; at `k=256`,
+  **~0.39 s → ~0.009 s (~40×)**.
+- **Multicore across configs (the main lever).** `run_trajectory` is a pure function of a
+  hash-seeded, immutable `SimConfig`, so the sweep is order-independent and embarrassingly
+  parallel. `run_sweep` uses joblib's process-based **loky** backend with
+  `inner_max_num_threads=1` (and the Makefile pins `*_NUM_THREADS=1`) to use all cores
+  without BLAS oversubscription. `--n-jobs -1` (default) uses every core; `--batch-size 1`
+  suits the small heavy full-scale grid. Measured on a 14-core box, the full scaled-`k`
+  sweep (`configs/default.yaml`, 3888 configs) runs in **~58 s parallel vs ~561 s serial
+  (9.7×)**.
+- **Opt-in within-config parallelism.** `lottery_chunks > 1` splits the per-slot lottery
+  across slot-chunks with independent `SeedSequence.spawn` children. After the sparse fix the
+  lottery is a small fraction of an epoch, so this rarely helps — it exists for a single
+  isolated config with an enormous `n_slots`. **Reproducibility caveat:** the chunked stream
+  differs from the serial stream and *changes with `n_chunks`*, so `lottery_chunks` must be a
+  pinned, recorded parameter, never derived from the core count.
+- **RNG reproducibility.** Every draw is a deterministic spawn off `SeedSequence(hash(config))`
+  — child 0 draws stake, child `e+1` drives epoch `e`, which spawns lottery/aux sub-streams.
+  Results are identical regardless of parallel scheduling order. The sparse sampler consumes
+  the RNG differently from the original dense one, so committed baselines here were
+  regenerated against the sparse sampler.
+
+## Faithfulness notes (from review)
+
+- `update_D(..., fixed_point=True)` mirrors the spec's integer `f`-truncation
+  (`int(f·1000)/1000 = 0.033`), reproducing the on-chain estimator's ~1% systematic
+  overestimate. Default `False` keeps the exact-`f`, analysis-faithful behaviour.
+- TSI counts *blocks* (wins), so at full uncle recovery the estimate equilibrates at the
+  **block-count ceiling** `-ln(1-f)/f ≈ 1.017`, not 1.0 — a deterministic overshoot, not
+  noise. Figures overlay this ceiling; `theory.block_count_ceiling` computes it.
+
+## Layout
+
+```
+src/tsi_sim/   constants config rng stake lottery latency blocktree uncles tsi epoch
+               engine metrics theory sweep verify  plotting/{style,figures,make_figures}
+scripts/       run_sweep.py  make_figures.py  verify.py   (thin shims; installed as
+                                                            tsi-sweep / tsi-figures / tsi-verify)
+configs/       smoke.yaml  default.yaml  fullscale.yaml
+tests/         test_{lottery,uncles,tsi_counting,blocktree,config,rng,stake,theory,
+                     latency,theory_convergence}.py
+```
diff --git a/tools/simulators/tsi/tsi-sim-mc/configs/default.yaml b/tools/simulators/tsi/tsi-sim-mc/configs/default.yaml
new file mode 100644
index 0000000..fddf3a8
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim-mc/configs/default.yaml
@@ -0,0 +1,16 @@
+# Full scaled-k parameter sweep. Mean accuracy / q_eff / convergence are k-invariant,
+# so k is scaled down for tractability; re-run headline + variance figures at full scale
+# with fullscale.yaml.
+# Sweep axes (cartesian product x replicates); every value is a list.
+n_nodes: [1000, 2000, 4000]                # number of nodes / stake holders
+stake_dist: [uniform, pareto]              # stake distribution (uniform = equal, pareto = heavy-tailed)
+latency: [0, 1, 2, 4, 6, 8, 10, 15, 20, 30] # L: network latency in slots; block visible at t+L
+max_uncles: [0, 1, 2, 3, 4, 8, 10]         # U: max uncle references per block (0 = baseline)
+uncle_strategy: [oldest, random]           # uncle selection: oldest-first fill vs random coin-flip
+replicates: 12                             # independent RNG replicates per grid cell
+base:                                      # per-run settings shared by every cell (not swept)
+  k: 256                                   # security parameter (scaled; T = 6*floor(256/f) = 46080 slots)
+  epochs: 75                               # epochs simulated per trajectory
+  f: 0.03333333333333333                  # slot activation coefficient (default 1/30); configurable
+  genesis_d_factor: 0.5                    # genesis D_est = factor x true total stake
+  pareto_shape: 1.16                       # Pareto (Lomax) tail index for the pareto stake distribution
diff --git a/tools/simulators/tsi/tsi-sim-mc/configs/fullscale.yaml b/tools/simulators/tsi/tsi-sim-mc/configs/fullscale.yaml
new file mode 100644
index 0000000..669cb94
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim-mc/configs/fullscale.yaml
@@ -0,0 +1,16 @@
+# Full-scale (true k=2160) confirmation runs. Same latency/uncle grid as default.yaml but
+# at the true security parameter (k not scaled down); heavier per config, so N and reps are
+# kept modest. Note: T = 6*floor(2160/f) = 388,800 slots per epoch.
+# Sweep axes (cartesian product x replicates); every value is a list.
+n_nodes: [1000]                             # number of nodes / stake holders
+stake_dist: [pareto]                        # heavy-tailed (realistic) stake distribution
+latency: [0, 1, 2, 4, 6, 8, 10, 15, 20, 30] # L: network latency in slots; block visible at t+L
+max_uncles: [0, 1, 2, 3, 4, 8, 10]          # U: max uncle references per block (0 = baseline)
+uncle_strategy: [oldest]                    # uncle selection: oldest-first fill
+replicates: 10                              # independent RNG replicates per grid cell
+base:                                       # per-run settings shared by every cell (not swept)
+  k: 2160                                   # true security parameter (not scaled)
+  epochs: 60                                # epochs simulated per trajectory
+  f: 0.03333333333333333                   # slot activation coefficient (default 1/30); configurable
+  genesis_d_factor: 0.001                   # genesis D_est = factor x true total stake
+  pareto_shape: 1.16                        # Pareto (Lomax) tail index for the pareto stake distribution
diff --git a/tools/simulators/tsi/tsi-sim-mc/configs/smoke.yaml b/tools/simulators/tsi/tsi-sim-mc/configs/smoke.yaml
new file mode 100644
index 0000000..15947fd
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim-mc/configs/smoke.yaml
@@ -0,0 +1,14 @@
+# Tiny end-to-end smoke grid (fast; for dev + CI-style checks).
+# Sweep axes (cartesian product x replicates); every value is a list.
+n_nodes: [1000]                  # number of nodes / stake holders
+stake_dist: [uniform, pareto]    # stake distribution (uniform = equal, pareto = heavy-tailed)
+latency: [0, 4, 8]               # L: network latency in slots; block visible to others at t+L
+max_uncles: [0, 1, 4]            # U: max uncle references per block (0 = baseline, no uncles)
+uncle_strategy: [oldest]         # uncle selection: oldest-first fill
+replicates: 4                    # independent RNG replicates per grid cell
+base:                            # per-run settings shared by every cell (not swept)
+  k: 16                          # security parameter (scaled; T = 6*floor(16/f) = 2880 slots)
+  epochs: 25                     # epochs simulated per trajectory
+  f: 0.03333333333333333        # slot activation coefficient (default 1/30); configurable
+  genesis_d_factor: 0.5          # genesis D_est = factor x true total stake
+  pareto_shape: 1.16             # Pareto (Lomax) tail index for the pareto stake distribution
diff --git a/tools/simulators/tsi/tsi-sim-mc/pyproject.toml b/tools/simulators/tsi/tsi-sim-mc/pyproject.toml
new file mode 100644
index 0000000..6b87fd1
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim-mc/pyproject.toml
@@ -0,0 +1,47 @@
+[build-system]
+requires = ["hatchling"]
+build-backend = "hatchling.build"
+
+[project]
+name = "tsi-sim"
+version = "0.1.0"
+description = "Monte-Carlo simulation of Cryptarchia Total Stake Inference with uncle references"
+readme = "README.md"
+requires-python = ">=3.11"
+license = { text = "MIT" }
+authors = [{ name = "Logos" }]
+dependencies = [
+    "numpy>=1.26",
+    "pandas>=2.1",
+    "pyarrow>=14",
+    "matplotlib>=3.8",
+    "scipy>=1.11",
+    "pyyaml>=6.0",
+    "tqdm>=4.66",
+    "joblib>=1.3",
+]
+
+[project.optional-dependencies]
+dev = ["pytest>=8.0", "pytest-xdist>=3.5", "ruff>=0.5", "mypy>=1.8"]
+accel = ["numba>=0.60"]   # optional; only used if a single-config tree build is a bottleneck
+
+[project.scripts]
+tsi-sweep = "tsi_sim.sweep:main"
+tsi-verify = "tsi_sim.verify:main"
+tsi-figures = "tsi_sim.plotting.make_figures:main"
+
+[tool.hatch.build.targets.wheel]
+packages = ["src/tsi_sim"]
+
+[tool.ruff]
+line-length = 100
+target-version = "py311"
+
+[tool.ruff.lint]
+select = ["E", "F", "I", "UP", "B", "NPY"]
+ignore = ["E741"]  # allow single-char names like L (latency), q, m — domain symbols
+
+[tool.pytest.ini_options]
+testpaths = ["tests"]
+markers = ["slow: long-running (full-scale k) checks"]
+addopts = "-q -m 'not slow'"   # fast subset by default; `make test-all` runs everything
diff --git a/tools/simulators/tsi/tsi-sim-mc/requirements-dev.txt b/tools/simulators/tsi/tsi-sim-mc/requirements-dev.txt
new file mode 100644
index 0000000..aefbcb6
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim-mc/requirements-dev.txt
@@ -0,0 +1 @@
+-e .[dev]
diff --git a/tools/simulators/tsi/tsi-sim-mc/requirements.txt b/tools/simulators/tsi/tsi-sim-mc/requirements.txt
new file mode 100644
index 0000000..d6e1198
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim-mc/requirements.txt
@@ -0,0 +1 @@
+-e .
diff --git a/tools/simulators/tsi/tsi-sim-mc/scripts/make_figures.py b/tools/simulators/tsi/tsi-sim-mc/scripts/make_figures.py
new file mode 100644
index 0000000..314a219
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim-mc/scripts/make_figures.py
@@ -0,0 +1,14 @@
+#!/usr/bin/env python
+"""Thin shim so `python scripts/make_figures.py` works without installing."""
+
+from __future__ import annotations
+
+import sys
+from pathlib import Path
+
+sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))
+
+from tsi_sim.plotting.make_figures import main  # noqa: E402
+
+if __name__ == "__main__":
+    main()
diff --git a/tools/simulators/tsi/tsi-sim-mc/scripts/run_sweep.py b/tools/simulators/tsi/tsi-sim-mc/scripts/run_sweep.py
new file mode 100644
index 0000000..77d98a0
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim-mc/scripts/run_sweep.py
@@ -0,0 +1,14 @@
+#!/usr/bin/env python
+"""Run a TSI parameter sweep from a YAML config -> results parquet."""
+
+from __future__ import annotations
+
+import sys
+from pathlib import Path
+
+sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))
+
+from tsi_sim.sweep import main  # noqa: E402
+
+if __name__ == "__main__":
+    main()
diff --git a/tools/simulators/tsi/tsi-sim-mc/scripts/verify.py b/tools/simulators/tsi/tsi-sim-mc/scripts/verify.py
new file mode 100644
index 0000000..126eda3
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim-mc/scripts/verify.py
@@ -0,0 +1,14 @@
+#!/usr/bin/env python
+"""Thin shim so `python scripts/verify.py` works without installing; see tsi_sim.verify."""
+
+from __future__ import annotations
+
+import sys
+from pathlib import Path
+
+sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))
+
+from tsi_sim.verify import main  # noqa: E402
+
+if __name__ == "__main__":
+    raise SystemExit(main())
diff --git a/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/__init__.py b/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/__init__.py
new file mode 100644
index 0000000..9816edf
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/__init__.py
@@ -0,0 +1,3 @@
+"""Cryptarchia Total Stake Inference simulator (uncle references)."""
+
+__version__ = "0.1.0"
diff --git a/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/blocktree.py b/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/blocktree.py
new file mode 100644
index 0000000..1f807d2
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/blocktree.py
@@ -0,0 +1,127 @@
+"""Block tree, latency-driven forks, and honest longest-chain fork choice.
+
+Blocks are stored in parallel arrays (id == index). A virtual genesis is block 0 at
+slot -1, height 0. Every real block is produced at an active slot by one winning node
+and points at the best tip *visible to that node at production time*, which is what makes
+network latency (and same-slot multi-winners) produce forks.
+
+Fork choice is honest longest-chain with a first-seen tie-break (prefer higher height,
+then earlier slot, then lower id) — no adversary is modelled, so the spec's density /
+deep-fork rules never engage.
+"""
+
+from __future__ import annotations
+
+import heapq
+from dataclasses import dataclass
+
+import numpy as np
+
+from .latency import LatencyModel
+
+GENESIS = 0
+
+
+@dataclass
+class BlockTree:
+    slot: np.ndarray          # int64, slot of each block (genesis = -1)
+    parent: np.ndarray        # int64, parent id (genesis = -1)
+    height: np.ndarray        # int64, chain height (genesis = 0)
+    leader: np.ndarray        # int64, producing node id (genesis = -1)
+    uncles: list[tuple[int, ...]]   # referenced uncle ids per block (filled later)
+
+    @property
+    def n_blocks(self) -> int:
+        return self.slot.shape[0]
+
+    def ancestors(self, block_id: int) -> list[int]:
+        """Ancestor chain of ``block_id`` from itself down to (excluding) genesis."""
+        out: list[int] = []
+        b = block_id
+        while b > GENESIS:
+            out.append(b)
+            b = int(self.parent[b])
+        return out
+
+    def canonical_chain(self) -> list[int]:
+        """Honest longest-chain: ancestors of the best tip over the whole tree.
+
+        Returns real block ids (genesis excluded), tip-first.
+        """
+        tip = self._best_over_all()
+        return self.ancestors(tip)
+
+    def _rank(self, bid: int) -> tuple[int, int, int]:
+        # Preference order for "better tip": higher height, earlier slot, lower id.
+        return (int(self.height[bid]), -int(self.slot[bid]), -bid)
+
+    def _best_over_all(self) -> int:
+        best = GENESIS
+        best_rank = self._rank(GENESIS)
+        for bid in range(1, self.n_blocks):
+            r = self._rank(bid)
+            if r > best_rank:
+                best_rank, best = r, bid
+        return best
+
+
+def build_tree(
+    active_slots: np.ndarray,
+    winners_per_slot: list[np.ndarray],
+    latency: LatencyModel,
+    rng: np.random.Generator,
+) -> BlockTree:
+    """Construct the block tree from grouped lottery winners under a latency model."""
+    # Preallocate with genesis in slot 0.
+    slot = [-1]
+    parent = [-1]
+    height = [0]
+    leader = [-1]
+
+    # global_best = best publicly-visible tip so far, as (height, slot, id).
+    def better(a: tuple[int, int, int], b: tuple[int, int, int]) -> tuple[int, int, int]:
+        # higher height, then earlier slot, then lower id
+        ah, as_, ai = a
+        bh, bs, bi = b
+        if ah != bh:
+            return a if ah > bh else b
+        if as_ != bs:
+            return a if as_ < bs else b
+        return a if ai < bi else b
+
+    global_best = (0, -1, GENESIS)
+    own_best: dict[int, tuple[int, int, int]] = {}
+    # min-heap of (visible_at, block_id) awaiting public visibility
+    pending: list[tuple[int, int]] = []
+
+    next_id = 1
+    for si in range(active_slots.shape[0]):
+        t = int(active_slots[si])
+        # advance visibility frontier to slot t
+        while pending and pending[0][0] <= t:
+            _, bid = heapq.heappop(pending)
+            cand = (height[bid], slot[bid], bid)
+            global_best = better(global_best, cand)
+        for v in winners_per_slot[si].tolist():
+            gb = global_best
+            ob = own_best.get(v, (0, -1, GENESIS))
+            chosen = better(gb, ob)
+            p_id = chosen[2]
+            h = chosen[0] + 1
+            bid = next_id
+            next_id += 1
+            slot.append(t)
+            parent.append(p_id)
+            height.append(h)
+            leader.append(v)
+            own_best[v] = (h, t, bid)
+            va = latency.visible_at(t, rng)
+            heapq.heappush(pending, (va, bid))
+
+    return BlockTree(
+        slot=np.asarray(slot, np.int64),
+        parent=np.asarray(parent, np.int64),
+        height=np.asarray(height, np.int64),
+        leader=np.asarray(leader, np.int64),
+        uncles=[() for _ in range(next_id)],
+    )
diff --git a/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/concurrency.py b/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/concurrency.py
new file mode 100644
index 0000000..1937dc8
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/concurrency.py
@@ -0,0 +1,62 @@
+"""Concurrent block-proposal analysis.
+
+Every lottery win is a block *proposal*. Two proposals produced within ``L`` slots of each
+other cannot see one another (a block becomes visible only after the network latency ``L``),
+so they are mutually concurrent — competing forks. Bucketing the timeline into
+non-overlapping windows of ``L`` slots and counting proposals per bucket gives a direct
+view of how many proposals are concurrent, and the busiest bucket is the peak number of
+concurrent proposals. This is the quantity that bounds how many uncles can appear, so it
+informs the ``MAX_UNCLES`` choice.
+"""
+
+from __future__ import annotations
+
+from dataclasses import replace
+
+import numpy as np
+
+from .config import SimConfig
+from .lottery import sample_wins, win_probs
+from .rng import seedseq_for
+from .stake import make_stake
+
+
+def window_counts(winner_slots: np.ndarray, epoch_len: int, bucket: int) -> np.ndarray:
+    """Proposals per non-overlapping ``bucket``-slot window over ``[0, epoch_len)``."""
+    if epoch_len <= 0:
+        return np.empty(0, np.int64)
+    bucket = max(int(bucket), 1)
+    n_windows = (epoch_len + bucket - 1) // bucket
+    if winner_slots.size == 0:
+        return np.zeros(n_windows, np.int64)
+    return np.bincount(winner_slots // bucket, minlength=n_windows)
+
+
+def proposal_slots(config: SimConfig, replicate: int = 0) -> np.ndarray:
+    """Simulate one epoch of block proposals at the *true* lottery difficulty (D=D_true).
+
+    Returns the sorted slots at which proposals (all lottery winners, including forks) occur.
+    Independent of the TSI trajectory — the proposal process depends only on stake, ``f``,
+    and ``epoch_len`` — so this is a cheap, self-contained re-simulation for the plots.
+    """
+    cfg = replace(config, replicate=replicate)
+    root = seedseq_for(cfg)
+    children = root.spawn(2)
+    stake = make_stake(cfg, np.random.default_rng(children[0]))
+    p_win = win_probs(stake, float(stake.sum()), cfg.f)
+    winner_slots, _ = sample_wins(p_win, cfg.epoch_len, np.random.default_rng(children[1]))
+    return winner_slots
+
+
+def concurrency_stats(config: SimConfig, replicate: int = 0) -> dict:
+    """Per-bucket proposal-count stats for one simulated epoch (bucket = ``max(L, 1)``)."""
+    ws = proposal_slots(config, replicate)
+    bucket = max(config.latency, 1)
+    counts = window_counts(ws, config.epoch_len, bucket)
+    return {
+        "bucket": bucket,
+        "counts": counts,
+        "max": int(counts.max()) if counts.size else 0,
+        "mean": float(counts.mean()) if counts.size else 0.0,
+        "p99": float(np.percentile(counts, 99)) if counts.size else 0.0,
+    }
diff --git a/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/config.py b/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/config.py
new file mode 100644
index 0000000..592e19e
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/config.py
@@ -0,0 +1,163 @@
+"""Configuration dataclasses for single runs and parameter sweeps."""
+
+from __future__ import annotations
+
+import itertools
+from dataclasses import dataclass, field, replace
+from typing import Any, Literal
+
+from . import constants
+
+StakeDist = Literal["uniform", "pareto"]
+UncleStrategy = Literal["oldest", "random"]
+
+
+@dataclass(frozen=True)
+class SimConfig:
+    """A single fully-specified simulation run (one grid cell, one replicate)."""
+
+    # --- network / stake ---
+    n_nodes: int = 1000
+    stake_dist: StakeDist = "uniform"
+    pareto_shape: float = 1.16          # Pareto (Lomax) tail index; ~80/20 by default
+    uniform_random: bool = False        # if True, draw i.i.d. uniform stakes; else equal
+    total_stake: float = 1.0e9          # FIXED across distributions for comparability
+
+    # --- network latency (slots) ---
+    latency: int = 0                    # L: block visible to others at t + L
+    latency_stochastic: bool = False    # if True, L is the mean of a stochastic model
+
+    # --- uncle references ---
+    uncle_window: int = constants.W_DEFAULT   # W
+    max_uncles: int = 0                       # U (0 = baseline, no uncles)
+    uncle_strategy: UncleStrategy = "oldest"
+    # Coin-flip inclusion prob for the "random" strategy. Only 0.5 reproduces the spec's
+    # unbiased coin (cryptarchia-v1-protocol.md); other values are a deliberate, non-spec
+    # sensitivity knob, not protocol behaviour.
+    uncle_random_p: float = 0.5
+
+    # --- consensus / TSI ---
+    f: float = constants.F
+    beta: float = constants.BETA_DEFAULT
+    k: int = 64                          # scaled by default; full scale = 2160
+    genesis_d_factor: float = 0.5        # genesis D = factor * true total stake
+    epochs: int = 40
+    # If True, mirror the spec's integer fixed-point f-truncation (f_p = int(f*1000)/1000),
+    # which the on-chain estimator uses; this reproduces its ~1% systematic overestimate.
+    # Default False keeps the analysis-faithful exact-f behaviour.
+    fixed_point: bool = False
+    per_node_dest: bool = False          # Phase-2 hook: per-node D_est (unused in reduced model)
+
+    # --- performance ---
+    # >1 parallelises the per-slot lottery across slot-chunks (opt-in; must be pinned and
+    # recorded because it changes the RNG stream — see lottery.sample_wins_chunked).
+    lottery_chunks: int = 1
+
+    # --- bookkeeping ---
+    replicate: int = 0
+    root_seed: int = 12345
+
+    def __post_init__(self) -> None:
+        # frozen dataclass: validation only (no attribute assignment)
+        if self.stake_dist not in ("uniform", "pareto"):
+            raise ValueError(f"stake_dist must be uniform|pareto, got {self.stake_dist!r}")
+        if self.uncle_strategy not in ("oldest", "random"):
+            raise ValueError(f"uncle_strategy must be oldest|random, got {self.uncle_strategy!r}")
+        checks = {
+            "n_nodes": self.n_nodes >= 1,
+            "k": self.k >= 1,
+            "epochs": self.epochs >= 1,
+            "latency": self.latency >= 0,
+            "max_uncles": self.max_uncles >= 0,
+            "uncle_window": self.uncle_window >= 1,
+            "lottery_chunks": self.lottery_chunks >= 1,
+            "uncle_random_p": 0.0 <= self.uncle_random_p <= 1.0,
+            "f": 0.0 < self.f < 1.0,
+            "beta": self.beta > 0.0,
+            "genesis_d_factor": self.genesis_d_factor > 0.0,
+            "pareto_shape": self.pareto_shape > 0.0,
+            "total_stake": self.total_stake > 0.0,
+        }
+        bad = [name for name, ok in checks.items() if not ok]
+        if bad:
+            raise ValueError(f"invalid SimConfig field(s): {bad}")
+
+    # derived geometry -------------------------------------------------------
+    @property
+    def epoch_len(self) -> int:
+        return constants.epoch_len(self.k, self.f)
+
+    @property
+    def period_T(self) -> int:
+        return constants.period_T(self.k, self.f)
+
+    def key(self) -> tuple:
+        """Hashable identity used to seed the RNG deterministically.
+
+        Must include EVERY field that affects the run (guarded by test_rng), otherwise two
+        distinct configs would share an RNG stream.
+        """
+        return (
+            self.n_nodes, self.stake_dist, self.pareto_shape, self.uniform_random,
+            self.total_stake, self.latency, self.latency_stochastic, self.uncle_window,
+            self.max_uncles, self.uncle_strategy, self.uncle_random_p, self.f, self.beta,
+            self.k, self.genesis_d_factor, self.epochs, self.fixed_point, self.per_node_dest,
+            self.lottery_chunks, self.replicate,
+        )
+
+
+@dataclass
+class SweepConfig:
+    """A cartesian grid of runs plus replicates, all sharing ``base`` settings."""
+
+    n_nodes: list[int] = field(default_factory=lambda: [1000])
+    stake_dist: list[StakeDist] = field(default_factory=lambda: ["uniform", "pareto"])
+    latency: list[int] = field(default_factory=lambda: [0, 1, 2, 4, 8])
+    max_uncles: list[int] = field(default_factory=lambda: [0, 1, 2, 3, 4])
+    uncle_strategy: list[UncleStrategy] = field(default_factory=lambda: ["oldest", "random"])
+    f: list[float] = field(default_factory=lambda: [constants.F])
+    replicates: int = 8
+    base: dict[str, Any] = field(default_factory=dict)
+
+    def expand(self) -> list[SimConfig]:
+        """Materialise every ``SimConfig`` in the grid × replicates."""
+        base = SimConfig(**self.base)
+        cells: list[SimConfig] = []
+        axes = itertools.product(
+            self.n_nodes, self.stake_dist, self.latency, self.max_uncles,
+            self.uncle_strategy, self.f,
+        )
+        for n, dist, lat, u, strat, fval in axes:
+            # U=0 is strategy-independent; keep only one strategy to avoid duplicate work.
+            if u == 0 and strat != self.uncle_strategy[0]:
+                continue
+            for rep in range(self.replicates):
+                cells.append(
+                    replace(
+                        base,
+                        n_nodes=n,
+                        stake_dist=dist,
+                        latency=lat,
+                        max_uncles=u,
+                        uncle_strategy=strat,
+                        f=fval,
+                        replicate=rep,
+                    )
+                )
+        return cells
+
+    @classmethod
+    def from_dict(cls, d: dict[str, Any]) -> SweepConfig:
+        d = dict(d)
+        base = d.pop("base", {})
+        known = {
+            "n_nodes", "stake_dist", "latency", "max_uncles",
+            "uncle_strategy", "f", "replicates",
+        }
+        unknown = set(d) - known
+        if unknown:
+            raise ValueError(
+                f"unknown sweep keys: {sorted(unknown)} (did you mean one of {sorted(known)}? "
+                "per-run settings belong under 'base:')"
+            )
+        return cls(base=base, **d)
diff --git a/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/constants.py b/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/constants.py
new file mode 100644
index 0000000..baba55b
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/constants.py
@@ -0,0 +1,39 @@
+"""Protocol constants and epoch/window geometry.
+
+All slot geometry derives from the pair ``(k, f)`` so a scaled-down ``k`` (used for
+parameter sweeps) automatically shrinks the epoch and measurement window. See
+``cryptarchia-v1-protocol.md`` and ``cryptarchia-total-stake-inference.md``.
+"""
+
+from __future__ import annotations
+
+# --- True protocol values (full scale) -------------------------------------
+K_TRUE = 2160          # security parameter (blocks)
+F = 1.0 / 30.0         # slot activation coefficient (fixed; never tuned)
+W_DEFAULT = 300        # uncle reference window w_u (slots)
+BETA_DEFAULT = 1.0     # TSI learning rate
+SLOT_SECONDS = 1       # slot length
+
+
+def floor_k_over_f(k: int, f: float = F) -> int:
+    """``floor(k / f)`` — the base quantum of the epoch schedule."""
+    return int(k / f)
+
+
+def epoch_len(k: int, f: float = F) -> int:
+    """Epoch length in slots: ``10 * floor(k/f)``."""
+    return 10 * floor_k_over_f(k, f)
+
+
+def period_T(k: int, f: float = F) -> int:
+    """TSI measurement window length ``T`` in slots: ``6 * floor(k/f)``.
+
+    This is the first ``6*floor(k/f)`` slots of the (previous) epoch over which the
+    block density is measured.
+    """
+    return 6 * floor_k_over_f(k, f)
+
+
+def expected_blocks_in_window(k: int, f: float = F) -> float:
+    """Expected honest-chain block count in the measurement window at equilibrium."""
+    return period_T(k, f) * f
diff --git a/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/engine.py b/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/engine.py
new file mode 100644
index 0000000..d08e34b
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/engine.py
@@ -0,0 +1,39 @@
+"""Multi-epoch trajectory driver for a single config."""
+
+from __future__ import annotations
+
+from typing import Any
+
+import numpy as np
+
+from . import tsi
+from .config import SimConfig
+from .epoch import simulate_epoch
+from .metrics import metric_row
+from .rng import seedseq_for
+from .stake import make_stake
+
+
+def run_trajectory(config: SimConfig) -> list[dict[str, Any]]:
+    """Run ``config.epochs`` epochs of TSI, returning one metric row per epoch.
+
+    Stake is drawn once (it is fixed; only ``D_est`` evolves). ``D_est`` starts at the
+    hardcoded genesis value ``genesis_d_factor * D_true`` and is updated each epoch from
+    the measured density. The RNG is a spawn hierarchy off the config's root SeedSequence:
+    child 0 draws the stake, child ``e+1`` drives epoch ``e`` — so every draw is a
+    deterministic, order-independent function of the config identity.
+    """
+    root = seedseq_for(config)
+    children = root.spawn(config.epochs + 1)
+    stake = make_stake(config, np.random.default_rng(children[0]))
+    d_true = float(stake.sum())
+    d_est = config.genesis_d_factor * d_true
+    T = config.period_T
+
+    rows: list[dict[str, Any]] = []
+    for epoch in range(config.epochs):
+        er = simulate_epoch(config, stake, d_est, children[epoch + 1])
+        d_next = tsi.update_D(d_est, er.m, T, config.f, config.beta, config.fixed_point)
+        rows.append(metric_row(config, epoch, d_est, d_next, d_true, er))
+        d_est = d_next
+    return rows
diff --git a/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/epoch.py b/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/epoch.py
new file mode 100644
index 0000000..37c97a4
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/epoch.py
@@ -0,0 +1,84 @@
+"""Single-epoch simulation: lottery -> block tree -> uncles -> density counting."""
+
+from __future__ import annotations
+
+from dataclasses import dataclass
+
+import numpy as np
+
+from . import concurrency, lottery, tsi
+from .blocktree import build_tree
+from .config import SimConfig
+from .latency import make_latency
+
+
+@dataclass
+class EpochResult:
+    m: int                 # TSI block count in window
+    q: float               # honest active-slot fraction (window)
+    q_eff: float           # uncle-recovered active-slot fraction (window)
+    n_active: int          # active slots in window
+    n_honest: int          # honest slots in window
+    n_recovered: int       # orphan slots recovered by uncles in window
+    total_winners_window: int   # total lottery wins in window (incl. multi-winner)
+    n_blocks: int          # real blocks produced this epoch
+    n_canonical: int       # canonical chain length
+    n_orphans: int         # orphaned blocks
+    max_concurrent: int    # most block proposals in any latency-sized (max(L,1)) bucket
+    mean_concurrent: float # mean proposals per latency-sized bucket
+
+
+def simulate_epoch(
+    config: SimConfig, stake: np.ndarray, d_est: float, epoch_ss: np.random.SeedSequence
+) -> EpochResult:
+    f = config.f
+    T = config.period_T
+    # independent sub-streams: one for the lottery, one for the tree/uncle auxiliary draws
+    lottery_ss, aux_ss = epoch_ss.spawn(2)
+    aux_rng = np.random.default_rng(aux_ss)
+
+    p_win = lottery.win_probs(stake, d_est, f)
+    if config.lottery_chunks > 1:
+        winner_slots, winner_nodes = lottery.sample_wins_chunked(
+            p_win, config.epoch_len, lottery_ss, config.lottery_chunks
+        )
+    else:
+        winner_slots, winner_nodes = lottery.sample_wins(
+            p_win, config.epoch_len, np.random.default_rng(lottery_ss)
+        )
+    active_slots, groups = lottery.group_by_slot(winner_slots, winner_nodes)
+
+    latency = make_latency(config)
+    tree = build_tree(active_slots, groups, latency, aux_rng)
+    canonical = tree.canonical_chain()
+    from .uncles import annotate_uncles
+
+    annotate_uncles(tree, canonical, config, aux_rng)
+
+    m = tsi.density_m(tree, canonical, T)
+    ref = tsi.referenced_uncle_ids(tree, canonical)
+    ss = tsi.slot_stats(tree, canonical, ref, active_slots, T)
+
+    total_winners_window = int((winner_slots < T).sum())
+    n_real = tree.n_blocks - 1
+    # Concurrent proposals: bucket the whole epoch into latency-sized windows and count
+    # block proposals (every winner is a proposal) per bucket. The max bucket is the peak
+    # number of mutually-concurrent proposals (they cannot see each other within L slots).
+    bucket = max(config.latency, 1)
+    counts = concurrency.window_counts(winner_slots, config.epoch_len, bucket)
+    max_concurrent = int(counts.max()) if counts.size else 0
+    mean_concurrent = float(counts.mean()) if counts.size else 0.0
+    return EpochResult(
+        m=m,
+        q=ss.q,
+        q_eff=ss.q_eff,
+        n_active=ss.n_active,
+        n_honest=ss.n_honest,
+        n_recovered=ss.n_recovered,
+        total_winners_window=total_winners_window,
+        n_blocks=n_real,
+        n_canonical=len(canonical),
+        n_orphans=n_real - len(canonical),
+        max_concurrent=max_concurrent,
+        mean_concurrent=mean_concurrent,
+    )
diff --git a/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/latency.py b/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/latency.py
new file mode 100644
index 0000000..213c7c1
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/latency.py
@@ -0,0 +1,54 @@
+"""Network latency models.
+
+Latency ``L`` is the number of slots between a block being produced and it becoming
+visible to the rest of the network. ``L`` is deliberately named to avoid clashing with
+``D`` (the stake estimate). A leader at slot ``t`` can only build on blocks whose
+``visible_at <= t`` (its own block is visible to itself immediately), which is what
+produces latency-induced forks.
+"""
+
+from __future__ import annotations
+
+from typing import Protocol
+
+import numpy as np
+
+
+class LatencyModel(Protocol):
+    def visible_at(self, produced_slot: int, rng: np.random.Generator) -> int:
+        """Slot at which a block produced at ``produced_slot`` becomes visible to others."""
+        ...
+
+
+class FixedSlotLatency:
+    """Deterministic integer-slot latency: visible to all others at ``t + L``."""
+
+    def __init__(self, latency: int) -> None:
+        self.latency = int(latency)
+
+    def visible_at(self, produced_slot: int, rng: np.random.Generator) -> int:
+        return produced_slot + self.latency
+
+
+class RealisticLatency:
+    """Stochastic latency with mean ``L`` slots (optional sensitivity model).
+
+    Rounds an exponential draw (mean ``L``) up to whole slots. A stand-in for the
+    reference notebook's blend/broadcast delay model; not used by the primary sweep.
+    """
+
+    def __init__(self, mean_latency: float) -> None:
+        self.mean_latency = float(mean_latency)
+
+    def visible_at(self, produced_slot: int, rng: np.random.Generator) -> int:
+        if self.mean_latency <= 0:
+            return produced_slot
+        draw = rng.exponential(self.mean_latency)
+        return produced_slot + int(np.ceil(draw))
+
+
+def make_latency(config) -> LatencyModel:  # noqa: ANN001 - avoid import cycle with config
+    """Build the latency model for a config."""
+    if config.latency_stochastic:
+        return RealisticLatency(config.latency)
+    return FixedSlotLatency(config.latency)
diff --git a/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/lottery.py b/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/lottery.py
new file mode 100644
index 0000000..a68641b
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/lottery.py
@@ -0,0 +1,114 @@
+"""Stake-weighted slot lottery (sparse sampler).
+
+Per node ``i`` and slot, an independent Bernoulli win with probability
+``phi_f(alpha_i) = 1 - (1 - f)^alpha_i`` where ``alpha_i = w_i / D_est``. Multiple winners
+in a slot are possible (a guaranteed fork).
+
+The number of slots a node wins is exactly ``Binomial(n_slots, p_i)``, and the won slots
+are a uniformly-random distinct subset — this is *distributionally identical* to drawing an
+independent Bernoulli(p_i) in every slot, but avoids materialising the dense
+``(n_nodes, n_slots)`` array (which was ~95% of the whole simulator's runtime). Winners are
+returned as sparse ``(winner_slots, winner_nodes)`` coordinates, sorted by slot.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+from joblib import Parallel, delayed
+
+
+def phi(f: float, alpha: np.ndarray | float) -> np.ndarray | float:
+    """Leader-lottery win probability ``1 - (1 - f)^alpha``."""
+    return 1.0 - (1.0 - f) ** alpha
+
+
+def win_probs(stake: np.ndarray, d_est: float, f: float) -> np.ndarray:
+    """Per-node win probability ``phi_f(w_i / D_est)``."""
+    return phi(f, stake / d_est)
+
+
+def _winners_from_counts(
+    counts: np.ndarray, offset: int, span: int, rng: np.random.Generator
+) -> tuple[list[np.ndarray], list[np.ndarray]]:
+    """For each winning node, sample ``counts[i]`` distinct slots in ``[offset, offset+span)``."""
+    nz = np.nonzero(counts)[0]
+    slot_parts: list[np.ndarray] = []
+    node_parts: list[np.ndarray] = []
+    for i in nz:
+        c = int(counts[i])
+        slots_i = rng.choice(span, size=c, replace=False).astype(np.int64) + offset
+        slot_parts.append(slots_i)
+        node_parts.append(np.full(c, i, np.int64))
+    return slot_parts, node_parts
+
+
+def _finalize(
+    slot_parts: list[np.ndarray], node_parts: list[np.ndarray]
+) -> tuple[np.ndarray, np.ndarray]:
+    if not slot_parts:
+        return np.empty(0, np.int64), np.empty(0, np.int64)
+    winner_slots = np.concatenate(slot_parts)
+    winner_nodes = np.concatenate(node_parts)
+    order = np.argsort(winner_slots, kind="stable")
+    return winner_slots[order], winner_nodes[order]
+
+
+def sample_wins(
+    p_win: np.ndarray, n_slots: int, rng: np.random.Generator, chunk: int = 8192
+) -> tuple[np.ndarray, np.ndarray]:
+    """Sample lottery wins over ``n_slots`` slots (sparse; ``chunk`` kept for API compat)."""
+    counts = rng.binomial(n_slots, p_win)
+    slot_parts, node_parts = _winners_from_counts(counts, 0, n_slots, rng)
+    return _finalize(slot_parts, node_parts)
+
+
+def sample_wins_chunked(
+    p_win: np.ndarray,
+    n_slots: int,
+    seedseq: np.random.SeedSequence,
+    n_chunks: int,
+    n_jobs: int = -1,
+) -> tuple[np.ndarray, np.ndarray]:
+    """Parallel sparse lottery: partition slots into ``n_chunks`` independent ranges.
+
+    Correct because ``Binomial(n_slots, p) = sum_c Binomial(L_c, p)`` and per-chunk distinct
+    subsets are independent. Deterministic given ``(seedseq, n_chunks)`` — but the exact
+    winner identities differ from the serial sampler and *change with* ``n_chunks``, so
+    ``n_chunks`` must be a pinned, recorded config parameter, never derived from core count.
+
+    Note: after the sparse rewrite the lottery is a small fraction of an epoch, so this has
+    little ROI versus across-config parallelism; it exists for the rare isolated config with
+    an enormous ``n_slots`` and no across-config work to fill cores.
+    """
+    if n_chunks <= 1:
+        return sample_wins(p_win, n_slots, np.random.default_rng(seedseq))
+    bounds = np.linspace(0, n_slots, n_chunks + 1).astype(np.int64)
+    children = seedseq.spawn(n_chunks)
+
+    def one_chunk(c: int) -> tuple[list[np.ndarray], list[np.ndarray]]:
+        lo, hi = int(bounds[c]), int(bounds[c + 1])
+        span = hi - lo
+        rng = np.random.default_rng(children[c])
+        counts = rng.binomial(span, p_win)
+        return _winners_from_counts(counts, lo, span, rng)
+
+    results = Parallel(n_jobs=n_jobs, prefer="threads")(
+        delayed(one_chunk)(c) for c in range(n_chunks)
+    )
+    slot_parts: list[np.ndarray] = []
+    node_parts: list[np.ndarray] = []
+    for sp, npar in results:
+        slot_parts.extend(sp)
+        node_parts.extend(npar)
+    return _finalize(slot_parts, node_parts)
+
+
+def group_by_slot(
+    winner_slots: np.ndarray, winner_nodes: np.ndarray
+) -> tuple[np.ndarray, list[np.ndarray]]:
+    """Group sorted winner coordinates into ``(active_slots, winners_per_active_slot)``."""
+    if winner_slots.size == 0:
+        return np.empty(0, np.int64), []
+    active_slots, starts = np.unique(winner_slots, return_index=True)
+    groups = np.split(winner_nodes, starts[1:])
+    return active_slots, groups
diff --git a/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/metrics.py b/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/metrics.py
new file mode 100644
index 0000000..0fe9b7f
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/metrics.py
@@ -0,0 +1,76 @@
+"""Per-epoch metric rows and equilibrium summaries."""
+
+from __future__ import annotations
+
+from typing import Any
+
+import numpy as np
+
+from .config import SimConfig
+from .epoch import EpochResult
+
+# Config fields recorded on every row for grouping/plotting.
+_CONFIG_FIELDS = (
+    "n_nodes", "stake_dist", "pareto_shape", "latency", "uncle_window",
+    "max_uncles", "uncle_strategy", "uncle_random_p", "f", "beta", "k",
+    "genesis_d_factor", "epochs", "fixed_point", "replicate",
+)
+
+
+def metric_row(
+    config: SimConfig,
+    epoch: int,
+    d_in: float,
+    d_out: float,
+    d_true: float,
+    er: EpochResult,
+) -> dict[str, Any]:
+    row: dict[str, Any] = {field: getattr(config, field) for field in _CONFIG_FIELDS}
+    row.update(
+        epoch=epoch,
+        d_in=d_in,
+        d_out=d_out,
+        d_true=d_true,
+        ratio=d_out / d_true,
+        m=er.m,
+        measured_density=er.m / config.period_T,
+        q=er.q,
+        q_eff=er.q_eff,
+        n_active=er.n_active,
+        n_honest=er.n_honest,
+        n_recovered=er.n_recovered,
+        total_winners_window=er.total_winners_window,
+        n_blocks=er.n_blocks,
+        n_canonical=er.n_canonical,
+        n_orphans=er.n_orphans,
+        orphan_rate=er.n_orphans / er.n_blocks if er.n_blocks else float("nan"),
+        max_concurrent=er.max_concurrent,
+        mean_concurrent=er.mean_concurrent,
+    )
+    return row
+
+
+def equilibrium_stats(ratios: np.ndarray, burn_in: int) -> dict[str, float]:
+    """Mean/variance of the stake ratio after ``burn_in`` epochs."""
+    ratios = np.asarray(ratios, dtype=float)
+    if ratios.size == 0:
+        return {"mean_ratio": float("nan"), "var_ratio": float("nan"), "std_ratio": float("nan")}
+    tail = ratios[burn_in:]
+    if tail.size == 0:
+        tail = ratios[-1:]
+    return {
+        "mean_ratio": float(np.mean(tail)),
+        "var_ratio": float(np.var(tail)),
+        "std_ratio": float(np.std(tail)),
+    }
+
+
+def epochs_to_within(ratios: np.ndarray, target: float, eps: float) -> int:
+    """First epoch index after which ``|ratio - target| <= eps`` holds for the rest."""
+    within = np.abs(np.asarray(ratios, dtype=float) - target) <= eps
+    n = within.size
+    if n == 0:
+        return 0
+    # first index i such that within[i:] are all True == 1 + last index that is False
+    false_idx = np.flatnonzero(~within)
+    return int(false_idx[-1] + 1) if false_idx.size else 0
diff --git a/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/plotting/__init__.py b/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/plotting/__init__.py
new file mode 100644
index 0000000..527b68e
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/plotting/__init__.py
@@ -0,0 +1 @@
+"""Academic-quality figure generation."""
diff --git a/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/plotting/figures.py b/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/plotting/figures.py
new file mode 100644
index 0000000..baced35
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/plotting/figures.py
@@ -0,0 +1,339 @@
+"""Figure builders. Each takes the results DataFrame and returns a Matplotlib Figure.
+
+Rows are per (config, epoch). We summarise each config to its *equilibrium* by averaging
+the tail epochs (after a burn-in), keeping replicates so we can draw percentile bands.
+"""
+
+from __future__ import annotations
+
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+
+from .. import concurrency, theory
+from ..config import SimConfig
+from . import style
+
+CONFIG_COLS = ["n_nodes", "stake_dist", "latency", "max_uncles", "uncle_strategy", "k"]
+
+
+def _ceiling_line(ax, f: float) -> None:
+    """Overlay the intrinsic block-count overshoot ceiling ``-ln(1-f)/f`` (~1.017)."""
+    c = float(theory.block_count_ceiling(f))
+    ax.axhline(c, color="0.55", lw=0.9, ls=":", zorder=0)
+    ax.text(0.99, c, f" full-recovery ceiling {c:.3f}", transform=ax.get_yaxis_transform(),
+            ha="right", va="bottom", fontsize=7, color="0.4")
+
+
+def equilibrium(df: pd.DataFrame, burn_frac: float = 0.5) -> pd.DataFrame:
+    """Per-(config, replicate) equilibrium means over the tail epochs."""
+    cutoff = df["epochs"] * burn_frac
+    tail = df[df["epoch"] >= cutoff]
+    keys = [*CONFIG_COLS, "replicate"]
+    return (
+        tail.groupby(keys, as_index=False)
+        .agg(
+            ratio=("ratio", "mean"),
+            q=("q", "mean"),
+            q_eff=("q_eff", "mean"),
+            var_ratio=("ratio", "var"),
+            orphan_rate=("orphan_rate", "mean"),
+            n_active=("n_active", "mean"),
+        )
+    )
+
+
+def _series_over_replicates(eq: pd.DataFrame, xcol: str, ycol: str, xvals) -> np.ndarray:
+    """Build a ``(n_replicates, len(xvals))`` matrix of ``ycol`` for band plots."""
+    reps = sorted(eq["replicate"].unique())
+    mat = np.full((len(reps), len(xvals)), np.nan)
+    for i, r in enumerate(reps):
+        sub = eq[eq["replicate"] == r].set_index(xcol)[ycol]
+        for j, x in enumerate(xvals):
+            if x in sub.index:
+                mat[i, j] = sub.loc[x]
+    return mat
+
+
+def _provenance(df: pd.DataFrame) -> str:
+    k = [int(x) for x in sorted(df["k"].unique())]
+    n = [int(x) for x in sorted(df["n_nodes"].unique())]
+    reps = int(df["replicate"].nunique())
+    return f"tsi-sim  |  k={k}  N={n}  reps={reps}  f={float(df['f'].iloc[0]):.4g}"
+
+
+# --- Figure 1: accuracy vs U per latency ------------------------------------
+def accuracy_vs_u(df: pd.DataFrame, stake_dist: str, strategy: str = "oldest") -> plt.Figure:
+    style.apply_style()
+    eq = equilibrium(df)
+    eq = eq[(eq["stake_dist"] == stake_dist) & (eq["n_nodes"] == df["n_nodes"].max())]
+    eq = eq[(eq["uncle_strategy"] == strategy) | (eq["max_uncles"] == 0)]
+    latencies = sorted(eq["latency"].unique())
+    uvals = sorted(eq["max_uncles"].unique())
+
+    fig, ax = plt.subplots()
+    for i, lat in enumerate(latencies):
+        sub = eq[eq["latency"] == lat]
+        mat = _series_over_replicates(sub, "max_uncles", "ratio", uvals)
+        style.band_plot(ax, uvals, mat, color=style.color_for(i), label=f"L={lat}")
+    ax.axhline(1.0, color="0.4", lw=1.0, ls="--", zorder=0)
+    _ceiling_line(ax, float(df["f"].iloc[0]))
+    ax.set_xlabel("max uncles per block  $U$")
+    ax.set_ylabel(r"inferred / true stake  $\langle \hat D / D_{\mathrm{true}} \rangle$")
+    ax.set_title(f"TSI accuracy vs uncle cap ({stake_dist} stake)")
+    ax.set_xticks(uvals)
+    ax.legend(title="latency (slots)", ncol=2)
+    return fig
+
+
+# --- Figure 2: q_eff vs U per latency ---------------------------------------
+def qeff_vs_u(df: pd.DataFrame, stake_dist: str, strategy: str = "oldest") -> plt.Figure:
+    style.apply_style()
+    eq = equilibrium(df)
+    eq = eq[(eq["stake_dist"] == stake_dist) & (eq["n_nodes"] == df["n_nodes"].max())]
+    eq = eq[(eq["uncle_strategy"] == strategy) | (eq["max_uncles"] == 0)]
+    latencies = sorted(eq["latency"].unique())
+    uvals = sorted(eq["max_uncles"].unique())
+
+    fig, ax = plt.subplots()
+    for i, lat in enumerate(latencies):
+        sub = eq[eq["latency"] == lat]
+        mat = _series_over_replicates(sub, "max_uncles", "q_eff", uvals)
+        style.band_plot(ax, uvals, mat, color=style.color_for(i), label=f"L={lat}")
+        base_q = sub[sub["max_uncles"] == 0]["q"].mean()
+        if not np.isnan(base_q):
+            ax.axhline(base_q, color=style.color_for(i), lw=0.8, ls=":", alpha=0.6)
+    ax.axhline(1.0, color="0.4", lw=1.0, ls="--", zorder=0)
+    ax.set_xlabel("max uncles per block  $U$")
+    ax.set_ylabel(r"effective utilisation  $q_{\mathrm{eff}}$")
+    ax.set_title(f"Active-slot recovery vs uncle cap ({stake_dist} stake)")
+    ax.set_xticks(uvals)
+    ax.legend(title="latency (slots)", ncol=2)
+    return fig
+
+
+# --- Figure 3: convergence (ratio vs epoch) per U ---------------------------
+def convergence(df: pd.DataFrame, stake_dist: str, latency: int,
+                strategy: str = "oldest") -> plt.Figure:
+    style.apply_style()
+    sub = df[(df["stake_dist"] == stake_dist) & (df["latency"] == latency)
+             & (df["n_nodes"] == df["n_nodes"].max())]
+    sub = sub[(sub["uncle_strategy"] == strategy) | (sub["max_uncles"] == 0)]
+    uvals = sorted(sub["max_uncles"].unique())
+    epochs = sorted(sub["epoch"].unique())
+
+    fig, ax = plt.subplots()
+    for i, u in enumerate(uvals):
+        su = sub[sub["max_uncles"] == u]
+        mat = np.full((su["replicate"].nunique(), len(epochs)), np.nan)
+        for ri, r in enumerate(sorted(su["replicate"].unique())):
+            s = su[su["replicate"] == r].set_index("epoch")["ratio"]
+            for j, e in enumerate(epochs):
+                if e in s.index:
+                    mat[ri, j] = s.loc[e]
+        style.band_plot(ax, epochs, mat, color=style.color_for(i), label=f"U={u}")
+    ax.axhline(1.0, color="0.4", lw=1.0, ls="--", zorder=0)
+    if any(u > 0 for u in uvals):
+        _ceiling_line(ax, float(sub["f"].iloc[0]))
+    ax.set_xlabel("epoch")
+    ax.set_ylabel(r"$\hat D / D_{\mathrm{true}}$")
+    ax.set_title(f"Convergence ({stake_dist}, latency L={latency})")
+    ax.legend(title="uncle cap", ncol=2)
+    return fig
+
+
+# --- Figure 4: accuracy heatmap over latency x U ----------------------------
+def heatmap_accuracy(df: pd.DataFrame, stake_dist: str, n_nodes: int,
+                     strategy: str = "oldest") -> plt.Figure:
+    style.apply_style()
+    eq = equilibrium(df)
+    eq = eq[(eq["stake_dist"] == stake_dist) & (eq["n_nodes"] == n_nodes)]
+    eq = eq[(eq["uncle_strategy"] == strategy) | (eq["max_uncles"] == 0)]
+    piv = eq.groupby(["latency", "max_uncles"])["ratio"].mean().unstack("max_uncles")
+    latencies = piv.index.to_numpy()
+    uvals = piv.columns.to_numpy()
+    data = piv.to_numpy()
+
+    fig, ax = plt.subplots()
+    vmax = np.nanmax(np.abs(data - 1.0))
+    im = ax.imshow(data, origin="lower", aspect="auto", cmap=style.DIVERGING_CMAP,
+                   vmin=1 - vmax, vmax=1 + vmax)
+    ax.set_xticks(range(len(uvals)), uvals)
+    ax.set_yticks(range(len(latencies)), latencies)
+    ax.set_xlabel("max uncles per block  $U$")
+    ax.set_ylabel("network latency  $L$ (slots)")
+    ax.set_title(f"Accuracy $\\hat D/D_{{\\mathrm{{true}}}}$ ({stake_dist}, N={n_nodes})")
+    for yi in range(len(latencies)):
+        for xi in range(len(uvals)):
+            v = data[yi, xi]
+            if not np.isnan(v):
+                safe = 0.98 <= v <= 1.02
+                ax.text(xi, yi, f"{v:.2f}", ha="center", va="center", fontsize=7,
+                        color="black", fontweight="bold" if safe else "normal")
+    fig.colorbar(im, ax=ax, label=r"$\hat D / D_{\mathrm{true}}$")
+    ax.grid(False)
+    return fig
+
+
+# --- Figure 5: distribution comparison --------------------------------------
+def dist_comparison(df: pd.DataFrame, latency: int, strategy: str = "oldest") -> plt.Figure:
+    style.apply_style()
+    eq = equilibrium(df)
+    eq = eq[(eq["latency"] == latency) & (eq["n_nodes"] == df["n_nodes"].max())]
+    eq = eq[(eq["uncle_strategy"] == strategy) | (eq["max_uncles"] == 0)]
+    uvals = sorted(eq["max_uncles"].unique())
+    fig, ax = plt.subplots()
+    for i, dist in enumerate(sorted(eq["stake_dist"].unique())):
+        sub = eq[eq["stake_dist"] == dist]
+        mat = _series_over_replicates(sub, "max_uncles", "ratio", uvals)
+        style.band_plot(ax, uvals, mat, color=style.color_for(i), label=dist)
+    ax.axhline(1.0, color="0.4", lw=1.0, ls="--", zorder=0)
+    _ceiling_line(ax, float(df["f"].iloc[0]))
+    ax.set_xlabel("max uncles per block  $U$")
+    ax.set_ylabel(r"$\hat D / D_{\mathrm{true}}$")
+    ax.set_title(f"Stake-distribution comparison (latency L={latency})")
+    ax.set_xticks(uvals)
+    ax.legend(title="distribution")
+    return fig
+
+
+# --- Figure 6: variance vs U ------------------------------------------------
+def variance_vs_u(df: pd.DataFrame, stake_dist: str, latency: int,
+                  strategy: str = "oldest") -> plt.Figure:
+    style.apply_style()
+    eq = equilibrium(df)
+    eq = eq[(eq["stake_dist"] == stake_dist) & (eq["latency"] == latency)
+            & (eq["n_nodes"] == df["n_nodes"].max())]
+    eq = eq[(eq["uncle_strategy"] == strategy) | (eq["max_uncles"] == 0)]
+    uvals = sorted(eq["max_uncles"].unique())
+    # per-epoch tail variance (mean of each replicate's within-tail variance), matching the
+    # marginal Var[D/D_true] that theory.variance_ratio denotes -- NOT the variance across
+    # replicate means (which is ~T smaller).
+    var_by_u = eq.groupby("max_uncles")["var_ratio"].mean()
+    qeff_by_u = eq.groupby("max_uncles")["q_eff"].mean()
+    f = float(df["f"].iloc[0])
+    T = int(round(6 * float(df["k"].iloc[0]) / f))   # measurement window length
+    fig, ax = plt.subplots()
+    ax.plot(uvals, [var_by_u.get(u, np.nan) for u in uvals], "o-",
+            color=style.color_for(0), label="empirical per-epoch tail variance")
+    theo = [float(theory.variance_ratio(f, min(qeff_by_u.get(u, np.nan), 1.0), T)) for u in uvals]
+    ax.plot(uvals, theo, "s--", color=style.color_for(1), label=r"theory $(q_{\mathrm{eff}})$")
+    ax.set_xlabel("max uncles per block  $U$")
+    ax.set_ylabel(r"$\mathrm{Var}[\hat D / D_{\mathrm{true}}]$")
+    ax.set_title(f"Estimate variance vs uncle cap ({stake_dist}, L={latency})")
+    ax.set_xticks(uvals)
+    ax.legend()
+    return fig
+
+
+# --- Figure 7: strategy comparison ------------------------------------------
+def strategy_comparison(df: pd.DataFrame, stake_dist: str, latency: int) -> plt.Figure:
+    style.apply_style()
+    eq = equilibrium(df)
+    eq = eq[(eq["stake_dist"] == stake_dist) & (eq["latency"] == latency)
+            & (eq["n_nodes"] == df["n_nodes"].max()) & (eq["max_uncles"] > 0)]
+    uvals = sorted(eq["max_uncles"].unique())
+    fig, ax = plt.subplots()
+    for i, strat in enumerate(sorted(eq["uncle_strategy"].unique())):
+        sub = eq[eq["uncle_strategy"] == strat]
+        mat = _series_over_replicates(sub, "max_uncles", "ratio", uvals)
+        style.band_plot(ax, uvals, mat, color=style.color_for(i), label=strat)
+    ax.axhline(1.0, color="0.4", lw=1.0, ls="--", zorder=0)
+    ax.set_xlabel("max uncles per block  $U$")
+    ax.set_ylabel(r"$\hat D / D_{\mathrm{true}}$")
+    ax.set_title(f"Uncle-selection strategy ({stake_dist}, L={latency})")
+    ax.set_xticks(uvals)
+    ax.legend(title="strategy")
+    return fig
+
+
+# --- Figure 8: diagnostics (orphan rate vs latency) -------------------------
+def orphan_diagnostics(df: pd.DataFrame, stake_dist: str) -> plt.Figure:
+    style.apply_style()
+    eq = equilibrium(df)
+    eq = eq[(eq["stake_dist"] == stake_dist) & (eq["max_uncles"] == 0)
+            & (eq["n_nodes"] == df["n_nodes"].max())]
+    g = eq.groupby("latency")["orphan_rate"].agg(["mean", "std"])
+    fig, ax = plt.subplots()
+    ax.errorbar(g.index, g["mean"], yerr=g["std"], marker="o", color=style.color_for(1),
+                capsize=3)
+    ax.set_xlabel("network latency  $L$ (slots)")
+    ax.set_ylabel("orphan rate (orphans / blocks)")
+    ax.set_title(f"Fork/orphan rate vs latency ({stake_dist} stake)")
+    return fig
+
+
+def _config_from_df(df: pd.DataFrame, stake_dist: str, n_nodes: int, latency: int) -> SimConfig:
+    """Reconstruct a single-epoch SimConfig from the swept parameters for re-simulation."""
+    r = df.iloc[0]
+    return SimConfig(
+        n_nodes=int(n_nodes), stake_dist=stake_dist, latency=int(latency),
+        k=int(r["k"]), pareto_shape=float(r["pareto_shape"]), f=float(r["f"]), epochs=1,
+    )
+
+
+# --- Figure 9: peak concurrent block proposals vs latency -------------------
+def concurrency_vs_latency(df: pd.DataFrame, stake_dist: str, n_nodes: int,
+                           reps: int = 6) -> plt.Figure:
+    """Proposals per latency-sized bucket vs L: the max bucket is the peak concurrency."""
+    style.apply_style()
+    latencies = sorted(int(x) for x in df["latency"].unique())
+    f = float(df["f"].iloc[0])
+    maxc, p99c, meanc = [], [], []
+    for lat in latencies:
+        cfg = _config_from_df(df, stake_dist, n_nodes, lat)
+        stats = [concurrency.concurrency_stats(cfg, replicate=r) for r in range(reps)]
+        maxc.append(max(s["max"] for s in stats))
+        p99c.append(float(np.mean([s["p99"] for s in stats])))
+        meanc.append(float(np.mean([s["mean"] for s in stats])))
+
+    fig, ax = plt.subplots()
+    ax.plot(latencies, maxc, "o-", color=style.color_for(1), label="max (peak concurrency)")
+    ax.plot(latencies, p99c, "s--", color=style.color_for(0), label="99th percentile bucket")
+    ax.plot(latencies, meanc, "^:", color=style.color_for(2), label="mean per bucket")
+    # expected proposals per bucket ~ L * (-ln(1-f)); reference for the mean
+    ref = [max(lat, 1) * (-np.log(1 - f)) for lat in latencies]
+    ax.plot(latencies, ref, color="0.6", lw=0.9, ls="-", zorder=0,
+            label=r"expected mean $=L\,(-\ln(1-f))$")
+    ax.set_xlabel("network latency  $L$ (slots) = bucket size")
+    ax.set_ylabel("block proposals per $L$-slot bucket")
+    ax.set_title(f"Concurrent block proposals per latency-window ({stake_dist}, N={n_nodes})")
+    ax.legend()
+    return fig
+
+
+# --- Figure 10: proposals-per-bucket across time (small multiples) ----------
+def concurrency_timeseries(df: pd.DataFrame, stake_dist: str, n_nodes: int,
+                           window_slots: int = 3000) -> plt.Figure:
+    """Block proposals per latency-sized bucket across time, one panel per latency.
+
+    Each panel steps through the proposals-per-bucket over the first ``window_slots``; the
+    dashed line marks the peak concurrency observed over the *whole* epoch for that latency.
+    """
+    style.apply_style()
+    latencies = sorted(int(x) for x in df["latency"].unique() if x > 0)
+    if len(latencies) > 4:
+        idx = np.linspace(0, len(latencies) - 1, 4).round().astype(int)
+        latencies = [latencies[i] for i in idx]
+
+    fig, axes = plt.subplots(len(latencies), 1, sharex=True,
+                             figsize=(6.4, 1.5 * len(latencies) + 0.5))
+    axes = np.atleast_1d(axes)
+    for i, (lat, ax) in enumerate(zip(latencies, axes, strict=True)):
+        cfg = _config_from_df(df, stake_dist, n_nodes, lat)
+        ws = concurrency.proposal_slots(cfg, replicate=0)
+        epoch_max = int(concurrency.window_counts(ws, cfg.epoch_len, lat).max())
+        span = min(window_slots, cfg.epoch_len)
+        counts = concurrency.window_counts(ws[ws < span], span, lat)
+        centers = (np.arange(counts.size) + 0.5) * lat
+        ax.step(centers, counts, where="mid", color=style.color_for(i), lw=1.2)
+        ax.axhline(epoch_max, color="0.45", ls="--", lw=0.9)
+        ax.text(0.995, 0.92, f"L={lat}  ·  epoch peak = {epoch_max}", transform=ax.transAxes,
+                ha="right", va="top", fontsize=8)
+        ax.set_ylim(0, epoch_max + 1)
+        ax.set_ylabel("proposals")
+    axes[-1].set_xlabel("slot (time)")
+    axes[0].set_title(
+        f"Block proposals across time, bucketed by latency ({stake_dist}, N={n_nodes})"
+    )
+    return fig
diff --git a/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/plotting/make_figures.py b/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/plotting/make_figures.py
new file mode 100644
index 0000000..9d13d2e
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/plotting/make_figures.py
@@ -0,0 +1,76 @@
+"""Render academic figures from a results frame (importable + CLI ``tsi-figures``)."""
+
+from __future__ import annotations
+
+import argparse
+import datetime as _dt
+from pathlib import Path
+
+import pandas as pd
+
+from . import figures as F
+from . import style
+
+
+def timestamped_figdir(outdir: str | Path, label: str) -> Path:
+    """A fresh ``<outdir>/<YYYY-MM-DD_HHMMSS>_<label>`` folder (never overwrites)."""
+    ts = _dt.datetime.now().strftime("%Y-%m-%d_%H%M%S")
+    base = Path(outdir) / f"{ts}_{label}"
+    d, i = base, 2
+    while d.exists():
+        d = base.with_name(f"{base.name}_{i}")
+        i += 1
+    return d
+
+
+def render(df: pd.DataFrame, out: str | Path) -> int:
+    """Render every applicable figure from ``df`` into directory ``out``; return the count."""
+    out = Path(out)
+    prov = F._provenance(df)
+    dists = sorted(df["stake_dist"].unique())
+    latencies = sorted(df["latency"].unique())
+    lat_focus = latencies[len(latencies) // 2] if latencies else 0
+    written: list[Path] = []
+
+    for dist in dists:
+        written += style.save(F.accuracy_vs_u(df, dist), out / f"01_accuracy_vs_u_{dist}", prov)
+        written += style.save(F.qeff_vs_u(df, dist), out / f"02_qeff_vs_u_{dist}", prov)
+        written += style.save(F.convergence(df, dist, lat_focus),
+                              out / f"03_convergence_{dist}_L{lat_focus}", prov)
+        written += style.save(F.orphan_diagnostics(df, dist), out / f"08_orphans_{dist}", prov)
+        for n in sorted(df["n_nodes"].unique()):
+            n = int(n)
+            written += style.save(F.heatmap_accuracy(df, dist, n),
+                                  out / f"04_heatmap_{dist}_N{n}", prov)
+            written += style.save(F.concurrency_vs_latency(df, dist, n),
+                                  out / f"09_concurrency_vs_latency_{dist}_N{n}", prov)
+            written += style.save(F.concurrency_timeseries(df, dist, n),
+                                  out / f"10_concurrency_timeseries_{dist}_N{n}", prov)
+        written += style.save(F.variance_vs_u(df, dist, lat_focus),
+                              out / f"06_variance_vs_u_{dist}_L{lat_focus}", prov)
+        if df[(df["stake_dist"] == dist) & (df["max_uncles"] > 0)]["uncle_strategy"].nunique() > 1:
+            written += style.save(F.strategy_comparison(df, dist, lat_focus),
+                                  out / f"07_strategy_{dist}_L{lat_focus}", prov)
+
+    written += style.save(
+        F.dist_comparison(df, lat_focus), out / f"05_dist_comparison_L{lat_focus}", prov
+    )
+    return len(written)
+
+
+def main(argv: list[str] | None = None) -> None:
+    ap = argparse.ArgumentParser(description="Render TSI figures")
+    ap.add_argument("--results", required=True)
+    ap.add_argument("--outdir", default="figures",
+                    help="parent dir; a dated sub-folder is created (never overwrites)")
+    ap.add_argument("--label", default=None, help="run label (default: results file stem)")
+    ap.add_argument("--out", default=None, help="explicit output dir (skips the dated folder)")
+    args = ap.parse_args(argv)
+    label = args.label or Path(args.results).stem
+    out = Path(args.out) if args.out else timestamped_figdir(args.outdir, label)
+    n = render(pd.read_parquet(args.results), out)
+    print(f"wrote {n} figures -> {out}")
+
+
+if __name__ == "__main__":  # pragma: no cover
+    main()
diff --git a/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/plotting/style.py b/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/plotting/style.py
new file mode 100644
index 0000000..280193a
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/plotting/style.py
@@ -0,0 +1,92 @@
+"""Shared academic matplotlib theme, palette, and helpers.
+
+Palette: Okabe-Ito — the standard colorblind-safe qualitative set — for categorical
+series (the ``U`` lines); ``cividis`` (perceptually uniform, CVD-safe) for heatmaps.
+Figures are saved as vector PDF (primary) plus 300-dpi PNG preview.
+"""
+
+from __future__ import annotations
+
+from collections.abc import Sequence
+from pathlib import Path
+
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+import numpy as np
+
+# Okabe-Ito colorblind-safe qualitative palette
+OKABE_ITO = [
+    "#0072B2",  # blue
+    "#D55E00",  # vermillion
+    "#009E73",  # bluish green
+    "#CC79A7",  # reddish purple
+    "#E69F00",  # orange
+    "#56B4E9",  # sky blue
+    "#F0E442",  # yellow
+    "#000000",  # black
+]
+SEQUENTIAL_CMAP = "cividis"
+# Diverging map centered at ratio = 1 (accuracy heatmaps)
+DIVERGING_CMAP = "RdBu_r"
+
+
+def apply_style() -> None:
+    """Install the shared rcParams theme (idempotent)."""
+    mpl.rcParams.update(
+        {
+            "figure.dpi": 150,
+            "savefig.dpi": 300,
+            "savefig.bbox": "tight",
+            "figure.figsize": (6.4, 4.0),
+            "font.size": 10,
+            "font.family": "sans-serif",
+            "axes.titlesize": 11,
+            "axes.labelsize": 10,
+            "axes.spines.top": False,
+            "axes.spines.right": False,
+            "axes.grid": True,
+            "grid.alpha": 0.25,
+            "grid.linewidth": 0.6,
+            "legend.frameon": False,
+            "legend.fontsize": 8.5,
+            "lines.linewidth": 1.8,
+            "axes.prop_cycle": mpl.cycler(color=OKABE_ITO),
+            "mathtext.default": "regular",
+        }
+    )
+
+
+def color_for(index: int) -> str:
+    return OKABE_ITO[index % len(OKABE_ITO)]
+
+
+def band_plot(
+    ax: plt.Axes,
+    x: Sequence[float],
+    series: np.ndarray,
+    *,
+    color: str,
+    label: str | None = None,
+    percentiles: tuple[float, float] = (10, 90),
+) -> None:
+    """Plot the mean of ``series`` (shape ``(n_replicates, len(x))``) with a percentile band."""
+    x = np.asarray(x, dtype=float)
+    mean = np.nanmean(series, axis=0)
+    lo = np.nanpercentile(series, percentiles[0], axis=0)
+    hi = np.nanpercentile(series, percentiles[1], axis=0)
+    ax.plot(x, mean, color=color, label=label)
+    ax.fill_between(x, lo, hi, color=color, alpha=0.18, linewidth=0)
+
+
+def save(fig: plt.Figure, out_stem: str | Path, provenance: str | None = None) -> list[Path]:
+    """Save ``fig`` as a 300-dpi PNG. ``out_stem`` has no suffix. Returns written paths."""
+    out_stem = Path(out_stem)
+    out_stem.parent.mkdir(parents=True, exist_ok=True)
+    if provenance:
+        # leave headroom so the footnote never collides with the x-axis label
+        fig.subplots_adjust(bottom=0.16)
+        fig.text(0.99, 0.005, provenance, fontsize=6, alpha=0.5, va="bottom", ha="right")
+    p = out_stem.with_suffix(".png")
+    fig.savefig(p)
+    plt.close(fig)
+    return [p]
diff --git a/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/rng.py b/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/rng.py
new file mode 100644
index 0000000..96316da
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/rng.py
@@ -0,0 +1,32 @@
+"""Deterministic, order-independent RNG derivation.
+
+Each ``SimConfig`` maps to an independent root ``SeedSequence`` seeded from a hash of its
+identity plus the global root seed. This guarantees a config yields the same random stream
+regardless of the order in which parallel workers run it. The engine ``spawn``\\s children
+of this root — one per epoch, and independent sub-streams within an epoch — so every draw
+(including the optional parallel chunked lottery) is a deterministic function of the root.
+"""
+
+from __future__ import annotations
+
+import hashlib
+
+import numpy as np
+
+from .config import SimConfig
+
+
+def _entropy(config: SimConfig) -> int:
+    payload = repr((config.root_seed, config.key())).encode()
+    digest = hashlib.blake2b(payload, digest_size=16).digest()
+    return int.from_bytes(digest, "big")
+
+
+def seedseq_for(config: SimConfig) -> np.random.SeedSequence:
+    """Return the reproducible root ``SeedSequence`` for this exact config+replicate."""
+    return np.random.SeedSequence(_entropy(config))
+
+
+def rng_for(config: SimConfig) -> np.random.Generator:
+    """Return the reproducible ``Generator`` for this exact config+replicate."""
+    return np.random.default_rng(seedseq_for(config))
diff --git a/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/stake.py b/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/stake.py
new file mode 100644
index 0000000..a9e3a32
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/stake.py
@@ -0,0 +1,32 @@
+"""Stake distribution generation.
+
+Total stake is held FIXED across distributions (via renormalisation) so that accuracy
+comparisons between ``uniform`` and ``pareto`` isolate the *shape* effect on the lottery
+(winner multiplicity / forking), not a difference in aggregate stake.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+
+from .config import SimConfig
+
+
+def make_stake(config: SimConfig, rng: np.random.Generator) -> np.ndarray:
+    """Return an ``(n_nodes,)`` non-negative stake vector summing to ``total_stake``."""
+    n = config.n_nodes
+    if config.stake_dist == "uniform":
+        if config.uniform_random:
+            w = rng.random(n)
+        else:
+            w = np.ones(n)
+    elif config.stake_dist == "pareto":
+        # numpy.pareto draws Lomax = Pareto(shape) - 1, heavy-tailed for small shape.
+        w = rng.pareto(config.pareto_shape, n) + 1.0
+    else:  # pragma: no cover - guarded by Literal typing
+        raise ValueError(f"unknown stake_dist: {config.stake_dist}")
+
+    total = w.sum()
+    if total <= 0:  # pragma: no cover - degenerate
+        raise ValueError("stake vector summed to zero")
+    return w * (config.total_stake / total)
diff --git a/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/sweep.py b/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/sweep.py
new file mode 100644
index 0000000..1f64e3e
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/sweep.py
@@ -0,0 +1,115 @@
+"""Parameter-sweep expansion, parallel execution, and result persistence.
+
+Across-config parallelism is the main multicore lever: ``run_trajectory`` is a pure
+function of an immutable, hash-seeded ``SimConfig``, so results are bitwise
+order-independent and the grid is embarrassingly parallel. We use joblib's process-based
+**loky** backend (this workload is CPU-bound Python that holds the GIL, so threads would
+serialise), and pin each worker to a single BLAS thread to avoid oversubscription.
+"""
+
+from __future__ import annotations
+
+import argparse
+import datetime as _dt
+import os
+from pathlib import Path
+
+import pandas as pd
+import yaml
+from joblib import Parallel, delayed
+from tqdm import tqdm
+
+from .config import SweepConfig
+from .engine import run_trajectory
+
+# Keep numpy/BLAS single-threaded inside each worker process (belt-and-braces alongside
+# joblib's inner_max_num_threads); prevents N_workers x N_blas_threads oversubscription.
+for _var in ("OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS", "NUMEXPR_NUM_THREADS"):
+    os.environ.setdefault(_var, "1")
+
+
+def run_sweep(
+    sweep: SweepConfig, n_jobs: int = -1, progress: bool = True, batch_size: str | int = "auto"
+) -> pd.DataFrame:
+    """Expand the grid, run every config across cores (loky), return one big frame.
+
+    ``n_jobs=-1`` uses all logical cores. ``batch_size`` defaults to joblib "auto" (good for
+    many tiny scaled-k tasks); pass ``1`` for a small grid of heavy full-scale configs so
+    they load-balance rather than pre-batch.
+    """
+    configs = sweep.expand()
+    runner = Parallel(
+        n_jobs=n_jobs,
+        backend="loky",
+        inner_max_num_threads=1,
+        batch_size=batch_size,
+        return_as="generator",
+    )(delayed(run_trajectory)(c) for c in configs)
+    if progress:
+        runner = tqdm(runner, total=len(configs), desc="configs")
+    rows: list[dict] = []
+    for traj in runner:
+        rows.extend(traj)
+    return pd.DataFrame(rows)
+
+
+def persist(df: pd.DataFrame, path: str | Path) -> None:
+    path = Path(path)
+    path.parent.mkdir(parents=True, exist_ok=True)
+    if path.suffix == ".parquet":
+        df.to_parquet(path, index=False)
+    else:
+        df.to_csv(path, index=False)
+
+
+def load_sweep_yaml(path: str | Path) -> SweepConfig:
+    with open(path) as fh:
+        data = yaml.safe_load(fh)
+    return SweepConfig.from_dict(data)
+
+
+def new_run_dir(outdir: str | Path, label: str) -> Path:
+    """Create a fresh timestamped run directory so runs never overwrite each other."""
+    ts = _dt.datetime.now().strftime("%Y-%m-%d_%H%M%S")
+    base = Path(outdir) / f"{ts}_{label}"
+    run_dir, i = base, 2
+    while run_dir.exists():                 # avoid same-second collisions
+        run_dir = base.with_name(f"{base.name}_{i}")
+        i += 1
+    run_dir.mkdir(parents=True)
+    return run_dir
+
+
+def main(argv: list[str] | None = None) -> None:
+    parser = argparse.ArgumentParser(description="Run a TSI parameter sweep + figures")
+    parser.add_argument("--config", required=True, help="sweep YAML")
+    parser.add_argument("--outdir", default="runs",
+                        help="parent dir; a dated sub-folder is created per run")
+    parser.add_argument("--label", default=None, help="run label (default: config name)")
+    parser.add_argument("--n-jobs", type=int, default=-1, help="-1 = all logical cores")
+    parser.add_argument("--batch-size", default="auto",
+                        help="'auto' (default) or an int; use 1 for a small heavy grid")
+    parser.add_argument("--no-figures", action="store_true",
+                        help="skip auto figure generation")
+    args = parser.parse_args(argv)
+
+    batch_size = int(args.batch_size) if args.batch_size != "auto" else "auto"
+    label = args.label or Path(args.config).stem
+    run_dir = new_run_dir(args.outdir, label)
+
+    sweep = load_sweep_yaml(args.config)
+    df = run_sweep(sweep, n_jobs=args.n_jobs, batch_size=batch_size)
+    results_path = run_dir / "results.parquet"
+    persist(df, results_path)
+    n_cfg = df[["n_nodes", "stake_dist", "latency", "max_uncles", "uncle_strategy", "replicate"]]
+    print(f"wrote {len(df)} rows ({len(n_cfg.drop_duplicates())} configs) -> {results_path}")
+
+    if not args.no_figures:
+        from .plotting.make_figures import render
+        n_fig = render(df, run_dir / "figures")
+        print(f"wrote {n_fig} figures -> {run_dir / 'figures'}")
+    print(f"run dir: {run_dir}")
+
+
+if __name__ == "__main__":  # pragma: no cover
+    main()
diff --git a/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/theory.py b/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/theory.py
new file mode 100644
index 0000000..8250745
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/theory.py
@@ -0,0 +1,60 @@
+"""Closed-form TSI results from ``analysis-total-stake-inference.md``.
+
+Used both as figure overlays and as ground truth for the verification checks. ``q`` is
+the honest active-slot utilisation; with uncle references, substitute the effective
+``q_eff`` to predict the improved accuracy.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+
+ArrayLike = np.ndarray | float
+
+
+def expected_ratio(f: float, q: ArrayLike) -> ArrayLike:
+    """Equilibrium ``E[D_inf] / D_true = log(1-f) / log(1 - f/q)`` for ``q in (f, 1]``."""
+    q = np.asarray(q, dtype=float)
+    return np.log(1.0 - f) / np.log(1.0 - f / q)
+
+
+def block_count_ceiling(f: float) -> float:
+    """Equilibrium ratio at *full* uncle recovery for the equal-stake limit.
+
+    TSI counts blocks (all lottery wins, rate ``-ln(1-f)`` per slot in the small-stake
+    limit), whereas ``f`` is the *active-slot* rate. So even with every orphan recovered the
+    estimate equilibrates at ``-ln(1-f)/f`` (~1.017 for f=1/30), not 1.0. This is a
+    deterministic overshoot floor, not noise; concentrated (Pareto) stake gives a smaller
+    value because of the concavity of ``phi``.
+    """
+    return -np.log(1.0 - f) / f
+
+
+def fixed_point_bias(f: float, precision: int = 1000) -> float:
+    """Extra multiplicative bias from the spec's integer f-truncation ``int(f*P)/P``."""
+    f_p = int(f * precision) / precision
+    return f / f_p
+
+
+def variance_ratio(f: float, q: ArrayLike, T: int, beta: float = 1.0) -> ArrayLike:
+    """Equilibrium ``Var[D_inf / D_true]``."""
+    q = np.asarray(q, dtype=float)
+    er = expected_ratio(f, q)
+    return (beta / f) ** 2 * (q / T) * er**2 * (1.0 - f) * f
+
+
+def variance_bound(f: float, T: int, beta: float = 1.0) -> float:
+    """Upper bound on ``Var[D_inf / D_true]`` at ``q = 1`` (perfect network)."""
+    return (beta / f) ** 2 / T * (1.0 - f) * f
+
+
+def beta_stability_bound(f: float, q: ArrayLike) -> ArrayLike:
+    """Stability threshold: convergence requires ``beta < 2f/((q-f) log(1/(1-f/q)))``."""
+    q = np.asarray(q, dtype=float)
+    return 2.0 * f / ((q - f) * np.log(1.0 / (1.0 - f / q)))
+
+
+def optimal_beta(f: float, q: ArrayLike) -> ArrayLike:
+    """Convergence-optimal learning rate ``f/((q-f) log(1/(1-f/q)))``."""
+    q = np.asarray(q, dtype=float)
+    return f / ((q - f) * np.log(1.0 / (1.0 - f / q)))
diff --git a/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/tsi.py b/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/tsi.py
new file mode 100644
index 0000000..9659f64
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/tsi.py
@@ -0,0 +1,91 @@
+"""Total Stake Inference: density counting and the per-epoch estimate update.
+
+The estimate update counts *blocks* exactly as the spec's ``density_over_slots`` does:
+``m = honest-chain blocks in window + deduplicated referenced uncle blocks (by their own
+slot) in window``. We additionally report slot-based ``q`` / ``q_eff`` (honest and
+uncle-recovered active-slot fractions) for comparison against the closed-form theory,
+which is written in terms of active-slot utilisation. The two differ only at rare
+multi-winner slots.
+"""
+
+from __future__ import annotations
+
+from dataclasses import dataclass
+
+import numpy as np
+
+from .blocktree import BlockTree
+
+
+def referenced_uncle_ids(tree: BlockTree, canonical_ids: list[int]) -> set[int]:
+    """Deduplicated set of uncle ids referenced by the canonical chain."""
+    ref: set[int] = set()
+    for b in canonical_ids:
+        ref.update(tree.uncles[b])
+    return ref
+
+
+def _in_window(slot: int, T: int) -> bool:
+    return 0 <= slot < T
+
+
+def density_m(tree: BlockTree, canonical_ids: list[int], T: int) -> int:
+    """Block count ``m`` for the TSI update (honest blocks + deduped uncles, in window)."""
+    s = tree.slot[canonical_ids]
+    honest = int(((s >= 0) & (s < T)).sum())
+    ref = referenced_uncle_ids(tree, canonical_ids)
+    uncle = sum(1 for u in ref if _in_window(int(tree.slot[u]), T))
+    return honest + uncle
+
+
+PRECISION = 1000  # spec on-chain fixed-point scale (cryptarchia-total-stake-inference.md)
+
+
+def update_D(
+    d_prev: float, m: int, T: int, f: float, beta: float, fixed_point: bool = False
+) -> float:
+    """Spec TSI recursion: ``max(1, D_prev * (1 - beta*(f_eff - m/T)/f_eff))``.
+
+    With ``fixed_point=True`` the target rate ``f`` is quantised exactly as the on-chain
+    algorithm does (``f_p = int(f*PRECISION)/PRECISION`` = 0.033 for f=1/30), which drives
+    the estimate to a measured density of 0.033 instead of 1/30 — a ~1% systematic
+    overestimate the deployed estimator exhibits but the exact-``f`` float model omits.
+    The remaining integer divisions in the spec (``tse/PRECISION``) are negligible at
+    realistic stakes and are not reproduced.
+    """
+    f_eff = (int(f * PRECISION) / PRECISION) if fixed_point else f
+    measured_density = m / T
+    d_new = d_prev * (1.0 - beta * (f_eff - measured_density) / f_eff)
+    return max(d_new, 1.0)
+
+
+@dataclass
+class SlotStats:
+    n_active: int          # active slots (>=1 winner) in window
+    n_honest: int          # honest-chain-occupied slots in window
+    n_recovered: int       # orphan-only slots recovered via referenced uncles
+    q: float               # n_honest / n_active
+    q_eff: float           # (n_honest + n_recovered) / n_active
+
+
+def slot_stats(
+    tree: BlockTree,
+    canonical_ids: list[int],
+    ref_uncle_ids: set[int],
+    active_slots: np.ndarray,
+    T: int,
+) -> SlotStats:
+    """Slot-based utilisation stats used for theory overlays."""
+    active_in = active_slots[(active_slots >= 0) & (active_slots < T)]
+    n_active = int(active_in.size)
+    honest_slots = {int(tree.slot[b]) for b in canonical_ids if _in_window(int(tree.slot[b]), T)}
+    recovered: set[int] = set()
+    for u in ref_uncle_ids:
+        su = int(tree.slot[u])
+        if _in_window(su, T) and su not in honest_slots:
+            recovered.add(su)
+    n_honest = len(honest_slots)
+    n_rec = len(recovered)
+    q = n_honest / n_active if n_active else float("nan")
+    q_eff = (n_honest + n_rec) / n_active if n_active else float("nan")
+    return SlotStats(n_active=n_active, n_honest=n_honest, n_recovered=n_rec, q=q, q_eff=q_eff)
diff --git a/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/uncles.py b/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/uncles.py
new file mode 100644
index 0000000..2ea398e
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/uncles.py
@@ -0,0 +1,82 @@
+"""Proposer-local uncle selection.
+
+For each canonical block ``B`` (processed oldest-first so ancestors' references are
+known), candidates are orphan (non-canonical) blocks ``U`` with
+``0 < slot_B - slot_U <= W`` that have not already been referenced by an ancestor of
+``B``. Two strategies match the spec: deterministic oldest-first, and random (oldest-first
+order, a coin of probability ``uncle_random_p`` per candidate, capped at ``U``). The spec's
+coin is unbiased (``uncle_random_p = 0.5``, the default); other values are a non-spec
+sensitivity knob. Dedup across ancestors is enforced by threading a ``referenced`` set down
+the canonical chain.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+
+from .blocktree import BlockTree
+from .config import SimConfig
+
+
+def _orphans_sorted(tree: BlockTree, canonical_ids: list[int]) -> tuple[np.ndarray, np.ndarray]:
+    """Return orphan block ids sorted by (slot, id) and their slots."""
+    canonical = np.zeros(tree.n_blocks, dtype=bool)
+    canonical[canonical_ids] = True
+    all_real = np.arange(1, tree.n_blocks)
+    orphan_ids = all_real[~canonical[1:]]
+    orphan_slots = tree.slot[orphan_ids]
+    order = np.lexsort((orphan_ids, orphan_slots))  # by slot, then id
+    return orphan_ids[order], orphan_slots[order]
+
+
+def annotate_uncles(
+    tree: BlockTree, canonical_ids: list[int], config: SimConfig, rng: np.random.Generator
+) -> None:
+    """Fill ``tree.uncles[B]`` for every canonical block ``B`` per the selection rule."""
+    u_max = config.max_uncles
+    if u_max <= 0:
+        return
+    w = config.uncle_window
+    orphan_ids, orphan_slots = _orphans_sorted(tree, canonical_ids)
+    if orphan_ids.size == 0:
+        return
+
+    referenced: set[int] = set()
+    # oldest canonical block first
+    for b in reversed(canonical_ids):
+        sb = int(tree.slot[b])
+        lo = int(np.searchsorted(orphan_slots, sb - w, side="left"))   # slot_U >= sb - W
+        hi = int(np.searchsorted(orphan_slots, sb, side="left"))       # slot_U <  sb
+        if hi <= lo:
+            continue
+        window_ids = orphan_ids[lo:hi]  # already oldest-first
+        selected = _select(window_ids, referenced, config, rng)
+        if selected:
+            tree.uncles[b] = tuple(selected)
+            referenced.update(selected)
+
+
+def _select(
+    window_ids: np.ndarray, referenced: set[int], config: SimConfig, rng: np.random.Generator
+) -> list[int]:
+    u_max = config.max_uncles
+    out: list[int] = []
+    if config.uncle_strategy == "oldest":
+        for bid in window_ids.tolist():
+            if bid in referenced:
+                continue
+            out.append(bid)
+            if len(out) >= u_max:
+                break
+    elif config.uncle_strategy == "random":
+        p = config.uncle_random_p
+        for bid in window_ids.tolist():
+            if bid in referenced:
+                continue
+            if rng.random() < p:
+                out.append(bid)
+                if len(out) >= u_max:
+                    break
+    else:  # pragma: no cover - guarded by Literal
+        raise ValueError(config.uncle_strategy)
+    return out
diff --git a/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/verify.py b/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/verify.py
new file mode 100644
index 0000000..efd73b8
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim-mc/src/tsi_sim/verify.py
@@ -0,0 +1,116 @@
+"""Analytic sanity checks: simulator vs closed-form theory (``tsi-verify``).
+
+Replicate runs are evaluated across CPU cores. Exits non-zero if any check fails.
+"""
+
+from __future__ import annotations
+
+from dataclasses import replace
+
+import numpy as np
+import pandas as pd
+from joblib import Parallel, delayed
+
+from .config import SimConfig
+from .engine import run_trajectory
+from .epoch import simulate_epoch
+from .rng import seedseq_for
+from .stake import make_stake
+from .theory import block_count_ceiling, expected_ratio, fixed_point_bias
+
+F = 1.0 / 30.0
+K = 128            # scaled: T = 6*floor(128/f) = 23040 slots
+EPOCHS = 45
+REPS = 12
+BURN = 25
+
+
+def tail_mean(cfg: SimConfig, col: str, reps: int = REPS, n_jobs: int = -1) -> tuple[float, float]:
+    """Mean over ``reps`` replicates of each trajectory's post-burn-in tail mean of ``col``."""
+    def one(r: int) -> float:
+        df = pd.DataFrame(run_trajectory(replace(cfg, replicate=r)))
+        return float(df[col].iloc[BURN:].mean())
+
+    vals = np.array(Parallel(n_jobs=n_jobs, backend="loky", inner_max_num_threads=1)(
+        delayed(one)(r) for r in range(reps)
+    ))
+    return float(vals.mean()), float(vals.std() / np.sqrt(reps))
+
+
+def check(name: str, ok: bool, detail: str) -> bool:
+    print(f"[{'PASS' if ok else 'FAIL'}] {name}: {detail}")
+    return ok
+
+
+def main() -> int:
+    results = []
+
+    # 1. Active-slot rate ~= f when D_est = D_true, L=0, U=0.
+    cfg = SimConfig(n_nodes=2000, stake_dist="uniform", latency=0, max_uncles=0,
+                    k=K, epochs=1, genesis_d_factor=1.0)
+    ss = seedseq_for(cfg)
+    children = ss.spawn(2)
+    stake = make_stake(cfg, np.random.default_rng(children[0]))
+    er = simulate_epoch(cfg, stake, float(stake.sum()), children[1])
+    active_rate = er.n_active / cfg.period_T
+    results.append(check("active-slot rate ~ f (L=0,U=0,D=D_true)",
+                         abs(active_rate - F) / F < 0.05,
+                         f"active_rate={active_rate:.5f} f={F:.5f}"))
+
+    # 2. U=0 equilibrium ratio ~= expected_ratio(f, measured q).
+    cfg = SimConfig(n_nodes=1000, stake_dist="uniform", latency=4, max_uncles=0,
+                    k=K, epochs=EPOCHS, genesis_d_factor=0.5)
+    ratio, se = tail_mean(cfg, "ratio")
+    q, _ = tail_mean(cfg, "q")
+    pred = float(expected_ratio(F, q))
+    results.append(check("U=0 ratio ~ theory(q)",
+                         abs(ratio - pred) < 0.02 + 2 * se,
+                         f"sim={ratio:.4f}±{se:.4f} theory(q={q:.3f})={pred:.4f}"))
+
+    # 3. Underestimate at higher latency (q < 1 => ratio < 1), U=0.
+    cfg = SimConfig(n_nodes=1000, stake_dist="uniform", latency=10, max_uncles=0,
+                    k=K, epochs=EPOCHS, genesis_d_factor=0.5)
+    ratio, se = tail_mean(cfg, "ratio")
+    q, _ = tail_mean(cfg, "q")
+    results.append(check("U=0 underestimates true stake at latency",
+                         ratio < 0.98 and q < 0.98,
+                         f"ratio={ratio:.4f} q={q:.3f}"))
+
+    # 4. q_eff -> 1 and ratio -> block-count ceiling as U grows (uncles recover forks).
+    #    The residual |ratio-1| ~ 0.017 is the intrinsic winners-vs-active-slots floor
+    #    (density_m counts blocks), NOT a convergence failure -- see check 5.
+    base = dict(n_nodes=1000, stake_dist="uniform", latency=8, uncle_strategy="oldest",
+                k=K, epochs=EPOCHS, genesis_d_factor=0.5)
+    r0, _ = tail_mean(SimConfig(max_uncles=0, **base), "ratio")
+    r4, _ = tail_mean(SimConfig(max_uncles=4, **base), "ratio")
+    qe4, _ = tail_mean(SimConfig(max_uncles=4, **base), "q_eff")
+    results.append(check("uncles recover accuracy (q_eff->1, |ratio-1| shrinks)",
+                         qe4 > 0.99 and abs(r4 - 1) < abs(r0 - 1) and abs(r4 - 1) < 0.03,
+                         f"ratio U0={r0:.4f} -> U4={r4:.4f}; q_eff(U4)={qe4:.4f}"))
+
+    # 5. Full recovery equilibrates at the block-count ceiling -ln(1-f)/f (~1.017), not 1.
+    ceiling = float(block_count_ceiling(F))
+    r4_val, se4 = tail_mean(SimConfig(max_uncles=4, **base), "ratio")
+    results.append(check("full-recovery ratio ~ block-count ceiling",
+                         abs(r4_val - ceiling) < 0.015 + 2 * se4,
+                         f"ratio(U4)={r4_val:.4f} ceiling=-ln(1-f)/f={ceiling:.4f}"))
+
+    # 6. Fixed-point mode adds the spec's ~1% f-truncation overestimate.
+    fp_base = dict(n_nodes=1000, stake_dist="uniform", latency=0, max_uncles=0,
+                   k=K, epochs=EPOCHS, genesis_d_factor=1.0)
+    r_float, _ = tail_mean(SimConfig(fixed_point=False, **fp_base), "ratio")
+    r_fixed, _ = tail_mean(SimConfig(fixed_point=True, **fp_base), "ratio")
+    bias = float(fixed_point_bias(F))
+    results.append(check("fixed_point mode adds ~1% f-truncation bias",
+                         r_fixed > r_float and abs(r_fixed / r_float - bias) < 0.01,
+                         f"float={r_float:.4f} fixed={r_fixed:.4f} ratio={r_fixed/r_float:.4f} "
+                         f"expected~{bias:.4f}"))
+
+    print()
+    n_pass = sum(results)
+    print(f"{n_pass}/{len(results)} checks passed")
+    return 0 if n_pass == len(results) else 1
+
+
+if __name__ == "__main__":  # pragma: no cover
+    raise SystemExit(main())
diff --git a/tools/simulators/tsi/tsi-sim-mc/tests/test_blocktree.py b/tools/simulators/tsi/tsi-sim-mc/tests/test_blocktree.py
new file mode 100644
index 0000000..367581a
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim-mc/tests/test_blocktree.py
@@ -0,0 +1,54 @@
+import numpy as np
+
+from tsi_sim.blocktree import build_tree
+from tsi_sim.latency import FixedSlotLatency
+
+
+def _winners(*groups):
+    return [np.array(g, dtype=np.int64) for g in groups]
+
+
+def test_latency_induces_fork_and_longest_chain():
+    # L=2. slot0: node0; slot1: node1 (can't see block1 yet -> forks on genesis);
+    # slot3: node2 (sees both, builds on the earlier-slot tip -> block1).
+    active = np.array([0, 1, 3], dtype=np.int64)
+    winners = _winners([0], [1], [2])
+    tree = build_tree(active, winners, FixedSlotLatency(2), np.random.default_rng(0))
+
+    assert tree.n_blocks == 4                      # genesis + 3
+    assert tree.height.tolist() == [0, 1, 1, 2]
+    # block 3 built on block 1 (earlier slot wins the height-1 tie)
+    assert tree.parent[3] == 1
+    assert tree.canonical_chain() == [3, 1]        # tip-first
+    # block 2 is the orphan
+    canon = set(tree.canonical_chain())
+    orphans = [b for b in range(1, tree.n_blocks) if b not in canon]
+    assert orphans == [2]
+
+
+def test_same_slot_multiwinner_forks_even_at_zero_latency():
+    active = np.array([0, 1], dtype=np.int64)
+    winners = _winners([0, 1], [2])   # two winners in slot 0 -> guaranteed fork
+    tree = build_tree(active, winners, FixedSlotLatency(0), np.random.default_rng(0))
+    # blocks 1 and 2 are siblings at height 1 on genesis
+    assert tree.parent[1] == 0 and tree.parent[2] == 0
+    assert tree.height[1] == 1 and tree.height[2] == 1
+    # block 3 at slot 1 extends one of them (height 2)
+    assert tree.height[3] == 2
+    assert len(tree.canonical_chain()) == 2
+
+
+def test_self_extension_within_latency():
+    # A single node winning consecutive slots builds on its own block despite latency.
+    active = np.array([0, 1], dtype=np.int64)
+    winners = _winners([5], [5])
+    tree = build_tree(active, winners, FixedSlotLatency(10), np.random.default_rng(0))
+    assert tree.parent[2] == 1          # node 5 self-extends
+    assert tree.height.tolist() == [0, 1, 2]
+
+
+def test_ancestors():
+    active = np.array([0, 1, 2], dtype=np.int64)
+    winners = _winners([0], [0], [0])   # one node, clean chain
+    tree = build_tree(active, winners, FixedSlotLatency(0), np.random.default_rng(0))
+    assert tree.ancestors(3) == [3, 2, 1]
diff --git a/tools/simulators/tsi/tsi-sim-mc/tests/test_concurrency.py b/tools/simulators/tsi/tsi-sim-mc/tests/test_concurrency.py
new file mode 100644
index 0000000..9e6cc93
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim-mc/tests/test_concurrency.py
@@ -0,0 +1,40 @@
+import numpy as np
+
+from tsi_sim import concurrency
+from tsi_sim.config import SimConfig
+
+
+def test_window_counts_partitions_all_proposals():
+    ws = np.array([0, 1, 2, 5, 5, 9], dtype=np.int64)
+    counts = concurrency.window_counts(ws, epoch_len=10, bucket=5)
+    assert counts.tolist() == [3, 3]          # slots 0-4 -> 3, slots 5-9 -> 3
+    assert counts.sum() == ws.size
+
+
+def test_window_counts_bucket_one_is_per_slot():
+    ws = np.array([0, 0, 3], dtype=np.int64)
+    counts = concurrency.window_counts(ws, epoch_len=4, bucket=1)
+    assert counts.tolist() == [2, 0, 0, 1]    # slot 0 has 2 concurrent, slot 3 has 1
+
+
+def test_window_counts_empty():
+    assert concurrency.window_counts(np.empty(0, np.int64), 10, 2).tolist() == [0, 0, 0, 0, 0]
+
+
+def test_concurrency_stats_scale_with_latency():
+    # Bigger latency bucket => more proposals per bucket (max/mean grow roughly with L).
+    small = concurrency.concurrency_stats(SimConfig(n_nodes=1000, latency=2, k=32, epochs=1))
+    large = concurrency.concurrency_stats(SimConfig(n_nodes=1000, latency=20, k=32, epochs=1))
+    assert large["mean"] > small["mean"]
+    assert large["max"] >= small["max"]
+    # mean per bucket ~ bucket * (-ln(1-f))
+    expected = large["bucket"] * (-np.log(1 - 1 / 30))
+    assert abs(large["mean"] - expected) < 0.25 * expected
+
+
+def test_proposal_slots_reproducible_and_sorted():
+    cfg = SimConfig(n_nodes=500, latency=4, k=16, epochs=1)
+    a = concurrency.proposal_slots(cfg, replicate=0)
+    b = concurrency.proposal_slots(cfg, replicate=0)
+    np.testing.assert_array_equal(a, b)
+    assert np.all(np.diff(a) >= 0)
diff --git a/tools/simulators/tsi/tsi-sim-mc/tests/test_config.py b/tools/simulators/tsi/tsi-sim-mc/tests/test_config.py
new file mode 100644
index 0000000..5d693fa
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim-mc/tests/test_config.py
@@ -0,0 +1,81 @@
+import dataclasses
+
+import pytest
+
+from tsi_sim.config import SimConfig, SweepConfig
+
+
+def test_expand_cardinality_and_u0_strategy_dedup():
+    sweep = SweepConfig(
+        n_nodes=[1000, 2000],
+        stake_dist=["uniform"],
+        latency=[0, 4],
+        max_uncles=[0, 1, 2],
+        uncle_strategy=["oldest", "random"],
+        replicates=3,
+        base={"k": 16, "epochs": 5},
+    )
+    configs = sweep.expand()
+    # U=0 keeps only the first strategy; U>0 keeps both.
+    # per (n,dist,lat): U0 x1 strat + U1 x2 + U2 x2 = 5 strat-U combos, x3 reps = 15
+    # x 2 n_nodes x 1 dist x 2 lat = 60
+    assert len(configs) == 60
+    u0 = [c for c in configs if c.max_uncles == 0]
+    assert all(c.uncle_strategy == "oldest" for c in u0)
+    assert {c.replicate for c in configs} == {0, 1, 2}
+
+
+def test_base_propagation():
+    sweep = SweepConfig(n_nodes=[500], stake_dist=["pareto"], latency=[2], max_uncles=[0],
+                        uncle_strategy=["oldest"], replicates=1,
+                        base={"k": 32, "epochs": 7, "fixed_point": True})
+    (c,) = sweep.expand()
+    assert c.k == 32 and c.epochs == 7 and c.fixed_point is True and c.stake_dist == "pareto"
+
+
+@pytest.mark.parametrize("kwargs", [
+    {"k": 0}, {"epochs": 0}, {"n_nodes": 0}, {"latency": -1}, {"max_uncles": -1},
+    {"uncle_window": 0}, {"lottery_chunks": 0}, {"uncle_random_p": 1.5}, {"f": 0.0},
+    {"f": 1.0}, {"beta": 0.0}, {"genesis_d_factor": 0.0}, {"pareto_shape": 0.0},
+    {"stake_dist": "zipf"}, {"uncle_strategy": "newest"},
+])
+def test_validation_rejects_bad_fields(kwargs):
+    with pytest.raises(ValueError):
+        SimConfig(**kwargs)
+
+
+def test_unknown_sweep_key_rejected():
+    with pytest.raises(ValueError, match="unknown sweep keys"):
+        SweepConfig.from_dict({"latencies": [0, 1], "base": {}})  # typo: latencies vs latency
+
+
+def test_from_dict_roundtrip_ok():
+    sw = SweepConfig.from_dict({"latency": [0, 3], "max_uncles": [0, 2], "replicates": 2,
+                                "base": {"k": 8, "epochs": 4}})
+    assert sw.latency == [0, 3] and sw.replicates == 2 and sw.base["k"] == 8
+
+
+def test_key_covers_every_field():
+    # Guard against the silent shared-RNG bug: key() must reflect all run-affecting fields.
+    ignored = {"root_seed"}  # root_seed enters _entropy separately, not via key()
+    names = {f.name for f in dataclasses.fields(SimConfig)} - ignored
+    a = SimConfig()
+    for name in names:
+        cur = getattr(a, name)
+        alt = _perturb(cur)
+        b = dataclasses.replace(a, **{name: alt})
+        assert a.key() != b.key(), f"key() does not distinguish field {name!r}"
+
+
+def _perturb(v):
+    if isinstance(v, bool):
+        return not v
+    if isinstance(v, int):
+        return v + 1
+    if isinstance(v, float):
+        return v + 0.001
+    if v == "uniform":
+        return "pareto"
+    if v == "oldest":
+        return "random"
+    return v
diff --git a/tools/simulators/tsi/tsi-sim-mc/tests/test_latency.py b/tools/simulators/tsi/tsi-sim-mc/tests/test_latency.py
new file mode 100644
index 0000000..a2ff97e
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim-mc/tests/test_latency.py
@@ -0,0 +1,35 @@
+import numpy as np
+
+from tsi_sim.config import SimConfig
+from tsi_sim.engine import run_trajectory
+from tsi_sim.latency import FixedSlotLatency, RealisticLatency, make_latency
+
+
+def test_fixed_slot_latency():
+    lat = FixedSlotLatency(3)
+    assert lat.visible_at(10, np.random.default_rng(0)) == 13
+
+
+def test_realistic_latency_zero_mean_immediate():
+    lat = RealisticLatency(0.0)
+    assert lat.visible_at(5, np.random.default_rng(0)) == 5
+
+
+def test_realistic_latency_positive_delays():
+    lat = RealisticLatency(4.0)
+    rng = np.random.default_rng(0)
+    vals = [lat.visible_at(100, rng) for _ in range(200)]
+    assert all(v >= 101 for v in vals)       # strictly after production
+    assert np.mean([v - 100 for v in vals]) > 1.0
+
+
+def test_make_latency_dispatch():
+    assert isinstance(make_latency(SimConfig(latency=2)), FixedSlotLatency)
+    assert isinstance(make_latency(SimConfig(latency=2, latency_stochastic=True)), RealisticLatency)
+
+
+def test_stochastic_latency_integration_runs():
+    cfg = SimConfig(n_nodes=500, k=8, epochs=5, latency=6, latency_stochastic=True, max_uncles=2)
+    rows = run_trajectory(cfg)
+    assert len(rows) == 5
+    assert all(np.isfinite(r["ratio"]) for r in rows)
diff --git a/tools/simulators/tsi/tsi-sim-mc/tests/test_lottery.py b/tools/simulators/tsi/tsi-sim-mc/tests/test_lottery.py
new file mode 100644
index 0000000..cf22408
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim-mc/tests/test_lottery.py
@@ -0,0 +1,77 @@
+import numpy as np
+
+from tsi_sim import lottery
+
+
+def test_phi_bounds():
+    f = 1 / 30
+    assert lottery.phi(f, 0.0) == 0.0
+    # a single all-stake node (alpha=1) wins at exactly rate f
+    assert abs(lottery.phi(f, 1.0) - f) < 1e-12
+
+
+def test_win_probs_monotone_in_stake():
+    stake = np.array([1.0, 2.0, 3.0])
+    p = lottery.win_probs(stake, d_est=6.0, f=1 / 30)
+    assert np.all(np.diff(p) > 0)
+
+
+def test_sample_wins_sorted_and_rate():
+    rng = np.random.default_rng(0)
+    n, slots = 500, 4000
+    p = np.full(n, 0.001)
+    ws, wn = lottery.sample_wins(p, slots, rng, chunk=512)
+    assert np.all(np.diff(ws) >= 0)                     # sorted by slot
+    assert ws.shape == wn.shape
+    assert np.all((wn >= 0) & (wn < n))
+    # expected wins ~ n * slots * p
+    assert abs(ws.size - n * slots * 0.001) < 4 * np.sqrt(n * slots * 0.001)
+
+
+def test_group_by_slot():
+    ws = np.array([0, 0, 2, 5, 5, 5])
+    wn = np.array([3, 7, 1, 2, 4, 9])
+    active, groups = lottery.group_by_slot(ws, wn)
+    assert list(active) == [0, 2, 5]
+    assert [g.tolist() for g in groups] == [[3, 7], [1], [2, 4, 9]]
+
+
+def test_group_by_slot_empty():
+    active, groups = lottery.group_by_slot(np.empty(0, int), np.empty(0, int))
+    assert active.size == 0 and groups == []
+
+
+def test_sparse_per_node_distinct_slots():
+    # Each node wins any slot at most once (independent Bernoulli-per-slot invariant).
+    rng = np.random.default_rng(1)
+    p = np.full(300, 0.05)
+    ws, wn = lottery.sample_wins(p, 2000, rng)
+    for node in np.unique(wn):
+        slots = ws[wn == node]
+        assert slots.size == np.unique(slots).size       # no duplicate (node, slot)
+
+
+def test_sparse_preserves_multiwinner_slots():
+    # With high p, some slots have >1 distinct winner (guaranteed forks) — must be possible.
+    rng = np.random.default_rng(2)
+    p = np.full(50, 0.3)
+    ws, _ = lottery.sample_wins(p, 500, rng)
+    _, counts = np.unique(ws, return_counts=True)
+    assert counts.max() >= 2
+
+
+def test_chunked_matches_serial_distribution():
+    # Chunked sampler is deterministic given (seedseq, n_chunks) and statistically matches
+    # serial. NOTE: SeedSequence.spawn is stateful, so each call needs a FRESH SeedSequence.
+    p = np.full(400, 0.02)
+    a_s, a_n = lottery.sample_wins_chunked(p, 100000, np.random.SeedSequence(123), 4, n_jobs=1)
+    b_s, b_n = lottery.sample_wins_chunked(p, 100000, np.random.SeedSequence(123), 4, n_jobs=1)
+    np.testing.assert_array_equal(a_s, b_s)              # deterministic
+    np.testing.assert_array_equal(a_n, b_n)
+    assert np.all(np.diff(a_s) >= 0)                    # sorted
+    for node in np.unique(a_n):                         # per-node distinct slots preserved
+        s = a_s[a_n == node]
+        assert s.size == np.unique(s).size
+    serial_n = lottery.sample_wins(p, 100000, np.random.default_rng(7))[0].size
+    assert abs(a_s.size - serial_n) < 6 * np.sqrt(serial_n)  # same rate
+
diff --git a/tools/simulators/tsi/tsi-sim-mc/tests/test_rng.py b/tools/simulators/tsi/tsi-sim-mc/tests/test_rng.py
new file mode 100644
index 0000000..8903527
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim-mc/tests/test_rng.py
@@ -0,0 +1,31 @@
+import numpy as np
+
+from tsi_sim.config import SimConfig
+from tsi_sim.engine import run_trajectory
+from tsi_sim.rng import rng_for, seedseq_for
+
+
+def test_seedseq_and_rng_deterministic():
+    cfg = SimConfig(k=8, epochs=3)
+    a = np.random.default_rng(seedseq_for(cfg)).random(5)
+    b = rng_for(cfg).random(5)
+    np.testing.assert_array_equal(a, b)
+
+
+def test_distinct_configs_get_distinct_streams():
+    c0 = SimConfig(k=8, epochs=3, latency=0)
+    c1 = SimConfig(k=8, epochs=3, latency=1)
+    assert not np.array_equal(rng_for(c0).random(4), rng_for(c1).random(4))
+
+
+def test_trajectory_is_order_independent_and_reproducible():
+    cfg = SimConfig(n_nodes=500, k=8, epochs=6, latency=4, max_uncles=2)
+    r1 = run_trajectory(cfg)
+    r2 = run_trajectory(cfg)
+    assert [row["ratio"] for row in r1] == [row["ratio"] for row in r2]
+
+
+def test_replicates_differ():
+    a = run_trajectory(SimConfig(n_nodes=500, k=8, epochs=6, latency=4, replicate=0))
+    b = run_trajectory(SimConfig(n_nodes=500, k=8, epochs=6, latency=4, replicate=1))
+    assert a[-1]["ratio"] != b[-1]["ratio"]
diff --git a/tools/simulators/tsi/tsi-sim-mc/tests/test_stake.py b/tools/simulators/tsi/tsi-sim-mc/tests/test_stake.py
new file mode 100644
index 0000000..6fca075
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim-mc/tests/test_stake.py
@@ -0,0 +1,40 @@
+import numpy as np
+
+from tsi_sim.config import SimConfig
+from tsi_sim.stake import make_stake
+
+
+def _rng():
+    return np.random.default_rng(0)
+
+
+def test_uniform_equal_and_sum():
+    cfg = SimConfig(n_nodes=100, stake_dist="uniform", total_stake=1e9)
+    w = make_stake(cfg, _rng())
+    assert w.shape == (100,)
+    assert np.allclose(w, w[0])                       # equal weights
+    assert abs(w.sum() - 1e9) < 1e-3
+
+
+def test_uniform_random_varies_but_sums():
+    cfg = SimConfig(n_nodes=200, stake_dist="uniform", uniform_random=True, total_stake=5e8)
+    w = make_stake(cfg, _rng())
+    assert w.std() > 0
+    assert abs(w.sum() - 5e8) < 1e-2
+    assert np.all(w >= 0)
+
+
+def test_pareto_sum_fixed_and_heavier_tailed():
+    n = 5000
+    uni = make_stake(SimConfig(n_nodes=n, stake_dist="uniform", total_stake=1e9), _rng())
+    par = make_stake(SimConfig(n_nodes=n, stake_dist="pareto", pareto_shape=1.16,
+                               total_stake=1e9), _rng())
+    assert abs(par.sum() - 1e9) < 1.0
+    assert par.max() > uni.max() * 5        # heavy tail: richest holds far more
+    assert np.all(par >= 0)
+
+
+def test_total_stake_fixed_across_distributions():
+    uni = make_stake(SimConfig(n_nodes=1000, stake_dist="uniform", total_stake=7e8), _rng())
+    par = make_stake(SimConfig(n_nodes=1000, stake_dist="pareto", total_stake=7e8), _rng())
+    assert abs(uni.sum() - par.sum()) < 1.0   # comparability guarantee
diff --git a/tools/simulators/tsi/tsi-sim-mc/tests/test_theory.py b/tools/simulators/tsi/tsi-sim-mc/tests/test_theory.py
new file mode 100644
index 0000000..8141e8b
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim-mc/tests/test_theory.py
@@ -0,0 +1,40 @@
+import numpy as np
+
+from tsi_sim import theory
+
+F = 1 / 30
+T = 10000
+
+
+def test_expected_ratio_unbiased_at_q1():
+    assert abs(float(theory.expected_ratio(F, 1.0)) - 1.0) < 1e-12
+
+
+def test_expected_ratio_monotone_in_q():
+    qs = np.linspace(0.5, 1.0, 20)
+    er = theory.expected_ratio(F, qs)
+    assert np.all(np.diff(er) > 0)          # accuracy improves as q -> 1
+    assert np.all(er <= 1.0 + 1e-12)        # always an underestimate
+
+
+def test_variance_bound_matches_at_q1():
+    v = float(theory.variance_ratio(F, 1.0, T))
+    assert abs(v - theory.variance_bound(F, T)) < 1e-15
+
+
+def test_optimal_beta_is_half_stability_bound():
+    for q in (0.7, 0.85, 0.95):
+        opt = float(theory.optimal_beta(F, q))
+        bound = float(theory.beta_stability_bound(F, q))
+        assert abs(opt - bound / 2) < 1e-12
+
+
+def test_block_count_ceiling_above_one():
+    c = theory.block_count_ceiling(F)
+    assert 1.015 < c < 1.02                 # -ln(1-1/30)/(1/30) ~ 1.01705
+
+
+def test_fixed_point_bias_about_one_percent():
+    b = theory.fixed_point_bias(F)
+    assert abs(b - (F / (33 / 1000))) < 1e-12
+    assert 1.005 < b < 1.02
diff --git a/tools/simulators/tsi/tsi-sim-mc/tests/test_theory_convergence.py b/tools/simulators/tsi/tsi-sim-mc/tests/test_theory_convergence.py
new file mode 100644
index 0000000..5490d96
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim-mc/tests/test_theory_convergence.py
@@ -0,0 +1,49 @@
+"""End-to-end statistical checks against closed-form theory (scaled k)."""
+
+import numpy as np
+import pandas as pd
+import pytest
+
+from tsi_sim.config import SimConfig
+from tsi_sim.engine import run_trajectory
+from tsi_sim.theory import expected_ratio
+
+F = 1 / 30
+
+
+def _tail_mean(col, **cfg):
+    reps = cfg.pop("reps", 6)
+    burn = cfg["epochs"] // 2
+    vals = []
+    for r in range(reps):
+        df = pd.DataFrame(run_trajectory(SimConfig(replicate=r, **cfg)))
+        vals.append(df[col].iloc[burn:].mean())
+    return float(np.mean(vals))
+
+
+def test_baseline_exact_without_forks():
+    # No latency, no uncles: active-slot rate == f, so the estimate is unbiased.
+    ratio = _tail_mean("ratio", n_nodes=1000, stake_dist="uniform", latency=0,
+                       max_uncles=0, k=64, epochs=35, genesis_d_factor=0.5, reps=6)
+    assert abs(ratio - 1.0) < 0.02
+
+
+@pytest.mark.slow
+def test_u0_matches_expected_ratio():
+    cfg = dict(n_nodes=1000, stake_dist="uniform", latency=6, max_uncles=0,
+               k=96, epochs=45, genesis_d_factor=0.5)
+    ratio = _tail_mean("ratio", reps=8, **cfg)
+    q = _tail_mean("q", reps=8, **cfg)
+    assert abs(ratio - float(expected_ratio(F, q))) < 0.02
+
+
+@pytest.mark.slow
+def test_uncles_recover_accuracy():
+    common = dict(n_nodes=1000, stake_dist="uniform", latency=8, uncle_strategy="oldest",
+                  k=96, epochs=45, genesis_d_factor=0.5)
+    r0 = _tail_mean("ratio", max_uncles=0, reps=8, **common)
+    r4 = _tail_mean("ratio", max_uncles=4, reps=8, **common)
+    qe4 = _tail_mean("q_eff", max_uncles=4, reps=8, **common)
+    assert qe4 > 0.99
+    assert abs(r4 - 1) < abs(r0 - 1)     # uncles reduce the error
+    assert abs(r4 - 1) < 0.03            # residual is the small block-count overshoot
diff --git a/tools/simulators/tsi/tsi-sim-mc/tests/test_tsi_counting.py b/tools/simulators/tsi/tsi-sim-mc/tests/test_tsi_counting.py
new file mode 100644
index 0000000..cc71cf0
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim-mc/tests/test_tsi_counting.py
@@ -0,0 +1,87 @@
+import numpy as np
+
+from tsi_sim.blocktree import BlockTree
+from tsi_sim.tsi import density_m, referenced_uncle_ids, slot_stats, update_D
+
+
+def make_tree(slots, parents, heights, uncles):
+    n = len(slots)
+    return BlockTree(
+        slot=np.array(slots, np.int64),
+        parent=np.array(parents, np.int64),
+        height=np.array(heights, np.int64),
+        leader=np.zeros(n, np.int64),
+        uncles=uncles,
+    )
+
+
+def test_density_counts_honest_plus_deduped_uncles_in_window():
+    # canonical 1(slot0),4(slot3); orphans 2(slot1),3(slot2). block4 refs uncles 2 and 3.
+    tree = make_tree(
+        slots=[-1, 0, 1, 2, 3],
+        parents=[-1, 0, 0, 0, 1],
+        heights=[0, 1, 1, 1, 2],
+        uncles=[(), (), (), (), (2, 3)],
+    )
+    canonical = [4, 1]  # tip-first
+    # window T=10 includes all slots
+    assert density_m(tree, canonical, T=10) == 4      # 2 honest (slots 0,3) + 2 uncles (1,2)
+    # window T=2 excludes slots 2,3 -> honest {slot0}=1, uncle slot1=1 (slot2 excluded)
+    assert density_m(tree, canonical, T=2) == 2
+
+
+def test_uncle_counted_by_own_slot_and_deduped():
+    tree = make_tree(
+        slots=[-1, 0, 5, 1],
+        parents=[-1, 0, 1, 0],
+        heights=[0, 1, 2, 1],
+        uncles=[(), (), (3,), ()],   # block2 (slot5) references orphan 3 (slot1)
+    )
+    canonical = [2, 1]
+    assert referenced_uncle_ids(tree, canonical) == {3}
+    # uncle counted by its OWN slot (1), so window T=2 includes it
+    assert density_m(tree, canonical, T=2) == 2       # honest slot0 + uncle slot1
+    # window that excludes the uncle's own slot
+    assert density_m(tree, canonical, T=1) == 1       # only honest slot0
+
+
+def test_slot_stats_q_and_qeff():
+    # active slots 0,1,2 in window; honest occupies 0,2; orphan at slot1 recovered by uncle
+    tree = make_tree(
+        slots=[-1, 0, 1, 2],
+        parents=[-1, 0, 0, 1],
+        heights=[0, 1, 1, 2],
+        uncles=[(), (), (), (2,)],   # block3 refs orphan 2 (slot1)
+    )
+    canonical = [3, 1]
+    active = np.array([0, 1, 2], np.int64)
+    ref = referenced_uncle_ids(tree, canonical)
+    ss = slot_stats(tree, canonical, ref, active, T=10)
+    assert ss.n_active == 3
+    assert ss.n_honest == 2            # slots 0 and 2
+    assert ss.n_recovered == 1         # slot 1 recovered via uncle
+    assert abs(ss.q - 2 / 3) < 1e-9
+    assert abs(ss.q_eff - 1.0) < 1e-9
+
+
+def test_update_D_fixed_point():
+    f, T = 1 / 30, 3000
+    m = int(round(T * f))              # measured density == f -> D unchanged
+    assert abs(update_D(1000.0, m, T, f, beta=1.0) - 1000.0) < 1e-6
+    # measured below f -> estimate drops
+    assert update_D(1000.0, m - 20, T, f, 1.0) < 1000.0
+    # clamp at 1
+    assert update_D(1.0, 0, T, f, 1.0) == 1.0
+
+
+def test_update_D_fixed_point_mode_targets_truncated_f():
+    # In fixed-point mode the target rate is int(f*1000)/1000 = 0.033, not 1/30.
+    f, T = 1 / 30, 10000
+    m = 340
+    # For the same measured density, fixed-point (lower target 0.033) raises the estimate
+    # more than exact-f, i.e. it is systematically higher.
+    assert (update_D(1000.0, m, T, f, 1.0, fixed_point=True)
+            > update_D(1000.0, m, T, f, 1.0, fixed_point=False))
+    # A density of exactly 0.033 is the fixed-point fixed point (estimate unchanged).
+    m_trunc = int(round(0.033 * T))    # 330
+    assert abs(update_D(1000.0, m_trunc, T, f, 1.0, fixed_point=True) - 1000.0) < 1e-6
diff --git a/tools/simulators/tsi/tsi-sim-mc/tests/test_uncles.py b/tools/simulators/tsi/tsi-sim-mc/tests/test_uncles.py
new file mode 100644
index 0000000..7e8b244
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim-mc/tests/test_uncles.py
@@ -0,0 +1,110 @@
+import numpy as np
+
+from tsi_sim.blocktree import BlockTree
+from tsi_sim.config import SimConfig
+from tsi_sim.uncles import annotate_uncles
+
+
+def make_tree(slots, parents, heights, leaders):
+    n = len(slots)
+    return BlockTree(
+        slot=np.array(slots, np.int64),
+        parent=np.array(parents, np.int64),
+        height=np.array(heights, np.int64),
+        leader=np.array(leaders, np.int64),
+        uncles=[() for _ in range(n)],
+    )
+
+
+def _canonical_and_orphan_tree():
+    # genesis(0); canonical chain 1(slot0)->3(slot3)->4(slot5); orphan 2(slot1)
+    tree = make_tree(
+        slots=[-1, 0, 1, 3, 5],
+        parents=[-1, 0, 0, 1, 3],
+        heights=[0, 1, 1, 2, 3],
+        leaders=[-1, 0, 1, 2, 3],
+    )
+    canonical = [4, 3, 1]  # tip-first
+    return tree, canonical
+
+
+def test_oldest_selection_and_window():
+    tree, canonical = _canonical_and_orphan_tree()
+    cfg = SimConfig(max_uncles=1, uncle_window=300, uncle_strategy="oldest")
+    annotate_uncles(tree, canonical, cfg, np.random.default_rng(0))
+    # orphan 2 (slot1) is within window of block 3 (slot3) -> referenced there
+    referenced = {u for b in canonical for u in tree.uncles[b]}
+    assert referenced == {2}
+
+
+def test_no_uncles_when_u_zero():
+    tree, canonical = _canonical_and_orphan_tree()
+    cfg = SimConfig(max_uncles=0)
+    annotate_uncles(tree, canonical, cfg, np.random.default_rng(0))
+    assert all(tree.uncles[b] == () for b in canonical)
+
+
+def test_window_excludes_out_of_range_orphan():
+    tree, canonical = _canonical_and_orphan_tree()
+    cfg = SimConfig(max_uncles=1, uncle_window=1, uncle_strategy="oldest")
+    annotate_uncles(tree, canonical, cfg, np.random.default_rng(0))
+    # orphan 2 at slot1; nearest canonical after it is block3 at slot3 -> gap 2 > W=1
+    referenced = {u for b in canonical for u in tree.uncles[b]}
+    assert referenced == set()
+
+
+def _wide_orphan_tree():
+    # canonical 1(0)->6(6); orphans 2,3,4,5 at slots 1,2,3,4 (all within window of block6)
+    tree = make_tree(
+        slots=[-1, 0, 1, 2, 3, 4, 6],
+        parents=[-1, 0, 0, 0, 0, 0, 1],
+        heights=[0, 1, 1, 1, 1, 1, 2],
+        leaders=[-1, 0, 1, 2, 3, 4, 0],
+    )
+    return tree, [6, 1]  # tip-first
+
+
+def test_random_strategy_deterministic_and_capped():
+    import numpy as np
+    tree_a, canon = _wide_orphan_tree()
+    tree_b, _ = _wide_orphan_tree()
+    cfg = SimConfig(max_uncles=2, uncle_window=300, uncle_strategy="random", uncle_random_p=0.5)
+    annotate_uncles(tree_a, canon, cfg, np.random.default_rng(7))
+    annotate_uncles(tree_b, canon, cfg, np.random.default_rng(7))
+    assert tree_a.uncles == tree_b.uncles                 # same seed -> identical
+    total = sum(len(tree_a.uncles[b]) for b in canon)
+    assert total <= cfg.max_uncles                        # capped
+
+
+def test_random_p_one_matches_oldest():
+    import numpy as np
+    tree_r, canon = _wide_orphan_tree()
+    tree_o, _ = _wide_orphan_tree()
+    annotate_uncles(tree_r, canon, SimConfig(max_uncles=2, uncle_strategy="random",
+                                             uncle_random_p=1.0), np.random.default_rng(1))
+    annotate_uncles(tree_o, canon, SimConfig(max_uncles=2, uncle_strategy="oldest"),
+                    np.random.default_rng(1))
+    assert tree_r.uncles == tree_o.uncles     # p=1 deterministically takes oldest-first
+
+
+def test_random_p_zero_selects_nothing():
+    import numpy as np
+    tree, canon = _wide_orphan_tree()
+    annotate_uncles(tree, canon, SimConfig(max_uncles=4, uncle_strategy="random",
+                                           uncle_random_p=0.0), np.random.default_rng(1))
+    assert all(tree.uncles[b] == () for b in canon)
+
+
+def test_dedup_across_ancestors():
+    # Two canonical blocks both within window of the single orphan: only one references it.
+    tree = make_tree(
+        slots=[-1, 0, 1, 2, 3],
+        parents=[-1, 0, 0, 1, 3],
+        heights=[0, 1, 1, 2, 3],
+        leaders=[-1, 0, 9, 2, 3],
+    )
+    canonical = [4, 3, 1]
+    cfg = SimConfig(max_uncles=4, uncle_window=300, uncle_strategy="oldest")
+    annotate_uncles(tree, canonical, cfg, np.random.default_rng(0))
+    counts = sum(len(tree.uncles[b]) for b in canonical)
+    assert counts == 1  # orphan 2 referenced exactly once despite two eligible blocks