"""n39 — the self-correcting vacuum, exact quantum test.
Z2 gauge-Higgs ring: 4 matter qubits (sites) + 4 link qubits. dim 256.
Gauss stabilizers G_i = sigma^x_i tau^x_{i-1,i} tau^x_{i,i+1}.
BMV-style state: dressed particle pair, second particle in a position
superposition (site 2 vs site 3), string included -> G=+1 in BOTH branches.
GATES (pre-registered):
 1 syndrome monitoring D[G_i]: coherence EXACTLY flat (rate 0), even with
   a gauge-invariant Hamiltonian on ([H,G]=0).
 2 field monitoring (tau^x on links, sigma^x on sites): coherence dies at
   2*gamma*(# differing X-eigenvalues) = 6*gamma for this pair. Predicted
   C(t) = exp(-6 gamma t), measured against it.
Corollary (stated): a syndrome-reading vacuum PROTECTS mass-position
superpositions -> BMV-positive, no DP floor; a field-reading vacuum is
book two (BMV-null, DP floor). The 2030s experiment picks the vacuum."""
import numpy as np

nq=8   # order: sigma_1..4 (sites), tau_12, tau_23, tau_34, tau_41 (links)
I2=np.eye(2); X=np.array([[0,1],[1,0]]); Z=np.array([[1,0],[0,-1]])
def op(k,M):
    out=np.array([[1.0]])
    for j in range(nq): out=np.kron(out,M if j==k else I2)
    return out
SX=[op(i,X) for i in range(4)]; SZ=[op(i,Z) for i in range(4)]
TX=[op(4+e,X) for e in range(4)]; TZ=[op(4+e,Z) for e in range(4)]
# Gauss: G_i = sx_i tx_{i-1,i} tx_{i,i+1}; links: e=i connects i,i+1 (mod 4)
G=[SX[i]@TX[(i-1)%4]@TX[i] for i in range(4)]

# vacuum |+>^8 : all X=+1
plus=np.ones(2)/np.sqrt(2)
psi0=np.array([1.0])
for _ in range(nq): psi0=np.kron(psi0,plus)
# branch a: particles at sites 1,2 (string on link 12=e index 0... site1-site2 is e=0? links: e=i connects i,i+1: site index 0..3; particles at site0,site1 string on link e0)
Sa=SZ[0]@TZ[0]@SZ[1]
Sb=SZ[0]@TZ[0]@TZ[1]@SZ[2]
pa=Sa@psi0; pb=Sb@psi0
for g in G:
    assert abs(pa@g@pa-1)<1e-12 and abs(pb@g@pb-1)<1e-12
print("both branches satisfy Gauss: G_i = +1 (dressed, strings included)")
psi=(pa+pb)/np.linalg.norm(pa+pb)
rho0=np.outer(psi,psi)

J,hf=0.4,0.3
H=-J*sum(SZ[i]@TZ[i]@SZ[(i+1)%4] for i in range(4))-hf*sum(TX)
assert max(np.abs(H@g-g@H).max() for g in G)<1e-12

def lind(rho,Ls,H=None):
    d=np.zeros_like(rho)
    if H is not None: d=-1j*(H@rho-rho@H)
    for L in Ls:
        LL=L.conj().T@L
        d=d+L@rho@L.conj().T-0.5*(LL@rho+rho@LL)
    return d
# NOTE: first version wrote L rho L - rho (missing gamma in the recoil term);
# it decayed everything at e^{-2t} and faked a gate-1 failure. Caught because
# the analytic prediction (D[G]rho = 0 on the code space) is exact.
def run(Ls,useH,T=4.0,dt=0.005,gam=0.5):
    rho=rho0.copy().astype(complex); a=pa.copy().astype(complex); b=pb.copy().astype(complex)
    Lg=[np.sqrt(gam)*L for L in Ls]; out=[]
    for s in range(int(T/dt)):
        Hc=H if useH else None
        k1=lind(rho,[np.sqrt(gam)*l for l in Ls],Hc)
        rho=rho+dt*k1  # Euler ok: linear, small dt; jumps unitary Hermitian
        if useH:
            a=a-1j*dt*(H@a); b=b-1j*dt*(H@b)
            a/=np.linalg.norm(a); b/=np.linalg.norm(b)
        if s%80==0: out.append(abs(np.conj(a)@rho@b))
    return np.array(out)

gam=0.5
print("\nGATE 1 — syndrome monitoring D[G_i], H off then on:")
c1=run(G,False); c2=run(G,True)
print(f"  H off: C(t)/C(0) final = {c1[-1]/c1[0]:.6f}   (prediction 1.000000)")
print(f"  H on : C(t)/C(0) final = {c2[-1]/c2[0]:.6f}   (prediction 1.000000)")
print("\nGATE 2 — field monitoring (all tau^x + sigma^x):")
c3=run(TX+SX,False)
ts=np.arange(len(c3))*0.005*80
pred=np.exp(-6*gam*ts)
err=np.abs(c3/c3[0]-pred).max()
print(f"  measured C(t)/C(0) vs exp(-6 gamma t): max dev = {err:.2e}")
print(f"  (rate = 2 gamma x 3 differing X-eigenvalues, as predicted)")
print("\nverdict: the syndrome vacuum protects the mass superposition;")
print("the field vacuum kills it. BMV distinguishes the two vacua.")