fix: stop D19/A11 Mpemba false-positives on symmetric initial states (#77)

Canonical amplitude damping with the default rho_0 = I/2 (and the excited
state |1><1|) previously classified as A11 F4 CONFIRMED PUBLICATION_GRADE
— a symmetry-protected zero slowest-mode overlap was read as an anomalous
Mpemba effect (issue #68). Three coupled fixes:

- overlap_c1 (D19) now returns the true biorthogonal coefficient
  |<l1,rho0>|/|<l1,r1>| instead of |<l1,rho0>|/||l1|| (only the zero test
  was meaningful before); denominator floored at EPS_DIV.
- Non-triviality guard is_trivial_overlap, wired through diagnose() via a
  new rho_steady_state argument to compute_mpemba_layer. Using the
  eigenPROJECTORS of rho_ss (robust to a degenerate / maximally mixed
  steady state, where an eigenvector basis is arbitrary), a rho_0 whose
  active blocks are disjoint from the slowest mode's cannot populate it —
  trivially fast relaxation, not Mpemba. MpembaResult gains trivial_overlap.
- A single initial state is CANDIDATE-grade, not PUBLICATION_GRADE: A11
  confidence capped at 0.70 since confirmation needs a reference family the
  single-state pipeline does not provide. The README dephasing example
  (rho_ss = I/2, no protecting symmetry sector) stays A11 but as
  CANDIDATE/CONFIRMATION, no longer self-certifying.

The classifier A11 branch now gates on the guarded is_mpemba_candidate flag.

Tests: biorthogonal-coefficient check, guard truth table, a fine-tuned
qutrit true-positive oracle, and canonical false-positive regressions in
tests/test_mpemba.py; updated A11 fixtures in tests/test_classification.py.
Anchors unchanged (21 green); pytest 427 passed; ruff + mypy clean.

Refs #68.


Claude-Session: https://claude.ai/code/session_01StgHHKB6YGxxjTKdWjcCtq

Co-authored-by: Claude <noreply@anthropic.com>

Author: marcohost33-maker

CHANGELOG.md

@@ -35,6 +35,31 @@ adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
   behaviour change; the semantics are now explicit in code.
 
 ### Fixed
+- **D19/A11 Mpemba detector no longer false-positives on trivially symmetric
+  initial states** (issue #68). Three coupled defects were corrected:
+  - `diagnostics/mpemba.py::overlap_c1` normalised the slowest-mode overlap by
+    `||l_1||`, so only its *zero test* was meaningful. It now returns the true
+    **biorthogonal** expansion coefficient `|<l_1, rho_0>| / |<l_1, r_1>|`
+    (denominator floored at `EPS_DIV`), the physical weight of the slow mode.
+  - A **non-triviality guard** (`is_trivial_overlap`, wired through `diagnose()`
+    via the new `rho_steady_state` argument of `compute_mpemba_layer`) now
+    distinguishes a *symmetry-protected* zero overlap from an anomalous skip. In
+    the sector decomposition set by the **eigenprojectors of `rho_ss`** (robust
+    to a degenerate / maximally mixed steady state, unlike an eigenvector basis),
+    a `rho_0` whose active blocks are disjoint from the slowest mode's can never
+    populate it — trivially fast relaxation, not Mpemba. Amplitude damping with
+    the default `rho_0 = I/2` (and the excited state `|1><1|`) previously
+    returned `A11 F4 CONFIRMED PUBLICATION_GRADE`; both now fall through the A11
+    branch. `MpembaResult` gains a `trivial_overlap` field.
+  - A single initial state skipping the slowest mode is now **CANDIDATE-grade,
+    not PUBLICATION_GRADE**: `A11` confidence is capped at 0.70 because
+    confirmation needs a reference family (e.g. thermal states across
+    temperatures) the single-state pipeline does not provide. The README
+    dephasing example (`rho_ss = I/2`, no protecting symmetry sector) stays
+    `A11` but as `CANDIDATE`/`CONFIRMATION`, no longer self-certifying.
+  - New/updated coverage in `tests/test_mpemba.py` (biorthogonal coefficient,
+    guard truth table, a fine-tuned qutrit true-positive oracle, canonical
+    false-positive regressions) and `tests/test_classification.py`.
 - **CI: pinned `ruff==0.15.20` and reformatted the `_diagnostics` import
   block.** Ruff's isort `I001` heuristic for `lines-after-imports` before a
   module-level constant changed between 0.15.14 and 0.15.20, flagging a

src/liouscope/_diagnostics.py

@@ -132,7 +132,9 @@ def diagnose(
     )
     mpemba: MpembaResult | None = None
     if include_mpemba:
-        mpemba = compute_mpemba_layer(L_super, rho_initial)
+        mpemba = compute_mpemba_layer(
+            L_super, rho_initial, rho_steady_state=rho_steady_state
+        )
 
     classification = classify_mechanism(
         spectral=spectral,

src/liouscope/_types.py

@@ -116,6 +116,7 @@ class MpembaResult:
     overlap_c1: float                     # D19 slowest-mode overlap
     is_mpemba_candidate: bool
     expansion_alpha: float                # D20 Phi_n scaling exponent
+    trivial_overlap: bool = False         # symmetry-protected zero overlap (not Mpemba)
 
 
 @dataclass(frozen=True, slots=True, kw_only=True)

src/liouscope/diagnostics/classification.py

@@ -78,6 +78,8 @@ def _gather_evidence(
     if mpemba is not None:
         ev["mpemba_overlap_c1"] = float(mpemba.overlap_c1)
         ev["mpemba_expansion_alpha"] = float(mpemba.expansion_alpha)
+        ev["mpemba_is_candidate"] = float(mpemba.is_mpemba_candidate)
+        ev["mpemba_trivial_overlap"] = float(mpemba.trivial_overlap)
     return ev
 
 
@@ -87,8 +89,11 @@ def _pick_a_class(
     relaxation: RelaxationResult,
 ) -> tuple[str, str]:
     """Return ``(a_class, f_family)`` based on evidence priorities."""
-    # F4 Mpemba check first (high salience for current literature risk)
-    if ev.get("mpemba_overlap_c1", 1.0) < 1.0e-4:
+    # F4 Mpemba check first (high salience for current literature risk). The
+    # candidate flag already folds in the non-triviality guard (issue #68), so a
+    # symmetry-protected zero overlap (diagonal rho_0 vs a coherence slow mode)
+    # no longer reaches A11 -- only a genuine, fine-tuned skip does.
+    if ev.get("mpemba_is_candidate", 0.0) > 0.5:
         return "A11", "F4"
     # F5 phantom relaxation
     if ev["pseudospectral_radius"] > 2.0 * ev.get("gap_to_gns_ratio", 1.0) and ev["henrici_eta"] > 1.0:
@@ -152,7 +157,13 @@ def _confidence(ev: dict[str, float], a_class: str) -> float:
     elif a_class == "A10":
         score = 0.70
     elif a_class == "A11" and ev.get("mpemba_overlap_c1", 1.0) < 1.0e-5:
-        score = 0.90
+        # A single initial state skipping the slowest mode is a *candidate*, not
+        # a confirmed anomalous Mpemba effect: confirmation needs a reference
+        # family (e.g. thermal states at different temperatures), which the
+        # single-state pipeline does not provide. Capped below the
+        # PUBLICATION_GRADE threshold so a lone overlap cannot self-certify
+        # (issue #68); genuine cases still surface as A11 CANDIDATE.
+        score = 0.70
     elif a_class == "A12":
         score = 0.20
     return float(min(max(score, 0.0), 1.0))

src/liouscope/diagnostics/mpemba.py

@@ -1,9 +1,28 @@
 """Mpemba layer: D19 slowest-mode overlap, D20 expansion-coefficient scaling.
 
 D19 -- :func:`overlap_c1`. The Mpemba-overlap test (Carollo 2021). Decompose
-the initial state in left-eigenvector basis of L and test whether the
+the initial state in the left-eigenvector basis of L and test whether the
 coefficient on the slowest mode (lambda_1 != 0) vanishes. If ``|c_1| << eps``,
-the system relaxes anomalously fast (Mpemba candidate).
+the state relaxes anomalously fast (strong-Mpemba *candidate*).
+
+The coefficient reported is the **biorthogonal** expansion coefficient
+``c_1 = <l_1, rho_0> / <l_1, r_1>`` -- the physically meaningful weight of the
+slowest right-mode ``r_1`` in the decomposition of ``rho_0`` (l_1, r_1 the left
+and right eigenvectors of the slowest non-zero mode). Normalising by
+``<l_1, r_1>`` rather than ``||l_1||`` (the previous behaviour, for which only
+the zero test was meaningful) makes ``c_1`` the actual expansion magnitude.
+
+**Non-triviality guard (LIOU-CLS-002, issue #68).** A vanishing ``c_1`` is only
+an *anomalous* Mpemba skip when a state that could generically excite the
+slowest mode is fine-tuned so it does not. A ``rho_0`` that **commutes with the
+steady state** ``rho_ss`` lies in the classical (population) sector of
+``rho_ss``; if the slowest mode is a *coherence* in the ``rho_ss`` eigenbasis,
+such a ``rho_0`` has ``c_1 = 0`` by symmetry -- it never populates that mode --
+which is trivially fast relaxation, not a Mpemba effect. The default
+``rho_0 = I/d`` and any diagonal state fall in this class. :func:`is_trivial_overlap`
+detects it, and :func:`compute_mpemba_layer` withholds the candidate flag so a
+canonical amplitude-damping / dephasing system cannot earn an anomalous-Mpemba
+label from its default inputs.
 
 D20 -- :func:`expansion_alpha`. Scaling exponent of overlap coefficients
 ``|c_n|`` against the index n. Polynomial scaling ``Phi_n ~ exp(alpha L)``
@@ -15,32 +34,119 @@
 import numpy as np
 import scipy.linalg as sla
 
-from .._consts import EPS_GAP
+from .._consts import EPS_DIV, EPS_GAP, EPS_HERMITICITY
 from .._types import MpembaResult
 from ..numerics.kronecker import vec
 
 
+def _slowest_mode(
+    L_super: np.ndarray, *, atol: float = EPS_GAP
+) -> tuple[np.ndarray, np.ndarray] | None:
+    """Return ``(l_slow, r_slow)`` for the slowest non-zero mode, or ``None``.
+
+    The slowest mode is the non-zero eigenvalue with the largest (least
+    negative) real part. Returns ``None`` when every mode is a zero mode.
+    """
+    L_super = np.asarray(L_super)
+    eigvals, vl, vr = sla.eig(L_super, left=True, right=True)
+    nonzero_mask = np.abs(eigvals) > atol
+    if not np.any(nonzero_mask):
+        return None
+    eigvals_nz = eigvals[nonzero_mask]
+    order = np.argsort(-np.real(eigvals_nz))
+    slowest_idx = np.where(nonzero_mask)[0][order[0]]
+    return vl[:, slowest_idx], vr[:, slowest_idx]
+
+
 def overlap_c1(
     L_super: np.ndarray,
     rho_initial: np.ndarray,
     *,
     atol: float = EPS_GAP,
 ) -> float:
-    """Magnitude of the slowest-mode overlap ``|c_1|``."""
-    L_super = np.asarray(L_super)
-    eigvals, vl, vr = sla.eig(L_super, left=True, right=True)
-    # Filter the zero mode
-    nonzero_mask = np.abs(eigvals) > atol
-    eigvals_nz = eigvals[nonzero_mask]
-    if eigvals_nz.size == 0:
+    """Magnitude of the biorthogonal slowest-mode overlap ``|c_1|``.
+
+    ``c_1 = <l_1, vec(rho_0)> / <l_1, r_1>``; the denominator is floored at
+    :data:`liouscope._consts.EPS_DIV` so a genuinely defective (``<l_1, r_1> = 0``)
+    slow mode yields a large finite value rather than a spurious NaN.
+    """
+    mode = _slowest_mode(L_super, atol=atol)
+    if mode is None:
         return 0.0
-    # Slowest non-zero mode = largest real part (least negative)
-    order = np.argsort(-np.real(eigvals_nz))
-    slowest_idx = np.where(nonzero_mask)[0][order[0]]
-    l_slow = vl[:, slowest_idx]
+    l_slow, r_slow = mode
     rho_vec0 = vec(np.asarray(rho_initial))
-    norm = max(np.linalg.norm(l_slow), 1.0e-12)
-    return float(abs(np.vdot(l_slow, rho_vec0)) / norm)
+    denom = abs(np.vdot(l_slow, r_slow))
+    return float(abs(np.vdot(l_slow, rho_vec0)) / max(denom, EPS_DIV))
+
+
+def _steady_state_eigenprojectors(
+    rho_steady_state: np.ndarray, *, tol: float = EPS_HERMITICITY
+) -> list[np.ndarray]:
+    """Orthogonal projectors onto the distinct eigenspaces of ``rho_ss``.
+
+    Grouping degenerate eigenvalues into a single projector makes the sector
+    decomposition well-defined even when ``rho_ss`` is degenerate or maximally
+    mixed -- the regime where an individual-eigenvector basis is arbitrary.
+    """
+    rho_steady_state = np.asarray(rho_steady_state, dtype=complex)
+    w, u = np.linalg.eigh(rho_steady_state)
+    projectors: list[np.ndarray] = []
+    i = 0
+    n = w.size
+    while i < n:
+        j = i + 1
+        while j < n and abs(w[j] - w[i]) <= tol:
+            j += 1
+        block = u[:, i:j]
+        projectors.append(block @ block.conj().T)
+        i = j
+    return projectors
+
+
+def is_trivial_overlap(
+    L_super: np.ndarray,
+    rho_initial: np.ndarray,
+    rho_steady_state: np.ndarray,
+    *,
+    atol: float = EPS_GAP,
+    block_tol: float = 1.0e-9,
+) -> bool:
+    """Is a vanishing slowest-mode overlap symmetry-protected (not Mpemba)?
+
+    Returns ``True`` when the zero overlap is *structural*: in the sector
+    decomposition set by the eigenprojectors ``{P_a}`` of ``rho_ss``, the slowest
+    mode ``l_1`` and ``rho_0`` occupy **disjoint** blocks ``P_a (.) P_b``. Then
+    ``rho_0`` can never populate the slowest mode -- its fast relaxation is
+    trivial, not a fine-tuned Mpemba skip (e.g. a diagonal ``rho_0`` vs a
+    coherence slowest mode of an amplitude-damping channel).
+
+    Using eigen*projectors* (not individual eigenvectors) keeps the test robust
+    when ``rho_ss`` is degenerate. A maximally mixed ``rho_ss`` collapses to a
+    single projector (one block), so nothing is structurally decoupled and the
+    guard does not fire -- such a state has no protecting symmetry sector and is
+    handled by the (demoted) confidence tier instead.
+    """
+    L_super = np.asarray(L_super)
+    rho_initial = np.asarray(rho_initial, dtype=complex)
+    rho_steady_state = np.asarray(rho_steady_state, dtype=complex)
+
+    mode = _slowest_mode(L_super, atol=atol)
+    if mode is None:
+        return False
+    l_slow, _ = mode
+    d = rho_steady_state.shape[0]
+    l_mat = l_slow.reshape(d, d, order="F")
+
+    projectors = _steady_state_eigenprojectors(rho_steady_state)
+    l_scale = max(float(np.linalg.norm(l_mat)), 1.0)
+    rho_scale = max(float(np.linalg.norm(rho_initial)), 1.0)
+    for p_a in projectors:
+        for p_b in projectors:
+            l_block = float(np.linalg.norm(p_a @ l_mat @ p_b))
+            rho_block = float(np.linalg.norm(p_a @ rho_initial @ p_b))
+            if l_block > block_tol * l_scale and rho_block > block_tol * rho_scale:
+                return False  # a shared block -> rho_0 can reach the slow mode
+    return True
 
 
 def expansion_alpha(
@@ -56,7 +162,7 @@ def expansion_alpha(
     Returns ``alpha`` (slope). A flat distribution gives ``alpha`` near 0.
     """
     L_super = np.asarray(L_super)
-    eigvals, vl, _ = sla.eig(L_super, left=True, right=True)
+    eigvals, vl, vr = sla.eig(L_super, left=True, right=True)
     mask = np.abs(eigvals) > atol
     eigvals_nz = eigvals[mask]
     if eigvals_nz.size < 2:
@@ -67,8 +173,8 @@ def expansion_alpha(
     cs: list[float] = []
     for idx in nz_indices[: min(n_modes, len(nz_indices))]:
         l_n = vl[:, idx]
-        norm = max(np.linalg.norm(l_n), 1.0e-12)
-        cs.append(float(abs(np.vdot(l_n, rho_vec0)) / norm))
+        denom = abs(np.vdot(l_n, vr[:, idx]))
+        cs.append(float(abs(np.vdot(l_n, rho_vec0)) / max(denom, EPS_DIV)))
     if len(cs) < 2:
         return 0.0
     log_cs = np.log(np.clip(np.asarray(cs), 1.0e-30, None))
@@ -81,13 +187,26 @@ def compute_mpemba_layer(
     L_super: np.ndarray,
     rho_initial: np.ndarray,
     *,
+    rho_steady_state: np.ndarray | None = None,
     overlap_threshold: float = 1.0e-4,
 ) -> MpembaResult:
-    """Run D19, D20 and flag Mpemba candidacy."""
+    """Run D19, D20 and flag Mpemba candidacy.
+
+    ``is_mpemba_candidate`` is ``True`` only when ``|c_1| < overlap_threshold``
+    **and** the vanishing overlap is not symmetry-protected. Triviality can only
+    be assessed when ``rho_steady_state`` is supplied (the standalone default is
+    ``None``, which preserves the raw overlap test for direct callers).
+    """
     c1 = overlap_c1(L_super, rho_initial)
     alpha = expansion_alpha(L_super, rho_initial)
+    trivial = (
+        is_trivial_overlap(L_super, rho_initial, rho_steady_state)
+        if rho_steady_state is not None
+        else False
+    )
     return MpembaResult(
         overlap_c1=c1,
-        is_mpemba_candidate=c1 < overlap_threshold,
+        is_mpemba_candidate=bool(c1 < overlap_threshold and not trivial),
         expansion_alpha=alpha,
+        trivial_overlap=trivial,
     )

tests/test_classification.py

@@ -203,10 +203,16 @@ def _lep(**kw) -> LepResult:
 def _mpemba(**kw) -> MpembaResult:
     base = {
         "overlap_c1": 0.5,
-        "is_mpemba_candidate": False,
         "expansion_alpha": 1.0,
     }
     base.update(kw)
+    # Synthetic evidence carries no steady state, so mirror the real
+    # (non-trivial) candidate rule: a below-threshold overlap is a candidate
+    # unless the caller explicitly marks it trivial or overrides the flag.
+    base.setdefault(
+        "is_mpemba_candidate",
+        bool(base["overlap_c1"] < 1.0e-4 and not kw.get("trivial_overlap", False)),
+    )
     return MpembaResult(**base)
 
 
@@ -241,11 +247,13 @@ def test_branch_default_fallback_a12():
 
 
 def test_branch_a11_mpemba():
+    # issue #68: a single-state Mpemba overlap is CANDIDATE-grade, never
+    # PUBLICATION_GRADE-confirmed (confirmation needs a reference family).
     cls = _classify(mpemba=_mpemba(overlap_c1=1.0e-6))
     assert (cls.a_class, cls.f_family) == ("A11", "F4")
-    assert cls.confidence == pytest.approx(0.90)
-    assert cls.verdict == VERDICT_CONFIRMED
-    assert cls.tier == TIER_PUBLICATION
+    assert cls.confidence == pytest.approx(0.70)
+    assert cls.verdict == VERDICT_CANDIDATE
+    assert cls.tier == TIER_CONFIRMATION
 
 
 def test_branch_a11_mpemba_weak_confidence():
@@ -257,6 +265,15 @@ def test_branch_a11_mpemba_weak_confidence():
     assert cls.verdict == VERDICT_NOT_EXCLUDED
 
 
+def test_branch_a11_trivial_overlap_not_flagged():
+    # issue #68: a symmetry-protected zero overlap (trivial_overlap=True) must
+    # NOT reach A11 even though overlap_c1 is below the raw trigger.
+    cls = _classify(mpemba=_mpemba(overlap_c1=1.0e-6, trivial_overlap=True))
+    assert cls.a_class != "A11"
+    assert cls.evidence["mpemba_trivial_overlap"] == 1.0
+    assert cls.evidence["mpemba_is_candidate"] == 0.0
+
+
 def test_branch_a10_phantom_relaxation():
     cls = _classify(
         resolvent=_resolvent(pseudospectral_radius=5.0),