diff --git a/algebra/backend.py b/algebra/backend.py
index 92941e61..1d590c2d 100644
--- a/algebra/backend.py
+++ b/algebra/backend.py
@@ -55,20 +55,39 @@ def _f32_1d32(x: np.ndarray) -> np.ndarray:
     )
 
 
-def _is_f32_workload(*arrays: np.ndarray) -> bool:
-    """True when all arrays are float32 (Rust f32 kernel is parity-safe).
+def _f64_1d32(x: np.ndarray) -> np.ndarray:
+    """Contiguous f64 (32,) for core_rs PyReadonlyArray1<f64> bindings."""
+    return np.ascontiguousarray(
+        np.asarray(x, dtype=np.float64).reshape(-1)[:32], dtype=np.float64
+    )
 
-    float64 wave residual pins require Python SOT (or future f64 Rust GP).
-    Forcing f64→f32 would break 1e-9 chiral / leakage pins (ADR-0241).
-    """
+
+def _is_f32_workload(*arrays: np.ndarray) -> bool:
+    """True when all arrays are float32 (Rust f32 kernel is parity-safe)."""
     return all(np.asarray(a).dtype == np.float32 for a in arrays)
 
 
-def geometric_product(A: np.ndarray, B: np.ndarray) -> np.ndarray:
-    """Cl(4,1) geometric product via Rust f32 when enabled, else Python.
+def _is_f64_workload(*arrays: np.ndarray) -> bool:
+    """True when all arrays are float64.
 
-    float64 inputs always use the pure-Python product (semantic SOT for
-    wave-field residual math). float32 field-graph workloads get Rust.
+    The Rust f64 kernel (``core_rs.geometric_product_f64``) is a term-for-term
+    mirror of the pure-Python f64 kernel — same scatter order, no FMA — so it is
+    **bit-identical**, not merely close. That is what lets f64 workloads take the
+    Rust path without moving the 1e-9 chiral / leakage residual pins (ADR-0241);
+    the D9 parity suite gates that bit-identity. An older ``core_rs`` build
+    without the function raises ``AttributeError`` and falls through to Python.
+    """
+    return all(np.asarray(a).dtype == np.float64 for a in arrays)
+
+
+def geometric_product(A: np.ndarray, B: np.ndarray) -> np.ndarray:
+    """Cl(4,1) geometric product via Rust when enabled, else pure Python.
+
+    float32 field-graph workloads and float64 wave-field workloads both get the
+    Rust kernel when ``CORE_BACKEND=rust`` and ``core_rs`` is present — the f64
+    path is bit-identical to Python (see :func:`_is_f64_workload`), so it never
+    perturbs residual math; it is a speed swap, not a numeric one. Any other
+    dtype, or an older/absent ``core_rs``, uses the pure-Python product.
     """
     if _RUST and _is_f32_workload(A, B):
         try:
@@ -78,6 +97,14 @@ def geometric_product(A: np.ndarray, B: np.ndarray) -> np.ndarray:
             )
         except (AttributeError, TypeError, ValueError, Exception):
             pass
+    if _RUST and _is_f64_workload(A, B):
+        try:
+            return np.asarray(
+                _rs.geometric_product_f64(_f64_1d32(A), _f64_1d32(B)),
+                dtype=np.float64,
+            )
+        except (AttributeError, TypeError, ValueError, Exception):
+            pass
     from algebra.cl41 import geometric_product as _gp
     return _gp(A, B)
 
diff --git a/core-rs/src/lib.rs b/core-rs/src/lib.rs
index d07e5270..975a8163 100644
--- a/core-rs/src/lib.rs
+++ b/core-rs/src/lib.rs
@@ -20,6 +20,7 @@ pub mod vault;
 pub mod versor;
 
 use cga::cga_inner_raw;
+use cl41::geometric_product_f64 as cl41_geometric_product_f64;
 use cl41::geometric_product_raw;
 use diffusion::{graph_diffusion_step, unitize_f32};
 use versor::{
@@ -45,6 +46,23 @@ fn geometric_product(
     f32_array_to_numpy(py, &result)
 }
 
+/// Full geometric product in Cl(4,1), f64.  Bit-identical to the pure-Python
+/// f64 kernel: the Rust scalar path (`cl41::geometric_product_f64`) mirrors it
+/// term-for-term — same i-major/j-minor scatter order, same left-associative
+/// `sign * ai * bj`, no FMA contraction — so enabling it never perturbs the
+/// f64 wave-field residual math (ADR-0244 §2.6; parity gated by the D9 suite).
+#[pyfunction]
+fn geometric_product_f64(
+    py: Python<'_>,
+    a: numpy::PyReadonlyArray1<'_, f64>,
+    b: numpy::PyReadonlyArray1<'_, f64>,
+) -> PyResult<PyObject> {
+    let a_slice = read_f64_cl41_mv(&a)?;
+    let b_slice = read_f64_cl41_mv(&b)?;
+    let result = cl41_geometric_product_f64(a_slice, b_slice);
+    f64_array_to_numpy(py, &result)
+}
+
 /// Sandwich product V*F*reverse(V).
 #[pyfunction]
 fn versor_apply(
@@ -360,6 +378,7 @@ fn f64_array_to_numpy(py: Python<'_>, data: &[f64; 32]) -> PyResult<PyObject> {
 #[pymodule]
 fn core_rs(m: &Bound<'_, PyModule>) -> PyResult<()> {
     m.add_function(wrap_pyfunction!(geometric_product, m)?)?;
+    m.add_function(wrap_pyfunction!(geometric_product_f64, m)?)?;
     m.add_function(wrap_pyfunction!(versor_apply, m)?)?;
     m.add_function(wrap_pyfunction!(versor_apply_with_closure, m)?)?;
     m.add_function(wrap_pyfunction!(versor_apply_with_closure_f64, m)?)?;
diff --git a/core/physics/cognitive_lifecycle.py b/core/physics/cognitive_lifecycle.py
index c16daee5..0e9caeb2 100644
--- a/core/physics/cognitive_lifecycle.py
+++ b/core/physics/cognitive_lifecycle.py
@@ -57,6 +57,7 @@ lazily via the ``core.physics`` barrel; enforced by
 
 from __future__ import annotations
 
+import functools
 import hashlib
 import json
 from dataclasses import dataclass, field
@@ -505,6 +506,27 @@ def _spectral_gap(evals: np.ndarray, tol: float) -> tuple[float, float, float]:
     return lam0, gap, energy_tol
 
 
+@functools.lru_cache(maxsize=128)
+def _cached_eigh(hamiltonian_id: str, matrix_bytes: bytes) -> tuple[np.ndarray, np.ndarray]:
+    """Memoized symmetric eigendecomposition (ADR-0244 §2.8 / directive M2).
+
+    ``ProblemHamiltonian`` is frozen and content-addressed, so a fresh LAPACK
+    ``eigh`` on an identical matrix (repeated active-turn / biography checks) is
+    wasted AMX compute. Keyed on the immutable ``hamiltonian_id`` *and* the raw
+    matrix bytes (collision-resistant: the id already content-addresses the
+    matrix; the bytes make a same-id/different-bytes hit impossible). The
+    returned arrays are frozen read-only so a cache hit cannot be mutated by a
+    caller — every hit yields bit-identical ``(evals, evecs)``.
+    """
+    matrix = np.frombuffer(matrix_bytes, dtype=np.float64).reshape(N_COMPONENTS, N_COMPONENTS)
+    evals, evecs = np.linalg.eigh(matrix)
+    evals = np.ascontiguousarray(evals)
+    evecs = np.ascontiguousarray(evecs)
+    evals.setflags(write=False)
+    evecs.setflags(write=False)
+    return evals, evecs
+
+
 def relax_to_ground(
     psi0: np.ndarray,
     hamiltonian: ProblemHamiltonian,
@@ -563,7 +585,7 @@ def relax_to_ground(
             return diag * v
 
     else:
-        evals_full, evecs_full = np.linalg.eigh(H_mat)
+        evals_full, evecs_full = _cached_eigh(hamiltonian.hamiltonian_id, H_mat.tobytes())
         lam0, gap, energy_tol = _spectral_gap(evals_full, tol_f)
         propagator = evecs_full @ np.diag(np.exp(-(evals_full - lam0) * dt_f)) @ evecs_full.T
 
diff --git a/tests/test_adr_0244_mechanical_sympathy.py b/tests/test_adr_0244_mechanical_sympathy.py
new file mode 100644
index 00000000..bc6adb5b
--- /dev/null
+++ b/tests/test_adr_0244_mechanical_sympathy.py
@@ -0,0 +1,69 @@
+"""ADR-0244 §2.8 mechanical-sympathy pins (cohesion directive Mandate 2).
+
+The eigendecomposition inside ``relax_to_ground`` for a non-diagonal, frozen,
+content-addressed ``ProblemHamiltonian`` must be memoized — a fresh LAPACK
+``eigh`` on an identical matrix is wasted AMX compute. The cache must return
+bit-identical, read-only ``(evals, evecs)`` so a hit can never be mutated into a
+different result, and a cached decomposition must equal a fresh one.
+
+(Mandate 1 — the Rust f64 ``geometric_product`` fast-path — is pinned by the
+bit-identical parity suite in ``tests/test_geometric_product_f64_parity.py``.)
+"""
+
+from __future__ import annotations
+
+import numpy as np
+
+from algebra.rotor import make_rotor_from_angle
+from core.physics.cognitive_lifecycle import (
+    _cached_eigh,
+    compile_quadratic_well,
+    relax_to_ground,
+)
+
+
+def _dense_well():
+    # c·(I − TTᵀ) is non-diagonal → exercises the eigh (not the diagonal) branch.
+    target = np.ascontiguousarray(make_rotor_from_angle(0.9, 7), dtype=np.float64)
+    well = compile_quadratic_well(target, curvature=1.0)
+    assert well.is_diagonal is False
+    return well
+
+
+def test_cached_eigh_matches_fresh_eigh() -> None:
+    well = _dense_well()
+    evals, evecs = _cached_eigh(well.hamiltonian_id, well.matrix.tobytes())
+    fresh_evals, fresh_evecs = np.linalg.eigh(well.matrix)
+    # Same matrix bytes → same LAPACK call → bit-identical.
+    assert np.array_equal(evals, fresh_evals)
+    assert np.array_equal(evecs, fresh_evecs)
+
+
+def test_cached_eigh_returns_readonly_arrays() -> None:
+    well = _dense_well()
+    evals, evecs = _cached_eigh(well.hamiltonian_id, well.matrix.tobytes())
+    assert evals.flags.writeable is False
+    assert evecs.flags.writeable is False
+
+
+def test_cached_eigh_hit_returns_identical_objects() -> None:
+    # Distinct target → distinct hamiltonian_id/bytes → clean miss then hit.
+    target = np.ascontiguousarray(make_rotor_from_angle(1.27, 8), dtype=np.float64)
+    well = compile_quadratic_well(target, curvature=1.3)
+    key = (well.hamiltonian_id, well.matrix.tobytes())
+    a_evals, a_evecs = _cached_eigh(*key)
+    b_evals, b_evecs = _cached_eigh(*key)
+    # A cache hit hands back the very same frozen objects (no recompute).
+    assert a_evals is b_evals
+    assert a_evecs is b_evecs
+
+
+def test_relaxation_is_deterministic_through_the_cache() -> None:
+    well = _dense_well()
+    start = np.ascontiguousarray(make_rotor_from_angle(0.2, 6), dtype=np.float64)
+    start = start / float(np.linalg.norm(start))
+    r1 = relax_to_ground(start, well)
+    r2 = relax_to_ground(start, well)
+    assert r1.certificate.certificate_id == r2.certificate.certificate_id
+    assert r1.certificate.converged is True
+    assert np.array_equal(r1.psi_steady, r2.psi_steady)
diff --git a/tests/test_geometric_product_f64_parity.py b/tests/test_geometric_product_f64_parity.py
index 9b3f5f5f..91d47cfd 100644
--- a/tests/test_geometric_product_f64_parity.py
+++ b/tests/test_geometric_product_f64_parity.py
@@ -187,11 +187,13 @@ def test_f64_backend_matches_python_sot_in_subprocess(seed: int) -> None:
 
 @pytest.mark.skipif(not _RUST_AVAILABLE, reason="core_rs extension not built")
 @pytest.mark.parametrize("seed", [0xD911, 0xD912])
-def test_f64_with_rust_opt_in_still_python_sot(seed: int) -> None:
-    """CORE_BACKEND=rust must not f32-truncate f64 GP (D9 honesty pin).
+def test_f64_with_rust_opt_in_matches_python_sot_bit_for_bit(seed: int) -> None:
+    """CORE_BACKEND=rust routes f64 GP to the Rust f64 kernel (ADR-0244 §2.6),
+    which is bit-identical to Python SOT — so the hex still matches exactly.
 
-    Rust exposes f32 geometric_product only; f64 remains Python SOT until a
-    future geometric_product_f64 PyO3 export is wired and parity-gated.
+    This is the D9 honesty pin, now stronger: not only no f32-truncation, but
+    no 1-ULP f64 divergence either. On an older ``core_rs`` build without the
+    export, backend falls through to Python and the hex still matches.
     """
     rs = _run_f64_backend("rust", seed)
     assert rs["using_rust"] is True
@@ -199,10 +201,31 @@ def test_f64_with_rust_opt_in_still_python_sot(seed: int) -> None:
     assert rs["dtype"] == "float64"
 
 
-def test_backend_source_documents_f64_python_sot() -> None:
-    """Structural pin: dispatch comments + _is_f32_workload gate remain."""
+@pytest.mark.skipif(not _RUST_AVAILABLE, reason="core_rs extension not built")
+def test_rust_f64_gp_is_bit_identical_to_python_n10000() -> None:
+    """Acceptance criterion 1 (ADR-0244 §2.6 / directive M1): the Rust f64
+    ``geometric_product`` equals the pure-Python f64 kernel **bit-for-bit** over
+    a large random panel — not tol-matched. A single ULP divergence would move
+    the f64 wave-field residual bytes and break I-02 replay under
+    ``CORE_BACKEND=rust``; this fails closed on the first mismatch.
+    """
+    if not hasattr(core_rs, "geometric_product_f64"):
+        pytest.skip("core_rs build predates geometric_product_f64 export")
+    rng = _rng(0xB17DE)
+    for _ in range(10_000):
+        a = np.ascontiguousarray(rng.standard_normal(N_COMPONENTS), dtype=np.float64)
+        b = np.ascontiguousarray(rng.standard_normal(N_COMPONENTS), dtype=np.float64)
+        rust = np.asarray(core_rs.geometric_product_f64(a, b), dtype=np.float64)
+        py = gp_py(a, b)
+        assert rust.tobytes() == py.tobytes()
+
+
+def test_backend_source_documents_f64_rust_bit_identical() -> None:
+    """Structural pin: both dtype gates exist and the f64 Rust path is
+    documented bit-identical (a speed swap, not a numeric one)."""
     src = (REPO / "algebra" / "backend.py").read_text(encoding="utf-8")
     assert "_is_f32_workload" in src
-    assert "float64 wave residual" in src or "float64" in src
-    # geometric_product only calls Rust under f32 workload gate
+    assert "_is_f64_workload" in src
     assert "if _RUST and _is_f32_workload(A, B):" in src
+    assert "if _RUST and _is_f64_workload(A, B):" in src
+    assert "bit-identical" in src