"""Fock-metric dissipativity test for the adapted heat generator (Remark r:stab-fock, Section 16). Question. Let K_h be the gauged generator gamma_h^{-1}(L_N + L_H(U))gamma_h of Section 18 (s:cons-heat), one step being exp(delta K_h), delta = h^4, with the onsite half-weight d_quad outside. Is Re _Fock <= C ||u||_Fock^2 uniformly in h ? Units. Physical normal coordinates U, physical heat time; one step has time delta = W^{-2} = h^4 for the generator (1/2) Delta. The quadratic normal data are those of Section 16 (s:tube-oscillator): r_sigma = E/2 + sigma i s_E, H = 1 - r^2, c = (W^2 H + H^2/4)^{1/2}, mu = W^2/(W^2 + H/2 + c), z = -log mu, lambda = exp(-i arg(c)/2). After the rotation U = lambda x (x real) the vacuum Psi = exp(-|c| x^2/2) is real, so the Fock metric L^2(Psi^2) on Wick functions is the flat metric on the raw functions Psi u. Part A: ordinary rank-one model, normal (u_+, u_-), one tangent direction y on a circle, row-geometry metric g_S = g_0 (1 - i w)^2 with w = (U_- - U_+)/(2 s). Variant: one complex normal coordinate with a U(1) connection D = d_y + A(y) J and a non-equivariant metric. Part B: supersymmetric normal factor: one coordinate (x, y | alpha, beta), STr U^2 = x^2 + y^2 + 2 alpha beta, Wick metric of the calibration. Part C: the odd stabilizer charge in the Fock metric. Run: python3 fock_dissipativity.py (about 35 minutes on 4 cores). """ import numpy as np import scipy.sparse as sp import scipy.sparse.linalg as spla from scipy.linalg import expm, eigh, sqrtm np.set_printoptions(precision=4, suppress=True) # ---------------------------------------------------------------------------- # Quadratic normal data # ---------------------------------------------------------------------------- def root_data(E, h, sigma): sE = np.sqrt(1 - E ** 2 / 4) r = E / 2 + sigma * 1j * sE H = 1 - r ** 2 W = h ** -2 c = np.sqrt(W ** 2 * H + H ** 2 / 4) # principal branch, Re c > 0 mu = W ** 2 / (W ** 2 + H / 2 + c) z = -np.log(mu) lam = np.exp(-1j * np.angle(c) / 2) return dict(sE=sE, H=H, c=c, a=abs(c), mu=mu, z=z, lam=lam) def top_herm(K): """Numerical abscissa: largest eigenvalue of (K + K^*)/2 (flat metric).""" Hm = (K + K.conj().T) / 2 if sp.issparse(Hm): return spla.eigsh(Hm.tocsc(), k=1, which="LA", tol=1e-8)[0][0].real return np.linalg.eigvalsh(Hm)[-1] def fd1d(n, R): """Dirichlet grid on (-R, R): points, second and first difference.""" dx = 2 * R / (n + 1) x = -R + dx * np.arange(1, n + 1) e = np.ones(n) D2 = sp.diags([e[:-1], -2 * e, e[:-1]], [-1, 0, 1]) / dx ** 2 D1 = sp.diags([-e[:-1], e[:-1]], [-1, 1]) / (2 * dx) return x, D2.tocsr(), D1.tocsr() # ---------------------------------------------------------------------------- # Part A: ordinary rank-one joint generator # ---------------------------------------------------------------------------- def gauged1d(d, x, D2, D1): """Psi^{-1} (1/2 lam^{-2} d^2) Psi with Psi = exp(-a x^2/2), a = |c|.""" a = d["a"]; X = sp.diags(x); I = sp.identity(len(x)) return 0.5 * d["lam"] ** -2 * (D2 - 2 * a * X @ D1 + a ** 2 * X @ X - a * I) def partA(E=0.8, r0=0.1, hs=(0.04, 0.02, 0.01, 0.005, 0.0025), n1=1600, n2=60, kmax=4): print("=" * 78) print(f"Part A: ordinary rank-one model, E = {E}, enlarged tube |x_sigma| < R = 2 r0, r0 = {r0}") R = 2 * r0 # (i) tangent mode k = 0 carries the top of both real parts (Re m > 0), and there # the generator separates; 1-d Dirichlet grids with n1 points. x, D2, D1 = fd1d(n1, R) print(" (i) k = 0, both normal coordinates: h | Fock abscissa | flat abscissa of the gauged" " operator | potential top 1/2 cos(arg c)(|c|^2R^2-|c|) summed") rows = [] for h in hs: fock = flat = pot = 0.0 for sg in (+1, -1): d = root_data(E, h, sg) fock += top_herm((0.5 * d["lam"] ** -2 * D2).tocsr()) flat += top_herm(gauged1d(d, x, D2, D1).tocsr()) pot += 0.5 * np.cos(np.angle(d["c"])) * (d["a"] ** 2 * R ** 2 - d["a"]) rows.append((h, fock, flat, pot)) print(f" h={h:7.4f} Fock {fock:11.4e} flat-gauged {flat:11.4e} potential {pot:11.4e}" f" flat*h^4/r0^2 = {flat*h**4/r0**2:.3f}") hh = np.array([r[0] for r in rows]); fl = np.array([r[2] for r in rows]) sl = np.polyfit(np.log(hh[-3:]), np.log(fl[-3:]), 1)[0] print(f" slope of log(flat-gauged) against log h over the last three h: {sl:.3f} (r0^2 h^-4: -4)") h = hs[2] print(f" r0-dependence at h = {h}:") for rr in (0.05, 0.1, 0.2): xr, D2r, D1r = fd1d(n1, 2 * rr) flat = sum(top_herm(gauged1d(root_data(E, h, sg), xr, D2r, D1r).tocsr()) for sg in (+1, -1)) fock = sum(top_herm((0.5 * root_data(E, h, sg)["lam"] ** -2 * D2r).tocsr()) for sg in (+1, -1)) print(f" r0={rr:5.2f}: Fock {fock:11.4e} flat-gauged {flat:11.4e} flat*h^4/r0^2 = {flat*h**4/rr**2:.3f}") # (ii) coupled 2-d problem for tangent modes k >= 1 (row-geometry metric). x, D2, D1 = fd1d(n2, R) I = sp.identity(n2, format="csr") xp, xm = np.meshgrid(x, x, indexing="ij") print(f" (ii) coupled normal plane, tangent modes 0..{kmax}, grid {n2}^2:") for h in (0.04, 0.01): dp, dm = root_data(E, h, +1), root_data(E, h, -1) s = dp["sE"]; g0 = 8 * s ** 2 w = (dm["lam"] * xm - dp["lam"] * xp) / (2 * s) m = (1 - 1j * w) ** -2 assert m.real.min() > 0 LN = 0.5 * dp["lam"] ** -2 * sp.kron(D2, I) + 0.5 * dm["lam"] ** -2 * sp.kron(I, D2) LNg = sp.kron(gauged1d(dp, x, D2, D1), I) + sp.kron(I, gauged1d(dm, x, D2, D1)) out = [] for k in range(1, kmax + 1): LH = -(k ** 2) / (2 * g0) * sp.diags(m.ravel()) out.append((k, top_herm((LN + LH).tocsr()), top_herm((LNg + LH).tocsr()))) print(f" h={h}: " + "; ".join(f"k={k}: Fock {a:.3e}, flat-gauged {b:.3e}" for k, a, b in out)) h = hs[-1] dp, dm = root_data(E, h, +1), root_data(E, h, -1) s = dp["sE"]; g0 = 8 * s ** 2 m = (1 - 1j * (dm["lam"] * xm - dp["lam"] * xp) / (2 * s)) ** -2 print(f" (iii) Fock-metric decomposition at h = {h} (closed forms on the tube):") print(f" OU part 1/2 lam^-2 (d^2 - 2|c| x d) abscissa 0, attained on the vacuum") print(f" gauge potential 1/2 lam^-2 (|c|^2 x^2 - |c|) abscissa {rows[-1][3]:.3e} (r0^2 h^-4)") print(f" OU + potential abscissa <= 0 (ground-state transform)") print(f" horizontal at U = 0 abscissa 0") print(f" coupling -(k^2/2g0)(m(U) - 1) abscissa (k^2/2g0) * {(1-m.real).max():.3e};" f" relative to the horizontal part: max(1 - Re m) = {(1-m.real).max():.3e} < 1") # (iv) one full step chi d exp(delta K) d chi in the Fock metric. h = 0.04 dp, dm = root_data(E, h, +1), root_data(E, h, -1) s = dp["sE"]; g0 = 8 * s ** 2 m = (1 - 1j * (dm["lam"] * xm - dp["lam"] * xp) / (2 * s)) ** -2 LN = 0.5 * dp["lam"] ** -2 * sp.kron(D2, I) + 0.5 * dm["lam"] ** -2 * sp.kron(I, D2) d = np.exp(-(dp["H"] * (dp["lam"] * xp) ** 2 + dm["H"] * (dm["lam"] * xm) ** 2) / 4).ravel() chi = (np.maximum(abs(xp), abs(xm)) <= r0).ravel().astype(float) delta = h ** 4 norms = [] for k in (0, 2): K = (LN - (k ** 2) / (2 * g0) * sp.diags(m.ravel())).tocsc() KH = K.conj().T.tocsc() step = lambda v: chi * d * spla.expm_multiply(delta * K, d * chi * v) stepH = lambda v: chi * np.conj(d) * spla.expm_multiply(delta * KH, np.conj(d) * chi * v) op = spla.LinearOperator((n2 * n2,) * 2, matvec=lambda v: stepH(step(v)), dtype=complex) norms.append((k, np.sqrt(spla.eigsh(op, k=1, which="LA", tol=1e-10)[0][0].real))) print(f" (iv) ||chi d exp(delta K) d chi||_Fock at h = {h}: " + ", ".join(f"k={k}: {v:.8f}" for k, v in norms) + f" (sup|d| = {abs(d).max():.6f})") return rows def partA_connection(E=0.8, r0=0.1, h=0.03, n=41, K=4, A0=0.7, A1=0.9): """One complex normal coordinate (x, y) with U(1) rotation J, a tangent circle, D = d_t + A(t) J, A(t) = A0 + A1 cos t, and the non-equivariant metric factor m(U) = (1 - i x lam/(2 s))^{-2}. [J, m] != 0 and [Delta_N, m] != 0.""" print("=" * 78) print("Part A': compact U(1) connection on a 2-dim normal plane, one tangent circle") R = 2 * r0 x, D2, D1 = fd1d(n, R) I = sp.identity(n, format="csr") X = sp.diags(x) d = root_data(E, h, +1) s = d["sE"]; g0 = 8 * s ** 2 lap = sp.kron(D2, I) + sp.kron(I, D2) J = sp.kron(X, D1) - sp.kron(D1, X) # x d_y - y d_x, skew ks = np.arange(-K, K + 1) nt = len(ks) It = sp.identity(nt, format="csr") Dt = sp.diags(1j * ks) Aop = A0 * It + 0.5 * A1 * (sp.eye(nt, k=1) + sp.eye(nt, k=-1)) D = sp.kron(Dt, sp.identity(n * n)) + sp.kron(Aop, J) xx, yy = np.meshgrid(x, x, indexing="ij") m = (1 - 1j * d["lam"] * xx / (2 * s)) ** -2 M = sp.kron(It, sp.diags(m.ravel())) LN = 0.5 * d["lam"] ** -2 * sp.kron(It, lap) LH = (D @ M @ D) / (2 * g0) skew = abs((D + D.conj().T)).max() print(f" ||D + D^*||_max = {skew:.1e} (D skew), ||[J, m]|| > 0: " f"{abs((J @ sp.diags(m.ravel()) - sp.diags(m.ravel()) @ J)).max():.2e}") print(f" Fock abscissa of L_N + L_H : {top_herm((LN + LH).tocsr()):.4e}" f" (Dirichlet ground state of the normal heat: " f"{-np.cos(np.angle(d['lam']**-2))*(np.pi/(2*R))**2:.4e})") # Re L_H = -(1/2g0) D^* Re(m) D <= 0 since D is skew; check the identity. HermLH = ((LH + LH.conj().T) / 2).tocsr() ident = HermLH + (D.conj().T @ sp.kron(It, sp.diags(m.real.ravel())) @ D).tocsr() / (2 * g0) print(f" ||Re L_H + (1/2g0) D^* Re(m) D||_max = {abs(ident).max():.1e}; min Re m = {m.real.min():.4f} > 0") # ---------------------------------------------------------------------------- # Part B: supersymmetric normal factor, Wick metric of Section 16 (s:tube-oscillator) # ---------------------------------------------------------------------------- # Grassmann algebra on alpha, beta: raw basis (1, alpha, beta, alpha beta). mul_a = np.zeros((4, 4)); mul_a[1, 0] = 1; mul_a[3, 2] = 1 mul_b = np.zeros((4, 4)); mul_b[2, 0] = 1; mul_b[3, 1] = -1 der_a = np.zeros((4, 4)); der_a[0, 1] = 1; der_a[2, 3] = 1 der_b = np.zeros((4, 4)); der_b[0, 2] = 1; der_b[1, 3] = -1 AB = mul_a @ mul_b # multiplication by alpha beta LapF = 2 * der_a @ der_b # Delta_F for STr = ... + 2 alpha beta I4 = np.eye(4) def fermion_metric(a): """Gram matrix (raw coordinates) of the Wick metric: the raw vectors Psi_F * {1, sqrt(2a) alpha, sqrt(2a) beta, 2a alpha beta + 1}, with Psi_F = exp(-a alpha beta), are orthonormal (E alpha beta = -1/(2a)).""" PsiF = I4 - a * AB one = np.array([1, 0, 0, 0.]) B = np.column_stack([PsiF @ one, PsiF @ (np.sqrt(2 * a) * mul_a @ one), PsiF @ (np.sqrt(2 * a) * mul_b @ one), PsiF @ (2 * a * AB @ one + one)]) Binv = np.linalg.inv(B) return Binv.conj().T @ Binv def metric_top(Xop, G): """Abscissa of Xop in the metric _G = u^* G v.""" S = np.real_if_close(sqrtm(G)) Si = np.linalg.inv(S) Y = S @ Xop @ Si return np.linalg.eigvalsh((Y + Y.conj().T) / 2)[-1] def metric_norm(Xop, G): S = np.real_if_close(sqrtm(G)) return np.linalg.norm(S @ Xop @ np.linalg.inv(S), 2) def hermite_ops(nmax): """Annihilation b on Hermite functions 0..nmax (flat metric).""" b = np.diag(np.sqrt(np.arange(1, nmax + 1)), 1) return b def partB(E=0.8, hs=(0.08, 0.04, 0.02, 0.01), r0=0.1): print("=" * 78) print(f"Part B: supersymmetric normal factor, E = {E}") for h in hs: d = root_data(E, h, +1) a, lam, H, mu, z = d["a"], d["lam"], d["H"], d["mu"], d["z"] t = h ** 4 G = fermion_metric(a) heatF = I4 + t * lam ** -2 * der_a @ der_b # e^{(t/2) lam^-2 Delta_F} dF = I4 - 0.5 * H * lam ** 2 * AB # exp(-H lam^2 alpha beta / 2) TF = dF @ heatF @ dF ev = np.sort_complex(np.linalg.eigvals(TF)) # normality of T_F in the Wick metric S = np.real_if_close(sqrtm(G)); Y = S @ TF @ np.linalg.inv(S) normal_defect = np.linalg.norm(Y @ Y.conj().T - Y.conj().T @ Y) genF = 0.5 * lam ** -2 * LapF # per unit time absF = metric_top(genF, G) print(f" h={h:5.3f}: eig T_F = {np.round(ev, 6)}; mu = {mu:.6f}, 1/mu = {1/mu:.6f};" f" ||[T,T*]||_Wick = {normal_defect:.1e}") print(f" odd heat abscissa (per unit delta) = {absF:.4e}; |c|/2 = {a/2:.4e};" f" ratio to h^-2 = {absF*h**2:.4f}") print(f" per step: ||heat_F||-1 = {metric_norm(heatF, G)-1:.3e}, ||d_F||-1 = " f"{metric_norm(dF, G)-1:.3e}, ||T_F|| = {metric_norm(TF, G):.8f}, e^(Re z) = {np.exp(z.real):.8f}") # Full SUSY normal heat on Hermite(x) x Hermite(y) x fermion, Wick metric. print(" frozen SUSY normal generator in the Wick metric, bosons truncated at a tube of" f" radius 2 r0, r0 = {r0}:") for h in hs[:3]: d = root_data(E, h, +1) a, lam = d["a"], d["lam"] nmax = int(a * (2 * r0) ** 2 / 2) # |x| <= 2 r0 <-> n <~ a R^2 / 2 nmax = min(nmax, 400) b = hermite_ops(nmax + 2) P = slice(0, nmax + 1) D2 = (a / 2) * ((b - b.T) @ (b - b.T))[P, P] # raw d^2 in Hermite basis, scale a bos_top = np.linalg.eigvalsh(0.5 * np.cos(np.angle(lam ** -2)) * D2)[-1] G = fermion_metric(a) absF = metric_top(0.5 * lam ** -2 * LapF, G) # product metric: abscissa = 2 * bosonic(1 coordinate) + fermionic print(f" h={h:5.3f}: nmax={nmax:4d} bosonic part {2*bos_top:.3e} odd part {absF:.3e}" f" total {2*bos_top+absF:.4e} (|c|/2 = {a/2:.4e})") # Invariant sector: u = phi(STr U^2), annihilated by every charge and body-constant, # so L_{H,h} u = 0 exactly and = . print(" invariant test vectors phi(S) = exp(-S/eps), S = x^2 + y^2 + 2 alpha beta," " eps = kappa/|c|:") for h in (0.04, 0.02, 0.01): d = root_data(E, h, +1) a, lam = d["a"], d["lam"] best = None for kap in np.linspace(2.0, 12.0, 201): eps = kap / a # raw f = A + B alpha beta, A = exp(-rho/eps), B = -(2/eps) A, rho = x^2 + y^2. # Wick norm^2 = (1/2)(||A||^2 + a^-2 ||B||^2); flat 2-d integrals in closed form. nA = np.pi * eps / 2; nB = (2 / eps) ** 2 * nA norm2 = 0.5 * (nA + nB / a ** 2) # quadratic form of (1/2) Delta_N in the Wick metric (rotation factor lam^-2): # Q = (1/2) + (1/2) a^-2 gA = np.pi * 1.0 # ||grad A||^2 = pi (any eps) in 2-d gB = (2 / eps) ** 2 * gA Q = 0.5 * (-0.5 * gA - (-(2 / eps)) * nA) + 0.5 / a ** 2 * (-0.5 * gB) val = (lam ** -2 * Q).real / norm2 if best is None or val > best[0]: best = (val, kap) print(f" h={h:5.3f}: max_kappa Re/||u||^2 = {best[0]:.4e} at kappa = {best[1]:.2f};" f" / |c| = {best[0]/a:.4f}; * h^2 = {best[0]*h**2:.4f}") # Metric independence among product metrics: spectral radius of T_F. print(" random metrics G on the odd factor: ||d_F||_G^2 ||heat_F||_G >= rho(T_F) = e^(Re z):") h = 0.02 d = root_data(E, h, +1) a, lam, H, z = d["a"], d["lam"], d["H"], d["z"] t = h ** 4 heatF = I4 + t * lam ** -2 * der_a @ der_b dF = I4 - 0.5 * H * lam ** 2 * AB rng = np.random.default_rng(1) worst = np.inf for _ in range(2000): Z = rng.standard_normal((4, 4)) + 1j * rng.standard_normal((4, 4)) Z = Z @ np.diag(np.exp(rng.uniform(-6, 6, 4))) Gr = Z.conj().T @ Z S = sqrtm(Gr); Si = np.linalg.inv(S) val = np.linalg.norm(S @ dF @ Si, 2) ** 2 * np.linalg.norm(S @ heatF @ Si, 2) worst = min(worst, val) print(f" h = {h}: min over 2000 random metrics = {worst:.8f}; e^(Re z) = {np.exp(z.real):.8f};" f" in the Wick metric = {metric_norm(dF, fermion_metric(a))**2*metric_norm(heatF, fermion_metric(a)):.8f}") # ---------------------------------------------------------------------------- # Part C: the odd charge q = sqrt2 (f_b^* a - f_a b^*) in the Fock metric # ---------------------------------------------------------------------------- def partC(nmax=40): print("=" * 78) print("Part C: odd stabilizer charge q (Section 16) in the Fock metric") nb = nmax + 1 a = sp.diags(np.sqrt(np.arange(1, nb)), 1, format="csr") Ib = sp.identity(nb, format="csr") # two fermion modes (alpha, beta) by Jordan-Wigner f = sp.csr_matrix(np.array([[0, 1], [0, 0.]])); Z2 = sp.diags([1, -1.]); I2 = sp.identity(2) fa = sp.kron(f, I2); fb = sp.kron(Z2, f) I4s = sp.identity(4) A = sp.kron(sp.kron(a, Ib), I4s); Bm = sp.kron(sp.kron(Ib, a), I4s) Fa = sp.kron(sp.kron(Ib, Ib), fa); Fb = sp.kron(sp.kron(Ib, Ib), fb) q = (np.sqrt(2) * (Fb.T @ A - Fa @ Bm.T)).tocsr() nvals = (np.add.outer(np.add.outer(np.arange(nb), np.arange(nb)).ravel(), np.array([0, 1, 1, 2.]))).ravel() N = sp.diags(nvals) keep = np.where(nvals <= nmax - 2)[0] sub = lambda M: M.tocsr()[keep][:, keep] print(f" ||q^2|| = {abs(sub(q @ q)).max():.1e}; ||qq^*+q^*q-2N|| = " f"{abs(sub(q @ q.T + q.T @ q - 2 * N)).max():.1e} on N <= {nmax-2}") qk = sub(q); nk = nvals[keep] for nn in (4, 16, 36): sel = np.where(nk <= nn)[0] Hq = ((qk + qk.T) / 2)[sel][:, sel] top = spla.eigsh(Hq.tocsc(), k=1, which="LA")[0][0] print(f" N <= {nn:3d}: abscissa of q = {top:.3f}; sqrt(N) = {np.sqrt(nn):.3f}") if __name__ == "__main__": partA() partA_connection() partB() partC()