# PLAN.md — Schwinger-Model Majorization Poset Extension for `tessera`

Hand this file to Claude Code as the working spec. Each phase has a
concrete deliverable, a file layout, and an acceptance test. Do not
skip acceptance tests; every phase gates the next.

The hypothesis this plan is built to test, plus falsification criteria,
limitations, and bibliographic context, are in
[quantum-methodology.md](quantum-methodology.md). The user-facing API
documentation for what's built so far is in [quantum.md](quantum.md).

## 0. Goal and scope

Extend `tessera` with a C++ subsystem that:

1. Represents the 1+1D Kogut–Susskind Schwinger-model Hamiltonian as
   a Matrix Product Operator (MPO) on a Jordan–Wigner'd spin chain.
2. Finds its ground state by DMRG.
3. Evolves it in real time by TDVP after a local quench creating a
   $q\bar{q}$ pair at separation $d$.
4. At each time slice, computes Schmidt spectra $\lambda_{[i,j]}(t)$
   for all contiguous-interval cuts $[i,j]$, builds a majorization
   poset, and compares it against (a) the Lieb–Robinson cone and
   (b) a causet causal order when the chain is embedded in a causet.

## 1. Architectural principle — minimize Python/C++ crossings

The entire simulation loop stays in C++. Python is a result viewer,
not a driver.

```
   C++ side:  build MPO  ->  DMRG  ->  TDVP loop  ->  spectra  ->  poset
                                                                    |
   Py side:                                                      export
                                                                    |
                                                    plot, save, interop
```

Only these cross the barrier:

- **In** (Py → C++): scalar config (N, m, g, Δt, T, cut family, seed).
- **Out** (C++ → Py): final/periodic snapshots of observables and
  the poset as plain data (`std::vector<double>`, adjacency lists).

Do not expose MPS tensors or MPOs to Python. Do not call Python
callbacks inside the time loop.

## 2. Dependencies

| Purpose | Library | Why |
|---|---|---|
| MPS/MPO, DMRG, TDVP | **ITensor (C++)** v3 | Mature, header-heavy C++, built-in DMRG+TDVP, handles quantum-number conservation |
| Linear algebra (small ops outside ITensor) | Eigen 3 | already likely in tessera's toolchain |
| Graph / poset | Boost.Graph | Hasse diagram, transitive reduction |
| Python bindings | pybind11 | already in tessera |
| Tests | Catch2 or doctest | whichever tessera already uses |

Add ITensor as a git submodule under `third_party/itensor` and wire
via CMake `add_subdirectory`. ITensor depends on BLAS/LAPACK — reuse
what tessera already has. Do not pull ITensor through a package manager;
pin a known-good commit.

## 3. File layout

```
tessera/
  include/tessera/quantum/
    schwinger_model.hpp        # MPO construction
    dmrg_runner.hpp            # ground-state driver
    tdvp_runner.hpp            # real-time driver
    schmidt.hpp                # spectrum extraction
    majorization.hpp           # partial order, poset
    causal_compare.hpp         # LR / tessera comparison
  src/quantum/
    schwinger_model.cpp
    dmrg_runner.cpp
    tdvp_runner.cpp
    schmidt.cpp
    majorization.cpp
    causal_compare.cpp
  tests/quantum/
    test_schwinger_spectrum.cpp   # Phase 1 acceptance
    test_majorization.cpp         # Phase 3 acceptance
    test_tdvp_string.cpp          # Phase 4 acceptance
    test_causal_compare.cpp       # Phase 5 acceptance
  python/tessera/quantum/
    __init__.py
    _bindings.cpp              # pybind11, only exports final data
  examples/quantum/
    run_schwinger.cpp
    qqbar_quench.cpp
```

## 4. Physics spec (fix once, reference everywhere)

Open-BC staggered-fermion Schwinger model, Jordan–Wigner'd,
Gauss-law-eliminated. One qubit per site, $N$ sites total.

$$H = H_\text{hop} + H_m + H_E$$

$$H_\text{hop} = \frac{1}{4a}\sum_{n=1}^{N-1}\bigl(X_n X_{n+1} + Y_n Y_{n+1}\bigr)$$

$$H_m = \frac{m}{2}\sum_{n=1}^{N}(-1)^{n}\,Z_n$$

$$H_E = \frac{g^2 a}{2}\sum_{n=1}^{N-1} L_n^2, \qquad L_n = L_0 + \sum_{k=1}^{n}\left[\frac{1-Z_k}{2} - \frac{1-(-1)^k}{2}\right]$$

with background field $L_0$ (keep as a parameter, default 0). Expand
$L_n^2$; the result is a long-range $Z_m Z_{m'}$ Hamiltonian whose
MPO bond dimension grows with $N$. ITensor's `AutoMPO` handles this
natively — use it, do not roll a hand-built MPO.

Reference for all numerics: Bañuls, Cichy, Cirac, Jansen,
*JHEP* **11**, 158 (2013).

## 5. Phased implementation

### Phase 0 — scaffolding (≤ 2 days)

- Add ITensor submodule; CMake build verifies with ITensor's own
  `sample/dmrg.cc` compiled against the tessera toolchain.
- Create the directory layout above with empty stubs.
- CI runs on each PR.

**Acceptance**: `make test_itensor_hello` returns the ground state
energy of the Heisenberg chain for $N=20$ to within $10^{-6}$ of the
ITensor reference value.

### Phase 1 — Schwinger MPO (2–3 days)

- Implement `schwinger_model.hpp`: a class that takes $(N, a, m, g, L_0)$
  and returns an `ITensor::MPO`.
- Use `AutoMPO` to avoid manually tracking the $L_n^2$ expansion.
- Conserve $U(1)$ total charge; quantum numbers on every index.

**Acceptance** (`test_schwinger_spectrum.cpp`):
DMRG ground-state energy per site for $N=20$, $a=1$, $m/g \in \{0.0,
0.125, 0.25\}$ matches Bañuls et al. 2013 Table 1 to $< 10^{-3}$.
Bond dimension 100 is enough for this check.

### Phase 2 — DMRG runner (2 days)

- `dmrg_runner.hpp`: thin wrapper around `ITensor::dmrg`.
  Inputs: `SchwingerMPO`, sweep schedule, max bond dim, noise, Krylov
  dim. Outputs: `MPS` ground state, converged energy, final bond dim.
- Expose `computeGroundState(config)` through pybind11 returning
  only $(E_0, \text{bond\_dim}, \text{truncation\_err})$.

**Acceptance**: re-runs Phase 1 with the wrapper; numerics unchanged.

### Phase 3 — Schmidt spectra and majorization poset (3 days)

- `schmidt.hpp`: given an MPS in right canonical form and an interval
  $[i,j]$, return $\lambda_{[i,j]}$ via SVD at the appropriate bond
  (ITensor exposes Schmidt values directly after orthogonalization).
- For contiguous cuts this is $O(N^2)$ spectra; store as a
  `std::vector<std::vector<double>>`.
- `majorization.hpp`:

  ```cpp
  // lambda <_maj mu iff, with both sorted descending, zero-padded,
  //   sum_{i=1..k} lambda_i  <=  sum_{i=1..k} mu_i   for all k,
  //   and total probability equal (trivially 1 here).
  bool majorizes(const std::vector<double>& mu,
                 const std::vector<double>& lambda,
                 double tol = 1e-12);
  ```

- Build a DAG whose nodes are intervals `[i,j]` and edges represent
  the majorization relation; apply Boost.Graph transitive reduction
  to obtain the Hasse diagram.
- Store as `tessera::Poset` (define this type; reuse tessera's existing
  causet adjacency structure if the API fits).

**Acceptance** (`test_majorization.cpp`):

1. Product state $|0\rangle^{\otimes N}$: every spectrum is
   $(1, 0, \ldots)$; Hasse diagram has no edges.
2. $N$-qubit GHZ: every single-site spectrum is $(\tfrac{1}{2},
   \tfrac{1}{2})$; pairwise majorization relations all equal (no
   strict edges).
3. Bell state across the center cut of $N=2$: spectrum
   $(\tfrac{1}{2}, \tfrac{1}{2})$; compared against product state
   $(\tfrac{1}{2}, \tfrac{1}{2}) \prec (1,0)$ correctly.

### Phase 4 — TDVP with $q\bar{q}$ quench (4 days)

- `tdvp_runner.hpp`: wraps ITensor's 2-site TDVP. Inputs: MPS, MPO,
  $\Delta t$, $T$, max bond dim, observables callback (pure C++).
- Initial state: DMRG ground state, then apply $\sigma^+_{i_0}
  \sigma^-_{i_0 + d}$ (normalized) to create a $q\bar{q}$ pair.
  Gauss's law requires appropriate string of $L$ shifts between them
  — handle via the electric-field string operator; reference Buyens
  et al. 2014 for the exact construction.
- At each recorded timestep, invoke observables callback that records:
  - $\langle L_n \rangle_t$ (electric-field profile)
  - $\langle Z_n \rangle_t$ (charge density)
  - $\lambda_{[i,j]}(t)$ for all contiguous cuts
  - full majorization poset

**Acceptance** (`test_tdvp_string.cpp`):

- Heavy-quark limit $m/g \gg 1$, $d=4$, $T = d\cdot a$ at $a=1$, $g=1$:
  $\langle L_n \rangle_t$ is approximately uniform at value $\pm 1$
  between $i_0$ and $i_0 + d$ (flux tube) and zero outside, to within
  $0.05$ after $t = T/2$.
- Energy conservation: $|\Delta E|/|E_0| < 10^{-3}$ over the run.

### Phase 5 — causal-order comparison (3 days)

- `causal_compare.hpp`: builds three partial orders on the label set
  $\{(\text{cut}, t)\}$:
  1. $\prec_\text{maj}$: from Phase 3 Hasse diagram, per time slice.
  2. $\prec_\text{LR}$: from the Lieb–Robinson velocity $v_\text{LR}$
     extracted by fitting OTOC ($\langle [\sigma^+_i(t), \sigma^-_j]^\dagger
     [\sigma^+_i(t), \sigma^-_j]\rangle$) front propagation.
  3. $\prec_\text{tessera}$: if the chain is embedded in a causet, inherit
     from tessera's existing causal order.
- Compute agreement statistics: Kendall-$\tau$, fraction of discordant
  pairs, graph edit distance on Hasse diagrams.
- Output: a single struct
  `struct CausalComparisonReport { double tau_maj_lr; double tau_maj_cs; ... }`

**Acceptance** (`test_causal_compare.cpp`):

- On a regular (total-order) causet chain: $\prec_\text{maj}$ agrees
  with $\prec_\text{LR}$ above some $v_E \leq v_\text{LR}$; both agree
  with $\prec_\text{tessera}$ since the causet is trivially totally
  ordered per time slice.
- Agreement rate quoted with uncertainty from bootstrap over Trotter
  seeds.

### Phase 6 — tessera integration (1–2 weeks, optional for v1)

**Status (as of 2026-04-25)**: data-extraction layer landed; MPO
rebuild on the extracted chain is the remaining piece.

Done:

- `tessera::Poset::fromSpacetime(Spacetime const&)` — inherits a Hasse
  cover Poset from the timelike-edge subgraph of any
  `tessera::Spacetime`. Uses transitive closure + reduction; preserves
  metric semantics by orienting earliest-time → latest-time and
  filtering on `Edge::getSquaredLength() < 0`.
- `tessera::quantum::extractCausetChain(Spacetime const&)` — packages
  the antichain layering, flat lattice ↔ Spacetime ID mapping,
  adjacent-slice timelike-edge hopping pairs, and the inherited
  partial-order Poset into a single `CausetChain`.
- Python bindings: `tessera.quantum.Poset.fromSpacetime(st)` and
  `tessera.quantum.extractCausetChain(st)` work on any
  `tessera.Spacetime` instance (CDT-built, hand-crafted, or future
  `tessera.CausalSet`).
- C++ acceptance: `tests/quantum/test_poset_from_spacetime.cpp` (7
  cases) and `tests/quantum/test_causet_chain.cpp` (5 cases). Python
  acceptance: `tests/quantum/test_phase6_causet_chain_python.py`.

Remaining for v1 close-out (deferred):

- Rebuild the Schwinger MPO on the `CausetChain` so
  `H_hop` follows the extracted hopping pairs rather than uniform
  nearest-neighbour. For trivial chains (one vertex per slice) this
  is just `params.N = chain.nSites` and the existing
  `buildSchwingerMpo` runs as-is.
- For non-trivial antichains (causet branches), the chain layout
  produces hopping pairs with stride > 1. Two paths:
    1. Stay on a 1D MPS, accept the long-range hopping cost in the
       MPO bond dimension.
    2. Switch to ITensor's tree-tensor-network support; tree TN
       handles branching antichains natively.
  Both paths are open work; pick after profiling the stride-1 case
  on a small CDT-derived chain.
- Replace ≼_cs in `computeCausalComparison` with the inherited
  `partialOrder` on the cut/time label set, then re-run the Phase 5
  invariants (`lrVsCs.kendallTau == 1.0` should still hold when
  the chain is totally ordered per slice; it may relax for branching
  causets, which is the Phase 6 hypothesis).

### Phase 7 — MERA / KAK / EML (research, not v1)

Do not attempt until Phases 1–5 produce publication-quality numbers.
Replace MPS with binary MERA; KAK-decompose disentanglers; EML-
parameterize KAK angles for symbolic regression over couplings.
Reference: Evenbly–Vidal 2009 for MERA algorithms; Hauru–Van Damme–
Haegeman 2021 for Riemannian optimization; Odrzywołek 2026 for EML.

## 6. Python layer (minimal)

`python/tessera/quantum/_bindings.cpp` exposes exactly:

```cpp
py::class_<QuantumConfig>(m, "QuantumConfig")
    .def(py::init<>())
    .def_readwrite("N", &QuantumConfig::N)
    .def_readwrite("a", &QuantumConfig::a)
    .def_readwrite("m", &QuantumConfig::m)
    .def_readwrite("g", &QuantumConfig::g)
    .def_readwrite("L0", &QuantumConfig::L0)
    .def_readwrite("maxBondDim", &QuantumConfig::maxBondDim)
    .def_readwrite("dt", &QuantumConfig::dt)
    .def_readwrite("T", &QuantumConfig::T);

py::class_<GroundStateResult>(m, "GroundStateResult")
    .def_readonly("energy", &GroundStateResult::energy)
    .def_readonly("bondDim", &GroundStateResult::bondDim)
    .def_readonly("truncationErr", &GroundStateResult::truncationErr);

py::class_<QuenchResult>(m, "QuenchResult")
    .def_readonly("times", &QuenchResult::times)
    .def_readonly("lProfile", &QuenchResult::lProfile)  // [t][n]
    .def_readonly("zProfile", &QuenchResult::zProfile)
    .def_readonly("posets", &QuenchResult::posets);       // one per t

m.def("computeGroundState", &computeGroundState);
m.def("runQqbarQuench", &runQqbarQuench);
```

No MPS, MPO, or ITensor type crosses the barrier.

## 7. Conventions and traps

- **Sign of $(-1)^n$ in $H_m$**: staggered convention, site index starts
  at $n=1$. Bañuls et al. use the same; match their sign.
- **Gauss's law**: after elimination, the long-range $Z_m Z_{m'}$ term
  has coefficients that depend on the parity of sites between $m$ and
  $m'$. Verify by independent sum on $N=4$ before trusting the MPO.
- **Truncation in TDVP**: 2-site TDVP grows bond dim; 1-site TDVP
  conserves it. Start with 2-site for the quench, switch to 1-site
  once bond dim saturates your budget.
- **Schmidt spectrum zeros**: pad to equal length before majorization
  comparison; otherwise partial-sum test is miscounted.
- **Boost.Graph transitive reduction**: not in Boost's public API
  directly; implement via DFS over comparability graph if needed.
- **OTOC computation**: expensive. Compute only at the final time and
  at one or two intermediate times for $v_\text{LR}$ estimation; do
  not run every step.

## 8. First week targets

- Day 1–2: Phase 0 done, ITensor builds inside tessera's CMake tree.
- Day 3–4: Phase 1 MPO constructed; matches Bañuls 2013 values.
- Day 5–6: Phase 2 wrapper; Phase 3 majorization infrastructure on
  trivial and GHZ test states.
- End of week 1: DMRG + majorization pipeline end-to-end on $N=20$
  Schwinger ground state.

## 9. Open questions to defer

- Whether to quantum-number-decompose spectra before majorization
  (majorization within each $U(1)$ sector vs. global). Default: global
  for v1; revisit.
- Whether to normalize spectra per cut (Schmidt coefficients squared
  sum to 1 automatically for normalized states; sanity-check this
  after TDVP truncation).
- Whether the "cut family" should include non-contiguous subsets
  eventually. Default: contiguous only, per the majorization
  well-definedness argument.

## 10. Non-goals for v1

- GPU acceleration of the tensor contractions. ITensor has partial
  CUDA support; ignore until CPU performance is a bottleneck.
- Non-abelian gauge groups ($SU(2)$, $SU(3)$). Schwinger is $U(1)$;
  that is the full v1 scope.
- 2+1D. Chain only for v1.
- Replacing ITensor. Do not write a bespoke MPS library.