Locality, all the way up to the spectrum

Paper: Subsystem-Resolved Spectral Theory for Quantum Many-Body Hamiltonians
Author: MD Nahidul Hasan Sabit
arXiv: 2604.21929 [quant-ph], 2026

A single dial turns a Hamiltonian into a family of Hamiltonians indexed by lattice subsets, and the recent paper of Sabit shows the spectra inherit two flavours of locality the operators always had: they are stable under exponential truncation, and they're approximately additive across disjoint blocks. We re-derive the bounds in a few lines, then verify them numerically on a 9-qubit XXZ chain with long-range couplings — including the satisfying finite-range corner where additivity becomes exact.

Why Python

I picked Python — specifically NumPy plus SciPy — for one concrete reason: the work is dense Hermitian linear algebra on Hilbert spaces of dimension up to $2^9 = 512$, and the bookkeeping is combinatorial (subsets, diameters, supports) rather than performance-critical. NumPy's eigvalsh lands on LAPACK's symmetric tridiagonal driver, which is roughly two orders of magnitude faster than anything I'd write by hand, and the kron / eye plumbing makes "build a Pauli string on the lattice" a one-liner. The interesting code is the index gymnastics, not the matrix multiply, and Python is exactly the right level for that. A compiled language would have bought me nothing here — the eigensolver is the bottleneck, and it's already in Fortran.

(One environment note: on the host I ran this on, OpenBLAS defaulted to a no-thread spin loop. Setting OPENBLAS_NUM_THREADS=2 cut a 256×256 eigvalsh from ~3 s to ~8 ms. If your eigvalsh mysteriously hangs on dim-256 problems, that's the first knob.)

The idea, in plain terms

Take a quantum spin system on a lattice $\Lambda$, with Hamiltonian decomposed as a sum over subsets,

where $\Phi(X)$ is supported on the sites in $X$ and the family $\{\Phi(X)\}$ has some exponential-decay control on long-range pieces. Standard physics. Now stop thinking about the operator $H$ and start thinking about a family of operators indexed by sub-regions. For each $S \subseteq \Lambda$ define

the part of $H$ that lives entirely inside $S$, together with its spectrum. This is the subsystem Hamiltonian and the subsystem spectrum. Everything in the paper is about how this map $S \mapsto \mathcal{S}(S)$ behaves.

Two questions immediately matter:

1. Locality. If I delete every term in $H_S$ whose support has diameter larger than $r$ — call the result $H_{S,r}$ — how badly does the spectrum change? You'd hope the change is exponentially small in $r$, because the deleted terms are themselves exponentially small in their diameter.
2. Disjointness. If $S = S_1 \sqcup S_2$ with the two pieces separated by lattice distance $D$, can I read off $\mathcal{S}(S_1 \cup S_2)$ from $\mathcal{S}(S_1)$ and $\mathcal{S}(S_2)$ alone? You'd hope so up to an error that's exponentially small in $D$, because the cross-terms straddling the gap are exponentially small.

The paper answers both, with explicit constants, and gives the punchline I find aesthetically pleasing: in the finite-range case (every $\Phi(X)$ has diameter $\le R$) the additivity result is not approximate but exact as soon as $D > R$. Locality of the operator algebra lifts all the way to the level of spectra, viewed as compact subsets of the real line under the Hausdorff metric.

The math

Write $\mathrm{diam}(X)$ for the diameter of $X$ in the lattice metric, and define the Sabit $\mu$-norm

This is finite as long as the $\Phi(X)$ decay faster than $e^{-\mu\,\mathrm{diam}(X)}$. It's the right object: it bundles the operator norms of all interactions through a vertex $x$, weighted by an exponential that's hostile toward long-range terms.

Truncation in operator norm. For $r \ge 0$, set

Each term $\Phi(X)$ in the residue contributes some operator norm $\|\Phi(X)\|$. Bound the sum naively, but pull a fictitious $e^{\mu\,\mathrm{diam}(X)} \cdot e^{-\mu\,\mathrm{diam}(X)}$ inside, and sum over a vertex $x$ that $X$ must contain:

The first inequality is the triangle inequality. The second is "every $X$ touches at least one $x \in S$, so we may overcount by summing over $x \in S$." The factor $e^{-\mu r}$ is the diameter cutoff. That's it.

Truncation in the spectrum. $H_S$ and $H_{S,r}$ are Hermitian. Weyl's inequality (the Hermitian instance of Bauer–Fike) says the eigenvalues of $A + E$ deviate from those of $A$ by at most $\|E\|$ each. The Hausdorff distance between $\sigma(H_S)$ and $\sigma(H_{S,r})$, viewed as multisets in $\mathbb{R}$, is at most $\|H_S - H_{S,r}\|$, hence

This is the result that makes me sit up. The bound on a set is just as clean as the bound on the operator. We pay no extra in dimension.

Disjointness. For $S_1, S_2$ disjoint with $D = d(S_1, S_2)$, decompose

Every cross-term has $\mathrm{diam}(X) \ge D$, because $X$ contains a point in each of $S_1, S_2$. The same accounting trick gives

But the spectrum of a sum of operators acting on different tensor factors is the Minkowski sum of the individual spectra: $\sigma(A\otimes I + I\otimes B) = \sigma(A) + \sigma(B) := \{a+b : a\in\sigma(A), b\in\sigma(B)\}$. So combining Weyl's inequality with the cross-term norm,

In the finite-range case where every $\Phi(X)$ is zero unless $\mathrm{diam}(X) \le R$, all cross-terms vanish identically once $D > R$, and the inequality becomes equality (with both sides equal to zero). That's the corner I most wanted to see numerically, because the bound's predictive power is checkable to machine precision.

Where this sits relative to prior work

The locality of operator algebras under finite-range or exponentially decaying interactions is old: Lieb–Robinson bounds (1972) gave the first speed-of-light statement, Hastings and others gave exponential clustering of correlations, and the toolkit of "quasi-locality" — that operators evolved for short times can be approximated by operators on local neighbourhoods — is well-established in the rigorous-physics literature.

Sabit's contribution, as I read it, is to lift those operator-level statements to the level of spectra as sets. It's a small conceptual move with a satisfying consequence: every classical result about exponential decay of perturbations gets a spectral cousin via Weyl's inequality, packaged in the Hausdorff metric. The disjoint-additivity statement in particular is something I hadn't seen written down so cleanly before — it's the spectral counterpart of "energies of non-interacting subsystems add," now made quantitative with an exponential error in the inter-block distance.

The model

I implemented the bounds for an XXZ chain on $\Lambda = \{0, 1, \ldots, N-1\}$, $N=9$, with one-body Z fields and two-body couplings that decay exponentially in distance:

with $J=1$, $\lambda=1$, $\Delta=0.5$. I want $\|\Phi\|_\mu < \infty$ for some $\mu > 0$, which means $\mu < \lambda$. I picked $\mu = 0.5$, comfortably inside the convergent regime. For 1-D, $\mathrm{diam}(\{i,j\}) = j-i$.

The full Hilbert space is $2^9 = 512$ dimensional — a doddle for np.linalg.eigvalsh. The interesting subsystem $H_S$ for the additivity test acts on $|S_1|+|S_2| = 6$ qubits, dimension 64.

subspectra.py — the framework

The module is small enough to read end-to-end. It encodes interactions as lists of Pauli terms tagged with their support, builds dense operators by Kronecker products, and computes everything else (operator norms, Hausdorff distances, Minkowski sums) on top.

"""
subspectra.py
=============

Subsystem-resolved spectral theory for finite-range and exponentially-decaying
quantum many-body Hamiltonians, after Sabit (arXiv:2604.21929).

The model lives on a 1-D chain Λ = {0, ..., N-1} of qubits.  An "interaction"
Φ is given as a list of Pauli terms (subset, pauli_string, coefficient); the
Hamiltonian is

    H = Σ_{X ⊆ Λ} Φ(X)

where Φ(X) is the sum of all interaction terms whose support is exactly X.

For a subset S ⊆ Λ we build

    H_S      = Σ_{X ⊆ S} Φ(X)
    H_{S,r}  = Σ_{X ⊆ S, diam(X) ≤ r} Φ(X)

and compute their spectra σ(H_S), σ(H_{S,r}) on the |S|-qubit Hilbert space.

The μ-norm of Φ is

    ||Φ||_μ = sup_{x ∈ Λ} Σ_{X ∋ x} ||Φ(X)|| · e^{μ diam(X)}

where ||·|| is the operator norm.  No Hausdorff metric package is needed —
the spectra are finite sets of reals so d_H(A,B) is two nested mins/maxes.
"""

from __future__ import annotations
import numpy as np
from dataclasses import dataclass, field
from typing import Iterable, Sequence

# ---------------------------------------------------------------------------
# Pauli ops
# ---------------------------------------------------------------------------

I2 = np.eye(2, dtype=complex)
PX = np.array([[0, 1], [1, 0]], dtype=complex)
PY = np.array([[0, -1j], [1j, 0]], dtype=complex)
PZ = np.array([[1, 0], [0, -1]], dtype=complex)
PAULI = {"I": I2, "X": PX, "Y": PY, "Z": PZ}


def kron_list(mats):
    out = mats[0]
    for m in mats[1:]:
        out = np.kron(out, m)
    return out


# ---------------------------------------------------------------------------
# Interactions
# ---------------------------------------------------------------------------


@dataclass(frozen=True)
class Term:
    """A single Pauli term: support is `sites` (sorted tuple), Pauli letters
    in `pauli` are aligned to `sites` in order, and `coeff` is the prefactor.
    """

    sites: tuple
    pauli: str
    coeff: complex


@dataclass
class Interaction:
    terms: list = field(default_factory=list)

    # ---- builders ---------------------------------------------------------

    def add(self, sites: Iterable[int], pauli: str, coeff: complex) -> None:
        sites = tuple(sorted(sites))
        assert len(sites) == len(pauli), (sites, pauli)
        self.terms.append(Term(sites, pauli, complex(coeff)))

    # ---- queries ----------------------------------------------------------

    def supports(self) -> set:
        return {t.sites for t in self.terms}

    def restrict(self, S: Iterable[int]) -> "Interaction":
        S = frozenset(S)
        return Interaction([t for t in self.terms if set(t.sites) <= S])

    def truncate_diam(self, r: int) -> "Interaction":
        out = []
        for t in self.terms:
            d = 0 if len(t.sites) <= 1 else t.sites[-1] - t.sites[0]
            if d <= r:
                out.append(t)
        return Interaction(out)

    # ---- operator construction -------------------------------------------

    def operator(self, sites_order: Sequence[int]) -> np.ndarray:
        """Return the dense operator on the Hilbert space of `sites_order`
        (a sorted list).  All terms must be supported within sites_order."""
        n = len(sites_order)
        dim = 1 << n
        H = np.zeros((dim, dim), dtype=complex)
        idx = {s: i for i, s in enumerate(sites_order)}
        for t in self.terms:
            ops = [I2] * n
            for s, ch in zip(t.sites, t.pauli):
                ops[idx[s]] = PAULI[ch]
            H += t.coeff * kron_list(ops)
        return H

    # ---- per-support norms ----------------------------------------------

    def support_norms(self) -> dict:
        """For each support X, return ||Φ(X)|| (operator norm)."""
        groups: dict = {}
        for t in self.terms:
            groups.setdefault(t.sites, []).append(t)
        out = {}
        for X, ts in groups.items():
            k = len(X)
            op = np.zeros((1 << k, 1 << k), dtype=complex)
            for t in ts:
                op += t.coeff * kron_list([PAULI[ch] for ch in t.pauli])
            out[X] = float(np.linalg.norm(op, ord=2))
        return out

    def mu_norm(self, mu: float, lattice: Sequence[int]) -> float:
        """sup_{x ∈ lattice} Σ_{X ∋ x} ||Φ(X)|| · e^{μ diam(X)}."""
        sn = self.support_norms()
        best = 0.0
        for x in lattice:
            s = 0.0
            for X, n in sn.items():
                if x in X:
                    d = 0 if len(X) <= 1 else X[-1] - X[0]
                    s += n * np.exp(mu * d)
            best = max(best, s)
        return best


# ---------------------------------------------------------------------------
# Spectra and Hausdorff distance
# ---------------------------------------------------------------------------


def spectrum(H: np.ndarray) -> np.ndarray:
    """Eigenvalues of a Hermitian operator, sorted ascending."""
    H = 0.5 * (H + H.conj().T)  # symmetrise to suppress numerical drift
    return np.linalg.eigvalsh(H)


def hausdorff(A: np.ndarray, B: np.ndarray) -> float:
    """One-sided sweeps in both directions; A and B are 1-D arrays of reals."""
    A = np.atleast_1d(A)
    B = np.atleast_1d(B)
    dab = np.max(np.min(np.abs(A[:, None] - B[None, :]), axis=1))
    dba = np.max(np.min(np.abs(B[:, None] - A[None, :]), axis=1))
    return float(max(dab, dba))


def minkowski_sum(A: np.ndarray, B: np.ndarray) -> np.ndarray:
    """Returns the *set* {a+b : a ∈ A, b ∈ B}.  We don't dedupe — Hausdorff
    works on multisets too."""
    return (A[:, None] + B[None, :]).reshape(-1)


# ---------------------------------------------------------------------------
# Demo: 1-D XXZ chain with exponentially-decaying long-range pair couplings
# ---------------------------------------------------------------------------


def build_xxz_long_range(
    N: int,
    *,
    decay: float = 1.0,
    delta: float = 0.5,
    coupling: float = 1.0,
    field_amp: float = 0.4,
    seed: int = 7,
) -> Interaction:
    """
    Φ({i}) = h_i Z_i with h_i ~ U[-field_amp, field_amp]
    Φ({i,j}) = J e^{-decay |i-j|} (X_i X_j + Y_i Y_j + Δ Z_i Z_j) for i<j
    """
    rng = np.random.default_rng(seed)
    inter = Interaction()
    for i in range(N):
        h = float(rng.uniform(-field_amp, field_amp))
        inter.add([i], "Z", h)
    for i in range(N):
        for j in range(i + 1, N):
            J = coupling * np.exp(-decay * (j - i))
            inter.add([i, j], "XX", J)
            inter.add([i, j], "YY", J)
            inter.add([i, j], "ZZ", J * delta)
    return inter

A few design notes worth calling out:

• Terms, not operators. I store interactions as a flat list of (sites, pauli, coeff) records. That makes restrict(S) and truncate_diam(r) one-line filter operations, and it lets support_norms group by sites cheaply. The operator is materialised on demand, only when I need the spectrum or an operator norm.
• Hausdorff on multisets. Eigenvalue lists are the natural multiset, and the Hausdorff distance formula (two nested $\max\min$s) doesn't care about multiplicities — it sees the support set. So I never deduplicate, and I never sort outside eigvalsh (which already returns sorted).
• Symmetrise before diagonalising. 0.5 * (H + H.conj().T) adds machine-noise resilience: it kills antisymmetric drift below $10^{-15}$ that otherwise leaks into eigvalsh as imaginary eigenvalues. Cheap insurance.

Per-support norms — the raw material

Before checking the bounds, let's look at what $\|\Phi(X)\|$ actually is for our model. The XXZ pair operator has eigenvalues $\{2+\Delta,\ -\Delta,\ -\Delta,\ -2+\Delta\}$, so $\|XX+YY+\Delta\,ZZ\| = 2+|\Delta|$. With $\Delta=0.5$ and the $e^{-\lambda(j-i)}$ prefactor, we expect each diameter-$d$ pair to come in at $(2+\Delta)\,e^{-d} \approx 2.5\,e^{-d}$:

X              | diam | ||Φ(X)||
----------------------------------------
[0]            |  0   | 0.1001
[1]            |  0   | 0.3178
[2]            |  0   | 0.2205
...
[0, 1]         |  1   | 0.9197
[1, 2]         |  1   | 0.9197
...
[0, 2]         |  2   | 0.3383
[1, 3]         |  2   | 0.3383
...

Sanity: $2.5\,e^{-1} \approx 0.9197$ ✓ and $2.5\,e^{-2} \approx 0.3383$ ✓. The $\mu$-norm sums these (with the $e^{\mu d}$ weight) and takes the worst vertex. For the seed-7 chain at $\mu = 0.5$ I get

||Φ||_μ = 6.830095

Plug $\lambda=1$, $\mu=0.5$, $\Delta=0.5$ into the geometric-sum estimate $\max|h_i| + (2+|\Delta|)\sum_{d=1}^{\infty} 2 e^{-(\lambda-\mu)d}$ and you land in the same ballpark. The vertex $x=6$ is the supremum; I print all terms through it in the support-norms table when debugging.

run_demo.py — the verification

"""
Demonstration: verify the three claims of arXiv:2604.21929 numerically on a
small XXZ chain with long-range exponentially-decaying couplings.

  1. Operator-norm truncation:
        ||H_S - H_{S,r}|| ≤ |S| e^{-μr} ||Φ||_μ
  2. Spectral truncation (Hausdorff):
        d_H(σ(H_S), σ(H_{S,r})) ≤ |S| e^{-μr} ||Φ||_μ
  3. Disjoint additivity (Hausdorff):
        d_H(σ(H_{S1∪S2}), σ(H_{S1}) + σ(H_{S2})) ≤ (|S1|+|S2|) e^{-μD} ||Φ||_μ

We choose μ < decay so ||Φ||_μ is finite.
"""
from __future__ import annotations
import sys
sys.path.insert(0, "/labs-output")

import numpy as np
from subspectra import (
    Interaction, build_xxz_long_range,
    spectrum, hausdorff, minkowski_sum,
)


def opnorm(M: np.ndarray) -> float:
    return float(np.linalg.norm(M, ord=2))


def main():
    N = 9
    decay = 1.0
    delta = 0.5
    mu = 0.5

    inter = build_xxz_long_range(N, decay=decay, delta=delta,
                                 coupling=1.0, field_amp=0.4, seed=7)

    Lambda = list(range(N))
    mu_norm = inter.mu_norm(mu, Lambda)
    print(f"Lattice Λ = {Lambda}")
    print(f"μ = {mu}, decay λ = {decay}  (need μ < λ for ||Φ||_μ < ∞)")
    print(f"||Φ||_μ = {mu_norm:.6f}")
    print()

    # =====================================================================
    # 1. + 2.   Truncation of H_S as r grows, S = entire chain.
    # =====================================================================
    S = list(range(N))
    H_S = inter.restrict(S).operator(S)
    sp_S = spectrum(H_S)
    print(f"=== H_S, S = full chain (|S| = {len(S)}, dim = {1<<len(S)}) ===")
    print(f"σ(H_S):  min = {sp_S[0]:+.5f}, max = {sp_S[-1]:+.5f}, "
          f"#eigs = {len(sp_S)}")
    print()

    print("  r | ||H_S - H_{S,r}|| | d_H spec |  bound   |  norm/bound  spec/bound")
    print(" ---+--------------------+----------+----------+-----------------------")
    for r in range(0, N):
        inter_trunc = inter.restrict(S).truncate_diam(r)
        H_Sr = inter_trunc.operator(S)
        diff = opnorm(H_S - H_Sr)
        sp_r = spectrum(H_Sr)
        dh = hausdorff(sp_S, sp_r)
        bound = len(S) * np.exp(-mu * r) * mu_norm
        print(f"  {r} | {diff:18.6f} | {dh:8.6f} | {bound:8.4f} | "
              f"{diff/bound:8.4f}    {dh/bound:8.4f}")
    print()

    # =====================================================================
    # 3.  Approximate additivity for disjoint blocks.
    # =====================================================================
    print("=== Disjoint additivity ===")
    print("  D | ||H_{S1∪S2} - (H_S1+H_S2)|| | d_H spec | bound  | n/b   s/b")
    print(" ---+-----------------------------+----------+--------+-----------")
    S1 = (0, 1, 2)
    block = 3
    for D in range(1, N):
        S2_start = S1[-1] + D
        S2 = tuple(range(S2_start, S2_start + block))
        if S2[-1] >= N:
            break
        d_actual = min(abs(a - b) for a in S1 for b in S2)
        assert d_actual == D, (S1, S2, d_actual, D)
        union = tuple(sorted(set(S1) | set(S2)))
        sites_union = list(union)

        H_union = inter.restrict(sites_union).operator(sites_union)

        H_S1_local = inter.restrict(S1).operator(list(S1))
        H_S2_local = inter.restrict(S2).operator(list(S2))
        I1 = np.eye(H_S1_local.shape[0])
        I2 = np.eye(H_S2_local.shape[0])
        H_sep = np.kron(H_S1_local, I2) + np.kron(I1, H_S2_local)

        cross_norm = opnorm(H_union - H_sep)
        sp_union = spectrum(H_union)
        sp_sum = minkowski_sum(spectrum(H_S1_local), spectrum(H_S2_local))
        dh = hausdorff(sp_union, sp_sum)
        bound = (len(S1) + len(S2)) * np.exp(-mu * D) * mu_norm
        print(f"  {D} | {cross_norm:27.6f} | {dh:8.6f} | {bound:6.3f} | "
              f"{cross_norm/bound:5.3f}  {dh/bound:5.3f}")
    print()

    # =====================================================================
    # 4.  Finite-range sanity — additivity is exact once D > range.
    # =====================================================================
    print("=== Finite-range sanity (range = 2, expect exact at D > 2) ===")
    fr = Interaction()
    rng = np.random.default_rng(11)
    for i in range(N):
        fr.add([i], "Z", float(rng.uniform(-0.4, 0.4)))
    for i in range(N - 1):
        fr.add([i, i + 1], "XX", 1.0)
        fr.add([i, i + 1], "YY", 1.0)
    for i in range(N - 2):
        fr.add([i, i + 2], "XX", 0.3)
        fr.add([i, i + 2], "YY", 0.3)
    print("  D | cross-norm | d_H spec")
    print(" ---+------------+----------")
    for D in (1, 2, 3, 4, 5):
        S1 = (0, 1, 2)
        S2 = tuple(range(2 + D, 2 + D + 3))
        if S2[-1] >= N:
            continue
        d_actual = min(abs(a - b) for a in S1 for b in S2)
        union = tuple(sorted(set(S1) | set(S2)))
        H_union = fr.restrict(union).operator(list(union))
        H_S1_local = fr.restrict(S1).operator(list(S1))
        H_S2_local = fr.restrict(S2).operator(list(S2))
        I1 = np.eye(H_S1_local.shape[0])
        I2 = np.eye(H_S2_local.shape[0])
        H_sep = np.kron(H_S1_local, I2) + np.kron(I1, H_S2_local)
        cross = opnorm(H_union - H_sep)
        sp_union = spectrum(H_union)
        sp_sum = minkowski_sum(spectrum(H_S1_local), spectrum(H_S2_local))
        dh = hausdorff(sp_union, sp_sum)
        print(f"  {d_actual} | {cross:10.3e} | {dh:8.3e}")


if __name__ == "__main__":
    main()

A subtle point in the additivity test: I build $H_{S_1\cup S_2}$ on the 6-qubit Hilbert space sorted(S1 ∪ S2), and I build the "separated" reference as kron(H_S1_local, I_S2) + kron(I_S1, H_S2_local). These match site-for-site only because $S_1 \prec S_2$ on the 1-D lattice (every site of $S_1$ is to the left of every site of $S_2$), so the natural Kronecker order coincides with the sorted-site order on $S_1 \cup S_2$. If the blocks interleaved on the line, I'd need an explicit qubit permutation; I avoided that case on purpose to keep the demonstration honest about what it's checking.

Results: truncation

Running the demo on $N=9$, seed 7:

Lattice Λ = [0, 1, 2, 3, 4, 5, 6, 7, 8]
μ = 0.5, decay λ = 1.0  (need μ < λ for ||Φ||_μ < ∞)
||Φ||_μ = 6.830095

=== H_S, S = full chain (|S| = 9, dim = 512) ===
σ(H_S):  min = -4.65197, max = +5.01280, #eigs = 512

  r | ||H_S - H_{S,r}|| | d_H spec |  bound   |  norm/bound  spec/bound
 ---+--------------------+----------+----------+-----------------------
  0 |           4.665320 | 2.805436 |  61.4709 |   0.0759      0.0456
  1 |           1.586976 | 1.224597 |  37.2840 |   0.0426      0.0328
  2 |           0.533267 | 0.384355 |  22.6139 |   0.0236      0.0170
  3 |           0.203629 | 0.111510 |  13.7160 |   0.0148      0.0081
  4 |           0.068664 | 0.031257 |   8.3192 |   0.0083      0.0038
  5 |           0.018887 | 0.008261 |   5.0458 |   0.0037      0.0016
  6 |           0.004607 | 0.002404 |   3.0605 |   0.0015      0.0008
  7 |           0.000839 | 0.000431 |   1.8563 |   0.0005      0.0002
  8 |           0.000000 | 0.000000 |   1.1259 |   0.0000      0.0000

A few things to read off this:

• The bound $|S|\,e^{-\mu r}\,\|\Phi\|_\mu$ is comfortably loose — about 13× larger than the actual operator-norm truncation error at $r=0$ and the gap shrinks slowly (the norm/bound column). That's expected: the bound treats every term as worst-case, but most of the residue cancels in operator norm. A tighter bound would need cancellation-aware reasoning, which the paper doesn't attempt and which depends on the specific interaction.
• The spectral Hausdorff distance is always smaller than the operator-norm distance, never larger. That's Weyl: deleting a small operator perturbation moves each eigenvalue by at most $\|E\|$, but typically less. In our data the spectral distance is roughly half the operator-norm distance.
• At $r=8$ everything goes exactly to zero. With $|S|=9$ and a 1-D chain, every $X \subseteq S$ has $\mathrm{diam}(X) \le 8$, so $H_{S,8} = H_S$ identically. The bound, of course, doesn't know that — it still permits $1.13$. This is a useful reminder that "the bound never beats reality" is a one-way street.

The decay ratio between consecutive rows in the truncation table is $\exp(-\mu) = \exp(-0.5) \approx 0.607$, which the bound respects exactly. The actual residue decays faster (ratios around $0.3$–$0.4$) because for our model the support norms decay at rate $\lambda=1$, not $\mu=0.5$ — the $\mu$-norm is conservative on purpose.

Results: spectral additivity

=== Disjoint additivity ===
  D | ||H_{S1∪S2} - (H_S1+H_S2)|| | d_H spec | bound  | n/b   s/b
 ---+-----------------------------+----------+--------+-----------
  1 |                    1.278037 | 0.622607 | 24.856 | 0.051  0.025
  2 |                    0.470163 | 0.117428 | 15.076 | 0.031  0.008
  3 |                    0.172963 | 0.056251 |  9.144 | 0.019  0.006
  4 |                    0.063630 | 0.022828 |  5.546 | 0.011  0.004

Same story, more dramatic. The cross-term operator norm decays at rate $\exp(-\lambda) = 0.368$, the bound at rate $\exp(-\mu) = 0.607$, and the spectral Hausdorff distance comes in at a quarter to a tenth of the operator norm. The Minkowski-sum spectrum approximates the joint spectrum to better than 0.06 already at $D=4$, against a worst-case bound of 5.5.

To make this concrete, here's the actual eigenvalue ladder at $D=3$ for the lowest eight states:

σ(H_S1)         σ(H_S2)         σ(H_union)        σ(S1)+σ(S2) (sorted)
------          ------          ----------        ----------
-1.2646         -1.7191         -2.9711           -2.9837
-1.0110         -0.9008         -2.7466           -2.7301
-0.6505         -0.7271         -2.3664           -2.3696
-0.2028         +0.0562         -2.1703           -2.1654
-0.0238         +0.2755         -2.0064           -1.9918
+0.8411         +0.5956         -1.9114           -1.9220
+1.0739         +1.0674         -1.8956           -1.9118
+1.2376         +1.3524         -1.7448           -1.7429

d_H(union, sum) = 0.0563

Reading down the third and fourth columns: every joint eigenvalue is shadowed by some Minkowski-sum eigenvalue within ~0.02, and vice versa. The two columns are not element-by-element equal — the Minkowski-sum multiset is much larger ($8\times8=64$ entries vs 64 entries on a 6-qubit space; here they happen to coincide in size because $|S_1|=|S_2|=3$) and the union spectrum's ordering doesn't have to line up with the sorted Minkowski sum entry-by-entry. What matters is the set distance, not pointwise agreement.

Results: finite-range, exact at the threshold

The finite-range sanity check is the most satisfying number in the run. I built a separate Hamiltonian with only nearest-neighbour and next-nearest-neighbour two-body terms — interaction range $R=2$ — and re-ran the additivity test:

=== Finite-range sanity (range = 2, expect exact at D > 2) ===
  D | cross-norm | d_H spec
 ---+------------+----------
  1 |  2.332e+00 | 1.638e+00
  2 |  6.000e-01 | 2.129e-01
  3 |  2.220e-16 | 1.510e-14
  4 |  2.220e-16 | 1.776e-14

At $D=3 > R=2$, every cross-term is structurally zero — no $\Phi(X)$ has $\mathrm{diam}(X) \ge 3$ — so the cross-term operator norm is double-precision machine epsilon ($2.2\times 10^{-16}$), and the spectral Hausdorff distance is $1.5\times 10^{-14}$, the slightly inflated number that comes from eigvalsh rounding errors propagating through the Minkowski sum. This is the result the paper's finite-range corollary predicts, and it lands precisely where it should: at the support-radius threshold, with no margin and no exponential tail.

What I didn't do, and what's interesting next

Two honest caveats:

• No higher-body interactions. My demo uses only one- and two-body terms. The framework handles arbitrary $|X|$, and the implementation is general (Term admits any Pauli string), but I didn't exercise three-site or higher terms. Nothing about the proof depends on that, but a stress test would be a fairer demo if you're putting this in front of someone who actually does many-body.
• No 2-D or higher dimension. I used $\mathrm{diam}(X) = \max(X) - \min(X)$, which is correct on a 1-D chain. To run the bounds on a square lattice you'd swap that for the graph-theoretic diameter of $X$ (or just $\max_{i,j\in X}\|i-j\|$ in some chosen metric). The rest of the code is dimension-agnostic. The exponential decay of the bound with $r$ would still hold, with the constant depending on lattice-ball volume growth.

A concrete next-step that would make the paper interesting in production: combine subsystem-spectral truncation with shift-and-invert eigensolvers for the full $H$, using the truncated subsystem spectra as preconditioners or eigenvalue trackers. The Hausdorff bound tells you how close the cheap (truncated) spectrum is to the true one, which is exactly the kind of guarantee a hybrid solver wants.

Reproducing this

The whole thing — subspectra.py + run_demo.py — runs in about three seconds on a 2-core 2.7 GHz Ivy Bridge container with NumPy and SciPy. No GPU, no MPI, no proprietary anything. Set OPENBLAS_NUM_THREADS=2 (or =1) and run:

$ export PYTHONPATH=$(pwd)
$ python3 run_demo.py
Lattice Λ = [0, 1, 2, 3, 4, 5, 6, 7, 8]
μ = 0.5, decay λ = 1.0  (need μ < λ for ||Φ||_μ < ∞)
||Φ||_μ = 6.830095
...

Bump $N$ to 12 and you're at $4096\times 4096$ Hermitian eigenproblems — still fine on a laptop, takes ~10 s. Past 14, dense diagonalisation gets uncomfortable; you'd want sparse eigensolvers (scipy.sparse.linalg.eigsh) and to stop materialising operators on the full Hilbert space, but the bookkeeping in subspectra.py doesn't change.

What this paper changes

If you've worked with Lieb–Robinson bounds before, the operator-norm truncation result will look familiar — it's morally the same exponential-decay bookkeeping, applied to a slightly different object (subsystem Hamiltonian instead of time-evolved local operator). The novelty is twofold:

1. A clean spectral metric. Hausdorff distance of finite spectra plus Weyl's inequality is an obvious tool, but writing down approximate-additivity as a Hausdorff bound on Minkowski sums is the first place I've seen the move made explicitly. It opens the door to spectral analogues of every quasi-locality result in the literature.

2. Finite-range exactness. The corollary that additivity is exact for $D > R$ in finite-range models is, frankly, a nice piece of bookkeeping more than a deep result — but it's the kind of bookkeeping that makes a framework usable. A reader who reaches for this paper to bound the spectral error of a tensor-network or DMRG approximation gets a clean go/no-go criterion: if the cross-distance exceeds the interaction range, the spectra simply add.

The paper's title undersells this last point. "Subsystem-resolved spectral theory" sounds like a generic framework; what's there is a specific, easy-to-use spectral perturbation calculus for many-body Hamiltonians.

Citation

MD Nahidul Hasan Sabit, Subsystem-Resolved Spectral Theory for Quantum Many-Body Hamiltonians, arXiv preprint arXiv:2604.21929 [quant-ph], 2026. https://arxiv.org/abs/2604.21929

@misc{sabit2026subsystem,
  title         = {Subsystem-Resolved Spectral Theory for Quantum Many-Body Hamiltonians},
  author        = {Sabit, MD Nahidul Hasan},
  year          = {2026},
  eprint        = {2604.21929},
  archivePrefix = {arXiv},
  primaryClass  = {quant-ph},
  url           = {https://arxiv.org/abs/2604.21929},
}

Background reading

• Lieb, E. H., and Robinson, D. W. The finite group velocity of quantum spin systems. Comm. Math. Phys. 28 (1972) — the original locality bound for spin systems.
• Hastings, M. B., and Koma, T. Spectral gap and exponential decay of correlations. Comm. Math. Phys. 265 (2006) — the cleanest exponential-clustering theorem in the gapped regime.
• Bauer, F. L., and Fike, C. T. Norms and exclusion theorems. Numer. Math. 2 (1960) — for the Hermitian case, this is Weyl, but in the form any numerical analyst will recognise.

Artefacts

The two Python files (subspectra.py, run_demo.py, aux_tables.py) and the captured demo output (demo_output.txt, aux_output.txt) are bundled in the post's tarball. Total ~12 kB of code, ~3 kB of output. Everything else in this post is reproducible from the source files alone.