"""Checks for the path blocks: the anchored cube, the coercivity bound 224 on group-valued fields, the coarse curvature and reconstruction bounds, and the thin-path ratio. Parts: (A) anchored cube: the Laplacian of {0,1}^3 with the vertex 0 deleted has characteristic polynomial (x-3)(x^2-6x+2)(x-2)^2(x-4)^2 (exact), so its eigenvalues are 3, 3 +- sqrt 7, and 2, 4 twice; the least, 3 - sqrt 7, exceeds 1/3. (B) SU(2) fields in the tree gauge on the periodic lattice of side 4: D_v^2 <= 224 Q, Q_1(V) <= 4 Q, Q(R(v)) = 4 Q_1(v), Q <= 8 Q_1 + 48 D_v^2 (floating point, random fields at five scales). (C) thin paths: Q / Q_k = s sin^2(pi/(2s)) for s = 2, 4, ..., 64, from the s/2 changed links with six plaquettes each and the six coarse plaquettes of the one nontrivial coarse link. (D) slice windings: three covariant continuous maps of eight slice links (u_1 and alternating products u_a u_b^-1 u_c ...) satisfy F(g u g^-1) = g F(u) g^-1 and F(g u) = g F(u), and the degrees of their coordinate loops add up to one; the normalised quaternion mean of the eight links, which treats them alike, takes values 2 apart at fields arbitrarily close to the field with u_t = -I on four links and I on the others. (E) the same for SU(3) on its maximal torus: for covariant continuous maps of six slice links the integer matrices W_t on H_1 of the torus add up to the identity. (F) the cut-off average V# = V^p exp(chi/16 sum_x log(V^p^-1 P_x)) on SU(2), rho = 1, c = pi/3: covariance, differential at I equal to the block average, the bounds 896, 56 c^2, 448 + 896 c^2, 8 + 2688 c^2 on random fields at five scales, slice windings 1 at the path and 0 elsewhere, and the coarse link e^{X/8} of a single small link e^X. (G) torus-valued fields (random and smooth, under random gauge transformations) on the lattice of side 8, two steps with the tree gauge restored between them: the cut-off is inactive, the coarse plaquette angles are the means of the sixteen 2x2 fluxes, the L1 norm drops by at least 4, the L2 norm does not increase, the maximum grows by at most 4, and Q_2 <= 1.001 Q when eta_0 A <= (rho/12)^2. (H) SU(2) fields of small noncommuting curvature (random and smooth, under random SU(2) gauge transformations) on the lattice of side 8, two steps: the cut-off is inactive, the coarse plaquette logarithm is the mean of the sixteen 2x2 square logarithms up to 3.2 r_P^2 with r_P <= 28.26 eta, and the growth bounds of sqrt E, A and eta of the nonabelian step hold. (I) large fields under the SU(2) Wilson measure: on the torus of side 8 a plaquette lies in exactly four unit cubes, a cube has 24 plaquettes and 80 neighbours, Z_{y,mu} has 24 cubes and a cube lies in Z_{y,mu} for at most eight coarse links; the constants C_0, 4 8 4 24 144 <= 4.5e5 and 8 b_0 = 11/(3 pi^2) < 4/3; the link-angle bound 1/(24 pi^3 beta^{3/2}) and the one-link tail p = v_* e^{-3 beta} by quadrature, with the cap probability decreasing in the distance of the centre and increasing in kappa; the staple sum -4I of the frustrated exterior on a lattice and the conditional probability of p_0 being bad, even for rho = 2.5; and on random SU(2) fields with localised large links, chi = 1 at every coarse link whose 24 cubes are good for rho = 1/12. Run: python3 path_blocks.py. Needs numpy.""" import itertools, math import numpy as np ok = [] def check(name, cond): ok.append(bool(cond)); print(('ok ' if cond else 'FAIL'), name) # (A) V = [a for a in itertools.product(range(2), repeat=3) if any(a)] L0 = np.zeros((7, 7), dtype=int) for i, a in enumerate(V): L0[i, i] = 3 for j, b in enumerate(V): if sum(x != y for x, y in zip(a, b)) == 1: L0[i, j] = -1 ev = np.sort(np.linalg.eigvalsh(L0)) want = np.sort([3, 3 + math.sqrt(7), 3 - math.sqrt(7), 2, 2, 4, 4]) check('(A) spectrum of the grounded cube Laplacian is 3, 3 +- sqrt7, 2, 2, 4, 4', np.allclose(ev, want)) # exact: det(L0 - 3) in integers equals the value of (t^3-7t)(t^2-6t+8)^2 at t=0... use # the integer characteristic polynomial via Faddeev-LeVerrier def charpoly(M): n = len(M); M = [[int(x) for x in r] for r in M] from fractions import Fraction as F c = [F(1)]; Mk = [[F(0)] * n for _ in range(n)] I = [[F(int(i == j)) for j in range(n)] for i in range(n)] A = [[F(x) for x in r] for r in M] for k in range(1, n + 1): Mk = [[sum(A[i][l] * Mk[l][j] for l in range(n)) + c[-1] * I[i][j] for j in range(n)] for i in range(n)] AM = [[sum(A[i][l] * Mk[l][j] for l in range(n)) for j in range(n)] for i in range(n)] c.append(-sum(AM[i][i] for i in range(n)) / k) return [int(x) for x in c] # expected (x-3)(x^2-6x+2)(x-2)^2(x-4)^2 def mul(p, q): r = [0] * (len(p) + len(q) - 1) for i, a in enumerate(p): for j, b in enumerate(q): r[i + j] += a * b return r exp = [1] for f in ([1, -3], [1, -6, 2], [1, -2], [1, -2], [1, -4], [1, -4]): exp = mul(exp, f) check('(A) exact characteristic polynomial (x-3)(x^2-6x+2)(x-2)^2(x-4)^2', charpoly(L0) == exp) check('(A) 3 - sqrt7 > 1/3, i.e. 64 > 63 after squaring 8/3 > sqrt7', 64 > 9 * 7) # (B) rng = np.random.default_rng(7) S = 4; d = 4; E = np.eye(d, dtype=int) sig = np.array([[[0, 1], [1, 0]], [[0, -1j], [1j, 0]], [[1, 0], [0, -1]]]) def su2(a): n = np.linalg.norm(a) if n == 0: return np.eye(2, dtype=complex) return math.cos(n) * np.eye(2) + math.sin(n) / n * 1j * np.tensordot(a, sig, 1) def tree(o, mu): return o[mu] == 0 and all(o[n] == 0 for n in range(mu + 1, d)) sites = list(itertools.product(range(S), repeat=d)) sh = lambda x, mu, k=1: tuple((np.array(x) + k * E[mu]) % S) def hs2(A): return float(np.sum(np.abs(A) ** 2)) def energies(U, Sz=S, step=1): Q = 0.0 pts = [x for x in itertools.product(range(0, Sz, step), repeat=d)] for x in pts: for mu in range(d): for nu in range(mu + 1, d): P = U[(x, mu)] @ U[(tuple((np.array(x) + step * E[mu]) % Sz), nu)] @ \ U[(tuple((np.array(x) + step * E[nu]) % Sz), mu)].conj().T @ U[(x, nu)].conj().T Q += hs2(P - np.eye(2)) return Q worst = [0, 0, 0, 0] for sc in [0.05, 0.3, 1.0, 2.0, 6.0]: for rep in range(2): U = {} for x in sites: o = [c % 2 for c in x] for mu in range(d): U[(x, mu)] = np.eye(2, dtype=complex) if tree(o, mu) else su2(sc * rng.normal(size=3)) Q = energies(U) v = {(y, mu): U[(sh(y, mu), mu)] for y in itertools.product(range(0, S, 2), repeat=d) for mu in range(d)} R, D = {}, 0.0 for x in sites: o = [c % 2 for c in x]; y = tuple(c // 2 * 2 for c in x) for mu in range(d): R[(x, mu)] = v[(y, mu)] if o[mu] == 1 else np.eye(2, dtype=complex) D += hs2(U[(x, mu)] - R[(x, mu)]) # coarse plaquettes: holonomy of 2x2 squares at even sites (coarse torus side 2 counts each twice per plane) Vc = {} for y in itertools.product(range(0, S, 2), repeat=d): for mu in range(d): Vc[(y, mu)] = U[(y, mu)] @ U[(sh(y, mu), mu)] Q1 = energies(Vc, S, 2) QR = energies(R) worst[0] = max(worst[0], D / Q); worst[1] = max(worst[1], Q1 / Q) worst[2] = max(worst[2], abs(QR - 4 * Q1) / max(1, Q1)); worst[3] = max(worst[3], Q / (8 * Q1 + 48 * D)) print(' max D^2/Q = %.4f, max Q_1/Q = %.4f, max |Q(R)-4Q_1| = %.2e, max Q/(8Q_1+48D^2) = %.4f' % tuple(worst)) check('(B) D_v^2 <= 224 Q on random SU(2) fields', worst[0] <= 224) check('(B) Q_1(V) <= 4 Q', worst[1] <= 4) check('(B) Q(R(v)) = 4 Q_1(v)', worst[2] < 1e-9) check('(B) Q <= 8 Q_1 + 48 D_v^2', worst[3] <= 1) # (C) T = np.diag([1j * math.pi / 2, -1j * math.pi / 2]) def expm(Hd): return np.diag(np.exp(np.diag(Hd))) good = True for k in range(1, 7): s = 2 ** k; n = 3 * s ch = expm(2 * T / s) # plaquettes containing a changed link (j,0,0,0), direction 0, j odd < s: six each, distinct Qf = 0.0 for j in range(1, s, 2): for nu in range(1, 4): # plaquette at (j,0,..) and at (j,0,..)-e_nu in the plane {0,nu} Qf += 2 * hs2(ch - np.eye(2)) Qk = 6 * hs2(expm(T) - np.eye(2)) r = Qf / Qk; want = s * math.sin(math.pi / (2 * s)) ** 2 good &= abs(r - want) < 1e-12 and r <= math.pi ** 2 / (4 * s) + 1e-15 print(' s = %3d: Q/Q_k = %.6f, s sin^2(pi/2s) = %.6f' % (s, r, want)) check('(C) thin-path ratio s sin^2(pi/(2s)) <= pi^2/(4s) for s = 2..64', good) # (D) X = np.array([[1j, 0], [0, -1j]]) def ex(th): return np.diag([np.exp(1j * th), np.exp(-1j * th)]) def angle(A): return math.atan2(A[0, 0].imag, A[0, 0].real) def degree(F, t, n=8, N=400): tot, prev = 0.0, None for k in range(N + 1): u = [np.eye(2, dtype=complex)] * n; u = list(u); u[t] = ex(2 * math.pi * k / N) a = angle(F(u)) if prev is not None: d = a - prev; d -= 2 * math.pi * round(d / (2 * math.pi)); tot += d prev = a return round(tot / (2 * math.pi)) inv = lambda A: A.conj().T maps = {'path map u_1': lambda u: u[0], 'u_1 u_2^-1 u_3': lambda u: u[0] @ inv(u[1]) @ u[2], 'u_4 u_7^-1 u_2 u_5^-1 u_8': lambda u: u[3] @ inv(u[6]) @ u[1] @ inv(u[4]) @ u[7]} rng = np.random.default_rng(7) def haar(): q = rng.normal(size=4); q /= np.linalg.norm(q) return np.array([[q[0] + 1j * q[1], q[2] + 1j * q[3]], [-q[2] + 1j * q[3], q[0] - 1j * q[1]]]) good = True for F in maps.values(): for _ in range(20): u = [haar() for _ in range(8)]; g = haar() good &= np.allclose(F([g @ a @ inv(g) for a in u]), g @ F(u) @ inv(g)) good &= np.allclose(F([g @ a for a in u]), g @ F(u)) check('(D) the three maps are covariant on slice fields', good) good = True for name, F in maps.items(): w = [degree(F, t) for t in range(8)] print(' %-30s w = %s' % (name, w)); good &= sum(w) == 1 check('(D) slice windings add up to one for covariant continuous maps', good) def qmean(u): q = sum(u) / len(u); r = math.sqrt(abs(np.linalg.det(q))) return q / r eps = 1e-6 Y = np.array([[0, 1], [-1, 0]], dtype=complex) def near(Xd): return [ex(math.pi) @ (np.eye(2) * math.cos(eps) + Xd * math.sin(eps))] * 4 + [np.eye(2, dtype=complex)] * 4 jump = np.linalg.norm(qmean(near(X)) - qmean(near(Y))) print(' quaternion mean at fields within %g of (-I)^4 I^4 in two directions differs by %.4f' % (eps, jump)) check('(D) the symmetric quaternion mean jumps by 2 near (-I,...,-I,I,...,I)', abs(jump - 2) < 1e-4) # (E) def tor3(a, b): return np.diag([np.exp(1j * a), np.exp(1j * b), np.exp(-1j * (a + b))]) def wind(vals): tot = 0.0 for v0, v1 in zip(vals, vals[1:]): d = math.atan2((v1 / v0).imag, (v1 / v0).real); tot += d return round(tot / (2 * math.pi)) def wmatrix(F, t, n=6, N=400): W = np.zeros((2, 2), dtype=int) for j in range(2): d0, d1 = [], [] for k in range(N + 1): u = [np.eye(3, dtype=complex) for _ in range(n)] th = 2 * math.pi * k / N u[t] = tor3(th, 0) if j == 0 else tor3(0, th) R = F(u); d0.append(R[0, 0]); d1.append(R[1, 1]) W[0, j], W[1, j] = wind(d0), wind(d1) return W def haar3(): z = rng.normal(size=(3, 3)) + 1j * rng.normal(size=(3, 3)) q, r = np.linalg.qr(z); q = q @ np.diag(np.diag(r) / abs(np.diag(r))) return q / np.linalg.det(q) ** (1 / 3) maps3 = {'u_1': lambda u: u[0], 'u_1 u_2^-1 u_3': lambda u: u[0] @ inv(u[1]) @ u[2], 'u_6 u_4^-1 u_5 u_2^-1 u_3': lambda u: u[5] @ inv(u[3]) @ u[4] @ inv(u[1]) @ u[2]} good = True for F in maps3.values(): for _ in range(10): u = [haar3() for _ in range(6)]; g = haar3() good &= np.allclose(F([g @ a @ inv(g) for a in u]), g @ F(u) @ inv(g)) good &= np.allclose(F([g @ a for a in u]), g @ F(u)) check('(E) the three SU(3) maps are covariant on slice fields', good) good = True for name, F in maps3.items(): S = sum(wmatrix(F, t) for t in range(6)) print(' %-26s sum of W_t = %s' % (name, S.tolist())); good &= (S == np.eye(2, dtype=int)).all() check('(E) for SU(3) the torus matrices W_t add up to the identity', good) # (F) Sd = 4; I2 = np.eye(2, dtype=complex); c = math.pi / 3; rho = 1.0 shf = lambda x, mu, k=1: tuple((np.array(x) + k * E[mu]) % Sd) def slog(g): # principal logarithm of SU(2) near I n = math.acos(max(-1.0, min(1.0, g[0, 0].real / 2 + g[1, 1].real / 2))) A = (g - g.conj().T) / 2 return A if n < 1e-12 else n / math.sin(n) * A def sexp(Z): n = math.sqrt(hs2(Z) / 2) return I2 * math.cos(n) + (Z if n < 1e-12 else math.sin(n) / n * Z) def phi(t): return 1.0 if t <= 0.5 else (0.0 if t >= 1 else 2 * (1 - t)) blk = list(itertools.product(range(2), repeat=d)) def sharp(U, y, mu): P = [U.get((tuple((np.array(y) + o) % Sd), mu), I2) @ U.get((shf(tuple((np.array(y) + o) % Sd), mu), mu), I2) for o in blk] v = [inv(P[0]) @ p for p in P] m = max(math.sqrt(hs2(w - I2)) for w in v); chi = phi(m / rho) M = I2 if chi == 0 else sexp(chi * sum(slog(w) for w in v) / 16) return P[0] @ M, P[0], M, P evens = list(itertools.product(range(0, Sd, 2), repeat=d)) def field(sc): U = {} for x in sites: o = [q % 2 for q in x] for mu in range(d): U[(x, mu)] = I2 if tree(o, mu) else su2(sc * rng.normal(size=3)) return U # covariance under the gauge transformations constant on blocks good = True for sc in [0.05, 0.3, 2.0]: U = field(sc); h = {y: haar() for y in evens} hb = lambda x: h[tuple(q // 2 * 2 for q in x)] Uh = {(x, mu): hb(x) @ U[(x, mu)] @ inv(hb(shf(x, mu))) for (x, mu) in U} for y in evens: for mu in range(d): good &= np.allclose(sharp(Uh, y, mu)[0], h[y] @ sharp(U, y, mu)[0] @ inv(h[shf(y, mu, 2)])) check('(F) the cut-off average is covariant under block gauge transformations', good) # linearisation: the block average t = 1e-6; err = 0.0 for rep in range(3): a = {} for x in sites: o = [q % 2 for q in x] for mu in range(d): a[(x, mu)] = np.zeros(3) if tree(o, mu) else rng.normal(size=3) U = {k: su2(t * a[k]) for k in a} for y in evens: for mu in range(d): avg = sum(a[(tuple((np.array(y) + o) % Sd), mu)] + a[(shf(tuple((np.array(y) + o) % Sd), mu), mu)] for o in blk) / 16 lin = (sharp(U, y, mu)[0] - I2) / t err = max(err, np.abs(lin - 1j * np.tensordot(avg, sig, 1)).max()) print(' max |(V#(e^{t phi}) - I)/t - V(phi)| at t = 1e-6: %.2e' % err) check('(F) the differential of the cut-off average at I is the block average', err < 1e-4) # bounds for all fields worst = [0, 0, 0, 0] for sc in [0.05, 0.3, 1.0, 2.0, 6.0]: for rep in range(2): U = field(sc); Q = energies(U, Sd) Vs, sP, sM = {}, 0.0, 0.0 for y in evens: for mu in range(d): V, P0, M, P = sharp(U, y, mu); Vs[(y, mu)] = V sP += sum(hs2(p - P0) for p in P); sM += hs2(M - I2) D = 0.0 for x in sites: o = [q % 2 for q in x]; y = tuple(q // 2 * 2 for q in x) for mu in range(d): D += hs2(U[(x, mu)] - (Vs[(y, mu)] if o[mu] == 1 else I2)) Q1 = energies(Vs, Sd, 2) worst = [max(worst[0], sP / Q), max(worst[1], sM / Q), max(worst[2], D / Q), max(worst[3], Q1 / Q)] print(' max sum d(P_x,P_2y)^2/Q = %.3f, sum d(M,I)^2/Q = %.4f, D^2/Q = %.3f, Q_1/Q = %.3f' % tuple(worst)) check('(F) sum of d(P_x, V^p)^2 <= 896 Q', worst[0] <= 896) check('(F) sum of d(V#, V^p)^2 <= 56 c^2 Q, c = pi/3', worst[1] <= 56 * c * c) check('(F) D_{V#}^2 <= (448 + 896 c^2) Q', worst[2] <= 448 + 896 * c * c) check('(F) Q_1(V#) <= (8 + 2688 c^2) Q', worst[3] <= 8 + 2688 * c * c) # slice windings: the identity at the path, zero elsewhere trans = [a for a in itertools.product(range(2), repeat=3)] def Fsharp(u): U = {((1,) + a, 0): u[i] for i, a in enumerate(trans)} return sharp(U, (0, 0, 0, 0), 0)[0] w = [degree(Fsharp, t) for t in range(8)] print(' slice windings of the cut-off average: %s' % w) check('(F) slice windings of the cut-off average are 1 at the path and 0 elsewhere', w == [1] + [0] * 7) # one link e^X on the straight path or at another transverse position gives e^{X/8} good = True for a in trans: for th in [0.05, 0.15]: X = th * np.array([[1j, 0], [0, -1j]]) V = sharp({((1,) + a, 0): sexp(X)}, (0, 0, 0, 0), 0)[0] good &= np.allclose(V, sexp(X / 8)) check('(F) one small link e^X at any transverse position gives the coarse link e^{X/8}', good) # (G) def tor(t): return np.diag([np.exp(1j * t), np.exp(-1j * t)]) def ang(g): return math.atan2(g[0, 0].imag, g[0, 0].real) def shz(x, mu, Sz, k=1): return tuple((np.array(x) + k * E[mu]) % Sz) def gfix(U, Sz): h = {} for y in itertools.product(range(0, Sz, 2), repeat=d): for o in sorted(blk, key=sum): x = tuple((np.array(y) + o) % Sz) if not any(o): h[x] = I2; continue nu = max(i for i in range(d) if o[i]); p = list(o); p[nu] = 0 xp = tuple((np.array(y) + p) % Sz); h[x] = h[xp] @ U[(xp, nu)] return {(x, mu): h[x] @ U[(x, mu)] @ inv(h[shz(x, mu, Sz)]) for (x, mu) in U} def sharpz(U, y, mu, Sz): P = [] for o in blk: x = tuple((np.array(y) + o) % Sz); P.append(U[(x, mu)] @ U[(shz(x, mu, Sz), mu)]) v = [inv(P[0]) @ p for p in P] m = max(math.sqrt(hs2(w - I2)) for w in v); chi = phi(m / rho) M = I2 if chi == 0 else sexp(chi * sum(slog(w) for w in v) / 16) return P[0] @ M, chi def plaq(U, x, mu, nu, Sz): return U[(x, mu)] @ U[(shz(x, mu, Sz), nu)] @ inv(U[(shz(x, nu, Sz), mu)]) @ inv(U[(x, nu)]) def angles(U, Sz): return {(x, mu, nu): ang(plaq(U, x, mu, nu, Sz)) for x in itertools.product(range(Sz), repeat=d) for mu in range(d) for nu in range(mu + 1, d)} def step(U, Sz): U = gfix(U, Sz); W, cmin = {}, 1.0 for y in itertools.product(range(0, Sz, 2), repeat=d): for mu in range(d): W[(tuple(q // 2 for q in y), mu)], c = sharpz(U, y, mu, Sz); cmin = min(cmin, c) return W, cmin Sz = 8; eta_G = rho / 12 good_formula, good_tree, rows = True, True, [] cph = rng.uniform(0, 2 * math.pi, size=(d, d)) for eps in [1e-5, 3e-3, 0.004]: U = {} g = {x: tor(rng.uniform(0, 2 * math.pi)) for x in itertools.product(range(Sz), repeat=d)} for x in itertools.product(range(Sz), repeat=d): for mu in range(d): a = eps * rng.normal() if eps != 0.004 else \ eps * sum(math.sin(2 * math.pi * x[n] / Sz + cph[mu, n]) for n in range(d) if n != mu) U[(x, mu)] = g[x] @ tor(a) @ inv(g[shz(x, mu, Sz)]) Th = angles(U, Sz); Q0 = sum(8 * math.sin(t / 2) ** 2 for t in Th.values()) A = sum(math.sqrt(2) * abs(t) for t in Th.values()); eta0 = max(math.sqrt(2) * abs(t) for t in Th.values()) cur, cs = U, Sz for j in range(2): Ug = gfix(cur, cs) good_tree &= all(np.allclose(Ug[(tuple((np.array(y) + o) % cs), mu)], I2) for y in itertools.product(range(0, cs, 2), repeat=d) for o in blk for mu in range(d) if tree(list(o), mu)) Tc = angles(cur, cs); mx = max(math.sqrt(2) * abs(t) for t in Tc.values()) nxt, cmin = step(cur, cs); ns = cs // 2; Tn = angles(nxt, ns) err = 0.0 if mx <= eta_G: for (y, mu, nu), t in Tn.items(): f = 0.0 for o in blk: x = tuple((2 * np.array(y) + o) % cs) for a, b in [(0, 0), (1, 0), (0, 1), (1, 1)]: z = tuple((np.array(x) + a * E[mu] + b * E[nu]) % cs); f += Tc[(z, mu, nu)] err = max(err, abs(t - f / 16)) good_formula &= err < 1e-9 and cmin == 1.0 L1c = sum(abs(t) for t in Tc.values()); L1n = sum(abs(t) for t in Tn.values()) L2c = sum(t * t for t in Tc.values()); L2n = sum(t * t for t in Tn.values()) mxn = max(abs(t) for t in Tn.values()) rows.append((eps, j, mx <= eta_G, cmin, err, L1n / L1c, L2n / L2c, mxn / max(abs(t) for t in Tc.values()))) cur, cs = nxt, ns Qk = sum(8 * math.sin(t / 2) ** 2 for t in angles(cur, cs).values()) rows.append(('Q', eps, eta0 * A <= eta_G ** 2, Qk / Q0)) for r in rows: if r[0] == 'Q': print(' eps = %g: eta0*A <= eta_G^2: %s, Q_2/Q = %.4f' % r[1:]) else: print(' eps = %g, step %d: hypothesis %s, min chi %.3f, formula err %.1e, L1 ratio %.4f, L2 ratio %.4f, sup ratio %.3f' % r) check('(G) the tree gauge fixing puts every tree link at I', good_tree) check('(G) torus fields: coarse angles are the means of the 2x2 fluxes, cut-off inactive', good_formula) check('(G) L1 ratio <= 1/4, L2 ratio <= 1, sup ratio <= 4 at every step', all(r[5] <= 0.25 + 1e-12 and r[6] <= 1 + 1e-12 and r[7] <= 4 + 1e-12 for r in rows if r[0] != 'Q' and r[2])) check('(G) Q_2(V#^2) <= 1.001 Q when eta0 * A <= eta_G^2', all(r[3] <= 1.001 for r in rows if r[0] == 'Q' and r[2]) and any(r[2] for r in rows if r[0] == 'Q')) # (H) def nrm(Z): return math.sqrt(hs2(Z)) def logs(U, Sz): return {(x, mu, nu): slog(plaq(U, x, mu, nu, Sz)) for x in itertools.product(range(Sz), repeat=d) for mu in range(d) for nu in range(mu + 1, d)} def stats(L): v = [nrm(Z) for Z in L.values()] return sum(t * t for t in v), sum(v), max(v) def pathp(U, x, mu, Sz): return U[(x, mu)] @ U[(shz(x, mu, Sz), mu)] def nonab_step(cur, cs): Ug = gfix(cur, cs); nxt, cmin = step(cur, cs); ns = cs // 2 Ln = logs(nxt, ns); worst, rmax = 0.0, 0.0 for (yc, mu, nu), Th in Ln.items(): y = 2 * np.array(yc) corner = [(y, mu), (y + 2 * E[mu], nu), (y + 2 * E[nu], mu), (y, nu)] K, Dmax = [], [0.0] * 4 for o in blk: Pl = [pathp(Ug, tuple((c + o) % cs), m, cs) for c, m in corner] K.append(slog(Pl[0] @ Pl[1] @ inv(Pl[2]) @ inv(Pl[3]))) for i, (c, m) in enumerate(corner): P0 = pathp(Ug, tuple(c % cs), m, cs) Dmax[i] = max(Dmax[i], nrm(slog(inv(P0) @ Pl[i]))) Fbar = sum(K) / 16; r = sum(Dmax) + nrm(K[0]) rmax = max(rmax, r); worst = max(worst, nrm(Th - Fbar) / (3.2 * r * r)) return nxt, ns, cmin, worst, rmax good_lem, good_steps, rows = True, True, [] for kind, eps in [('random', 4e-5), ('random', 7e-5), ('smooth', 1e-4)]: U = {} g = {x: haar() for x in itertools.product(range(Sz), repeat=d)} for x in itertools.product(range(Sz), repeat=d): for mu in range(d): a = eps * rng.normal(size=3) if kind == 'random' else \ eps * np.array([sum(math.sin(2 * math.pi * x[n] / Sz + cph[mu, n] + k) for n in range(d) if n != mu) for k in range(3)]) U[(x, mu)] = g[x] @ su2(a) @ inv(g[shz(x, mu, Sz)]) cur, cs = U, Sz for j in range(2): E0, A0, e0 = stats(logs(cur, cs)) nxt, ns, cmin, worst, rmax = nonab_step(cur, cs) E1, A1, e1 = stats(logs(nxt, ns)) ok_h = e0 <= rho / 1000 good_lem &= (not ok_h) or (cmin == 1.0 and worst <= 1 and rmax <= 28.26 * e0) good_steps &= (not ok_h) or (math.sqrt(E1) <= (1 + 6030 * e0) * math.sqrt(E0) and A1 <= A0 / 4 * (1 + 9147 * e0) and e1 <= 4 * e0 * (1 + 647 * e0)) rows.append((kind, eps, j, e0, ok_h, cmin, worst, math.sqrt(E1 / E0), 4 * A1 / A0, e1 / e0)) cur, cs = nxt, ns for r in rows: print(' %s eps = %g, step %d: eta = %.2e (<= rho/1000: %s), min chi %.0f, max err/(3.2 r^2) = %.3f, ' 'sqrt(E1/E0) = %.4f, 4 A1/A0 = %.4f, eta1/eta0 = %.3f' % r) check('(H) nonabelian small fields: cut-off inactive and |log V#_P - Fbar_P| <= 3.2 r_P^2, r_P <= 28.26 eta', good_lem) check('(H) sqrt E grows by at most 1 + 6030 eta, A by (1/4)(1 + 9147 eta), eta by 4(1 + 647 eta)', good_steps) check('(H) the nonabelian error term is present (max err/(3.2 r^2) > 0)', max(r[6] for r in rows) > 1e-4) # (I) def quad(f, a, b, n=20000): t = np.linspace(a, b, n + 1); w = np.full(n + 1, 2.0); w[0] = w[-1] = 1; w[1:-1:2] = 4 return float(np.sum(w * f(t)) * (b - a) / (3 * n)) Lt = 8 cubes = list(itertools.product(range(Lt), repeat=d)) def cube_plaqs(z): out = set() for w in itertools.product(range(2), repeat=d): for mu in range(d): for nu in range(mu + 1, d): if w[mu] == 0 and w[nu] == 0: out.add((tuple((np.array(z) + w) % Lt), mu, nu)) return out cp = {z: cube_plaqs(z) for z in cubes} inc = {} for z, P in cp.items(): for p in P: inc[p] = inc.get(p, 0) + 1 check('(I) a unit cube has 24 plaquettes and a plaquette lies in exactly 4 unit cubes', all(len(P) == 24 for P in cp.values()) and set(inc.values()) == {4} and len(inc) == 6 * Lt ** 4) nb = len([w for w in itertools.product((-1, 0, 1), repeat=d) if any(w)]) check('(I) a unit cube has 80 neighbours at l-infinity distance one', nb == 80 and 80 ** 2 <= 6400) def Zset(y, mu): return {tuple((2 * np.array(y) + t * E[mu] + np.array(w)) % Lt) for t in range(3) for w in itertools.product(range(2), repeat=d) if w[mu] == 0} Zs = {(y, mu): Zset(y, mu) for y in itertools.product(range(Lt // 2), repeat=d) for mu in range(d)} cnt = {} for k, Z in Zs.items(): for z in Z: cnt[z] = cnt.get(z, 0) + 1 check('(I) |Z_{y,mu}| = 24 and a cube lies in Z_{y,mu} for at most 8 coarse links', all(len(Z) == 24 for Z in Zs.values()) and max(cnt.values()) <= 8) C0 = 0.5 + 2 / 3 * math.log(24 * math.pi ** 3) b0 = 11 * 2 / (48 * math.pi ** 2) print(' C_0 = %.4f, 4*8*4*24*144 = %d, 8 b_0 = %.4f' % (C0, 4 * 8 * 4 * 24 * 144, 8 * b0)) check('(I) 4*8*4*24*144 <= 4.5e5 and the one-loop slope 8 b_0 = 11/(3 pi^2) is below 4/3', 4 * 8 * 4 * 24 * 144 <= 4.5e5 and abs(8 * b0 - 11 / (3 * math.pi ** 2)) < 1e-12 and 8 * b0 < 4 / 3) good = True for beta in [1, 3, 10, 100, 1e4]: r = 1 / (4 * math.sqrt(beta)) pr = 2 / math.pi * quad(lambda t: np.sin(t) ** 2, 0, r) good &= pr >= 1 / (24 * math.pi ** 3 * beta ** 1.5) check('(I) a Haar link has angle <= (4 sqrt beta)^-1 with probability >= 1/(24 pi^3 beta^{3/2})', good) vstar = 2 / math.pi * quad(lambda t: np.sin(t) ** 2, math.acos(0.75), math.pi / 3) good = True; rows = [] for beta in [0.5, 1, 2, 5, 10]: for rho_ in [1.0, 0.5, 1 / 12]: t0 = math.acos(1 - rho_ ** 2 / 4) for kap in np.linspace(0, 6 * beta, 7): m = kap # stabilise the exponentials num = quad(lambda t: np.exp(kap * np.cos(t) - m) * np.sin(t) ** 2, t0, math.pi) den = quad(lambda t: np.exp(kap * np.cos(t) - m) * np.sin(t) ** 2, 0, math.pi) good &= num / den >= vstar * math.exp(-3 * beta) * (1 - 1e-9) rows.append((beta, vstar * math.exp(-3 * beta))) print(' v_* = %.5f; p = %s' % (vstar, ', '.join('%.3e at beta %g' % (q, b) for b, q in rows))) check('(I) the angle exceeds arccos(1 - rho^2/4) with probability >= v_* e^{-3 beta} for kappa <= 6 beta', good) # cap about I under the density e^{kappa h.u}: decreasing in the angle phi of h, increasing in kappa def capprob(phi_, kap, th0, n=241): th = np.linspace(0, math.pi, n); al = np.linspace(0, math.pi, n) T, A = np.meshgrid(th, al, indexing='ij') f = np.exp(kap * (math.cos(phi_) * np.cos(T) + math.sin(phi_) * np.sin(T) * np.cos(A) - 1)) * np.sin(T) ** 2 * np.sin(A) W = np.outer(np.where(th <= th0, 1.0, 0.0), np.ones(n)) return float(np.trapezoid(np.trapezoid(f * W, al, axis=1), th) / np.trapezoid(np.trapezoid(f, al, axis=1), th)) good = True for kap in [1.0, 4.0, 12.0]: vals = [capprob(ph, kap, 0.6) for ph in np.linspace(0, math.pi, 9)] good &= all(vals[i] >= vals[i + 1] - 1e-6 for i in range(len(vals) - 1)) kv = [capprob(0.0, k, 0.6) for k in [0.0, 1.0, 4.0, 12.0]] good &= all(kv[i] <= kv[i + 1] + 1e-6 for i in range(3)) check('(I) the probability of a cap about I is largest at h = I and increases with kappa', good) # frustrated exterior on a lattice of side 4 S4 = 4; x0 = (1, 1, 1, 1) sh4 = lambda x, mu, k=1: tuple((np.array(x) + k * E[mu]) % S4) Uf = {(x, mu): I2.copy() for x in itertools.product(range(S4), repeat=d) for mu in range(d)} for nu, sgn in [(2, 1), (3, 1), (1, -1), (2, -1), (3, -1)]: Uf[(sh4(x0, nu, sgn), 0)] = -I2 def pl4(U, x, mu, nu): return U[(x, mu)] @ U[(sh4(x, mu), nu)] @ inv(U[(sh4(x, nu), mu)]) @ inv(U[(x, nu)]) good = True for rep in range(5): Ue = haar(); Uf[(x0, 0)] = Ue tot = 0.0 for nu in range(1, d): tot += pl4(Uf, x0, 0, nu).trace().real / 2 + pl4(Uf, sh4(x0, nu, -1), 0, nu).trace().real / 2 good &= abs(tot - (-2 * Ue.trace().real)) < 1e-10 and np.allclose(pl4(Uf, x0, 0, 1), Ue) check('(I) frustrated exterior: the six plaquettes of e add to -2 Tr U_e and p_0 has holonomy U_e', good) fr = [] for beta in [0.25, 0.5, 1, 2, 10]: t0 = math.acos(1 - 2.5 ** 2 / 4) num = quad(lambda t: np.exp(-4 * beta * (np.cos(t) + 1)) * np.sin(t) ** 2, t0, math.pi) den = quad(lambda t: np.exp(-4 * beta * (np.cos(t) + 1)) * np.sin(t) ** 2, 0, math.pi) fr.append(num / den) print(' frustrated exterior, rho = 2.5: P(p_0 bad) = %s at beta = 0.25, 0.5, 1, 2, 10' % ', '.join('%.6f' % v for v in fr)) check('(I) under the frustrated exterior p_0 is bad with probability increasing to 1', all(fr[i] <= fr[i + 1] + 1e-12 for i in range(4)) and fr[-1] > 0.999) # chi = 1 off the bad cubes rho_b = 1 / 12; good, nact, nfree, nnear = True, 0, 0, 0 for rep in range(4): g = {x: haar() for x in itertools.product(range(Sz), repeat=d)} centres = [tuple(rng.integers(0, Sz, size=d)) for _ in range(2)] U = {} for x in itertools.product(range(Sz), repeat=d): for mu in range(d): big = any(max(min((x[n] - c[n]) % Sz, (c[n] - x[n]) % Sz) for n in range(d)) == 0 for c in centres) a = rng.normal(size=3) * (0.8 if big else 1.5e-3) U[(x, mu)] = g[x] @ su2(a) @ inv(g[shz(x, mu, Sz)]) bad = set() for x in itertools.product(range(Sz), repeat=d): for mu in range(d): for nu in range(mu + 1, d): if nrm(plaq(U, x, mu, nu, Sz) - I2) > rho_b: for w in itertools.product(range(2), repeat=d): if w[mu] == 0 and w[nu] == 0: bad.add(tuple((np.array(x) - np.array(w)) % Sz)) Ug = gfix(U, Sz) for y in itertools.product(range(0, Sz, 2), repeat=d): for mu in range(d): Z = {tuple((np.array(y) + t * E[mu] + np.array(w)) % Sz) for t in range(3) for w in itertools.product(range(2), repeat=d) if w[mu] == 0} chi = sharpz(Ug, y, mu, Sz)[1] if Z & bad: nnear += 1; nact += chi < 1 else: nfree += 1; good &= chi == 1.0 print(' coarse links with a bad cube in Z: %d (chi < 1 at %d); without: %d' % (nnear, nact, nfree)) check('(I) chi = 1 at every coarse link whose 24 cubes are good (rho = 1/12)', good and nfree > 0 and nnear > 0) print('%d checks passed, %d failed' % (sum(ok), len(ok) - sum(ok)))