diff --git a/tools/simulators/blend/pd/src/pd/config.py b/tools/simulators/blend/pd/src/pd/config.py
index 7c663b8..faaa2c6 100644
--- a/tools/simulators/blend/pd/src/pd/config.py
+++ b/tools/simulators/blend/pd/src/pd/config.py
@@ -34,6 +34,13 @@ class SimConfig:
     unresponsive_frac: float = 0.0      # ratio of nodes that do NOT relay any messages (swept)
     churn_mode: ChurnMode = "uniform"   # "uniform" = independent nodes; "regional" = whole regions
     redundancy: int = 1                 # copies per emission via R independent cascades (swept)
+
+    # --- cover traffic (1 slot = 1 second) ---
+    cover_rate_mult: float = 1.0        # emissions/slot network-wide; 1.0 = one per second
+    block_interval_slots: int = 30      # slots per block, i.e. f = 1/30
+    slots_per_epoch: int = 648_000      # epoch length, for the per-node emission quota
+    stake_inference_ratio: float = 1.0  # D_hat/D from the consensus study; scales the stake ceiling
+    traffic_window_slots: int = 600     # simulated timeline length (seconds) per cell
     n_rounds: int = 200                 # random-sender rounds per topology
     transport_jitter_mean_ms: float = 5.0
     processing_lags_ms: tuple[float, ...] = (10.0, 50.0, 100.0)
@@ -66,6 +73,15 @@ class SimConfig:
             raise ValueError(f"need 0 <= unresponsive_frac < 1, got {self.unresponsive_frac}")
         if self.redundancy < 1:
             raise ValueError(f"redundancy must be >= 1, got {self.redundancy}")
+        if self.cover_rate_mult < 0.0:
+            raise ValueError(f"cover_rate_mult must be >= 0, got {self.cover_rate_mult}")
+        if self.block_interval_slots < 2:
+            raise ValueError(f"block_interval_slots must be >= 2, got {self.block_interval_slots}")
+        if self.slots_per_epoch < 1 or self.traffic_window_slots < 1:
+            raise ValueError("slots_per_epoch and traffic_window_slots must be >= 1")
+        if self.stake_inference_ratio <= 0.0:
+            raise ValueError(
+                f"stake_inference_ratio (D_hat/D) must be > 0, got {self.stake_inference_ratio}")
         if self.n_regions < 1:
             raise ValueError(f"n_regions must be >= 1, got {self.n_regions}")
         if self.n_regions > 1:
@@ -114,7 +130,9 @@ class SimConfig:
         return (
             self.n_nodes, self.degree, self.n_regions, self.region_locality,
             self.blend_hops, self.max_blend_delay,
-            self.unresponsive_frac, self.churn_mode, self.redundancy, self.n_rounds,
+            self.unresponsive_frac, self.churn_mode, self.redundancy,
+            self.cover_rate_mult, self.block_interval_slots, self.slots_per_epoch,
+            self.stake_inference_ratio, self.traffic_window_slots, self.n_rounds,
             self.transport_jitter_mean_ms, self.processing_lags_ms, self.processing_lag_probs,
             self.link_latency_dist, self.link_latency_mean_ms, self.coverage_pcts,
             self.f_adv, self.adversary_mode, self.n_placements, self.worstcase_max_n,
diff --git a/tools/simulators/blend/pd/src/pd/quota.py b/tools/simulators/blend/pd/src/pd/quota.py
new file mode 100644
index 0000000..0ed222c
--- /dev/null
+++ b/tools/simulators/blend/pd/src/pd/quota.py
@@ -0,0 +1,120 @@
+"""Emission budget: how much stake a node can hold and still look like everyone else.
+
+Cover traffic gives every node the **same number of emissions per epoch**. A node emits in a slot
+either because a cover slot came up or because it won the block lottery, and a block cancels the
+next scheduled cover, so the emission *count* carries no information about who produces blocks --
+that uniformity is the anonymity property cover traffic buys.
+
+It only holds while a node's block proposals fit inside its quota. The quota is
+
+    q = cover_rate_mult / n_nodes          emissions per node per slot
+
+and the Cryptarchia lottery gives a node of **inferred** relative stake ``alpha = sigma / D_hat``
+a per-slot win probability ``phi(alpha) = 1 - (1-f)**alpha``. Requiring ``phi(alpha) <= q``:
+
+    alpha_max = ln(1 - q) / ln(1 - f)
+
+Note the stake here is relative to the *inferred* total ``D_hat``, not the true total ``D``, because
+that is the denominator the lottery threshold is derived from. Converting to true relative stake
+``s = sigma / D`` therefore needs the estimator's accuracy ratio: since every node's win rate scales
+with ``D/D_hat``, the network's block rate is ``f * D/D_hat`` and
+
+    s_max = (D_hat / D) * alpha_max
+
+so an estimator that runs low tightens the true-stake ceiling in exact proportion. ``D_hat/D`` is an
+input here (``stake_inference_ratio``), measured by the consensus-side study, not assumed to be 1.
+
+``alpha_max`` is the *mean* bind. Wins are a Bernoulli process, so a node sitting exactly at it
+overruns its quota in about half of all epochs; ``max_alpha_for_confidence`` gives the ceiling that
+keeps the quota with a stated probability.
+"""
+
+from __future__ import annotations
+
+import math
+
+
+def emission_quota_per_slot(n_nodes: int, cover_rate_mult: float = 1.0) -> float:
+    """Emissions allowed per node per slot. The default rate puts one emission per slot on the
+    whole network, so each node gets ``1/n_nodes``."""
+    if n_nodes < 1:
+        raise ValueError("n_nodes must be >= 1")
+    return cover_rate_mult / n_nodes
+
+
+def win_prob(alpha: float, f: float) -> float:
+    """Cryptarchia per-slot block-lottery probability for inferred relative stake ``alpha``."""
+    if not (0.0 < f < 1.0):
+        raise ValueError("need 0 < f < 1")
+    return 1.0 - (1.0 - f) ** alpha
+
+
+def alpha_max(n_nodes: int, f: float, cover_rate_mult: float = 1.0) -> float:
+    """Largest **inferred** relative stake whose expected block rate fits the emission quota.
+
+    Exact bind ``ln(1-q)/ln(1-f)``; the familiar ``q/f`` is a small-q approximation that runs
+    about 1.7 % high at f = 1/30 and so overstates the tolerable stake.
+    """
+    q = emission_quota_per_slot(n_nodes, cover_rate_mult)
+    if q >= 1.0:
+        return math.inf                      # quota exceeds one emission per slot: never binds
+    return math.log1p(-q) / math.log1p(-f)
+
+
+def s_max_true(n_nodes: int, f: float, stake_inference_ratio: float = 1.0,
+               cover_rate_mult: float = 1.0) -> float:
+    """The ceiling expressed in **true** relative stake, ``(D_hat/D) * alpha_max``."""
+    if stake_inference_ratio <= 0.0:
+        raise ValueError("stake_inference_ratio must be > 0")
+    return stake_inference_ratio * alpha_max(n_nodes, f, cover_rate_mult)
+
+
+def expected_blocks_per_epoch(alpha: float, f: float, slots_per_epoch: int) -> float:
+    """Expected proposals won by a node of inferred relative stake ``alpha`` over one epoch."""
+    return win_prob(alpha, f) * slots_per_epoch
+
+
+def quota_per_epoch(n_nodes: int, slots_per_epoch: int, cover_rate_mult: float = 1.0) -> float:
+    """Emissions a node is allowed in one epoch."""
+    return emission_quota_per_slot(n_nodes, cover_rate_mult) * slots_per_epoch
+
+
+def quota_exceedance_prob(alpha: float, f: float, n_nodes: int, slots_per_epoch: int,
+                          cover_rate_mult: float = 1.0) -> float:
+    """P(a node of inferred stake ``alpha`` wins more proposals in an epoch than its quota).
+
+    Wins are Binomial(slots_per_epoch, phi(alpha)); the Poisson limit is used, which is accurate
+    here because phi is tiny and the epoch is long.
+    """
+    lam = expected_blocks_per_epoch(alpha, f, slots_per_epoch)
+    quota = quota_per_epoch(n_nodes, slots_per_epoch, cover_rate_mult)
+    k = math.floor(quota)
+    # P(X > k) for X ~ Poisson(lam), summed up from 0 (k is small in every regime of interest)
+    if lam <= 0.0:
+        return 0.0
+    term = math.exp(-lam)
+    cdf = term
+    for i in range(1, k + 1):
+        term *= lam / i
+        cdf += term
+    return max(0.0, min(1.0, 1.0 - cdf))
+
+
+def max_alpha_for_confidence(f: float, n_nodes: int, slots_per_epoch: int,
+                             confidence: float = 0.99,
+                             cover_rate_mult: float = 1.0) -> float:
+    """Largest inferred stake that keeps inside the quota with probability ``confidence``.
+
+    Always below :func:`alpha_max`, because a node sitting on the mean bind overruns half the time.
+    """
+    if not (0.0 < confidence < 1.0):
+        raise ValueError("need 0 < confidence < 1")
+    lo, hi = 0.0, alpha_max(n_nodes, f, cover_rate_mult)
+    if not math.isfinite(hi):
+        return hi
+    tol = 1.0 - confidence
+    for _ in range(80):                      # bisection on a monotone exceedance probability
+        mid = 0.5 * (lo + hi)
+        p = quota_exceedance_prob(mid, f, n_nodes, slots_per_epoch, cover_rate_mult)
+        lo, hi = (mid, hi) if p <= tol else (lo, mid)
+    return lo
diff --git a/tools/simulators/blend/pd/src/pd/traffic.py b/tools/simulators/blend/pd/src/pd/traffic.py
new file mode 100644
index 0000000..84c0168
--- /dev/null
+++ b/tools/simulators/blend/pd/src/pd/traffic.py
@@ -0,0 +1,245 @@
+"""Cover traffic on a timeline: blending, mixing, and what the network sees between blocks.
+
+The rest of the simulator samples *independent* rounds, drawing each relay's hold from
+``mixclock.mix_wait`` -- the stationary residual to that relay's next release. That is exact for one
+message meeting an independent clock, but it has no notion of time, so two messages can never meet
+at the same relay. Everything this module measures is defined by them meeting, so here each node
+owns **one free-running clock** shared by every message that passes through it, and messages are
+played out on a real timeline.
+
+Two quantities, easy to confuse:
+
+* **mixing** -- how many messages a relay is *holding* at once. A relay holds a message from its
+  arrival until that relay's next tick.
+* **blending** -- how many messages a relay has *seen* between two consecutive releases. Every
+  broadcast reaches every node, so a relay sees the whole network's broadcast stream; an observer
+  watching it release cannot tell which of those it forwarded. This, not the held count, is the
+  anonymity set, and it grows with the broadcast rate and with the release interval:
+  ``blending ~ rate * interval`` where the interval of a ``Uniform{0..M}`` clock averages ``M/2``.
+
+Emissions are the union of cover traffic and block proposals. Each node emits at rate
+``cover_rate_mult / n_nodes`` per slot; winning the block lottery consumes the next scheduled cover,
+so every node emits the same number of times whether or not it produces blocks (see ``quota`` for
+the stake ceiling that keeps that true). Cover and block messages travel the same cascade over
+independently drawn paths, and both end in a broadcast, so they are indistinguishable in transit.
+"""
+
+from __future__ import annotations
+
+from dataclasses import dataclass, field
+
+import numpy as np
+from scipy.sparse.csgraph import dijkstra
+
+from .config import SimConfig
+from .graph import Graph
+from .mixclock import mix_wait
+
+
+class ReleaseClock:
+    """One node's free-running release clock: ticks separated by ``Uniform{0..M}`` whole seconds.
+
+    The first tick is the stationary residual from t=0, so a clock sampled once reproduces
+    ``mixclock.mix_wait``. Ticks are extended lazily, which keeps a million-node network cheap --
+    only the relays a message actually visits ever grow a clock.
+    """
+
+    __slots__ = ("_m", "_rng", "_ticks")
+
+    def __init__(self, max_blend_delay: int, rng: np.random.Generator) -> None:
+        self._m = int(max_blend_delay)
+        self._rng = rng
+        first = 0.0 if self._m <= 0 else float(mix_wait(rng, self._m, 1)[0]) / 1000.0
+        self._ticks: list[float] = [first]
+
+    def next_tick_at_or_after(self, t: float) -> float:
+        """Time of this clock's first tick at or after ``t`` (``t`` itself when M = 0)."""
+        if self._m <= 0:
+            return t
+        while self._ticks[-1] < t:
+            step = float(self._rng.integers(0, self._m + 1))
+            self._ticks.append(self._ticks[-1] + step)
+        for tick in self._ticks:                 # tick lists stay short (window / mean interval)
+            if tick >= t:
+                return tick
+        return self._ticks[-1]
+
+    def ticks_in(self, lo: float, hi: float) -> list[float]:
+        """Ticks in ``(lo, hi]`` -- the release opportunities used to bracket blending."""
+        self.next_tick_at_or_after(hi)
+        return [t for t in self._ticks if lo < t <= hi]
+
+
+@dataclass
+class Hold:
+    """One message waiting at one relay."""
+    node: int
+    arrived: float
+    released: float
+
+
+@dataclass
+class TrafficWindow:
+    """Raw event record of a simulated window; :func:`traffic_metrics` reduces it."""
+    emitted_cover: int = 0
+    emitted_block: int = 0
+    cancelled_cover: int = 0
+    broadcasts: list[float] = field(default_factory=list)       # when each message floods
+    holds: list[Hold] = field(default_factory=list)
+    window_seconds: float = 0.0
+    block_slots: list[float] = field(default_factory=list)      # when blocks were proposed
+    cover_slots: list[float] = field(default_factory=list)      # when cover was emitted
+    clocks: dict[int, ReleaseClock] = field(default_factory=dict)
+
+
+def _clock(clocks: dict[int, ReleaseClock], node: int, max_blend_delay: int,
+           rng: np.random.Generator) -> ReleaseClock:
+    c = clocks.get(node)
+    if c is None:
+        c = ReleaseClock(max_blend_delay, rng)
+        clocks[node] = c
+    return c
+
+
+def simulate_window(graph: Graph, config: SimConfig, rng: np.random.Generator,
+                    window_slots: int, max_blend_delay: int | None = None,
+                    blend_hops: int | None = None) -> TrafficWindow:
+    """Play ``window_slots`` seconds of network traffic and record every hold and broadcast.
+
+    Each slot, the network emits once per ``cover_rate_mult`` on average -- the emitter is uniform
+    because every node carries the same per-slot rate. A block proposal is drawn at the network
+    block rate and, per the quota rule, cancels that node's next cover emission.
+    """
+    n = graph.n
+    k = int(config.blend_hops if blend_hops is None else blend_hops)
+    m = int(config.max_blend_delay if max_blend_delay is None else max_blend_delay)
+    f = 1.0 / config.block_interval_slots
+    win = TrafficWindow(window_seconds=float(window_slots))
+    clocks = win.clocks
+    cancelled: set[int] = set()          # nodes owing a cancelled cover after a block proposal
+
+    for slot in range(window_slots):
+        t0 = float(slot)
+        n_emissions = rng.poisson(config.cover_rate_mult)      # network-wide emissions this slot
+        block_this_slot = rng.random() < f
+        senders = rng.integers(0, n, size=n_emissions).tolist() if n_emissions else []
+        if block_this_slot:
+            senders.append(int(rng.integers(0, n)))            # the proposer also emits
+        for i, sender in enumerate(senders):
+            is_block = block_this_slot and i == len(senders) - 1
+            if not is_block and sender in cancelled:
+                cancelled.discard(sender)                      # this cover is the one forfeited
+                win.cancelled_cover += 1
+                continue
+            if is_block:
+                cancelled.add(sender)
+                win.emitted_block += 1
+                win.block_slots.append(t0)
+            else:
+                win.emitted_cover += 1
+                win.cover_slots.append(t0)
+            relays = rng.choice(n - 1, size=k, replace=False)
+            relays[relays >= sender] += 1
+            sources = np.empty(k + 1, dtype=np.int64)
+            sources[0] = sender
+            sources[1:] = relays
+            data = (graph.base + rng.exponential(config.transport_jitter_mean_ms,
+                                                 size=graph.base.shape[0]) + graph.p[graph.src])
+            dist = dijkstra(graph.weighted_csr(data), directed=True, indices=sources)
+            t = t0
+            dropped = False
+            for hop in range(k):
+                leg = dist[hop, relays[hop]]
+                if not np.isfinite(leg):
+                    dropped = True
+                    break
+                arrived = t + float(leg) / 1000.0
+                released = _clock(clocks, int(relays[hop]), m, rng).next_tick_at_or_after(arrived)
+                win.holds.append(Hold(int(relays[hop]), arrived, released))
+                t = released
+            if not dropped:
+                win.broadcasts.append(t)                       # the final relay floods at t
+    return win
+
+
+def traffic_metrics(win: TrafficWindow, config: SimConfig,
+                    max_blend_delay: int | None = None) -> dict:
+    """Reduce a window to the reported quantities.
+
+    ``blending`` is evaluated per release: every broadcast in the network is seen by every node, so
+    the anonymity set of a release at time ``T`` is the number of broadcasts since that relay's
+    previous release.
+    """
+    m = int(config.max_blend_delay if max_blend_delay is None else max_blend_delay)
+    # cover emitted between consecutive block proposals, measured gap by gap
+    if len(win.block_slots) >= 2:
+        cover = np.sort(np.asarray(win.cover_slots))
+        blocks = np.sort(np.asarray(win.block_slots))
+        gaps = np.diff(np.searchsorted(cover, blocks, "right")).astype(float)
+        per_block = float(gaps.mean())
+    else:
+        per_block = float("nan")
+    out = {
+        "emitted_cover": win.emitted_cover,
+        "emitted_block": win.emitted_block,
+        "cancelled_cover": win.cancelled_cover,
+        "broadcasts_seen": len(win.broadcasts),
+        "cover_per_block_interval": per_block,
+        "hold_events": len(win.holds),
+    }
+    if not win.holds:
+        out.update(hold_seconds_mean=float("nan"), delayed_frac=0.0,
+                   queue_mean=0.0, queue_p90=0.0, queue_max=0.0,
+                   blending_mean=float("nan"), blending_p50=float("nan"),
+                   blending_p90=float("nan"), blending_max=float("nan"))
+        return out
+    hold_s = np.array([h.released - h.arrived for h in win.holds])
+    out["hold_seconds_mean"] = float(hold_s.mean())
+    out["delayed_frac"] = float(np.mean(hold_s > 0.0))
+
+    # mixing: concurrent holds at a relay, time-weighted over the window
+    by_node: dict[int, list[Hold]] = {}
+    for h in win.holds:
+        by_node.setdefault(h.node, []).append(h)
+    occupancy, peaks = [], []
+    for holds in by_node.values():
+        events = sorted([(h.arrived, 1) for h in holds] + [(h.released, -1) for h in holds])
+        cur = peak = 0
+        area = 0.0
+        prev = events[0][0]
+        for tstamp, delta in events:
+            area += cur * (tstamp - prev)
+            prev = tstamp
+            cur += delta
+            peak = max(peak, cur)
+        occupancy.append(area / win.window_seconds)
+        peaks.append(peak)
+    out["queue_mean"] = float(np.sum(occupancy) / max(len(by_node), 1))
+    out["queue_p90"] = float(np.percentile(peaks, 90))
+    out["queue_max"] = float(np.max(peaks))
+
+    # Blending: the anonymity set of a release. A timed-release relay would have flushed anything
+    # that arrived before its previous tick, so the message it forwards must have arrived in the
+    # last inter-tick interval -- and since every broadcast reaches every node, the candidates are
+    # all broadcasts in that interval. Uses each clock's real previous tick, not the mean interval.
+    casts = np.sort(np.asarray(win.broadcasts))
+    sets: list[int] = []
+    for node, holds in by_node.items():
+        clock = win.clocks.get(node)
+        for r in sorted({h.released for h in holds}):
+            if m <= 0:
+                prev = r                                  # no delay: the set is whatever is instant
+            elif clock is not None:
+                earlier = [t for t in clock.ticks_in(-1e18, r) if t < r]
+                prev = earlier[-1] if earlier else max(0.0, r - m / 2.0)
+            else:
+                prev = max(0.0, r - m / 2.0)
+            # strictly before the release: when this relay is the last hop it broadcasts at exactly
+            # `r`, and that is its own output, not a candidate input it saw.
+            sets.append(int(np.searchsorted(casts, r, "left")
+                            - np.searchsorted(casts, prev, "right")))
+    if sets:
+        arr = np.asarray(sets, dtype=float)
+        out.update(blending_mean=float(arr.mean()), blending_p50=float(np.percentile(arr, 50)),
+                   blending_p90=float(np.percentile(arr, 90)), blending_max=float(arr.max()))
+    return out
diff --git a/tools/simulators/blend/pd/tests/test_quota.py b/tools/simulators/blend/pd/tests/test_quota.py
new file mode 100644
index 0000000..4becf61
--- /dev/null
+++ b/tools/simulators/blend/pd/tests/test_quota.py
@@ -0,0 +1,85 @@
+"""The emission-quota stake ceiling: exact bind, the D_hat/D normalisation, and epoch compliance."""
+
+import math
+
+import numpy as np
+
+from pd.quota import (
+    alpha_max,
+    emission_quota_per_slot,
+    expected_blocks_per_epoch,
+    max_alpha_for_confidence,
+    quota_exceedance_prob,
+    quota_per_epoch,
+    s_max_true,
+    win_prob,
+)
+
+F = 1.0 / 30.0
+
+
+def test_quota_is_one_emission_per_slot_network_wide():
+    n = 20_000
+    assert emission_quota_per_slot(n) * n == 1.0            # whole network emits once per slot
+    assert emission_quota_per_slot(n, 4.0) * n == 4.0       # the multiplier scales it
+
+
+def test_alpha_max_is_where_the_win_rate_equals_the_quota():
+    n = 20_000
+    a = alpha_max(n, F)
+    assert abs(win_prob(a, F) - emission_quota_per_slot(n)) < 1e-15
+
+
+def test_alpha_max_is_below_the_q_over_f_approximation():
+    """q/f is a small-q expansion and errs optimistic, so the exact bind must be lower."""
+    for n in (1_000, 20_000, 10**6):
+        exact = alpha_max(n, F)
+        approx = emission_quota_per_slot(n) / F
+        assert exact < approx
+        assert abs(approx / exact - 1) < 0.02               # ~1.7% at f = 1/30
+
+
+def test_alpha_max_scales_inversely_with_network_size_and_with_cover_rate():
+    assert abs(alpha_max(20_000, F) / alpha_max(200_000, F) - 10.0) < 0.01
+    assert abs(alpha_max(20_000, F, 8.0) / alpha_max(20_000, F, 1.0) - 8.0) < 0.01
+
+
+def test_true_stake_ceiling_is_scaled_by_the_inference_ratio():
+    """The lottery uses sigma/D_hat, so the ceiling in TRUE stake carries the D_hat/D factor."""
+    n = 20_000
+    a = alpha_max(n, F)
+    assert s_max_true(n, F, 1.0) == a                       # accurate estimator: no correction
+    assert abs(s_max_true(n, F, 0.74) - 0.74 * a) < 1e-15   # deflated estimate tightens it
+    assert s_max_true(n, F, 0.64) < s_max_true(n, F, 0.74) < a
+
+
+def test_expected_blocks_equal_the_quota_at_alpha_max():
+    n, S = 20_000, 648_000
+    a = alpha_max(n, F)
+    assert abs(expected_blocks_per_epoch(a, F, S) - quota_per_epoch(n, S)) < 1e-6
+
+
+def test_a_node_at_the_mean_bind_overruns_about_half_the_time():
+    n, S = 20_000, 648_000
+    p = quota_exceedance_prob(alpha_max(n, F), F, n, S)
+    assert 0.35 < p < 0.65                                  # mean bind is a coin flip, as expected
+
+
+def test_confidence_ceiling_is_stricter_than_the_mean_bind():
+    n, S = 20_000, 648_000
+    safe = max_alpha_for_confidence(F, n, S, confidence=0.99)
+    assert safe < alpha_max(n, F)
+    assert quota_exceedance_prob(safe, F, n, S) <= 0.01 + 1e-9
+    assert 0.5 < safe / alpha_max(n, F) < 0.9               # Poisson noise eats real headroom
+
+
+def test_exceedance_matches_a_direct_simulation():
+    """Closed-form exceedance vs drawing epochs of block wins."""
+    n, S = 2_000, 20_000
+    a = alpha_max(n, F) * 0.8
+    closed = quota_exceedance_prob(a, F, n, S)
+    rng = np.random.default_rng(0)
+    quota = quota_per_epoch(n, S)
+    wins = rng.binomial(S, win_prob(a, F), size=20_000)
+    emp = float(np.mean(wins > math.floor(quota)))
+    assert abs(closed - emp) < max(0.01, 0.1 * closed)
diff --git a/tools/simulators/blend/pd/tests/test_traffic.py b/tools/simulators/blend/pd/tests/test_traffic.py
new file mode 100644
index 0000000..fc57f31
--- /dev/null
+++ b/tools/simulators/blend/pd/tests/test_traffic.py
@@ -0,0 +1,115 @@
+"""Cover traffic on a timeline: shared clocks, the emission quota, blending and mixing."""
+
+import numpy as np
+
+from pd.config import SimConfig
+from pd.graph import build_graph
+from pd.traffic import ReleaseClock, simulate_window, traffic_metrics
+
+
+def _win(n_nodes=2000, degree=8, hops=3, M=3, mult=1.0, slots=600, seed=0):
+    cfg = SimConfig(n_nodes=n_nodes, degree=degree, blend_hops=hops, max_blend_delay=M,
+                    cover_rate_mult=mult)
+    g = build_graph(cfg)
+    w = simulate_window(g, cfg, np.random.default_rng(seed), window_slots=slots)
+    return w, traffic_metrics(w, cfg)
+
+
+# --- the clock ----------------------------------------------------------------------------------
+
+def test_clock_ticks_are_monotonic_and_spaced_within_the_bound():
+    c = ReleaseClock(3, np.random.default_rng(0))
+    c.next_tick_at_or_after(100.0)
+    ticks = c._ticks
+    assert all(b >= a for a, b in zip(ticks, ticks[1:], strict=False))
+    gaps = [b - a for a, b in zip(ticks, ticks[1:], strict=False)]
+    assert all(0 <= g <= 3 + 1e-9 for g in gaps)
+
+
+def test_clock_with_zero_delay_releases_immediately():
+    c = ReleaseClock(0, np.random.default_rng(0))
+    for t in (0.0, 1.5, 99.0):
+        assert c.next_tick_at_or_after(t) == t
+
+
+def test_next_tick_is_at_or_after_the_request_and_is_stable():
+    c = ReleaseClock(3, np.random.default_rng(1))
+    for t in (0.3, 5.0, 5.0, 12.7):
+        assert c.next_tick_at_or_after(t) >= t
+    assert c.next_tick_at_or_after(5.0) == c.next_tick_at_or_after(5.0)   # idempotent
+
+
+def test_one_clock_is_shared_so_messages_batch_at_the_same_tick():
+    """Two messages arriving before the same tick leave together -- that is the mixing."""
+    c = ReleaseClock(3, np.random.default_rng(2))
+    t1 = c.next_tick_at_or_after(10.0)
+    t2 = c.next_tick_at_or_after(10.0 + 1e-6)
+    assert t1 == t2 or t2 >= t1
+
+
+def test_first_tick_reproduces_the_stationary_residual():
+    """A clock sampled once matches mixclock's residual, so single-message stats are unchanged."""
+    M = 5
+    firsts = [ReleaseClock(M, np.random.default_rng(s))._ticks[0] for s in range(4000)]
+    assert abs(float(np.mean(firsts)) - (2 * M + 1) / 6) < 0.06        # mean residual (2M+1)/6
+
+
+# --- emissions and the quota --------------------------------------------------------------------
+
+def test_block_proposals_cancel_a_later_cover_emission():
+    w, _ = _win(slots=1500, seed=3)
+    assert w.emitted_block > 0
+    assert w.cancelled_cover > 0
+    # every cancellation is owed to a block, and cannot exceed the blocks emitted
+    assert w.cancelled_cover <= w.emitted_block
+
+
+def test_cover_between_blocks_matches_the_rate_times_the_block_interval():
+    w, m = _win(slots=3000, seed=4)
+    rate = (w.emitted_cover + w.emitted_block) / w.window_seconds
+    assert abs(m["cover_per_block_interval"] - rate * 30) < 6          # ~30 at 1 msg/s
+
+
+# --- what the relays experience -------------------------------------------------------------------
+
+def test_mean_hold_is_the_renewal_residual():
+    for M in (3, 10):
+        _, m = _win(M=M, slots=900, seed=5)
+        assert abs(m["hold_seconds_mean"] - (2 * M + 1) / 6) < 0.25
+
+
+def test_blending_follows_the_size_biased_interval():
+    """Anonymity set = broadcasts seen in the last inter-tick gap. Gaps sampled at a release are
+    size-biased, so the mean is rate*(2M+1)/3 -- twice the mean hold, not rate*M/2."""
+    for M in (3, 10, 30):
+        w, m = _win(M=M, slots=1200, seed=1)
+        rate = (w.emitted_cover + w.emitted_block) / w.window_seconds
+        assert abs(m["blending_mean"] - rate * (2 * M + 1) / 3) < 0.12 * rate * (2 * M + 1) / 3
+
+
+def test_blending_grows_with_the_cover_rate_and_with_the_delay():
+    _, lo = _win(M=3, mult=1.0, slots=400, seed=2)
+    _, hi = _win(M=3, mult=8.0, slots=400, seed=2)
+    assert hi["blending_mean"] > 5 * lo["blending_mean"]               # ~linear in rate
+    _, slow = _win(M=30, mult=1.0, slots=400, seed=2)
+    assert slow["blending_mean"] > 4 * lo["blending_mean"]             # ~linear in delay
+
+
+def test_mixing_is_negligible_at_the_baseline_rate():
+    """One message per second over thousands of nodes: a relay essentially never holds two."""
+    _, m = _win(n_nodes=4000, mult=1.0, slots=600, seed=6)
+    assert m["queue_mean"] < 0.05
+    assert m["queue_max"] <= 3
+
+
+def test_mixing_grows_when_the_network_is_loaded():
+    _, lo = _win(n_nodes=500, mult=1.0, slots=400, seed=7)
+    _, hi = _win(n_nodes=500, mult=32.0, slots=400, seed=7)
+    assert hi["queue_max"] > lo["queue_max"]
+    assert hi["queue_mean"] > lo["queue_mean"]
+
+
+def test_every_hop_is_recorded_as_a_hold():
+    w, m = _win(hops=3, slots=300, seed=8)
+    delivered = len(w.broadcasts)
+    assert m["hold_events"] >= 3 * delivered                          # 3 relays per delivered msg