06-ilang-source / 04-profile
04profileverified
Determines for a class of hopping fields how the mutual cost of two held sources depends on their distance and on the dimension of the field graph, and what ties field places to body positions.
check_control_decay.py7 KB · raw
"""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><x| - J * sum_x (|x><x+1| + h.c.), coupled operator Phi_x = |x><empty| + h.c.
Second-order mutual cost from the single-object contract (generic form (4.3)):
M2 = -2 g1 g2 Re <phi|Phi_x1 R Phi_x2|phi>,
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><empty| + |empty><x| applied to u."""
w = np.zeros(n)
w[x] += u[0]
w[0] += u[x]
return w
theta = np.arccosh(Delta / (2 * J))
xi = 1 / theta
ds = np.arange(1, 161)
M2 = np.array([-2 * (phi @ Phi_apply(x1, R_apply(Phi_apply(x1 + d, phi)))) for d in ds])
form = -np.exp(-ds * theta) / (J * np.sinh(theta))
ok = np.abs(form) > 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")
check_exact_ring.py9 KB · raw
"""Check (ii) of 06-ilang-source/04-profile.
One field object on a ring of L places, two sources at distance d.
Exact mutual cost M = E(g1,g2) - E(g1,0) - E(0,g2) + E(0,0) from
(a) exact diagonalization of h = C1 + g1 P_x1 + g2 P_x2,
(b) the pole-cleared eigenvalue equation: F(eps) = delta * [(4.10) LHS - RHS], delta = eps + 2 t_f,
with the ring Green function (4.11); eps0(g) = smallest zero of F on (-inf, eps1), eps1 = -2t cos(2pi/L)
(F = det(eps - h) / prod_{k != 0}(eps - eps_k): no pole at eps = -2t, so q = 0 is included;
F is evaluated in the variable delta, q = 2 arcsin(sqrt(delta/4t)), q = i theta for delta < 0),
(c) for g1 = g2 the reduced equation (4.12), t sin q sin(qL/2) - g cos(qd/2) cos(q(L-d)/2) = 0,
against the second-order form (4.5) and the third-order term (4.8).
Plus unequal and mixed-sign sources, including roots at q = 0 (delta = 0).
Plus the exact finite-L bracket of (4.8) on the ring (Step 5).
Plus the mean over placements (4.17): Mbar = (1/L) sum_y M(0, y), y = 0..L-1 (y = 0: both sources on
place 0), so that E(0, d) - mean_y E(0, y) = M(d) - Mbar; (4.17) states that Mbar = O(g1^2 g2^2).
"""
from pathlib import Path
import numpy as np
from scipy.linalg import eigh
from scipy.optimize import brentq
OUT = Path(__file__).parent / "output"
OUT.mkdir(exist_ok=True)
LOG = []
def out(s=""):
print(s)
LOG.append(s)
t = 1.0
def ring(L):
return -t * (np.eye(L, k=1) + np.eye(L, k=-1) + np.eye(L, k=L - 1) + np.eye(L, k=-(L - 1)))
def E_exact(C1, sources):
h = C1.copy()
for x, g in sources:
h[x, x] += g
return eigh(h, eigvals_only=True, subset_by_index=[0, 0])[0]
def ell(d, L):
"""Step 3: Lambda^+(x, x+d) on the ring."""
return (L * L - 1) / (12 * L) - d * (L - d) / (2 * L)
def qofd(delta):
"""q with delta = eps + 2t = 4 t sin^2(q/2); q = i theta for delta < 0."""
return 2 * np.arcsin(np.sqrt(complex(delta / (4 * t))))
def F_cleared(delta, g1, g2, d, L):
"""delta * [(1 - g1 G0)(1 - g2 G0) - g1 g2 Gd^2] with the ring Green function (4.11), written as
delta (1 - (g1+g2) G0) + g1 g2 (G0 - Gd) delta (G0 + Gd); at delta = 0 the limit
-(g1 + g2)/L - 2 g1 g2 (l0 - ld)/(L t)."""
if delta == 0.0:
return -(g1 + g2) / L - 2 * g1 * g2 * (ell(0, L) - ell(d, L)) / (L * t)
q = qofd(delta)
den = t * np.sin(q) * np.sin(q * L / 2)
G0 = np.cos(q * L / 2) / (2 * den)
Gm = -np.sin(q * (L - d) / 2) * np.sin(q * d / 2) / den # G0 - Gd
Gp = np.cos(q * (L - d) / 2) * np.cos(q * d / 2) / den # G0 + Gd
return float(np.real(delta * (1 - (g1 + g2) * G0) + g1 * g2 * Gm * (delta * Gp)))
def F412(delta, g, d, L):
"""(4.12) for g1 = g2 = g: t sin q sin(qL/2) - g cos(qd/2) cos(q(L-d)/2);
one source (d = None): 2 t sin q sin(qL/2) - g cos(qL/2)."""
q = qofd(delta)
if d is None:
return float(np.real(2 * t * np.sin(q) * np.sin(q * L / 2) - g * np.cos(q * L / 2)))
return float(np.real(t * np.sin(q) * np.sin(q * L / 2) - g * np.cos(q * d / 2) * np.cos(q * (L - d) / 2)))
def smallest_zero(f, L, gneg):
"""eps = -2t + delta for the smallest zero delta of f on (-inf, delta1), delta1 = eps1 + 2t:
f < 0 below every eigenvalue; scan upward in delta."""
dmax = 2 * t * (1 - np.cos(2 * np.pi / L))
neg = -np.geomspace(gneg + 1e-3, 1e-18, 3000)
pos = np.concatenate([np.geomspace(1e-18, dmax / 2, 3000),
dmax - np.geomspace(dmax / 2, dmax * 1e-12, 3000)])
ds = np.concatenate([neg, [0.0], pos])
vals = [f(x) for x in ds]
assert vals[0] < 0
for a, b, fa, fb in zip(ds[:-1], ds[1:], vals[:-1], vals[1:]):
if fb == 0.0:
return -2 * t + b
if np.sign(fa) != np.sign(fb):
return -2 * t + brentq(f, a, b, xtol=1e-20, rtol=1e-15)
raise RuntimeError("no root below eps1")
def eps_cleared(g1, g2, d, L):
return smallest_zero(lambda x: F_cleared(x, g1, g2, d, L), L, max(0, -g1) + max(0, -g2))
out("Exact mutual cost on a ring (one field object, t_f = 1), gamma = g L / t_f, w = d(L-d)/L^2")
out("rel2 = (M_exact - M2)/M2, pred3 = M3/M2 from (4.8), rel3 = (M_exact - M2 - M3)/M2,")
out("pred3_L = -(gamma/5)(1-10w+20w^2)/(1-6w) (leading in 1/L); dev10 / dev12 = max deviation of the")
out("smallest zero of the pole-cleared (4.10) / of (4.12) from eigh (one- and two-source energies)")
for L in [64, 256]:
C1 = ring(L)
P0 = np.full((L, L), 1.0 / L) # projector on the lowest eigenvector of C1
Lp = np.linalg.inv(C1 + 2 * t * np.eye(L) + t * P0) - P0 / t # reduced resolvent R = L^+/t_f
Lp2 = Lp @ Lp
out(f"-- L = {L}")
out(" gamma d w M_exact M2 rel2 pred3 pred3_L rel3 dev10 dev12")
for gam in [0.03, 0.1, 0.3, 1.0, 3.0, -0.3, -1.0]:
g = gam * t / L
E00 = -2 * t
E10 = E_exact(C1, [(0, g)])
c10 = eps_cleared(g, 0.0, 0, L)
r10 = smallest_zero(lambda x: F412(x, g, None, L), L, max(0, -g))
for d in [1, L // 8, L // 4, L // 2]:
E11 = E_exact(C1, [(0, g), (d, g)])
c11 = eps_cleared(g, g, d, L)
r11 = smallest_zero(lambda x: F412(x, g, d, L), L, max(0, -2 * g))
Mx = E11 - 2 * E10 + E00
w = d * (L - d) / L**2
M2 = g * g / (t * L * L) * (d * (L - d) - (L * L - 1) / 6)
l0, ld, q0, qd = Lp[0, 0], Lp[0, d], Lp2[0, 0], Lp2[0, d] # R and R^2 entries (t_f = 1)
M3 = g * g * (2 * g) / L * (2 * l0 * ld + ld**2 - (q0 + 2 * qd) / L)
pL = -(gam / 5) * (1 - 10 * w + 20 * w * w) / (1 - 6 * w)
dev10 = max(abs(c11 - E11), abs(c10 - E10))
dev12 = max(abs(r11 - E11), abs(r10 - E10))
out(f" {gam:5.2f} {d:5d} {w:6.4f} {Mx: .5e} {M2: .5e} {(Mx-M2)/M2: .3e} {M3/M2: .3e} {pL: .3e} "
f"{(Mx-M2-M3)/M2: .2e} {dev10:.1e} {dev12:.1e}")
out("")
out("Unequal and mixed-sign sources (pole-cleared (4.10) against eigh):")
cases = [(64, 8, 0.3 / 64, 0.6 / 64), (64, 8, 1.0 / 64, 2.0 / 64), (64, 8, -1.0 / 64, -2.0 / 64),
(64, 8, 1.0 / 64, -0.5 / 64), (64, 20, 3.0 / 64, -1.0 / 64), (4, 2, 1.0, -0.5), (4, 1, 1.0, -0.5)]
# a case with the root exactly at q = 0: 1/g1 + 1/g2 = -2 (l0 - ld)/t (Step 6)
L, d, g1 = 64, 8, 1.0 / 64
cases.append((L, d, g1, 1.0 / (-2 * (ell(0, L) - ell(d, L)) / t - 1 / g1)))
for L, d, g1, g2 in cases:
C1 = ring(L)
ex = E_exact(C1, [(0, g1), (d, g2)])
cl = eps_cleared(g1, g2, d, L)
F0 = F_cleared(0.0, g1, g2, d, L)
Mx = ex - E_exact(C1, [(0, g1)]) - E_exact(C1, [(d, g2)]) + (-2 * t)
M2 = g1 * g2 / (t * L * L) * (d * (L - d) - (L * L - 1) / 6)
out(f" L={L:3d} d={d:3d} g1={g1: .5f} g2={g2: .5f}: E_eigh = {ex:.15f}, E_(4.10 cleared) = {cl:.15f}, "
f"F(q=0) = {F0: .2e}, (M_exact - M2)/M2 = {(Mx-M2)/M2: .3e}")
out("")
out("Exact finite-L bracket of (4.8) on the ring: B_matrix = 2 l0 ld + ld^2 - (l0^(2) + 2 ld^(2))/L against")
out("B_closed = (L^2/60)(1 - 10w + 20w^2) - 1/12 + w/3 + 1/(15 L^2); max_d |B_matrix - B_closed| / max_d |B_closed|")
for L in [3, 4, 5, 16, 64, 256]:
C1 = ring(L)
P0 = np.full((L, L), 1.0 / L)
Lp = np.linalg.inv(C1 + 2 * t * np.eye(L) + t * P0) - P0 / t
Lp2 = Lp @ Lp
ds = np.arange(L)
Bm = 2 * Lp[0, 0] * Lp[0, ds] + Lp[0, ds] ** 2 - (Lp2[0, 0] + 2 * Lp2[0, ds]) / L
w = ds * (L - ds) / L**2
Bc = L * L / 60 * (1 - 10 * w + 20 * w * w) - 1 / 12 + w / 3 + 1 / (15 * L * L)
out(f" L={L:4d}: {np.max(np.abs(Bm - Bc)) / np.max(np.abs(Bc)):.2e}")
out("")
out("Mean over placements (4.17), g1 = g2 = g, all y = 0..L-1 (y = 0: both sources on one place):")
out("Mbar = (1/L) sum_y M(0,y); inv = E(0,d) - mean_y E(0,y) = M(d) - Mbar (computed from total contracts only);")
out("res = max_d |inv - M2 - M3| / max_d |M2| (M2, M3 from (4.5), (4.8)); (4.17): Mbar has no term of order < g^4")
out(" L gamma Mbar*L^2/t_f Mbar*L^2/(t_f gamma^4) |Mbar|/max_d|M2| res")
for L in [64, 256]:
C1 = ring(L)
P0 = np.full((L, L), 1.0 / L)
Lp = np.linalg.inv(C1 + 2 * t * np.eye(L) + t * P0) - P0 / t
Lp2 = Lp @ Lp
ds = np.arange(L)
for gam in [0.1, 0.3, 1.0, -0.3, -1.0]:
g = gam * t / L
E00, E10 = -2 * t, E_exact(C1, [(0, g)])
Ms = np.array([E_exact(C1, [(0, g), (y, g)]) - 2 * E10 + E00 for y in ds])
Mbar = Ms.mean()
M2 = g * g / (t * L * L) * (ds * (L - ds) - (L * L - 1) / 6)
M3 = g * g * (2 * g) / L * (2 * Lp[0, 0] * Lp[0, ds] + Lp[0, ds] ** 2 - (Lp2[0, 0] + 2 * Lp2[0, ds]) / L)
inv = Ms - Mbar
res = np.max(np.abs(inv - M2 - M3)) / np.max(np.abs(M2))
out(f" {L:3d} {gam:6.2f} {Mbar*L*L/t: .4e} {Mbar*L*L/(t*gam**4): .4e} "
f"{abs(Mbar)/np.max(np.abs(M2)):.3e} {res:.2e}")
(OUT / "exact_ring.txt").write_text("\n".join(LOG) + "\n", encoding="utf-8")
check_remainder_bound.py5 KB · raw
"""Check of (4.18) and (4.19) of 06-ilang-source/04-profile (item 2: range of the second-order form).
One field object, t_f = 1, on a ring (L = 64), a square torus (16 x 16) and a cubic torus (8^3).
delta(g) := lowest eigenvalue of C1 + V minus eps0 = lowest eigenvalue of Lambda + V (exact diagonalization).
(a) residual of the pole-free secular equation (4.18),
|X| delta (c1 c2 - b^2 g1 g2) - (g1 c2 + g2 c1 - 2 b g1 g2), c_i = 1 + g_i s0, b = s_d,
s(delta) = [(Lambda - delta)^{-1} Q] entries,
at the exact delta, for couplings inside and far outside the regime of (4.19);
(b) the remainders r2 = max_x2 |M - M2| / calM and r3 = max_x2 |M - M2 - M3| / calM (x2 = x1 included),
calM = 2 |g1 g2| l0 / |X|, against the bounds 4 eta/(1 - eta) and 4 eta^2/(1 - eta), eta = 16 l0 max|g_i|.
x2 runs over one representative of every class of separations under the symmetry of the graph.
"""
from itertools import product
from pathlib import Path
import numpy as np
from scipy.linalg import eigh
OUT = Path(__file__).parent / "output"
OUT.mkdir(exist_ok=True)
LOG = []
def out(s=""):
print(s)
LOG.append(s)
def laplacian(L, D):
n = L**D
idx = np.arange(n).reshape((L,) * D)
A = np.zeros((n, n))
for ax in range(D):
nb = np.roll(idx, -1, axis=ax).ravel()
A[idx.ravel(), nb] += 1
A[nb, idx.ravel()] += 1
return 2 * D * np.eye(n) - A, idx
def lowest(Lam, V):
h = Lam.copy()
for x, g in V:
h[x, x] += g
return eigh(h, eigvals_only=True, subset_by_index=[0, 0])[0]
graphs = [("ring L=64", 64, 1), ("square torus L=16", 16, 2), ("cubic torus L=8", 8, 3)]
for name, L, D in graphs:
Lam, idx = laplacian(L, D)
n = L**D
w, U = eigh(Lam)
assert abs(w[0]) < 1e-10 and w[1] > 1e-6
Lp = (U[:, 1:] / w[1:]) @ U[:, 1:].T
Lp2 = (U[:, 1:] / w[1:] ** 2) @ U[:, 1:].T
l0, om1 = Lp[0, 0], w[1]
reps = [idx[tuple(c)] for c in product(range(L // 2 + 1), repeat=D) if list(c) == sorted(c)]
out(f"== {name}: |X| = {n}, l0 = Lambda^+(x,x) = {l0:.6f}, omega_1 = {om1:.6f}, l0 |X| omega_1 = {l0*n*om1:.3f}, "
f"{len(reps)} separations")
def S(delta, x2):
r = 1.0 / (w[1:] - delta)
return (U[0, 1:] * r) @ U[0, 1:], (U[0, 1:] * r) @ U[x2, 1:]
# (a) residual of (4.18) at the exact lowest eigenvalue
worst = 0.0
for g1, g2 in [(0.01, 0.02), (0.3, -0.2), (3.0, 5.0), (-4.0, 1.5), (-10.0, -10.0)]:
for x2 in reps[1:]:
dl = lowest(Lam, [(0, g1), (x2, g2)])
s0, sd = S(dl, x2)
c1, c2 = 1 + g1 * s0, 1 + g2 * s0
lhs, rhs = n * dl * (c1 * c2 - sd**2 * g1 * g2), g1 * c2 + g2 * c1 - 2 * sd * g1 * g2
scale = abs(n * dl) * (abs(c1 * c2) + abs(sd**2 * g1 * g2)) + abs(g1 * c2) + abs(g2 * c1) + abs(2 * sd * g1 * g2)
worst = max(worst, abs(lhs - rhs) / scale)
out(f" (4.18): max relative residual (|lhs - rhs| / sum of |terms|) at the exact lowest eigenvalue = {worst:.1e} "
f"(couplings from (0.01,0.02) to (-10,-10))")
# (b) remainders against the bound (4.19)
out(" eta g1 g2 r2 4eta/(1-eta) r3 4eta^2/(1-eta)")
for eta in [0.05, 0.25, 0.5, 0.9]:
gm = eta / (16 * l0)
for g1, g2 in [(gm, gm), (gm, -gm / 2), (-gm, -gm / 3)]:
d1 = lowest(Lam, [(0, g1)])
d2 = lowest(Lam, [(0, g2)])
calM = 2 * abs(g1 * g2) * l0 / n
r2 = r3 = 0.0
for x2 in reps:
V = [(0, g1 + g2)] if x2 == 0 else [(0, g1), (x2, g2)]
M = lowest(Lam, V) - d1 - d2
ld = Lp[0, x2]
M2 = -2 * g1 * g2 * ld / n
M3 = g1 * g2 * (g1 + g2) / n * (2 * l0 * ld + ld**2 - (Lp2[0, 0] + 2 * Lp2[0, x2]) / n)
r2 = max(r2, abs(M - M2) / calM)
r3 = max(r3, abs(M - M2 - M3) / calM)
out(f" {eta:4.2f} {g1: .3e} {g2: .3e} {r2:.2e} {4*eta/(1-eta):.2e} {r3:.2e} "
f"{4*eta**2/(1-eta):.2e} {'ok' if (r2 <= 4*eta/(1-eta) and r3 <= 4*eta**2/(1-eta)) else 'VIOLATED'}")
# the pole case of Step 6: ring L = 4, g = (1, -1/2), d = 2 -> delta = 0 is a root of (4.18)
Lam, idx = laplacian(4, 1)
w, U = eigh(Lam)
Lp = (U[:, 1:] / w[1:]) @ U[:, 1:].T
dl = lowest(Lam, [(0, 1.0), (2, -0.5)])
c1, c2, b = 1 + 1.0 * Lp[0, 0], 1 - 0.5 * Lp[0, 0], Lp[0, 2]
out("")
out(f"Ring L=4, g=(1,-1/2), d=2: exact delta = {dl:.2e}; right side of (4.18) at delta=0: "
f"{1.0*c2 - 0.5*c1 - 2*b*1.0*(-0.5):.2e}")
(OUT / "remainder_bound.txt").write_text("\n".join(LOG) + "\n", encoding="utf-8")
check_second_order.py6 KB · raw
"""Check (i) of 06-ilang-source/04-profile.
Second-order mutual cost M2 computed numerically from the single-object contract
C1 = -t_f * (adjacency of the graph) via (4.3):
M2 = -2 N_f g1 g2 phi(x1) phi(x2) R(x1, x2), R = reduced resolvent of C1 at the eigenvalue of phi,
compared with the ring closed form (4.5) and the leading torus forms (4.6), (4.7).
Also: lowest + highest state on even rings (sum over the two ends, (4.13)).
"""
from pathlib import Path
from math import lgamma, log, pi
import numpy as np
import scipy.sparse as sp
from scipy.sparse.linalg import cg, LinearOperator
OUT = Path(__file__).parent / "output"
OUT.mkdir(exist_ok=True)
LOG = []
def out(s=""):
print(s)
LOG.append(s)
def torus_contract(L, D, t=1.0):
"""C1 = -t * sum over nearest-neighbour edges of the L^D torus (L >= 3)."""
n = L**D
idx = np.arange(n).reshape((L,) * D)
r = np.concatenate([idx.ravel()] * D)
c = np.concatenate([np.roll(idx, -1, axis=ax).ravel() for ax in range(D)])
A = sp.coo_matrix((np.ones(len(r)), (r, c)), shape=(n, n)).tocsr()
return -t * (A + A.T), idx
def phi_and_R(C1, which):
"""Eigenvector phi (which = 0: lowest, -1: highest) and reduced resolvent R of C1 at its eigenvalue."""
e, v = np.linalg.eigh(C1)
j = which % len(e)
w = np.array([0.0 if k == j else 1.0 / (e[k] - e[j]) for k in range(len(e))])
return v[:, j], (v * w) @ v.T
def m2_dense(C1, x1, x2, which, Nf=1, g1=1.0, g2=1.0, pr=None):
"""(4.3): M2 = -2 N_f g1 g2 phi(x1) phi(x2) R(x1, x2)."""
phi, R = pr if pr is not None else phi_and_R(C1, which)
return -2 * Nf * g1 * g2 * phi[x1] * phi[x2] * R[x1, x2]
def resolvent_column_cg(C1, phi, eps, x2, t=1.0):
"""R(., x2) from (C1 - eps + t |phi><phi|) y = e_x2 - phi phi(x2) by conjugate gradients."""
n = C1.shape[0]
op = LinearOperator((n, n), matvec=lambda u: C1 @ u - eps * u + t * phi * (phi @ u))
b = -phi * phi[x2]
b[x2] += 1.0
y, info = cg(op, b, rtol=1e-13, atol=0.0, maxiter=20000)
assert info == 0
return y
# ---------------- ring: closed form (4.5) ----------------
out("== Ring: numerical (4.3) vs closed form (4.5); M2 * t_f * L^2 / (N_f g1 g2) vs d(L-d) - (L^2-1)/6 ==")
for L in [3, 4, 16, 101, 400]:
for t in [1.0, 2.0]:
C1 = -t * (np.eye(L, k=1) + np.eye(L, k=-1) + np.eye(L, k=L - 1) + np.eye(L, k=-(L - 1)))
if L == 3:
C1 = -t * (np.ones((3, 3)) - np.eye(3))
dev = 0.0
pr = phi_and_R(C1, 0)
for d in range(1, L):
M2 = m2_dense(C1, 0, d, 0, pr=pr)
closed = (d * (L - d) - (L * L - 1) / 6) / (t * L * L)
dev = max(dev, abs(M2 - closed) / ((L * L - 1) / (6 * t * L * L)))
out(f"L={L:4d} t_f={t}: max_d |M2 - (4.5)| / max|(4.5)| = {dev:.2e}")
out("Ring L=3, d=1, t_f=1: M2 = %.12f, 2/27 = %.12f" % (m2_dense(-(np.ones((3, 3)) - np.eye(3)), 0, 1, 0), 2 / 27))
out("")
out("== Even rings: lowest + highest state (sum over the two ends, (4.13)); odd ring: top degeneracy ==")
for L in [4, 16, 64, 15]:
C1 = -(np.eye(L, k=1) + np.eye(L, k=-1) + np.eye(L, k=L - 1) + np.eye(L, k=-(L - 1)))
if L % 2 == 0:
p0, p1 = phi_and_R(C1, 0), phi_and_R(C1, -1)
s = max(abs(m2_dense(C1, 0, d, 0, pr=p0) + m2_dense(C1, 0, d, -1, pr=p1)) for d in range(1, L))
out(f"L={L}: max_d |M2_lowest + M2_highest| = {s:.2e}")
else:
e = np.linalg.eigvalsh(C1)
out(f"L={L} (odd): two highest eigenvalues of C1 = {e[-1]:.12f}, {e[-2]:.12f} (degenerate)")
# ---------------- square and cubic tori ----------------
kappa2 = 2 * lgamma(0.25) - log(2) - 0.5 * log(pi) # ln(Gamma(1/4)^2 / (2 sqrt(pi)))
xi3 = -2.837297 / (4 * pi)
out("")
out(f"kappa_2 = {kappa2:.6f}, xi_3 = {xi3:.6f}")
def torus_table(D, Ls, rs, dirs, lead):
for L in Ls:
C1, idx = torus_contract(L, D)
n = L**D
phi = np.full(n, n**-0.5)
eps = phi @ (C1 @ phi)
res = np.linalg.norm(C1 @ phi - eps * phi)
x2 = 0
y = resolvent_column_cg(C1, phi, eps, x2)
out(f"-- D={D}, L={L}: |X|={n}, eps0={eps:.6f}, ||C1 phi - eps0 phi|| = {res:.1e}, phi>0 (Perron)")
out(" d vector |d| m_num m_lead m_num-m_lead")
maxdev = 0.0
for dv in dirs:
for r in rs:
d = tuple(r * a for a in dv)
if max(d) > L // 4:
continue
x1 = idx[d]
M2 = -2 * phi[x1] * phi[x2] * y[x1] # N_f = g1 = g2 = t_f = 1
m_num = M2 * n # M2 * t_f / (n g1 g2)
dist = float(np.linalg.norm(d))
m_lead = lead(dist, L)
out(f" {str(d):15s} {dist:6.2f} {m_num: .6f} {m_lead: .6f} {m_num - m_lead: .2e}")
if 4 <= dist <= L / 8:
maxdev = max(maxdev, abs(m_num - m_lead))
out(f" max |m_num - m_lead| over 4 <= |d| <= L/8: {maxdev:.2e}")
out("")
out("== Square torus: m = M2 t_f/(n g1 g2) vs leading form (4.6): m_lead = -(1/pi)[ln(L/|d|) - kappa_2] ==")
torus_table(2, [64, 128, 256], [2, 4, 8, 16, 32], [(1, 0), (1, 1), (2, 1)],
lambda r, L: -(log(L / r) - kappa2) / pi)
out("")
out("== Cubic torus: m vs leading form (4.7): m_lead = -2[1/(4 pi |d|) + xi_3/L] ==")
torus_table(3, [32, 48, 64], [2, 3, 4, 6, 8], [(1, 0, 0), (1, 1, 0), (1, 1, 1)],
lambda r, L: -2 * (1 / (4 * pi * r) + xi3 / L))
(OUT / "second_order.txt").write_text("\n".join(LOG) + "\n", encoding="utf-8")
output/control_decay.txt5 KB · raw
(a) Shared empty place, second order: open chain of L = 400 places, J = 1.0, x1 = 100, g1 = g2 = 1, N_f = 1
Delta/J = 2.02: phi = |empty> (overlap 1.000000000000000), gap = 0.0201; xi = 1/arccosh(Delta/2J) = 7.07695; fitted decay length = 7.07695; max |M2/(4.14) - 1| over d = 1..160 = 5.7e-13
d : 1 5 10 20 40 80
M2 : -6.12e+00 -3.48e+00 -1.72e+00 -4.18e-01 -2.48e-02 -8.69e-05
d^6 |M2| : 6.12e+00 5.44e+04 1.72e+06 2.67e+07 1.01e+08 2.28e+07
Delta/J = 2.1: phi = |empty> (overlap 1.000000000000000), gap = 0.1001; xi = 1/arccosh(Delta/2J) = 3.17536; fitted decay length = 3.17536; max |M2/(4.14) - 1| over d = 1..88 = 2.2e-15
d : 1 5 10 20 40 80
M2 : -2.28e+00 -6.47e-01 -1.34e-01 -5.74e-03 -1.06e-05 -3.57e-11
d^6 |M2| : 2.28e+00 1.01e+04 1.34e+05 3.68e+05 4.33e+04 9.37e+00
Delta/J = 2.5: phi = |empty> (overlap 1.000000000000000), gap = 0.5001; xi = 1/arccosh(Delta/2J) = 1.44270; fitted decay length = 1.44270; max |M2/(4.14) - 1| over d = 1..40 = 2.4e-15
d : 1 5 10 20 40 80
M2 : -6.67e-01 -4.17e-02 -1.30e-03 -1.27e-06 -1.21e-12 -1.10e-24
d^6 |M2| : 6.67e-01 6.51e+02 1.30e+03 8.14e+01 4.97e-03 2.89e-13
Delta/J = 3.0: phi = |empty> (overlap 1.000000000000000), gap = 1.0001; xi = 1/arccosh(Delta/2J) = 1.03904; fitted decay length = 1.03904; max |M2/(4.14) - 1| over d = 1..29 = 2.9e-15
d : 1 5 10 20 40 80
M2 : -3.42e-01 -7.27e-03 -5.91e-05 -3.91e-09 -1.71e-17 -3.26e-34
d^6 |M2| : 3.42e-01 1.14e+02 5.91e+01 2.50e-01 7.00e-08 8.55e-23
(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, L = 400, x1 = 100
limit = g1^2 g2^2 (A11 B22 + A22 B11) (fourth order, d -> infinity), A = K^{-1}, B = K^{-2}
eps0(0,0) = 0.0e+00
g d: 1 10 40 160 limit(d=160) M(160)/limit - 1
0.05 -8.487e-04 1.347e-06 1.493e-06 1.493e-06 1.5000e-06 -4.45e-03
0.10 -3.332e-03 2.302e-05 2.358e-05 2.358e-05 2.4000e-05 -1.76e-02
0.20 -1.239e-02 3.564e-04 3.583e-04 3.583e-04 3.8400e-04 -6.68e-02
(c) Cell field: 18 cells (places {0,1}^18, dimension 262144), J = 1, x1 = 2 (cells 0..17), g1 = g2 = g; exact M by sparse diagonalization (eigsh)
Delta/J = 3.0, g = 0.1 (|g1|+|g2| = 0.2 < (Delta-2J)/2 = 0.5): E0 = 0 (vacuum), xi = 1.03904; fitted decay length on 2 <= d <= 12: of M 1.03921, of M2 1.03902
d : 1 2 3 4 5 6 7 8 9 10 11 12 13 14
M : -3.36e-03 -1.28e-03 -4.91e-04 -1.88e-04 -7.18e-05 -2.74e-05 -1.05e-05 -4.00e-06 -1.53e-06 -5.84e-07 -2.23e-07 -8.51e-08 -3.24e-08 -1.22e-08
M/M2 : 0.98566 0.98767 0.98889 0.98953 0.98983 0.98997 0.99003 0.99005 0.99007 0.99007 0.99007 0.99008 0.99012 0.99037
M/M2 over 1 <= d <= 14: min 0.98566, max 0.99037; over 4 <= d <= 14: spread 8.4e-04
Delta/J = 3.0, g = 0.2 (|g1|+|g2| = 0.4 < (Delta-2J)/2 = 0.5): E0 = 0 (vacuum), xi = 1.03904; fitted decay length on 2 <= d <= 12: of M 1.03974, of M2 1.03902
d : 1 2 3 4 5 6 7 8 9 10 11 12 13 14
M : -1.29e-02 -4.96e-03 -1.90e-03 -7.28e-04 -2.79e-04 -1.06e-04 -4.07e-05 -1.55e-05 -5.93e-06 -2.27e-06 -8.66e-07 -3.31e-07 -1.26e-07 -4.73e-08
M/M2 : 0.94531 0.95254 0.95706 0.95941 0.96055 0.96106 0.96129 0.96140 0.96144 0.96146 0.96147 0.96149 0.96164 0.96257
M/M2 over 1 <= d <= 14: min 0.94531, max 0.96257; over 4 <= d <= 14: spread 3.2e-03
Delta/J = 2.5, g = 0.1 (|g1|+|g2| = 0.2 < (Delta-2J)/2 = 0.25): E0 = 0 (vacuum), xi = 1.44270; fitted decay length on 2 <= d <= 12: of M 1.44387, of M2 1.44223
d : 1 2 3 4 5 6 7 8 9 10 11 12 13 14
M : -6.31e-03 -3.17e-03 -1.59e-03 -7.96e-04 -3.99e-04 -2.00e-04 -9.99e-05 -4.99e-05 -2.50e-05 -1.25e-05 -6.24e-06 -3.11e-06 -1.54e-06 -7.33e-07
M/M2 : 0.96196 0.96543 0.96864 0.97092 0.97238 0.97327 0.97379 0.97409 0.97426 0.97436 0.97444 0.97458 0.97503 0.97655
M/M2 over 1 <= d <= 14: min 0.96196, max 0.97655; over 4 <= d <= 14: spread 5.6e-03
output/exact_ring.txt9 KB · raw
Exact mutual cost on a ring (one field object, t_f = 1), gamma = g L / t_f, w = d(L-d)/L^2
rel2 = (M_exact - M2)/M2, pred3 = M3/M2 from (4.8), rel3 = (M_exact - M2 - M3)/M2,
pred3_L = -(gamma/5)(1-10w+20w^2)/(1-6w) (leading in 1/L); dev10 / dev12 = max deviation of the
smallest zero of the pole-cleared (4.10) / of (4.12) from eigh (one- and two-source energies)
-- L = 64
gamma d w M_exact M2 rel2 pred3 pred3_L rel3 dev10 dev12
0.03 1 0.0154 -3.30464e-08 -3.32326e-08 -5.602e-03 -5.619e-03 -5.625e-03 1.63e-05 2.2e-16 2.2e-16
0.03 8 0.1094 -1.25477e-08 -1.25796e-08 -2.535e-03 -2.530e-03 -2.540e-03 -5.34e-06 2.2e-16 2.2e-16
0.03 16 0.1875 4.54887e-09 4.58658e-09 -8.221e-03 -8.249e-03 -8.250e-03 2.78e-05 2.2e-16 2.2e-16
0.03 32 0.2500 1.82646e-08 1.83195e-08 -2.994e-03 -2.999e-03 -3.000e-03 5.02e-06 4.4e-16 4.4e-16
0.10 1 0.0154 -3.62402e-07 -3.69251e-07 -1.855e-02 -1.873e-02 -1.875e-02 1.81e-04 4.4e-16 4.4e-16
0.10 8 0.1094 -1.38586e-07 -1.39773e-07 -8.489e-03 -8.432e-03 -8.466e-03 -5.74e-05 4.4e-16 4.4e-16
0.10 16 0.1875 4.95764e-08 5.09620e-08 -2.719e-02 -2.750e-02 -2.750e-02 3.06e-04 4.4e-16 4.4e-16
0.10 32 0.2500 2.01527e-07 2.03550e-07 -9.940e-03 -9.995e-03 -1.000e-02 5.53e-05 4.4e-16 4.4e-16
0.30 1 0.0154 -3.14194e-06 -3.32326e-06 -5.456e-02 -5.619e-02 -5.625e-02 1.63e-03 2.2e-16 2.2e-16
0.30 8 0.1094 -1.22553e-06 -1.25796e-06 -2.578e-02 -2.530e-02 -2.540e-02 -4.80e-04 2.2e-16 2.2e-16
0.30 16 0.1875 4.22078e-07 4.58658e-07 -7.975e-02 -8.249e-02 -8.250e-02 2.73e-03 2.2e-16 2.2e-16
0.30 32 0.2500 1.77793e-06 1.83195e-06 -2.949e-02 -2.999e-02 -3.000e-02 4.96e-04 4.4e-16 4.4e-16
1.00 1 0.0154 -3.06651e-05 -3.69251e-05 -1.695e-01 -1.873e-01 -1.875e-01 1.78e-02 4.4e-16 4.4e-16
1.00 8 0.1094 -1.27422e-05 -1.39773e-05 -8.836e-02 -8.432e-02 -8.466e-02 -4.04e-03 4.4e-16 4.4e-16
1.00 16 0.1875 3.84439e-06 5.09620e-06 -2.456e-01 -2.750e-01 -2.750e-01 2.93e-02 4.4e-16 4.4e-16
1.00 32 0.2500 1.84310e-05 2.03550e-05 -9.452e-02 -9.995e-02 -1.000e-01 5.43e-03 4.4e-16 4.4e-16
3.00 1 0.0154 -1.94134e-04 -3.32326e-04 -4.158e-01 -5.619e-01 -5.625e-01 1.46e-01 0.0e+00 0.0e+00
3.00 8 0.1094 -9.24888e-05 -1.25796e-04 -2.648e-01 -2.530e-01 -2.540e-01 -1.18e-02 2.2e-16 2.2e-16
3.00 16 0.1875 1.88729e-05 4.58658e-05 -5.885e-01 -8.249e-01 -8.250e-01 2.36e-01 2.2e-16 2.2e-16
3.00 32 0.2500 1.36755e-04 1.83195e-04 -2.535e-01 -2.999e-01 -3.000e-01 4.64e-02 2.2e-16 2.2e-16
-0.30 1 0.0154 -3.51540e-06 -3.32326e-06 5.782e-02 5.619e-02 5.625e-02 1.63e-03 4.4e-16 4.4e-16
-0.30 8 0.1094 -1.28903e-06 -1.25796e-06 2.470e-02 2.530e-02 2.540e-02 -5.93e-04 0.0e+00 0.0e+00
-0.30 16 0.1875 4.97778e-07 4.58658e-07 8.529e-02 8.249e-02 8.250e-02 2.81e-03 0.0e+00 0.0e+00
-0.30 32 0.2500 1.88780e-06 1.83195e-06 3.049e-02 2.999e-02 3.000e-02 5.01e-04 4.4e-16 4.4e-16
-1.00 1 0.0154 -4.45014e-05 -3.69251e-05 2.052e-01 1.873e-01 1.875e-01 1.79e-02 4.4e-16 4.4e-16
-1.00 8 0.1094 -1.50412e-05 -1.39773e-05 7.612e-02 8.432e-02 8.466e-02 -8.20e-03 4.4e-16 4.4e-16
-1.00 16 0.1875 6.66095e-06 5.09620e-06 3.070e-01 2.750e-01 2.750e-01 3.21e-02 4.4e-16 4.4e-16
-1.00 32 0.2500 2.25039e-05 2.03550e-05 1.056e-01 9.995e-02 1.000e-01 5.62e-03 4.4e-16 4.4e-16
-- L = 256
gamma d w M_exact M2 rel2 pred3 pred3_L rel3 dev10 dev12
0.03 1 0.0039 -2.22219e-09 -2.23535e-09 -5.886e-03 -5.906e-03 -5.906e-03 1.94e-05 4.4e-16 4.4e-16
0.03 32 0.1094 -7.84745e-10 -7.86746e-10 -2.543e-03 -2.539e-03 -2.540e-03 -4.26e-06 4.4e-16 4.4e-16
0.03 64 0.1875 2.83784e-10 2.86137e-10 -8.225e-03 -8.250e-03 -8.250e-03 2.46e-05 4.4e-16 4.4e-16
0.03 128 0.2500 1.14102e-09 1.14444e-09 -2.996e-03 -3.000e-03 -3.000e-03 4.24e-06 4.4e-16 4.4e-16
0.10 1 0.0039 -2.43535e-08 -2.48372e-08 -1.948e-02 -1.969e-02 -1.969e-02 2.11e-04 2.2e-16 2.2e-16
0.10 32 0.1094 -8.66715e-09 -8.74163e-09 -8.520e-03 -8.464e-03 -8.466e-03 -5.65e-05 2.2e-16 2.2e-16
0.10 64 0.1875 3.09285e-09 3.17930e-09 -2.719e-02 -2.750e-02 -2.750e-02 3.07e-04 2.2e-16 2.2e-16
0.10 128 0.2500 1.25896e-08 1.27160e-08 -9.944e-03 -1.000e-02 -1.000e-02 5.54e-05 2.2e-16 2.2e-16
0.30 1 0.0039 -2.10755e-07 -2.23535e-07 -5.717e-02 -5.906e-02 -5.906e-02 1.89e-03 2.2e-16 2.2e-16
0.30 32 0.1094 -7.66398e-08 -7.86746e-08 -2.586e-02 -2.539e-02 -2.540e-02 -4.73e-04 2.2e-16 2.2e-16
0.30 64 0.1875 2.63314e-08 2.86137e-08 -7.976e-02 -8.250e-02 -8.250e-02 2.73e-03 2.2e-16 2.2e-16
0.30 128 0.2500 1.11068e-07 1.14444e-07 -2.950e-02 -3.000e-02 -3.000e-02 4.97e-04 2.2e-16 2.2e-16
1.00 1 0.0039 -2.04534e-06 -2.48372e-06 -1.765e-01 -1.969e-01 -1.969e-01 2.04e-02 2.2e-16 2.2e-16
1.00 32 0.1094 -7.96704e-07 -8.74163e-07 -8.861e-02 -8.464e-02 -8.466e-02 -3.97e-03 4.4e-16 4.4e-16
1.00 64 0.1875 2.39830e-07 3.17930e-07 -2.457e-01 -2.750e-01 -2.750e-01 2.93e-02 2.2e-16 2.2e-16
1.00 128 0.2500 1.15137e-06 1.27160e-06 -9.456e-02 -1.000e-01 -1.000e-01 5.44e-03 2.2e-16 2.2e-16
3.00 1 0.0039 -1.28013e-05 -2.23535e-05 -4.273e-01 -5.906e-01 -5.906e-01 1.63e-01 4.4e-16 4.4e-16
3.00 32 0.1094 -5.78102e-06 -7.86746e-06 -2.652e-01 -2.539e-01 -2.540e-01 -1.13e-02 2.2e-16 2.2e-16
3.00 64 0.1875 1.17749e-06 2.86137e-06 -5.885e-01 -8.250e-01 -8.250e-01 2.37e-01 0.0e+00 0.0e+00
3.00 128 0.2500 8.54277e-06 1.14444e-05 -2.535e-01 -3.000e-01 -3.000e-01 4.64e-02 2.2e-16 2.2e-16
-0.30 1 0.0039 -2.37166e-07 -2.23535e-07 6.098e-02 5.906e-02 5.906e-02 1.92e-03 0.0e+00 0.0e+00
-0.30 32 0.1094 -8.06263e-08 -7.86746e-08 2.481e-02 2.539e-02 2.540e-02 -5.85e-04 4.4e-16 4.4e-16
-0.30 64 0.1875 3.10548e-08 2.86137e-08 8.531e-02 8.250e-02 8.250e-02 2.81e-03 0.0e+00 0.0e+00
-0.30 128 0.2500 1.17935e-07 1.14444e-07 3.050e-02 3.000e-02 3.000e-02 5.02e-04 4.4e-16 4.4e-16
-1.00 1 0.0039 -3.02602e-06 -2.48372e-06 2.183e-01 1.969e-01 1.969e-01 2.15e-02 4.4e-16 4.4e-16
-1.00 32 0.1094 -9.41059e-07 -8.74163e-07 7.653e-02 8.464e-02 8.466e-02 -8.11e-03 4.4e-16 4.4e-16
-1.00 64 0.1875 4.15579e-07 3.17930e-07 3.071e-01 2.750e-01 2.750e-01 3.21e-02 4.4e-16 4.4e-16
-1.00 128 0.2500 1.40593e-06 1.27160e-06 1.056e-01 1.000e-01 1.000e-01 5.64e-03 4.4e-16 4.4e-16
Unequal and mixed-sign sources (pole-cleared (4.10) against eigh):
L= 64 d= 8 g1= 0.00469 g2= 0.00937: E_eigh = -1.999791523936965, E_(4.10 cleared) = -1.999791523936965, F(q=0) = -2.25e-04, (M_exact - M2)/M2 = -3.896e-02
L= 64 d= 8 g1= 0.01562 g2= 0.03125: E_eigh = -1.999382359521463, E_(4.10 cleared) = -1.999382359521462, F(q=0) = -7.86e-04, (M_exact - M2)/M2 = -1.335e-01
L= 64 d= 8 g1=-0.01562 g2=-0.03125: E_eigh = -2.000878297773575, E_(4.10 cleared) = -2.000878297773575, F(q=0) = 6.79e-04, (M_exact - M2)/M2 = 1.042e-01
L= 64 d= 8 g1= 0.01562 g2=-0.00781: E_eigh = -1.999895398537231, E_(4.10 cleared) = -1.999895398537231, F(q=0) = -1.09e-04, (M_exact - M2)/M2 = -2.154e-02
L= 64 d= 20 g1= 0.04688 g2=-0.01562: E_eigh = -1.999715219474732, E_(4.10 cleared) = -1.999715219474732, F(q=0) = -3.31e-04, (M_exact - M2)/M2 = -1.396e-01
L= 4 d= 2 g1= 1.00000 g2=-0.50000: E_eigh = -2.000000000000000, E_(4.10 cleared) = -2.000000000000000, F(q=0) = 0.00e+00, (M_exact - M2)/M2 = -1.484e-01
L= 4 d= 1 g1= 1.00000 g2=-0.50000: E_eigh = -1.972818212999035, E_(4.10 cleared) = -1.972818212999035, F(q=0) = -3.12e-02, (M_exact - M2)/M2 = -1.849e-01
L= 64 d= 8 g1= 0.01562 g2=-0.01408: E_eigh = -2.000000000000000, E_(4.10 cleared) = -2.000000000000000, F(q=0) = 6.78e-21, (M_exact - M2)/M2 = -4.655e-03
Exact finite-L bracket of (4.8) on the ring: B_matrix = 2 l0 ld + ld^2 - (l0^(2) + 2 ld^(2))/L against
B_closed = (L^2/60)(1 - 10w + 20w^2) - 1/12 + w/3 + 1/(15 L^2); max_d |B_matrix - B_closed| / max_d |B_closed|
L= 3: 7.49e-16
L= 4: 1.48e-16
L= 5: 9.91e-16
L= 16: 8.49e-16
L= 64: 3.75e-15
L= 256: 2.26e-13
Mean over placements (4.17), g1 = g2 = g, all y = 0..L-1 (y = 0: both sources on one place):
Mbar = (1/L) sum_y M(0,y); inv = E(0,d) - mean_y E(0,y) = M(d) - Mbar (computed from total contracts only);
res = max_d |inv - M2 - M3| / max_d |M2| (M2, M3 from (4.5), (4.8)); (4.17): Mbar has no term of order < g^4
L gamma Mbar*L^2/t_f Mbar*L^2/(t_f gamma^4) |Mbar|/max_d|M2| res
64 0.10 1.2774e-08 1.2774e-04 7.666e-06 2.28e-04
64 0.30 9.6820e-07 1.1953e-04 6.456e-05 2.03e-03
64 1.00 9.4761e-05 9.4761e-05 5.687e-04 2.17e-02
64 -0.30 1.1827e-06 1.4601e-04 7.886e-05 2.09e-03
64 -1.00 1.8410e-04 1.8410e-04 1.105e-03 2.37e-02
256 0.10 1.2767e-08 1.2767e-04 7.660e-06 2.29e-04
256 0.30 9.6940e-07 1.1968e-04 6.463e-05 2.04e-03
256 1.00 9.4897e-05 9.4897e-05 5.694e-04 2.18e-02
256 -0.30 1.1840e-06 1.4617e-04 7.893e-05 2.10e-03
256 -1.00 1.8427e-04 1.8427e-04 1.106e-03 2.38e-02
output/remainder_bound.txt4 KB · raw
== ring L=64: |X| = 64, l0 = Lambda^+(x,x) = 5.332031, omega_1 = 0.009631, l0 |X| omega_1 = 3.286, 33 separations
(4.18): max relative residual (|lhs - rhs| / sum of |terms|) at the exact lowest eigenvalue = 7.6e-13 (couplings from (0.01,0.02) to (-10,-10))
eta g1 g2 r2 4eta/(1-eta) r3 4eta^2/(1-eta)
0.05 5.861e-04 5.861e-04 7.46e-03 2.11e-01 3.11e-05 1.05e-02 ok
0.05 5.861e-04 -2.930e-04 1.87e-03 2.11e-01 4.45e-06 1.05e-02 ok
0.05 -5.861e-04 -1.954e-04 5.01e-03 2.11e-01 1.44e-05 1.05e-02 ok
0.25 2.930e-03 2.930e-03 3.67e-02 1.33e+00 7.73e-04 3.33e-01 ok
0.25 2.930e-03 -1.465e-03 9.26e-03 1.33e+00 1.11e-04 3.33e-01 ok
0.25 -2.930e-03 -9.768e-04 2.53e-02 1.33e+00 3.60e-04 3.33e-01 ok
0.50 5.861e-03 5.861e-03 7.19e-02 4.00e+00 3.07e-03 2.00e+00 ok
0.50 5.861e-03 -2.930e-03 1.83e-02 4.00e+00 4.42e-04 2.00e+00 ok
0.50 -5.861e-03 -1.954e-03 5.14e-02 4.00e+00 1.45e-03 2.00e+00 ok
0.90 1.055e-02 1.055e-02 1.25e-01 3.60e+01 9.81e-03 3.24e+01 ok
0.90 1.055e-02 -5.275e-03 3.23e-02 3.60e+01 1.43e-03 3.24e+01 ok
0.90 -1.055e-02 -3.516e-03 9.46e-02 3.60e+01 4.71e-03 3.24e+01 ok
== square torus L=16: |X| = 256, l0 = Lambda^+(x,x) = 0.489924, omega_1 = 0.152241, l0 |X| omega_1 = 19.094, 45 separations
(4.18): max relative residual (|lhs - rhs| / sum of |terms|) at the exact lowest eigenvalue = 5.2e-12 (couplings from (0.01,0.02) to (-10,-10))
eta g1 g2 r2 4eta/(1-eta) r3 4eta^2/(1-eta)
0.05 6.379e-03 6.379e-03 9.15e-03 2.11e-01 6.45e-05 1.05e-02 ok
0.05 6.379e-03 -3.189e-03 2.29e-03 2.11e-01 9.25e-06 1.05e-02 ok
0.05 -6.379e-03 -2.126e-03 6.17e-03 2.11e-01 3.00e-05 1.05e-02 ok
0.25 3.189e-02 3.189e-02 4.45e-02 1.33e+00 1.57e-03 3.33e-01 ok
0.25 3.189e-02 -1.595e-02 1.13e-02 1.33e+00 2.29e-04 3.33e-01 ok
0.25 -3.189e-02 -1.063e-02 3.15e-02 1.33e+00 7.64e-04 3.33e-01 ok
0.50 6.379e-02 6.379e-02 8.60e-02 4.00e+00 6.10e-03 2.00e+00 ok
0.50 6.379e-02 -3.189e-02 2.21e-02 4.00e+00 9.02e-04 2.00e+00 ok
0.50 -6.379e-02 -2.126e-02 6.46e-02 4.00e+00 3.13e-03 2.00e+00 ok
0.90 1.148e-01 1.148e-01 1.47e-01 3.60e+01 1.89e-02 3.24e+01 ok
0.90 1.148e-01 -5.741e-02 3.86e-02 3.60e+01 2.86e-03 3.24e+01 ok
0.90 -1.148e-01 -3.827e-02 1.21e-01 3.60e+01 1.05e-02 3.24e+01 ok
== cubic torus L=8: |X| = 512, l0 = Lambda^+(x,x) = 0.224606, omega_1 = 0.585786, l0 |X| omega_1 = 67.364, 35 separations
(4.18): max relative residual (|lhs - rhs| / sum of |terms|) at the exact lowest eigenvalue = 2.7e-11 (couplings from (0.01,0.02) to (-10,-10))
eta g1 g2 r2 4eta/(1-eta) r3 4eta^2/(1-eta)
0.05 1.391e-02 1.391e-02 9.27e-03 2.11e-01 6.72e-05 1.05e-02 ok
0.05 1.391e-02 -6.957e-03 2.33e-03 2.11e-01 9.61e-06 1.05e-02 ok
0.05 -1.391e-02 -4.638e-03 6.26e-03 2.11e-01 3.13e-05 1.05e-02 ok
0.25 6.957e-02 6.957e-02 4.51e-02 1.33e+00 1.64e-03 3.33e-01 ok
0.25 6.957e-02 -3.478e-02 1.14e-02 1.33e+00 2.38e-04 3.33e-01 ok
0.25 -6.957e-02 -2.319e-02 3.19e-02 1.33e+00 7.96e-04 3.33e-01 ok
0.50 1.391e-01 1.391e-01 8.71e-02 4.00e+00 6.34e-03 2.00e+00 ok
0.50 1.391e-01 -6.957e-02 2.24e-02 4.00e+00 9.39e-04 2.00e+00 ok
0.50 -1.391e-01 -4.638e-02 6.55e-02 4.00e+00 3.26e-03 2.00e+00 ok
0.90 2.504e-01 2.504e-01 1.49e-01 3.60e+01 1.96e-02 3.24e+01 ok
0.90 2.504e-01 -1.252e-01 3.91e-02 3.60e+01 2.98e-03 3.24e+01 ok
0.90 -2.504e-01 -8.348e-02 1.23e-01 3.60e+01 1.10e-02 3.24e+01 ok
Ring L=4, g=(1,-1/2), d=2: exact delta = -1.74e-16; right side of (4.18) at delta=0: 5.55e-16
output/second_order.txt7 KB · raw
== Ring: numerical (4.3) vs closed form (4.5); M2 * t_f * L^2 / (N_f g1 g2) vs d(L-d) - (L^2-1)/6 ==
L= 3 t_f=1.0: max_d |M2 - (4.5)| / max|(4.5)| = 5.62e-16
L= 3 t_f=2.0: max_d |M2 - (4.5)| / max|(4.5)| = 5.62e-16
L= 4 t_f=1.0: max_d |M2 - (4.5)| / max|(4.5)| = 1.78e-16
L= 4 t_f=2.0: max_d |M2 - (4.5)| / max|(4.5)| = 1.78e-16
L= 16 t_f=1.0: max_d |M2 - (4.5)| / max|(4.5)| = 1.50e-15
L= 16 t_f=2.0: max_d |M2 - (4.5)| / max|(4.5)| = 1.50e-15
L= 101 t_f=1.0: max_d |M2 - (4.5)| / max|(4.5)| = 6.56e-14
L= 101 t_f=2.0: max_d |M2 - (4.5)| / max|(4.5)| = 6.56e-14
L= 400 t_f=1.0: max_d |M2 - (4.5)| / max|(4.5)| = 1.47e-12
L= 400 t_f=2.0: max_d |M2 - (4.5)| / max|(4.5)| = 1.47e-12
Ring L=3, d=1, t_f=1: M2 = 0.074074074074, 2/27 = 0.074074074074
== Even rings: lowest + highest state (sum over the two ends, (4.13)); odd ring: top degeneracy ==
L=4: max_d |M2_lowest + M2_highest| = 1.39e-17
L=16: max_d |M2_lowest + M2_highest| = 1.30e-15
L=64: max_d |M2_lowest + M2_highest| = 2.78e-15
L=15 (odd): two highest eigenvalues of C1 = 1.956295201468, 1.956295201468 (degenerate)
kappa_2 = 1.310533, xi_3 = -0.225785
== Square torus: m = M2 t_f/(n g1 g2) vs leading form (4.6): m_lead = -(1/pi)[ln(L/|d|) - kappa_2] ==
-- D=2, L=64: |X|=4096, eps0=-4.000000, ||C1 phi - eps0 phi|| = 0.0e+00, phi>0 (Perron)
d vector |d| m_num m_lead m_num-m_lead
(2, 0) 2.00 -0.695059 -0.686022 -9.04e-03
(4, 0) 4.00 -0.469300 -0.465387 -3.91e-03
(8, 0) 8.00 -0.253040 -0.244751 -8.29e-03
(16, 0) 16.00 -0.056436 -0.024116 -3.23e-02
(2, 2) 2.83 -0.573480 -0.575705 2.22e-03
(4, 4) 5.66 -0.358125 -0.355069 -3.06e-03
(8, 8) 11.31 -0.149590 -0.134433 -1.52e-02
(16, 16) 22.63 0.027659 0.086202 -5.85e-02
(4, 2) 4.47 -0.431846 -0.429872 -1.97e-03
(8, 4) 8.94 -0.218862 -0.209237 -9.63e-03
(16, 8) 17.89 -0.027187 0.011399 -3.86e-02
max |m_num - m_lead| over 4 <= |d| <= L/8: 8.29e-03
-- D=2, L=128: |X|=16384, eps0=-4.000000, ||C1 phi - eps0 phi|| = 0.0e+00, phi>0 (Perron)
d vector |d| m_num m_lead m_num-m_lead
(2, 0) 2.00 -0.915339 -0.906658 -8.68e-03
(4, 0) 4.00 -0.688478 -0.686022 -2.46e-03
(8, 0) 8.00 -0.467770 -0.465387 -2.38e-03
(16, 0) 16.00 -0.252726 -0.244751 -7.97e-03
(32, 0) 32.00 -0.056370 -0.024116 -3.23e-02
(2, 2) 2.83 -0.793395 -0.796340 2.95e-03
(4, 4) 5.66 -0.575856 -0.575705 -1.52e-04
(8, 8) 11.31 -0.358750 -0.355069 -3.68e-03
(16, 16) 22.63 -0.149758 -0.134433 -1.53e-02
(32, 32) 45.25 0.027599 0.086202 -5.86e-02
(4, 2) 4.47 -0.650663 -0.650508 -1.55e-04
(8, 4) 8.94 -0.432210 -0.429872 -2.34e-03
(16, 8) 17.89 -0.218948 -0.209237 -9.71e-03
(32, 16) 35.78 -0.027219 0.011399 -3.86e-02
max |m_num - m_lead| over 4 <= |d| <= L/8: 7.97e-03
-- D=2, L=256: |X|=65536, eps0=-4.000000, ||C1 phi - eps0 phi|| = 0.0e+00, phi>0 (Perron)
d vector |d| m_num m_lead m_num-m_lead
(2, 0) 2.00 -1.135885 -1.127294 -8.59e-03
(4, 0) 4.00 -0.908750 -0.906658 -2.09e-03
(8, 0) 8.00 -0.686940 -0.686022 -9.18e-04
(16, 0) 16.00 -0.467447 -0.465387 -2.06e-03
(32, 0) 32.00 -0.252650 -0.244751 -7.90e-03
(2, 2) 2.83 -1.013850 -1.016976 3.13e-03
(4, 4) 5.66 -0.795763 -0.796340 5.77e-04
(8, 8) 11.31 -0.576473 -0.575705 -7.68e-04
(16, 16) 22.63 -0.358907 -0.355069 -3.84e-03
(32, 32) 45.25 -0.149800 -0.134433 -1.54e-02
(4, 2) 4.47 -0.870843 -0.871144 3.00e-04
(8, 4) 8.94 -0.651019 -0.650508 -5.11e-04
(16, 8) 17.89 -0.432288 -0.429872 -2.42e-03
(32, 16) 35.78 -0.218969 -0.209237 -9.73e-03
(64, 32) 71.55 -0.027227 0.011399 -3.86e-02
max |m_num - m_lead| over 4 <= |d| <= L/8: 7.90e-03
== Cubic torus: m vs leading form (4.7): m_lead = -2[1/(4 pi |d|) + xi_3/L] ==
-- D=3, L=32: |X|=32768, eps0=-6.000000, ||C1 phi - eps0 phi|| = 0.0e+00, phi>0 (Perron)
d vector |d| m_num m_lead m_num-m_lead
(2, 0, 0) 2.00 -0.071711 -0.065466 -6.24e-03
(3, 0, 0) 3.00 -0.041074 -0.038940 -2.13e-03
(4, 0, 0) 4.00 -0.026635 -0.025677 -9.58e-04
(6, 0, 0) 6.00 -0.013006 -0.012414 -5.91e-04
(8, 0, 0) 8.00 -0.006579 -0.005783 -7.96e-04
(2, 2, 0) 2.83 -0.042082 -0.042158 7.58e-05
(3, 3, 0) 4.24 -0.023502 -0.023402 -1.01e-04
(4, 4, 0) 5.66 -0.014303 -0.014023 -2.79e-04
(6, 6, 0) 8.49 -0.005343 -0.004645 -6.98e-04
(8, 8, 0) 11.31 -0.001184 0.000044 -1.23e-03
(2, 2, 2) 3.46 -0.031314 -0.031833 5.18e-04
(3, 3, 3) 5.20 -0.016597 -0.016518 -7.90e-05
(4, 4, 4) 6.93 -0.009248 -0.008860 -3.88e-04
(6, 6, 6) 10.39 -0.002170 -0.001203 -9.67e-04
(8, 8, 8) 13.86 0.001009 0.002626 -1.62e-03
max |m_num - m_lead| over 4 <= |d| <= L/8: 9.58e-04
-- D=3, L=48: |X|=110592, eps0=-6.000000, ||C1 phi - eps0 phi|| = 4.3e-14, phi>0 (Perron)
d vector |d| m_num m_lead m_num-m_lead
(2, 0, 0) 2.00 -0.076384 -0.070170 -6.21e-03
(3, 0, 0) 3.00 -0.045711 -0.043644 -2.07e-03
(4, 0, 0) 4.00 -0.031220 -0.030381 -8.38e-04
(6, 0, 0) 6.00 -0.017433 -0.017118 -3.15e-04
(8, 0, 0) 8.00 -0.010770 -0.010487 -2.83e-04
(2, 2, 0) 2.83 -0.046727 -0.046862 1.35e-04
(3, 3, 0) 4.24 -0.028077 -0.028105 2.88e-05
(4, 4, 0) 5.66 -0.018780 -0.018727 -5.23e-05
(6, 6, 0) 8.49 -0.009549 -0.009349 -2.00e-04
(8, 8, 0) 11.31 -0.005031 -0.004660 -3.71e-04
(2, 2, 2) 3.46 -0.035932 -0.036536 6.05e-04
(3, 3, 3) 5.20 -0.021112 -0.021222 1.10e-04
(4, 4, 4) 6.93 -0.013626 -0.013564 -6.21e-05
(6, 6, 6) 10.39 -0.006195 -0.005907 -2.88e-04
(8, 8, 8) 13.86 -0.002602 -0.002078 -5.23e-04
max |m_num - m_lead| over 4 <= |d| <= L/8: 8.38e-04
-- D=3, L=64: |X|=262144, eps0=-6.000000, ||C1 phi - eps0 phi|| = 0.0e+00, phi>0 (Perron)
d vector |d| m_num m_lead m_num-m_lead
(2, 0, 0) 2.00 -0.078728 -0.072522 -6.21e-03
(3, 0, 0) 3.00 -0.048046 -0.045996 -2.05e-03
(4, 0, 0) 4.00 -0.033543 -0.032733 -8.10e-04
(6, 0, 0) 6.00 -0.019720 -0.019470 -2.50e-04
(8, 0, 0) 8.00 -0.013004 -0.012839 -1.65e-04
(2, 2, 0) 2.83 -0.049065 -0.049214 1.49e-04
(3, 3, 0) 4.24 -0.030397 -0.030457 6.04e-05
(4, 4, 0) 5.66 -0.021076 -0.021079 3.44e-06
(6, 6, 0) 8.49 -0.011777 -0.011701 -7.63e-05
(8, 8, 0) 11.31 -0.007166 -0.007012 -1.55e-04
(2, 2, 2) 3.46 -0.038262 -0.038888 6.26e-04
(3, 3, 3) 5.20 -0.023417 -0.023574 1.56e-04
(4, 4, 4) 6.93 -0.015897 -0.015916 1.97e-05
(6, 6, 6) 10.39 -0.008369 -0.008259 -1.10e-04
(8, 8, 8) 13.86 -0.004654 -0.004430 -2.23e-04
max |m_num - m_lead| over 4 <= |d| <= L/8: 8.10e-04