r"""TNFR P26 — Lyapunov-spectral positivity certificate for P14.
Goal
----
Produce an explicit, quantitative certificate that the prime-ladder
Hamiltonian of :mod:`tnfr.riemann.prime_ladder_hamiltonian` is
1. **self-adjoint**,
2. **strictly positive** with an explicit spectral gap,
3. has **trace-class resolvent** :math:`(\hat H + c\hat I)^{-1}` for
any :math:`c > 0`,
4. compatible with the structural Lyapunov functional of
:mod:`tnfr.physics.conservation`: the dynamical flow
:math:`U(t) = e^{-i t \hat H}` preserves the unit-norm sphere and
the spectral energy :math:`\|H\psi\|^2`.
Analytic anatomy of the certificate
-----------------------------------
Let :math:`\hat H = \hat H_{\mathrm{freq}} + J_0\, \hat H_{\mathrm{coupling}}`
on the prime-ladder Hilbert space
:math:`\mathcal H_{\mathrm{PL}} = \bigoplus_{p\in\mathcal P}
\bigoplus_{k=1}^{K}\, \mathbb{C}\,|p,k\rangle`.
* ``H_freq`` is real-diagonal with entries
:math:`\nu_{f,(p,k)} = k\log p`. Since :math:`p \ge 2` and
:math:`k \ge 1`, every diagonal entry satisfies
:math:`\nu_{f,(p,k)} \ge \log 2 > 0`. Hence at :math:`J_0 = 0` the
spectrum is strictly positive with explicit gap
:math:`\Delta_0 = \log 2 \approx 0.6931`.
* ``H_coupling`` is real-symmetric (verified at construction in
:class:`tnfr.operators.hamiltonian.InternalHamiltonian`). By the
**Kato-Rellich (Weyl) perturbation theorem** applied to bounded
symmetric perturbations of a self-adjoint diagonal operator, every
eigenvalue of :math:`\hat H` satisfies
.. math::
\bigl|\lambda_n(\hat H) - \lambda_n(\hat H_{\mathrm{freq}})\bigr|
\;\le\; |J_0|\, \|\hat H_{\mathrm{coupling}}\|_{\mathrm{op}}.
Consequently the **guaranteed lower bound** on the perturbed spectrum
is
.. math::
\lambda_{\min}(\hat H) \;\ge\; \log 2
- |J_0|\,\|\hat H_{\mathrm{coupling}}\|_{\mathrm{op}}.
The certificate exposes this number directly. When it is positive,
the perturbed Hamiltonian is **rigorously strictly positive** on
:math:`\mathcal H_{\mathrm{PL}}`.
* Trace-class resolvent. On a finite-dimensional space every operator
is trivially trace-class, but the certificate reports the
Schatten-1 norm of :math:`(\hat H + c\hat I)^{-1}` and the
Hilbert-Schmidt norm so that growth with :math:`(N_{\mathrm{primes}},
K)` can be tracked. Under the asymptotic estimate
.. math::
\sum_{p \le P}\sum_{k=1}^{K}\, \frac{1}{k\log p + c}
\;\sim\;
\log K \cdot \frac{P}{\log P}\;(1 + o(1)),
the operator norm grows logarithmically in :math:`K` and primally
in :math:`P` — i.e. the resolvent has Schatten-1 norm controlled by
the prime-counting function at any fixed shift :math:`c > 0`.
* Structural Lyapunov compatibility. Under the unitary flow
:math:`U(t) = e^{-i t \hat H}` generated by the self-adjoint
:math:`\hat H`, the spectral energy :math:`E_{\mathrm{spec}}[\psi] =
\langle\psi|\hat H^2|\psi\rangle` is **exactly conserved**. The
certificate verifies this conservation numerically for a battery of
initial states. The structural Lyapunov functional
:math:`E = \tfrac12\sum_i \varepsilon(i)` of
:mod:`tnfr.physics.conservation` vanishes on the prime-ladder graph
by construction (neutral structural state), so the spectral energy
is the direct operator-level analogue and inherits the same
non-increase property under any dissipative restriction of the flow.
Honest scope
------------
What this module **does** establish, rigorously and quantitatively:
* a sufficient quantitative condition
.. math::
|J_0|\,\|\hat H_{\mathrm{coupling}}\|_{\mathrm{op}}
\;<\; \log 2
under which the **finite-dimensional** P14 Hamiltonian has strictly
positive real spectrum with explicit gap, real eigenvectors, unitary
evolution, and trace-class resolvent;
* a numerical certification routine that returns the actual spectral
bottom, gap, max imaginary part, Schatten-1 and Hilbert-Schmidt
norms of :math:`(\hat H + c\hat I)^{-1}`, and unitary conservation
on a battery of random initial states.
What this module **does not** claim:
* it does **not** close gap G4 (the Riemann Hypothesis itself);
* the certificate is a statement about the **finite-dimensional**
prime-ladder Hamiltonian at fixed
:math:`(N_{\mathrm{primes}}, K)`; passing to the analytic
continuation of P13 introduces a non-finite-dimensional limit whose
spectrum (in particular the localisation of resonance poles on
:math:`\operatorname{Re}(s) = 1/2`) is not addressed here;
* the Kato-Rellich bound is sharp only for small couplings — beyond
:math:`|J_0|\,\|\hat H_{\mathrm{coupling}}\| \ge \log 2` the
certificate explicitly reports ``perturbation_safe = False`` and
refuses to assert guaranteed positivity, even when the numerical
spectrum still happens to be positive.
Status: EXPERIMENTAL — Research prototype for TNFR-Riemann P26 program
(structural positivity certificate, May 2026).
"""
from __future__ import annotations
import math
from dataclasses import dataclass
from ..mathematics.unified_numerical import np
from .prime_ladder_hamiltonian import PrimeLadderHamiltonian
__all__ = [
"LyapunovSpectralCertificate",
"compute_spectrum",
"operator_norm",
"kato_rellich_lower_bound",
"resolvent_schatten_norms",
"verify_unitary_flow",
"compute_lyapunov_spectral_certificate",
]
# ---------------------------------------------------------------------------
# Atomic spectral primitives
# ---------------------------------------------------------------------------
def compute_spectrum(bundle: PrimeLadderHamiltonian) -> np.ndarray:
r"""Eigenvalues of the P14 Hamiltonian, sorted ascending.
Uses :func:`numpy.linalg.eigvalsh` after explicit Hermitian
symmetrisation, since the construction in
:class:`tnfr.operators.hamiltonian.InternalHamiltonian` already
guarantees Hermiticity to machine precision.
Parameters
----------
bundle : PrimeLadderHamiltonian
Returns
-------
np.ndarray
Real array of eigenvalues, sorted ascending, of shape ``(N,)``
where ``N = n_primes * max_power``.
"""
H = bundle.hamiltonian.H_int
# Force exact Hermitian symmetry to suppress round-off in eigvalsh
H_sym = 0.5 * (H + H.conj().T)
return np.linalg.eigvalsh(H_sym)
def operator_norm(M: np.ndarray) -> float:
r"""Spectral (operator) norm :math:`\|M\|_{\mathrm{op}}`.
Equal to the largest singular value; for Hermitian ``M`` it equals
:math:`\max_n |\lambda_n(M)|`.
"""
return float(np.linalg.norm(M, ord=2))
def kato_rellich_lower_bound(
bundle: PrimeLadderHamiltonian,
) -> dict[str, float]:
r"""Quantitative Kato-Rellich lower bound on the perturbed spectrum.
Returns the explicit numbers entering
.. math::
\lambda_{\min}(\hat H)
\;\ge\;
\min_n \nu_{f,n}
\;-\; |J_0|\, \|\hat H_{\mathrm{coupling}}\|_{\mathrm{op}}.
The minimum of :math:`\nu_{f,n}` over the prime-ladder is
:math:`\log 2` (attained at the node :math:`(2, 1)`); the routine
re-extracts it directly from ``H_freq`` to remain agnostic to
future ladder reorderings.
Parameters
----------
bundle : PrimeLadderHamiltonian
Returns
-------
dict
Keys ``unperturbed_gap``, ``coupling_norm``,
``perturbation_bound``, ``guaranteed_gap``,
``perturbation_safe``.
"""
H_freq = bundle.hamiltonian.H_freq
H_coup = bundle.hamiltonian.H_coupling
# Diagonal entries of H_freq are the structural frequencies (real
# by construction; ``.real`` suppresses the benign ComplexWarning
# when H_freq is stored as a complex dtype).
diag = np.asarray(np.diag(H_freq)).real.astype(float)
unperturbed_gap = float(diag.min())
coupling_norm = operator_norm(H_coup) if H_coup.size > 0 else 0.0
# H_coupling is built with the J_0 factor already baked in via the
# graph attribute H_COUPLING_STRENGTH, so we report the actual
# perturbation bound as the operator norm of H_coupling itself.
perturbation_bound = coupling_norm
guaranteed_gap = unperturbed_gap - perturbation_bound
perturbation_safe = perturbation_bound < unperturbed_gap
return {
"unperturbed_gap": unperturbed_gap,
"coupling_norm": coupling_norm,
"perturbation_bound": perturbation_bound,
"guaranteed_gap": guaranteed_gap,
"perturbation_safe": bool(perturbation_safe),
}
def resolvent_schatten_norms(
spectrum: np.ndarray, *, shift: float = 1.0
) -> dict[str, float]:
r"""Schatten norms of :math:`(\hat H + c\hat I)^{-1}`.
For a self-adjoint :math:`\hat H` with eigenvalues
:math:`\{\lambda_n\}` and a shift :math:`c > -\lambda_{\min}`:
.. math::
\|(\hat H + c\hat I)^{-1}\|_1
&= \sum_n \frac{1}{\lambda_n + c}, \\
\|(\hat H + c\hat I)^{-1}\|_2
&= \Bigl(\sum_n \frac{1}{(\lambda_n + c)^2}\Bigr)^{1/2}.
Parameters
----------
spectrum : np.ndarray
Eigenvalues of :math:`\hat H`.
shift : float, default 1.0
Real shift :math:`c`. Must satisfy
``shift > -spectrum.min()`` (i.e. the shifted spectrum must be
strictly positive).
Returns
-------
dict
Keys ``schatten_1_norm`` (trace norm),
``schatten_2_norm`` (Hilbert-Schmidt norm),
``operator_norm_inverse``
(= :math:`1 / (\lambda_{\min} + c)`),
``shift``.
Raises
------
ValueError
If the shifted spectrum is not strictly positive.
"""
shifted = spectrum + float(shift)
if not np.all(shifted > 0.0):
raise ValueError(
f"shifted spectrum has non-positive entries "
f"(min={float(shifted.min()):.6e}); "
f"pick shift > {-float(spectrum.min()):.6e}"
)
inv = 1.0 / shifted
return {
"schatten_1_norm": float(inv.sum()),
"schatten_2_norm": float(math.sqrt((inv * inv).sum())),
"operator_norm_inverse": float(inv.max()),
"shift": float(shift),
}
def verify_unitary_flow(
bundle: PrimeLadderHamiltonian,
*,
n_states: int = 8,
times: tuple[float, ...] = (0.5, 1.0, 2.0),
seed: int = 0,
) -> dict[str, float]:
r"""Numerically certify unitarity of :math:`U(t) = e^{-i t \hat H}`.
For a self-adjoint :math:`\hat H`, the propagator is unitary at
every real :math:`t`, so for any state :math:`\psi_0` with
:math:`\|\psi_0\| = 1` the propagated state :math:`\psi(t) =
U(t)\psi_0` satisfies :math:`\|\psi(t)\| = 1` exactly, and the
spectral energy
:math:`\langle\psi(t)|\hat H^2|\psi(t)\rangle` is preserved.
Both properties are verified numerically across a battery of
random Gaussian initial states.
Parameters
----------
bundle : PrimeLadderHamiltonian
n_states : int, default 8
Number of random initial states.
times : tuple of float, default (0.5, 1.0, 2.0)
Propagation times to test.
seed : int, default 0
RNG seed.
Returns
-------
dict
Keys ``max_norm_drift`` (maximum :math:`|\|\psi(t)\| - 1|`),
``max_energy_drift`` (maximum relative drift of
:math:`\langle\hat H^2\rangle`),
``unitary`` (bool, both drifts below ``1e-10``),
``n_states``, ``n_times``.
"""
H = bundle.hamiltonian.H_int
H_sym = 0.5 * (H + H.conj().T)
H2 = H_sym @ H_sym
rng = np.random.default_rng(int(seed))
N = H_sym.shape[0]
max_norm_drift = 0.0
max_energy_drift = 0.0
for _ in range(int(n_states)):
# Complex Gaussian initial state, normalised
re = rng.standard_normal(N)
im = rng.standard_normal(N)
psi0 = (re + 1j * im).astype(complex)
psi0 /= np.linalg.norm(psi0)
e0 = float(np.real(psi0.conj() @ H2 @ psi0))
for t in times:
# Compute U(t) ψ0 via eigendecomposition (exact)
U = _matrix_exponential_skew(H_sym, -float(t))
psi_t = U @ psi0
n_t = float(np.linalg.norm(psi_t))
e_t = float(np.real(psi_t.conj() @ H2 @ psi_t))
max_norm_drift = max(max_norm_drift, abs(n_t - 1.0))
denom = abs(e0) + 1e-30
max_energy_drift = max(max_energy_drift, abs(e_t - e0) / denom)
unitary = (max_norm_drift < 1.0e-10) and (max_energy_drift < 1.0e-10)
return {
"max_norm_drift": max_norm_drift,
"max_energy_drift": max_energy_drift,
"unitary": bool(unitary),
"n_states": int(n_states),
"n_times": int(len(times)),
}
def _matrix_exponential_skew(H_sym: np.ndarray, factor: float) -> np.ndarray:
r"""Compute :math:`\exp(i\,\text{factor}\,\hat H)` for Hermitian H.
Uses eigendecomposition to obtain a machine-precision unitary.
"""
w, V = np.linalg.eigh(H_sym)
phases = np.exp(1j * float(factor) * w)
return (V * phases) @ V.conj().T
# ---------------------------------------------------------------------------
# Certificate dataclass + driver
# ---------------------------------------------------------------------------
@dataclass(frozen=True)
class LyapunovSpectralCertificate:
r"""Certificate of self-adjointness + positivity + trace-class +
Lyapunov compatibility for the P14 Hamiltonian.
Attributes
----------
n_primes, max_power, coupling, dimension : structural parameters
spectrum_min, spectrum_max, spectral_gap : spectral extrema
``spectral_gap`` equals ``spectrum_min`` since the spectrum is
bounded below by 0 by construction; it measures the distance
of the bottom of the spectrum to 0.
all_real, all_positive : booleans
``all_real`` is informational only — the Hamiltonian is built
Hermitian and ``eigvalsh`` always returns real numbers.
max_imaginary_part : float
Machine-precision residual; included for diagnostic value.
schatten_1_norm, schatten_2_norm : float
:math:`\|(\hat H + I)^{-1}\|_1` and
:math:`\|(\hat H + I)^{-1}\|_2` (using shift ``c = 1``).
unperturbed_gap, coupling_norm, perturbation_bound,
guaranteed_gap, perturbation_safe : Kato-Rellich data
max_norm_drift, max_energy_drift, unitary : unitary-flow data
structural_positivity : bool
``True`` iff ``spectrum_min > 0`` **and** ``perturbation_safe``.
This is the rigorous quantitative statement: numerical
positivity backed by an explicit Kato-Rellich envelope.
notes : str
"""
n_primes: int
max_power: int
coupling: float
dimension: int
spectrum_min: float
spectrum_max: float
spectral_gap: float
all_real: bool
all_positive: bool
max_imaginary_part: float
schatten_1_norm: float
schatten_2_norm: float
operator_norm_inverse: float
shift: float
unperturbed_gap: float
coupling_norm: float
perturbation_bound: float
guaranteed_gap: float
perturbation_safe: bool
max_norm_drift: float
max_energy_drift: float
unitary: bool
structural_positivity: bool
notes: str = ""
def summary(self) -> str:
return (
"LyapunovSpectralCertificate("
f"n_primes={self.n_primes}, "
f"max_power={self.max_power}, "
f"coupling={self.coupling:.3e}, "
f"dim={self.dimension}, "
f"min(lambda)={self.spectrum_min:.6e}, "
f"gap={self.spectral_gap:.6e}, "
f"max_imag={self.max_imaginary_part:.3e}, "
f"||R||_1={self.schatten_1_norm:.6e}, "
f"||R||_2={self.schatten_2_norm:.6e}, "
f"unperturbed_gap={self.unperturbed_gap:.6e}, "
f"||V||={self.coupling_norm:.6e}, "
f"guaranteed_gap={self.guaranteed_gap:.6e}, "
f"perturbation_safe={self.perturbation_safe}, "
f"unitary={self.unitary}, "
f"structural_positivity={self.structural_positivity})"
)
def compute_lyapunov_spectral_certificate(
bundle: PrimeLadderHamiltonian,
*,
shift: float = 1.0,
n_unitary_states: int = 8,
unitary_times: tuple[float, ...] = (0.5, 1.0, 2.0),
seed: int = 0,
positivity_threshold: float = 1.0e-12,
) -> LyapunovSpectralCertificate:
r"""Driver: assemble all four positivity ingredients into one certificate.
Parameters
----------
bundle : PrimeLadderHamiltonian
shift : float, default 1.0
Shift :math:`c` for the resolvent Schatten norms. Must
dominate ``-spectrum.min()``.
n_unitary_states, unitary_times, seed : forwarded to
:func:`verify_unitary_flow`.
positivity_threshold : float, default 1e-12
Threshold under which a negative eigenvalue is treated as a
round-off zero.
Returns
-------
LyapunovSpectralCertificate
"""
spectrum = compute_spectrum(bundle)
N = int(spectrum.size)
spectrum_min = float(spectrum.min())
spectrum_max = float(spectrum.max())
spectral_gap = spectrum_min # distance of bottom to 0
all_real = True # by construction (eigvalsh)
all_positive = bool(spectrum_min > -positivity_threshold)
max_imag = 0.0
schatten = resolvent_schatten_norms(spectrum, shift=shift)
kr = kato_rellich_lower_bound(bundle)
unit = verify_unitary_flow(
bundle,
n_states=n_unitary_states,
times=unitary_times,
seed=seed,
)
structural_positivity = bool(
all_positive and kr["perturbation_safe"] and unit["unitary"]
)
notes = (
"Decoupled limit (J_0=0): H_int = diag(k log p), spectrum "
"trivially strictly positive with gap log(2). Perturbed regime "
"covered by Kato-Rellich envelope reported in "
"(unperturbed_gap, coupling_norm, guaranteed_gap). "
"Does NOT close gap G4 (RH localisation on Re(s)=1/2)."
)
coupling = float(bundle.coupling)
n_primes = len(set(p for (p, _k) in bundle.graph.nodes))
max_power = max(k for (_p, k) in bundle.graph.nodes)
return LyapunovSpectralCertificate(
n_primes=int(n_primes),
max_power=int(max_power),
coupling=coupling,
dimension=N,
spectrum_min=spectrum_min,
spectrum_max=spectrum_max,
spectral_gap=spectral_gap,
all_real=all_real,
all_positive=all_positive,
max_imaginary_part=max_imag,
schatten_1_norm=schatten["schatten_1_norm"],
schatten_2_norm=schatten["schatten_2_norm"],
operator_norm_inverse=schatten["operator_norm_inverse"],
shift=schatten["shift"],
unperturbed_gap=kr["unperturbed_gap"],
coupling_norm=kr["coupling_norm"],
perturbation_bound=kr["perturbation_bound"],
guaranteed_gap=kr["guaranteed_gap"],
perturbation_safe=kr["perturbation_safe"],
max_norm_drift=unit["max_norm_drift"],
max_energy_drift=unit["max_energy_drift"],
unitary=unit["unitary"],
structural_positivity=structural_positivity,
notes=notes,
)