"""Check (iii) of 06-ilang-source/04-profile. (a) Field with one shared empty place, (4.14): one field object with places {empty} + an open chain of L places, own term C1 = Delta * sum_x |x>, compared with -g1 g2 exp(-d/xi) / (J sinh(1/xi)), cosh(1/xi) = Delta/(2J); decay length fitted from ln|M2|. (b) The same field beyond second order: exact M(d) by diagonalization, against its large-d limit g1^2 g2^2 (A11 B22 + A22 B11), A = K^{-1}, B = K^{-2}, K = C1 restricted to the chain (Step 7). (c) Cell field (4.20): one field object with places {0,1}^L (L cells of an open chain), own term Cf = Delta sum_y n_y - J sum_y (s+_y s-_{y+1} + h.c.), coupled operator Phi_x = sigma^x_x. Exact M(d) by sparse diagonalization, against M2 = -2 g1 g2 K^{-1}(x1, x2) (K: the L x L chain block), and the decay length fitted from ln|M| against xi. """ from pathlib import Path import numpy as np import scipy.sparse as sp from scipy.linalg import eigh, solve, eigvalsh from scipy.sparse.linalg import eigsh OUT = Path(__file__).parent / "output" OUT.mkdir(exist_ok=True) LOG = [] def out(s=""): print(s) LOG.append(s) # ---------------- (a) second order of the shared-empty-place field ---------------- L, J, x1 = 400, 1.0, 100 # chain of L places (indices 1..L), place 0 = empty out(f"(a) Shared empty place, second order: open chain of L = {L} places, J = {J}, x1 = {x1}, g1 = g2 = 1, N_f = 1") def shared_contract(Delta, n=L + 1): C = np.zeros((n, n)) for x in range(1, n): C[x, x] = Delta if x + 1 < n: C[x, x + 1] = C[x + 1, x] = -J return C for ratio in [2.02, 2.1, 2.5, 3.0]: Delta = ratio * J n = L + 1 C = shared_contract(Delta) e, v = eigh(C) phi, eps = v[:, 0], e[0] P0 = np.outer(phi, phi) A = C - eps * np.eye(n) + (e[1] - e[0]) * P0 # invertible; equals C - eps on phi-perp def R_apply(u): """Reduced resolvent of C at eps applied to u.""" return solve(A, u - phi * (phi @ u), assume_a="sym") def Phi_apply(x, u): """Coupled operator Phi_x = |x> 1e-12 * abs(form[0]) rel = np.max(np.abs(M2[ok] / form[ok] - 1)) slope = np.polyfit(ds[ok], np.log(np.abs(M2[ok])), 1)[0] out(f"Delta/J = {ratio}: phi = |empty> (overlap {abs(phi[0]):.15f}), gap = {e[1]-e[0]:.4f}; " f"xi = 1/arccosh(Delta/2J) = {xi:.5f}; fitted decay length = {-1/slope:.5f}; " f"max |M2/(4.14) - 1| over d = 1..{ds[ok][-1]} = {rel:.1e}") out(" d : " + " ".join(f"{d:>9d}" for d in [1, 5, 10, 20, 40, 80])) out(" M2 : " + " ".join(f"{M2[d-1]: .2e}" for d in [1, 5, 10, 20, 40, 80])) out(" d^6 |M2| : " + " ".join(f"{d**6*abs(M2[d-1]):.2e}" for d in [1, 5, 10, 20, 40, 80])) # ---------------- (b) the shared-empty-place field beyond second order ---------------- out("") out("(b) Shared empty place, exact M(d) = eps0(g,g) - eps0(g,0) - eps0(0,g) + eps0(0,0), N_f = 1, Delta/J = 3, " f"L = {L}, x1 = {x1}") out(" limit = g1^2 g2^2 (A11 B22 + A22 B11) (fourth order, d -> infinity), A = K^{-1}, B = K^{-2}") Delta = 3.0 * J C = shared_contract(Delta) K = C[1:, 1:] Ki = np.linalg.inv(K) K2 = Ki @ Ki def eps0(srcs): h = C.copy() for x, g in srcs: h[x, 0] += g h[0, x] += g return eigvalsh(h, subset_by_index=[0, 0])[0] E00 = eps0([]) out(f" eps0(0,0) = {E00:.1e}") out(" g d: 1 10 40 160 limit(d=160) M(160)/limit - 1") for g in [0.05, 0.1, 0.2]: E1 = eps0([(x1, g)]) Ms = {} for d in [1, 10, 40, 160]: Ms[d] = eps0([(x1, g), (x1 + d, g)]) - E1 - eps0([(x1 + d, g)]) + E00 a, b = x1 - 1, x1 + 160 - 1 lim = g**4 * (Ki[a, a] * K2[b, b] + Ki[b, b] * K2[a, a]) out(f" {g:4.2f} " + " ".join(f"{Ms[d]: .3e}" for d in [1, 10, 40, 160]) + f" {lim: .4e} {Ms[160]/lim - 1: .2e}") # ---------------- (c) the cell field ---------------- Lc = 18 out("") out(f"(c) Cell field: {Lc} cells (places {{0,1}}^{Lc}, dimension {2**Lc}), J = 1, x1 = 2 (cells 0..{Lc-1}), " "g1 = g2 = g; exact M by sparse diagonalization (eigsh)") def cell_contract(Delta): N = 1 << Lc s = np.arange(N) pop = sum((s >> y) & 1 for y in range(Lc)) rows, cols, vals = [s], [s], [Delta * pop.astype(float)] for y in range(Lc - 1): m = ((s >> y) & 1) != ((s >> (y + 1)) & 1) rows.append(s[m]) cols.append(s[m] ^ ((1 << y) | (1 << (y + 1)))) vals.append(-J * np.ones(m.sum())) return sp.csr_matrix((np.concatenate(vals), (np.concatenate(rows), np.concatenate(cols))), shape=(N, N)) def sigx(x): N = 1 << Lc s = np.arange(N) return sp.csr_matrix((np.ones(N), (s, s ^ (1 << x))), shape=(N, N)) SX = [sigx(x) for x in range(Lc)] v0 = np.zeros(1 << Lc) v0[0] = 1.0 def ground(H): return eigsh(H, k=1, which="SA", tol=0, v0=v0 + 1e-3)[0][0] xa = 2 ds = np.arange(1, 15) for ratio, g in [(3.0, 0.1), (3.0, 0.2), (2.5, 0.1)]: Delta = ratio * J Cf = cell_contract(Delta) Kc = Delta * np.eye(Lc) - J * (np.eye(Lc, k=1) + np.eye(Lc, k=-1)) Kci = np.linalg.inv(Kc) xi = 1 / np.arccosh(Delta / (2 * J)) assert np.abs(Cf @ v0).max() == 0.0 # the vacuum is an exact eigenvector with cost 0, the lowest (Step 7); E0 = 0.0 # eigsh does not resolve this isolated one-dimensional sector, so it is set E1 = ground(Cf + g * SX[xa]) Ms, M2s = [], [] for d in ds: xb = xa + d Ms.append(ground(Cf + g * SX[xa] + g * SX[xb]) - E1 - ground(Cf + g * SX[xb]) + E0) M2s.append(-2 * g * g * Kci[xa, xb]) Ms, M2s = np.array(Ms), np.array(M2s) fit = (ds >= 2) & (ds <= 12) xi_M = -1 / np.polyfit(ds[fit], np.log(np.abs(Ms[fit])), 1)[0] xi_M2 = -1 / np.polyfit(ds[fit], np.log(np.abs(M2s[fit])), 1)[0] ratio_r = Ms / M2s out(f"Delta/J = {ratio}, g = {g} (|g1|+|g2| = {2*g} < (Delta-2J)/2 = {(Delta-2*J)/2}): E0 = 0 (vacuum), xi = {xi:.5f}; " f"fitted decay length on 2 <= d <= 12: of M {xi_M:.5f}, of M2 {xi_M2:.5f}") out(" d : " + " ".join(f"{d:>9d}" for d in ds)) out(" M : " + " ".join(f"{m: .2e}" for m in Ms)) out(" M/M2 : " + " ".join(f"{r:9.5f}" for r in ratio_r)) out(f" M/M2 over 1 <= d <= 14: min {ratio_r.min():.5f}, max {ratio_r.max():.5f}; " f"over 4 <= d <= 14: spread {ratio_r[3:].max() - ratio_r[3:].min():.1e}") (OUT / "control_decay.txt").write_text("\n".join(LOG) + "\n", encoding="utf-8")