diff --git a/tools/simulators/blend/pd/src/pd/quota.py b/tools/simulators/blend/pd/src/pd/quota.py
index 0810fa3..108f9bf 100644
--- a/tools/simulators/blend/pd/src/pd/quota.py
+++ b/tools/simulators/blend/pd/src/pd/quota.py
@@ -34,6 +34,7 @@ from __future__ import annotations
 import math
 
 import numpy as np
+from scipy.stats import poisson
 
 
 def assign_stake(n_nodes: int, dist: str, rng: np.random.Generator,
@@ -179,19 +180,17 @@ def quota_exceedance_prob(alpha: float, f: float, n_nodes: int, slots_per_epoch:
 
     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.
+
+    Uses scipy's survival function rather than summing the series by hand: ``exp(-lam)`` underflows
+    to zero past ``lam ~ 745``, which silently collapses a hand-rolled CDF to 0 and reports every
+    node as exceeding. That regime is reached as soon as the cover rate is raised (the quota, and
+    with it the tolerable block count, grows in proportion).
     """
     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))
+    return float(poisson.sf(math.floor(quota), lam))
 
 
 def max_alpha_for_confidence(f: float, n_nodes: int, slots_per_epoch: int,
diff --git a/tools/simulators/blend/pd/tests/test_quota.py b/tools/simulators/blend/pd/tests/test_quota.py
index 29aa767..813ccd8 100644
--- a/tools/simulators/blend/pd/tests/test_quota.py
+++ b/tools/simulators/blend/pd/tests/test_quota.py
@@ -157,3 +157,22 @@ def test_more_cover_traffic_raises_the_ceiling():
     rich = simulate_epoch_emissions(s, F, n, S, np.random.default_rng(10), cover_rate_mult=32.0)
     assert rich["compliant_frac"] > lean["compliant_frac"]
     assert rich["max_compliant_stake"] > lean["max_compliant_stake"]
+
+
+def test_exceedance_survives_a_large_quota():
+    """Regression: exp(-lam) underflows past lam ~ 745, which silently collapsed a hand-rolled
+    Poisson CDF to 0 and reported every node as exceeding. Raising the cover rate reaches that
+    regime immediately, so the mean bind must still be a coin flip at every rate."""
+    n, S = 20_000, 648_000
+    for rate in (1.0, 16.0, 64.0, 256.0):
+        p = quota_exceedance_prob(alpha_max(n, F, rate), F, n, S, rate)
+        assert 0.4 < p < 0.6, (rate, p)
+
+
+def test_the_safe_ceiling_approaches_the_mean_bind_as_the_quota_grows():
+    """Poisson noise shrinks relative to the mean, so the headroom needed for confidence shrinks."""
+    n, S = 20_000, 648_000
+    ratios = [max_alpha_for_confidence(F, n, S, 0.99, r) / alpha_max(n, F, r)
+              for r in (1.0, 16.0, 256.0)]
+    assert all(b > a for a, b in zip(ratios, ratios[1:], strict=False))
+    assert ratios[-1] > 0.9