feat: add new scripts for S_m x SO(3) double collapse examples and common utilities

.gitignore

@@ -556,4 +556,7 @@ marimo/_lsp/
 __marimo__/
 
 # Streamlit
-.streamlit/secrets.toml
+.streamlit/secrets.toml
+
+
+paper/_added_diffs/

paper/references.bib

@@ -208,6 +208,52 @@
   url = {https://arxiv.org/abs/quant-ph/0406168}
 }
 
+@article{tothguhne2005stabilizer,
+  author = {T{\'o}th, G{\'e}za and G{\"u}hne, Otfried},
+  title = {Entanglement detection in the stabilizer formalism},
+  journal = {Physical Review A},
+  volume = {72},
+  number = {2},
+  pages = {022340},
+  year = {2005},
+  doi = {10.1103/PhysRevA.72.022340},
+  note = {arXiv:quant-ph/0501020}
+}
+
+@article{sharmachan2012spinadapted,
+  author = {Sharma, Sandeep and Chan, Garnet Kin-Lic},
+  title = {Spin-adapted density matrix renormalization group algorithms for quantum chemistry},
+  journal = {The Journal of Chemical Physics},
+  volume = {136},
+  number = {12},
+  pages = {124121},
+  year = {2012},
+  doi = {10.1063/1.3695642}
+}
+
+@article{menczer2024quarterpetaflops,
+  author = {Menczer, Andor and van Damme, Maarten and Rask, Alan and Huntington, Lee and Hammond, Jeff and Xantheas, Sotiris S. and Ganahl, Martin and Legeza, {\"O}rs},
+  title = {Parallel Implementation of the Density Matrix Renormalization Group Method Achieving a Quarter petaFLOPS Performance on a Single DGX-H100 GPU Node},
+  journal = {Journal of Chemical Theory and Computation},
+  volume = {20},
+  number = {19},
+  pages = {8397--8404},
+  year = {2024},
+  doi = {10.1021/acs.jctc.4c00903},
+  note = {arXiv:2407.07411}
+}
+
+@article{huberklepmagronvolcic2022werner,
+  author = {Huber, Felix and Klep, Igor and Magron, Victor and Vol{\v{c}}i{\v{c}}, Jurij},
+  title = {Dimension-Free Entanglement Detection in Multipartite Werner States},
+  journal = {Communications in Mathematical Physics},
+  volume = {396},
+  pages = {1051--1070},
+  year = {2022},
+  doi = {10.1007/s00220-022-04485-9},
+  note = {arXiv:2108.08720}
+}
+
 @article{heinEisertBriegel2004,
   author = {Hein, M. and Eisert, J. and Briegel, H. J.},
   title = {Multiparty entanglement in graph states},

paper/shadow_maps_symmetric_states.tex

@@ -527,7 +527,7 @@ numerically identical to the cut bound in \cite{aschauer2026a}, but now a statem
 \subsection{Axial $U(1)$ symmetry: an exact worked example}
 \label{sec:u1-example}
 
-The isotypic mechanism of Proposition~\ref{prop:block-diagonal} is easiest to see concretely for an abelian symmetry group, where it reduces to an ordinary charge-conservation selection rule. That a global $U(1)$ (or other abelian) symmetry splits the singular values of a bipartite correlation-type matrix into charge sectors is by now a recurring theme in the entanglement-detection literature under the name \emph{symmetry-resolved entanglement}: it has been used to sharpen partial-transpose-moment witnesses \cite{nevenetal2021}, to decompose the singular values of the realignment matrix itself into $U(1)$-charge blocks \cite{rathetal2023}, and, in a very recent and structurally close construction, to build symmetry-resolved bounds on the negativity directly from charge-sector-resolved realignment-matrix singular values \cite{tarabungahaug2025}. The mechanism below is the same charge-conservation selection rule applied instead to the bigraduated shadow map $\mathcal M_S(\rho)$, worked out completely and in closed form for one illustrative four-qubit family rather than as a general witness-improvement scheme.
+The isotypic mechanism of Proposition~\ref{prop:block-diagonal} is easiest to see concretely for an abelian symmetry group, where it reduces to an ordinary charge-conservation selection rule. That a global $U(1)$ (or other abelian) symmetry splits the singular values of a bipartite correlation-type matrix into charge sectors is by now a recurring theme in the entanglement-detection literature under the name \emph{symmetry-resolved entanglement}: it has been used to sharpen partial-transpose-moment witnesses \cite{nevenetal2021}, to decompose the singular values of the realignment matrix itself into $U(1)$-charge blocks \cite{rathetal2023}, and, in a very recent and structurally close construction, to build symmetry-resolved bounds on the negativity directly from charge-sector-resolved realignment-matrix singular values \cite{tarabungahaug2025}. That literature block-diagonalizes the reduced state $\rho_A$ itself (or a moment thereof) and works almost exclusively with abelian charges; Proposition~\ref{prop:block-diagonal} instead block-diagonalizes the linear \emph{response map} $\mathcal M_S(\rho)$, and, in Section~\ref{sec:full-su2-symmetry} below, does so for the non-abelian group $SU(2)$ via the Casimir operator rather than a conserved charge. The mechanism below is the same charge-conservation selection rule applied instead to the bigraduated shadow map $\mathcal M_S(\rho)$, worked out completely and in closed form for one illustrative four-qubit family rather than as a general witness-improvement scheme.
 
 \begin{remark}[Abelian symmetries as selection rules]
 \label{rem:abelian-selection-rule}
@@ -603,6 +603,9 @@ This example illustrates Corollary~\ref{cor:joint-refinement} concretely: the fu
 % ============================================================
 
 Remark~\ref{rem:abelian-selection-rule} already anticipated that the non-abelian case needs the Casimir operator, not a Cartan generator, and stops short of working this out. We do so here, for the strong hypothesis of Section~\ref{sec:mechanisms} specialized further: not merely party-local, but \emph{collective}, i.e.\ the same single-qubit unitary $U_g$ acting on every party of the full system $P$ simultaneously,
+
+The representation-theoretic content of what follows --- Casimir diagonalization in place of group averaging, and Wigner--Eckart reduction of a tensor to a Clebsch--Gordan coefficient times a smaller reduced operator --- is not new machinery on its own. It is standard practice under the name \emph{spin-adapted DMRG} in ab-initio quantum chemistry, where exploiting the non-abelian $SU(2)$ spin symmetry of the electronic Hamiltonian via exactly this Wigner--Eckart factorization \cite{sharmachan2012spinadapted} routinely enables density-matrix-renormalization-group calculations on active spaces with well over a hundred orbitals, e.g.\ CAS$(113,76)$ on GPU hardware \cite{menczer2024quarterpetaflops}. We are not aware of this toolbox having been carried through as systematically for the entanglement-witness setting of this note, where the reduced object of interest is a witness-relevant response map rather than a Hamiltonian or an MPS tensor; the two settings share the group theory but differ in what is being reduced and why. A related but distinct use of collective-unitary symmetry appears in \cite{huberklepmagronvolcic2022werner}, which constructs dimension-free entanglement witnesses for Werner states (invariant under the diagonal conjugate action of the full unitary group, not merely collective single-qubit rotation) via a semidefinite-programming hierarchy; the guiding idea --- exploit a fixed symmetry to obtain a structural rather than dimension-dependent result --- is the same in spirit, but the symmetry group, the object being reduced, and the resulting statement (existence of a witness, versus the block structure of a given response map) are all different.
+
 \begin{equation}
 g\ \longmapsto\ \bigotimes_{a\in P}U_g,
 \qquad
@@ -860,7 +863,7 @@ satisfying $3A_1+7A_3=1+14=15=\norm{M_{S\to S^c}(\ketbra QQ)}_*$ exactly, matchi
 
 \subsection{Exact degeneracy from stabilizer structure}
 
-The degeneracies recorded for $\GHZ_3$, the Smolin state, and all $38$ four-qubit graph states share no continuous symmetry of the kind used above. Their common origin is instead a discrete, combinatorial fact about \emph{Pauli-diagonal} states, requiring only elementary group theory over $\mathbb F_2$, and it is this mechanism --- not Proposition~\ref{prop:block-diagonal} --- that is responsible for every degeneracy reported in \cite{aschauer2026a}. That reduced density matrices of stabilizer states are maximally mixed on their support, with the entanglement across any cut given by the rank of the corresponding restriction of the stabilizer group, is a classical fact \cite{fattal2004stabilizer}; Lemma~\ref{lem:code-support} below is the flat-support statement underlying that result, and Lemma~\ref{lem:stabilizer-degeneracy} translates it directly into the singular-value structure of the shadow map itself, rather than into an entanglement entropy.
+The degeneracies recorded for $\GHZ_3$, the Smolin state, and all $38$ four-qubit graph states share no continuous symmetry of the kind used above. Their common origin is instead a discrete, combinatorial fact about \emph{Pauli-diagonal} states, requiring only elementary group theory over $\mathbb F_2$, and it is this mechanism --- not Proposition~\ref{prop:block-diagonal} --- that is responsible for every degeneracy reported in \cite{aschauer2026a}. That reduced density matrices of stabilizer states are maximally mixed on their support, with the entanglement across any cut given by the rank of the corresponding restriction of the stabilizer group, is a classical fact \cite{fattal2004stabilizer}; the same $\mathbb F_2$-linear stabilizer-restriction toolbox underlies the entanglement witnesses of Tóth and Gühne \cite{tothguhne2005stabilizer}, though their target object is a witness operator rather than the singular-value structure of a response map. Lemma~\ref{lem:code-support} below is the flat-support statement underlying that result, and Lemma~\ref{lem:stabilizer-degeneracy} translates it directly into the singular-value structure of the shadow map itself, rather than into an entanglement entropy or a witness.
 
 \paragraph{Setup.} Identify each single-party Pauli index with $\mathbb F_2^2$ via $I\mapsto(0,0)$, $X\mapsto(1,0)$, $Y\mapsto(1,1)$, $Z\mapsto(0,1)$, so that an $n$-party Pauli string $\sigma_{\vec i}$ corresponds to $\vec i\in\mathbb F_2^{2n}$, and string multiplication (up to phase) becomes addition. Let $H\le\mathbb F_2^{2n}$ be an isotropic subgroup (i.e.\ its elements pairwise commute as operators) not containing $-I$, and let
 \[

scripts/symmetric_states/_from_other_discussion/01_double_collapse_table.py

@@ -0,0 +1,101 @@
+"""
+01_double_collapse_table.py
+
+Reproduces and extends Example ex:dicke-network of shadow_maps_symmetric_states.tex
+(sec:sm-so3-combination): for the doubly-symmetric state
+
+    |Q_m> = N * sum_{k=0}^m (-1)^k |D_m^k>_S ox |D_m^{m-k}>_S^c
+
+(S,S^c each m qubits, S_m-symmetric on each side AND fully collectively
+SO(3)-invariant), the fully active block M_{S->S^c}(|Q_m><Q_m|) collapses,
+in two independent steps, from 3^m x 3^m down to O(m) scalar channels:
+
+  Step 1 (Prop. multinomial-collapse): 3^m x 3^m -> binom(m+2,2) x binom(m+2,2)
+          via the S_m orbit-type embedding.
+  Step 2 (Prop. harmonic-decomposition + Cor. projector-norm-formula):
+          binom(m+2,2) -> exactly floor(m/2)+1 SCALAR values A_j, since
+          Sym^m(R^3) decomposes multiplicity-free into spherical harmonics
+          H_j, j = m, m-2, ..., (0 or 1).
+
+This script runs the check for m = 3, 4, 5 and prints a summary table.
+Only the m=3 case appears in the current paper draft (with values A_1, A_3
+whose SIGN should be checked against this script -- see README).
+
+Run:  python3 01_double_collapse_table.py
+"""
+
+import numpy as np
+from shadow_su2_common import (
+    build_Q, verify_collective_invariance, reduced_correlation_matrix,
+    total_J2, isotypic_eigenbasis, isotypic_block_report
+)
+
+np.set_printoptions(precision=6, suppress=True)
+
+
+def run_for_m(m):
+    print("=" * 70)
+    print(f"m = {m}  (source/target cluster size), j0 = {m}/2 = {m/2}")
+    print("=" * 70)
+
+    Q = build_Q(m)
+
+    # sanity: full collective SO(3) invariance
+    overlaps = verify_collective_invariance(Q, n_legs=2*m, n_trials=3, seed=m)
+    print(f"collective invariance check |<Q|U^ox{2*m}|Q>| (should be 1.0): "
+          f"{[f'{o:.10f}' for o in overlaps]}")
+
+    # Step 1: multinomial collapse
+    Chat, types = reduced_correlation_matrix(Q, m)
+    nuclear_norm_reduced = np.linalg.svd(Chat, compute_uv=False).sum()
+    print(f"Step 1 (multinomial collapse): 3^{m}x3^{m} -> "
+          f"{len(types)}x{len(types)},  ||hat-C||_* = {nuclear_norm_reduced:.6f}")
+
+    # Step 2: SO(3) Casimir block-diagonalization within the reduced space
+    J2 = total_J2(m)
+    # restrict J2 (3^m x 3^m) to the binom(m+2,2)-dim symmetric subspace via
+    # the SAME orbit embedding used to build Chat (Sym^m(R^3) = image of that embedding)
+    from shadow_su2_common import build_orbit_embedding
+    U, _ = build_orbit_embedding(m)
+    J2_red = U.T @ J2 @ U
+
+    evals, evecs, groups = isotypic_eigenbasis(J2_red)
+    results = isotypic_block_report(Chat, groups, evecs, label=f"Step 2, m={m}")
+
+    # nuclear norm consistency check: sum_j (2j+1)|A_j| == ||hat-C||_*
+    total = 0.0
+    Aj_summary = {}
+    for j, (block, sv) in results.items():
+        Aj = np.diag(block).mean()
+        Aj_summary[j] = Aj
+        total += (2*j+1) * abs(Aj)
+    print(f"consistency: sum_j (2j+1)|A_j| = {total:.6f}  "
+          f"vs ||hat-C||_* = {nuclear_norm_reduced:.6f}  "
+          f"(match: {np.isclose(total, nuclear_norm_reduced)})")
+    return Aj_summary
+
+
+if __name__ == "__main__":
+    all_results = {}
+    for m in (3, 4, 5):
+        all_results[m] = run_for_m(m)
+        print()
+
+    print("=" * 70)
+    print("SUMMARY TABLE  (exact A_j values, m = 3,4,5)")
+    print("=" * 70)
+    for m, Aj in all_results.items():
+        j0 = m/2
+        items = ", ".join(f"A_{j:g}={v:+.6f}" for j, v in sorted(Aj.items(), reverse=True))
+        print(f"m={m} (j0={j0}):  {items}")
+
+    print()
+    print("Pattern check across m:")
+    print("  sign(A_j) should be (-1)^m for every j")
+    print("  top ratio A_m / A_{m-2} should equal 2m exactly")
+    for m, Aj in all_results.items():
+        js = sorted(Aj.keys(), reverse=True)
+        if len(js) >= 2:
+            top, second = js[0], js[1]
+            ratio = Aj[top]/Aj[second]
+            print(f"  m={m}: A_{top:g}/A_{second:g} = {ratio:.6f}  (expect {2*m})")

scripts/symmetric_states/_from_other_discussion/02_subsector_multiplicity_example.py

@@ -0,0 +1,80 @@
+"""
+02_subsector_multiplicity_example.py
+
+Extends Example ex:dicke-network beyond the FULLY ACTIVE block M_{S->S^c}.
+For m=4 (S = {A,B,C,D}, S^c = {E,F,G,H}), this looks at the SUB-sector block
+M_{V->T} with V = {A,B,C} subsetneq S (qubit D a "spectator") and
+T = {E,F,G} subsetneq S^c (qubit H a spectator).
+
+Unlike the fully active case, V (and T) here live in the FULL (unsymmetrized)
+(R^3)^{ox 3}, which has multiplicities m_j = 1,3,2,1 for j=0,1,2,3
+(Prop. branching-su2). In particular m_1 = 3, so the reduced block A_1 need
+NOT be a scalar -- this is the first concrete state in the note where a
+genuinely non-trivial (non-identity) reduced matrix A_j appears, rather than
+the always-scalar case forced by the multiplicity-free Sym^m(R^3) of
+Example ex:dicke-network.
+
+Finding (see README): A_3 is scalar (mult 1, as expected), A_1 is a genuine
+RANK-1 3x3 matrix (multiplicity 3, one nonzero singular value repeated 3x),
+and A_2, A_0 vanish identically for this particular state/cut -- an
+unexplained selection rule, flagged as an open question, not yet understood.
+
+Run:  python3 02_subsector_multiplicity_example.py
+"""
+
+import numpy as np
+from shadow_su2_common import build_Q, apply_leg, PAULIS, total_J2, isotypic_eigenbasis, isotypic_block_report
+from itertools import product
+
+np.set_printoptions(precision=4, suppress=True, linewidth=140)
+
+
+def main():
+    m = 4
+    Q = build_Q(m)   # legs 0,1,2,3 = A,B,C,D (source S); legs 4,5,6,7 = E,F,G,H (target S^c)
+
+    V_legs = [0, 1, 2]   # {A,B,C} subsetneq S;  D = leg 3 left inactive (spectator)
+    T_legs = [4, 5, 6]   # {E,F,G} subsetneq S^c; H = leg 7 left inactive (spectator)
+
+    # raw (unreduced) 27x27 block M_{V->T} -- no multinomial collapse here,
+    # since V,T are not the full active sector and the S_m symmetry does not
+    # act on a fixed 3-out-of-4 subset the same simple way (cf. Remark
+    # permutation-separating-example: only the residual S_3-within-V x S_3-within-T
+    # symmetry survives, and we do not exploit it here -- this is the raw block).
+    c = np.zeros((3, 3, 3, 3, 3, 3), dtype=complex)
+    for iA, iB, iC, iE, iF, iG in product(range(3), repeat=6):
+        ket = Q
+        for leg, ip in zip(V_legs, [iA, iB, iC]):
+            ket = apply_leg(ket, leg, PAULIS[ip])
+        for leg, ip in zip(T_legs, [iE, iF, iG]):
+            ket = apply_leg(ket, leg, PAULIS[ip])
+        c[iA, iB, iC, iE, iF, iG] = np.vdot(Q, ket).real
+
+    M = c.reshape(27, 27)
+    nuc = np.linalg.svd(M, compute_uv=False).sum()
+    print(f"V = {{A,B,C}} subsetneq S = {{A,B,C,D}},  T = {{E,F,G}} subsetneq S^c = {{E,F,G,H}}")
+    print(f"||M_{{V->T}}||_* = {nuc:.6f}\n")
+
+    # Casimir decomposition of the FULL (unsymmetrized) (R^3)^{ox3} on each side
+    J2 = total_J2(3)
+    evals, evecs, groups = isotypic_eigenbasis(J2)
+    print("multiplicities on full (R^3)^{ox3} (no S_3 symmetrization):",
+          {j: len(idx) for j, idx in groups.items()}, "  (expect {3:1, 2:2, 1:3, 0:1})\n")
+
+    results = isotypic_block_report(M, groups, evecs, label="isotypic blocks of M_{V->T}")
+
+    print("\nFull reduced matrices (rounded) and their singular values:")
+    total_check = 0.0
+    for j in sorted(results, reverse=True):
+        block, sv = results[j]
+        print(f"\nj={j}: block =\n{block}")
+        print(f"  singular values: {sv}")
+        total_check += (2*j+1) * 0  # placeholder, real check below using sv directly if block scalar
+    # nuclear-norm consistency: sum over ALL singular values of all blocks == ||M||_*
+    all_sv = np.concatenate([sv for (_, sv) in results.values()])
+    print(f"\nsum of all block singular values = {all_sv.sum():.6f}  vs ||M||_* = {nuc:.6f}  "
+          f"(match: {np.isclose(all_sv.sum(), nuc)})")
+
+
+if __name__ == "__main__":
+    main()

scripts/symmetric_states/_from_other_discussion/README.md

@@ -0,0 +1,72 @@
+# Scripts: S_m × SO(3) double collapse (fehlende Resultate im Entwurf)
+
+Diese Skripte reproduzieren und erweitern `Example ex:dicke-network` /
+`sec:sm-so3-combination` aus `shadow_maps_symmetric_states.tex`. Sie decken
+genau die Ergebnisse ab, die in der Team-Notiz vom [Datum der Sitzung] als
+"im Entwurf fehlend" markiert wurden.
+
+## Dateien
+
+- **`shadow_su2_common.py`** — gemeinsames Modul: Dicke-Zustände, der
+  doppelt-symmetrische Zustand `|Q_m>`, Pauli-Operatoren, die
+  Multinomial-Kollaps-Einbettung (Prop. `multinomial-collapse`), und die
+  exakten Casimir-Projektoren + Blockstruktur-Check (Remark
+  `casimir-projectors`, Cor. `projector-norm-formula`). Keine externen
+  Abhängigkeiten außer NumPy.
+
+- **`01_double_collapse_table.py`** — der volle aktive Block $M_{S\to S^c}$
+  für $m=3,4,5$. Reproduziert das bestehende $m=3$-Beispiel im Paper und
+  erweitert es auf $m=4,5$. Erzeugt die Tabelle der exakten $A_j$-Werte und
+  prüft zwei Muster:
+  - Vorzeichen $=(-1)^m$
+  - Top-Verhältnis $A_m/A_{m-2}=2m$ exakt
+
+  **Achtung:** Die Werte hier sind $A_1=-1/3$, $A_3=-2$ für $m=3$ —
+  **negativ**. Der aktuelle Paper-Entwurf gibt $A_1=+1/3$, $A_3=+2$ an
+  (positiv). Für dieselbe Zustandsdefinition sollte das Vorzeichen aber
+  eindeutig sein (kein Freiheitsgrad einer globalen Phase, da es sich um
+  einen Erwartungswert handelt). Bitte vor Übernahme ins Paper gegenprüfen —
+  betrifft nicht die im Text gezogene Konsequenz $3A_1+7A_3=15$, die nur
+  $|A_j|$ benutzt.
+
+- **`02_subsector_multiplicity_example.py`** — der Teilsektor-Fall
+  $V=\{A,B,C\}\subsetneq S=\{A,B,C,D\}$ (ein "Zuschauer"-Qubit $D$), analog
+  auf der Zielseite. Im Unterschied zu `01` lebt $V$ hier in der vollen,
+  nicht symmetrisierten $(\mathbb R^3)^{\otimes3}$ mit Multiplizitäten
+  $m_j=1,3,2,1$ für $j=0,1,2,3$. Ergebnis:
+  - $j=3$: skalar, $A_3=-0.8$ (Multiplizität 1, wie erwartet)
+  - $j=2$: **verschwindet identisch** (nicht nur klein — exakt Null,
+    unerklärt, siehe unten)
+  - $j=1$: **echte, nicht-skalare $3\times3$-Matrix**, Rang 1 (ein
+    Singulärwert $0{,}3$, dreifach über die Multiplizität)
+  - $j=0$: exakt Null
+
+  Das ist im aktuellen Entwurf komplett unbehandelt — bisher zeigt das
+  Paper nur den multiplizitätsfreien Fall $V=S$, bei dem $A_j$ zwangsläufig
+  skalar ist. Dieses Beispiel ist der erste konkrete Beleg für echte
+  Multiplizitätsraum-Struktur am doppelt-symmetrischen Zustand.
+
+## Offene Punkte (nicht in den Skripten gelöst)
+
+1. **Geschlossene Form für $A_j(m)$**: Nur numerisch gemustert (Tabelle in
+   `01`), nicht hergeleitet. Nächster Schritt wäre eine saubere
+   Clebsch-Gordan-Herleitung via `sympy` (z. B. über den "gestreckten
+   Zustand"-Trick), nicht durch Raten aus den drei Datenpunkten.
+2. **Warum verschwinden $j=2,0$ im Teilsektor-Beispiel (`02`)?** Exakt
+   Null, nicht nur klein — deutet auf eine zusätzliche Auswahlregel hin
+   (Parität? Eigenschaft der spezifischen $Q_4$-Konstruktion?). Nicht
+   untersucht.
+3. Das Vorzeichen-Problem oben (Punkt zu `01`) sollte geklärt werden, bevor
+   das Beispiel im Paper erweitert wird.
+
+## Ausführen
+
+```bash
+pip install -r requirements.txt
+python3 01_double_collapse_table.py
+python3 02_subsector_multiplicity_example.py
+```
+
+Beide Skripte sind eigenständig lauffähig (importieren nur
+`shadow_su2_common.py` aus demselben Verzeichnis) und laufen jeweils in
+wenigen Sekunden.

scripts/symmetric_states/_from_other_discussion/requirements.txt

@@ -0,0 +1 @@
+numpy>=1.20

scripts/symmetric_states/_from_other_discussion/shadow_su2_common.py

@@ -0,0 +1,241 @@
+"""
+shadow_su2_common.py
+
+Shared building blocks for the S_m x SO(3) double-collapse numerics
+(Section "Combining with permutation symmetry: multiplicity-free channels",
+sec:sm-so3-combination, in shadow_maps_symmetric_states.tex).
+
+Provides:
+  - Dicke states and the doubly-symmetric ("singlet-of-two-multiplets") state
+    |Q_m> = sum_k (-1)^k |D_m^k>_S ox |D_m^{m-k}>_S^c
+  - Pauli matrices and leg-wise operator application on a rank-2m tensor
+  - The S_m orbit-type embedding (multinomial collapse, Prop. multinomial-collapse)
+  - Exact SO(3) Casimir projectors via J_x,J_y,J_z on (R^3)^{ox k}
+    (Remark casimir-projectors), and the resulting isotypic block check
+    (Corollary projector-norm-formula)
+
+No external dependencies beyond numpy.
+"""
+
+import numpy as np
+from itertools import product
+from math import factorial
+
+# ---------------------------------------------------------------------
+# 1. Dicke states and the doubly-symmetric invariant state |Q_m>
+# ---------------------------------------------------------------------
+
+def dicke(m, k):
+    """|D_m^k>: equal superposition of all weight-k bitstrings on m qubits."""
+    dim = 2**m
+    psi = np.zeros(dim, dtype=complex)
+    n_terms = 0
+    for bits in range(dim):
+        if bin(bits).count("1") == k:
+            psi[bits] = 1.0
+            n_terms += 1
+    psi /= np.sqrt(n_terms)
+    return psi
+
+
+def build_Q(m):
+    """
+    |Q_m> = N * sum_{k=0}^m (-1)^k |D_m^k>_S ox |D_m^{m-k}>_S^c ,  S,S^c each m qubits.
+
+    This is (up to overall phase) the unique SO(3) singlet formed by coupling
+    the two spin-j0=m/2 permutation-symmetric multiplets on S and S^c to J=0
+    (Example ex:dicke-network for m=3; here for general m).
+
+    Returns the state reshaped as a rank-(2m) tensor: legs 0..m-1 are the
+    source cluster S, legs m..2m-1 are the target cluster S^c.
+    """
+    dimS = 2**m
+    Q = np.zeros(dimS*dimS, dtype=complex)
+    for k in range(m+1):
+        Dk = dicke(m, k)
+        Dmk = dicke(m, m-k)
+        Q += ((-1)**k) * np.kron(Dk, Dmk)
+    Q /= np.linalg.norm(Q)
+    return Q.reshape([2]*(2*m))
+
+
+def verify_collective_invariance(Q_flat_dim_legs, n_legs, n_trials=5, seed=0):
+    """
+    Sanity check: |<Q|U^{ox n_legs}|Q>| == 1 for random single-qubit U in SU(2),
+    confirming full collective SO(3) invariance of the state.
+    Q_flat_dim_legs: the state as a rank-n_legs tensor.
+    """
+    rng = np.random.default_rng(seed)
+    Q = Q_flat_dim_legs
+    results = []
+    for _ in range(n_trials):
+        v = rng.standard_normal(4)
+        v /= np.linalg.norm(v)
+        a, b, c, d = v
+        U = np.array([[a+1j*b, c+1j*d], [-c+1j*d, a-1j*b]], dtype=complex)
+        Qrot = Q
+        for leg in range(n_legs):
+            Qrot = apply_leg(Qrot, leg, U)
+        overlap = np.vdot(Q, Qrot)
+        results.append(abs(overlap))
+    return results
+
+
+# ---------------------------------------------------------------------
+# 2. Pauli matrices and leg-wise application
+# ---------------------------------------------------------------------
+
+X = np.array([[0, 1], [1, 0]], dtype=complex)
+Y = np.array([[0, -1j], [1j, 0]], dtype=complex)
+Z = np.array([[1, 0], [0, -1]], dtype=complex)
+PAULIS = [X, Y, Z]   # index 0,1,2 <-> x,y,z
+
+
+def apply_leg(psi, axis, P):
+    """Apply a 2x2 (or dxd) operator P on the given tensor leg of psi."""
+    out = np.tensordot(P, psi, axes=([1], [axis]))
+    return np.moveaxis(out, 0, axis)
+
+
+# ---------------------------------------------------------------------
+# 3. Multinomial (S_m orbit-type) collapse -- Prop. multinomial-collapse
+# ---------------------------------------------------------------------
+
+def orbit_types(m):
+    """All (a_x,a_y,a_z) with a_x+a_y+a_z=m -- the |T_m| = binom(m+2,2) orbit types."""
+    return [(a, b, m-a-b) for a in range(m+1) for b in range(m+1-a)]
+
+
+def _type_of(idx_tuple):
+    cnt = [0, 0, 0]
+    for i in idx_tuple:
+        cnt[i] += 1
+    return tuple(cnt)
+
+
+def multinomial_coeff(alpha):
+    m = sum(alpha)
+    denom = 1
+    for a in alpha:
+        denom *= factorial(a)
+    return factorial(m) // denom
+
+
+def build_orbit_embedding(m):
+    """
+    Orthonormal embedding U: 3^m -> binom(m+2,2), columns u_alpha
+    (normalized indicator vectors of each S_m orbit), as in the proof of
+    Prop. multinomial-collapse.
+    """
+    types = orbit_types(m)
+    idx_m = list(product(range(3), repeat=m))
+    U = np.zeros((3**m, len(types)))
+    for col, alpha in enumerate(types):
+        mask = np.array([1.0 if _type_of(t) == alpha else 0.0 for t in idx_m])
+        U[:, col] = mask / np.sqrt(mask.sum())
+    return U, types
+
+
+def reduced_correlation_matrix(Q, m):
+    """
+    Build the multinomial-collapsed reduced matrix hat-C (Eq. reduced-dicke-matrix)
+    for the FULLY ACTIVE block M_{S->S^c}(rho) of the m+m cluster state Q,
+    without ever materializing the full 3^m x 3^m matrix explicitly (uses one
+    orbit representative per row/column instead of a brute-force loop).
+    """
+    types = orbit_types(m)
+    ntypes = len(types)
+    Chat = np.zeros((ntypes, ntypes))
+    for a_col, alpha in enumerate(types):
+        rep_src = []
+        for i, cnt in enumerate(alpha):
+            rep_src += [i]*cnt
+        ket = Q
+        for leg in range(m):
+            ket = apply_leg(ket, leg, PAULIS[rep_src[leg]])
+        for b_col, beta in enumerate(types):
+            rep_tgt = []
+            for i, cnt in enumerate(beta):
+                rep_tgt += [i]*cnt
+            ket2 = ket
+            for leg in range(m):
+                ket2 = apply_leg(ket2, m+leg, PAULIS[rep_tgt[leg]])
+            val = np.vdot(Q, ket2).real
+            Chat[a_col, b_col] = val * np.sqrt(multinomial_coeff(alpha) * multinomial_coeff(beta))
+    return Chat, types
+
+
+# ---------------------------------------------------------------------
+# 4. Exact SO(3) Casimir projectors (Remark casimir-projectors) and the
+#    isotypic block check (Corollary projector-norm-formula)
+# ---------------------------------------------------------------------
+
+_EPS = np.zeros((3, 3, 3))
+_EPS[0, 1, 2] = _EPS[1, 2, 0] = _EPS[2, 0, 1] = 1
+_EPS[0, 2, 1] = _EPS[2, 1, 0] = _EPS[1, 0, 2] = -1
+J1_CARTESIAN = [-1j*_EPS[a] for a in range(3)]   # spin-1 generator, single leg
+
+
+def _kron_n(mats):
+    out = mats[0]
+    for M in mats[1:]:
+        out = np.kron(out, M)
+    return out
+
+
+def total_J2(k):
+    """J^2_tot = sum_a (sum_l J_a on leg l)^2 acting on (R^3)^{ox k}, as a 3^k x 3^k matrix."""
+    I3 = np.eye(3)
+    dim = 3**k
+    Jtot = [np.zeros((dim, dim), dtype=complex) for _ in range(3)]
+    for a in range(3):
+        for leg in range(k):
+            mats = [I3]*k
+            mats[leg] = J1_CARTESIAN[a]
+            Jtot[a] += _kron_n(mats)
+    return sum(Ja @ Ja for Ja in Jtot)
+
+
+def isotypic_eigenbasis(J2_matrix):
+    """
+    Diagonalize J^2 (real part), return (evals, evecs, groups) where groups
+    maps j -> list of eigenvector indices spanning that isotypic component.
+    Eigenvalues are j(j+1); j is recovered via j=(-1+sqrt(1+4*ev))/2.
+    """
+    J2_real = J2_matrix.real
+    evals, evecs = np.linalg.eigh(J2_real)
+    groups = {}
+    for i, ev in enumerate(evals):
+        j = round((-1 + np.sqrt(1 + 4*max(ev, 0))) / 2, 3)
+        groups.setdefault(j, []).append(i)
+    return evals, evecs, groups
+
+
+def isotypic_block_report(M, groups, evecs, label=""):
+    """
+    Rotate M into the J^2 eigenbasis and report, per isotype j:
+      - block dimension
+      - whether it is a scalar multiple of the identity (mean diag, spread, off-diag)
+      - max leakage into other isotypes (should vanish -- Cor. projector-norm-formula)
+    Returns dict j -> (block matrix, singular values).
+    """
+    M_rot = evecs.T @ M @ evecs
+    results = {}
+    if label:
+        print(f"--- {label} ---")
+    for j in sorted(groups, reverse=True):
+        idxs = groups[j]
+        block = M_rot[np.ix_(idxs, idxs)]
+        diag = np.diag(block)
+        offdiag_within = block - np.diag(diag)
+        offblock_max = 0.0
+        for j2 in groups:
+            if j2 == j:
+                continue
+            offblock_max = max(offblock_max, np.abs(M_rot[np.ix_(idxs, groups[j2])]).max())
+        sv = np.linalg.svd(block, compute_uv=False) if len(idxs) > 1 else np.abs(diag)
+        print(f"  j={j:.1f}  dim={len(idxs):2d}  mean(diag)={diag.mean(): .8f}  "
+              f"spread(diag)={diag.std():.2e}  |offdiag|max={np.abs(offdiag_within).max():.2e}  "
+              f"|leak to other j|max={offblock_max:.2e}")
+        results[j] = (block, sv)
+    return results