"""The flat kappa = 1 Bloch operator of the Kuhn lattice, rebuilt from Appendix B, with checks against the paper. Paper: "Linearized Einstein gravity from Lorentzian Regge calculus on the Kuhn lattice, and the divergence of its one-loop graviton limit", J. Washburn (2026), https://doi.org/10.5281/zenodo.23764491. Section, theorem and appendix numbers below are the paper's. What this program does * Assembles the Bloch matrix B(k) from the cell Hessian C/(16 sqrt3) printed in Appendix B, the vertex-move columns Gamma(k) of Section 4.2, the gauge-fixed operator A(k) = B(k) + Gamma(k) Gamma(k)^dagger, and its midpoint frame M(k). * Checks, and prints: B is Hermitian and satisfies the Bianchi identity B Gamma = 0; the kernel of B(0) is the ten constant strains; Proposition 6.1 at every axis corner pi e_c (v_c(d) = 2^{|d|-3} for d_c = 0, else 0; A v_c = 0, B v_c = 0, Gamma^dagger v_c = 0, |v_c|^2 = 31/16, a one-dimensional kernel, not a vertex move); the other half-period corners; Proposition 6.2 (det K = -1013375, one negative eigenvalue, the Schur complement (sqrt3/640) h^T K h near pi e_0); the numbers of Theorem 6.15 (q and h^T K h at record time, the doubler's time axis); and the link to Appendix F: on the paired covector u = (z+t, z+t, z-t, z-t) the cross sector is one-dimensional and M there equals -(2/sqrt3) G_x(U, V). * Checks the static field (Section 3.5, Theorem 5.4). A resting clock of proper mass M_p on a time edge has the action M_p sqrt(-l), so its source on that edge is tau = M_p/(2 sqrt12) at l = -12; the field equation is s H_0 h = -tau; the time-time readout is h_tt = h_(time edge)/12; and c_N = 8 pi G 2 M_p = 2 M_p/s with G = 1/(8 pi s). The ratio chi_tt q / c_N then equals (A^-1)_(time edge, time edge) q(u) / (-96 sqrt3), and it tends to 1 as the static wave u -> 0. Also checks the static resonance printed in Appendix B: the vector W at u = (-t*, t*, 0, 0), cos t* = 9/11, satisfies B(u) W = 0 and Gamma(u)^dagger W = 0, with its largest entry on the time edge. * Checks Proposition 4.5. Along the axis line z = (w, 1, 1, 1) and the null line z = (w, w, w, 1), w = e^{is}, the printed closed forms of the integer polynomial F1 satisfy det A(k) (z0 z1 z2 z3)^3 = const * F1(z) det(Gamma^dagger Gamma)^2 with one constant on both lines; A is singular at s = pi on the axis line (the corner pi e_0) and, on the null line, at the zeros of F1: cos s = 2/5 and the two roots of 3x^3 - 6x^2 - 12x + 16 with x = 2 cos s in [-2, 2]; and A is invertible on the null line for 0 < s <= arccos(3/5) at the sampled points. * With --static FILE.json, writes the data of the figure "The static field of a resting mass" on the gravity page: the ratio chi_tt q / c_N along three static directions, and where each direction meets the static resonance. * With --slice FILE.json, recomputes the data of the figure "Where the kappa = 1 operator vanishes" on the gravity page: the sign of the determinant of M in its two symmetry sectors on the plane k = a e_0 + b (1,1,1,1), and the grid points where it changes, for the graviton's E sector and the doubler's symmetric sector. Requires Python 3.9+ and numpy. Run: python3 kuhn_bloch.py [--static static.json] [--slice slice.json] License: MIT. Copyright (c) 2026 Recognition Physics Institute. """ import argparse import itertools import json import math import sys import numpy as np # Appendix B: Hessian of one simplex at the flat kappa = 1 lattice, d^2 S / d l_ij d l_kl = C / (16 sqrt3), edges ordered # 01, 02, 03, 04, 12, 13, 14, 23, 24, 34 along the simplex chain v_0 = x, v_j = x + e_a1 + ... + e_aj. C = np.array([ [2, -2, 1, 0, 0, 3, -2, -6, 2, 2], [-2, -4, 4, -1, 6, -6, 2, 3, -2, 2], [1, 4, 0, -2, -6, 0, 2, -3, 2, -2], [0, -1, -2, 2, 2, 2, -2, 2, -1, 0], [0, 6, -6, 2, -6, 6, -3, 0, 3, -6], [3, -6, 0, 2, 6, 0, 0, 6, -6, 3], [-2, 2, 2, -2, -3, 0, 0, -6, 4, 1], [-6, 3, -3, 2, 0, 6, -6, -6, 6, 0], [2, -2, 2, -1, 3, -6, 4, 6, -4, -2], [2, 2, -2, 0, -6, 3, 1, 0, -2, 2]], dtype=float) PAIRS = [(0, 1), (0, 2), (0, 3), (0, 4), (1, 2), (1, 3), (1, 4), (2, 3), (2, 4), (3, 4)] # The fifteen edge classes d in {0,1}^4 minus 0, indexed by bitmask: index = sum_a d_a 2^a - 1 (axes a = 0..3). CLASSES = np.array([[(m >> a) & 1 for a in range(4)] for m in range(1, 16)], dtype=float) K = np.array([[-2, 9, 9, 9], [9, 52, -3, -3], [9, -3, 52, -3], [9, -3, -3, 52]], dtype=float) # Proposition 6.2 def class_index(d): return int(sum(int(d[a]) << a for a in range(4))) - 1 def _terms(): """(class of f, class of f', C_ff'/(16 sqrt3), b_f' - b_f) over the 24 simplices of one cell.""" rows = [] for perm in itertools.permutations(range(4)): verts = [np.zeros(4, dtype=int)] for a in perm: v = verts[-1].copy() v[a] = 1 verts.append(v) edges = [(class_index(verts[j] - verts[i]), verts[i]) for i, j in PAIRS] for f, (df, bf) in enumerate(edges): for g, (dg, bg) in enumerate(edges): if C[f, g]: rows.append((df, dg, C[f, g] / (16 * math.sqrt(3)), bg - bf)) return (np.array([r[0] for r in rows]), np.array([r[1] for r in rows]), np.array([r[2] for r in rows]), np.array([r[3] for r in rows], dtype=float)) _TD, _TG, _TW, _TS = _terms() def B(k): """Bloch matrix of the flat Hessian (Section 4.2, Appendix B).""" out = np.zeros((15, 15), dtype=complex) np.add.at(out, (_TD, _TG), _TW * np.exp(1j * (_TS @ np.asarray(k, dtype=float)))) return out def Gamma(k): """Vertex-move columns: (Gamma(k) eta)_d = 2 (exp(i k.d) - 1) d.eta (Section 4.2).""" return 2.0 * (np.exp(1j * (CLASSES @ np.asarray(k, dtype=float))) - 1.0)[:, None] * CLASSES def A(k): """Gauge-fixed operator in the Feynman-type gauge.""" G = Gamma(k) return B(k) + G @ G.conj().T def M(k): """Midpoint frame Phi^dagger A Phi, real symmetric.""" phi = np.exp(0.5j * (CLASSES @ np.asarray(k, dtype=float))) return ((phi.conj()[:, None] * A(k)) * phi[None, :]).real def axis_permutation(perm): """Permutation of the edge classes induced by sending axis a to axis perm[a].""" P = np.zeros((15, 15)) for i, d in enumerate(CLASSES.astype(int)): e = np.zeros(4, dtype=int) for a in range(4): e[perm[a]] = d[a] P[class_index(e), i] = 1 return P def _range_basis(Q): w, V = np.linalg.eigh((Q + Q.T) / 2) return V[:, w > 0.5] TIME = 14 # class index of the time edge (1,1,1,1) def newton_ratio(u): """chi_tt(u) q(u) / c_N for a static covector u (sum u_a = 0), in the normalization of the module docstring.""" u = np.asarray(u, dtype=float) e = np.zeros(15, dtype=complex) e[TIME] = 1.0 h = np.linalg.solve(A(u), e) return float((h[TIME] * (u @ u)).real / (-96 * math.sqrt(3))) STATIC_DIRECTIONS = {'(-1,1,0,0)': (-1, 1, 0, 0), '(3,-1,-1,-1)': (3, -1, -1, -1), '(1,1,-1,-1)': (1, 1, -1, -1)} def resonance(direction, k_hi=3.1, n=600): """Smallest |k| > 0 along a static direction where the response changes sign through a pole (bisection on 1/ratio).""" d = np.asarray(direction, dtype=float); d /= np.linalg.norm(d) ks = np.linspace(1e-3, k_hi, n) inv = [1 / newton_ratio(k * d) for k in ks] for (k0, f0), (k1, f1) in zip(zip(ks, inv), zip(ks[1:], inv[1:])): if f0 > 0 >= f1: for _ in range(60): km = (k0 + k1) / 2 fm = 1 / newton_ratio(km * d) k0, k1 = (km, k1) if fm > 0 else (k0, km) return (k0 + k1) / 2 return None def static_data(path, step=0.005, k_max=2.7): out = {'quantity': 'chi_tt(k) q(k) / c_N along static directions; |k| = sqrt(sum u_a^2) for static u', 'directions': {}} for name, d in STATIC_DIRECTIONS.items(): v = np.asarray(d, dtype=float); v /= np.linalg.norm(v) ks = [step * i for i in range(1, int(round(k_max / step)) + 1)] out['directions'][name] = {'resonance_k': resonance(d), 'k': [round(k, 4) for k in ks], 'ratio': [round(newton_ratio(k * v), 6) for k in ks]} print(f' {name}: static resonance at |k| = {out["directions"][name]["resonance_k"]:.6f}') with open(path, 'w') as fh: json.dump(out, fh) print(f' wrote {path}') def report(name, ok, detail): print(f" [{'PASS' if ok else 'FAIL'}] {name}: {detail}") return ok def run_checks(): rng = np.random.default_rng(2026) ok = True print('Appendix B and Section 4.2: the operator') res_h = res_b = 0.0 for _ in range(20): k = rng.uniform(-math.pi, math.pi, 4) Bk = B(k) res_h = max(res_h, np.abs(Bk - Bk.conj().T).max()) res_b = max(res_b, np.abs(Bk @ Gamma(k)).max()) ok &= report('B(k) Hermitian at 20 random waves', res_h < 1e-12, f'max |B - B^dagger| = {res_h:.1e}') ok &= report('Bianchi identity B(k) Gamma(k) = 0', res_b < 1e-12, f'max residual = {res_b:.1e}') w = np.linalg.eigvalsh(B(np.zeros(4))) H = rng.normal(size=(4, 4)); H = H + H.T strain = np.array([d @ H @ d for d in CLASSES]) ok &= report('kernel of B(0) = the ten constant strains d -> d^T H d', int((np.abs(w) < 1e-10).sum()) == 10 and np.abs(B(np.zeros(4)) @ strain).max() < 1e-12, f'{int((np.abs(w) < 1e-10).sum())} zero eigenvalues') print('Proposition 6.1: the axis-corner zero modes') for c in range(4): k = np.zeros(4); k[c] = math.pi v = np.array([2.0 ** (d.sum() - 3) if d[c] == 0 else 0.0 for d in CLASSES]) G = Gamma(k) dimker = int((np.abs(np.linalg.eigvalsh(A(k))) < 1e-10).sum()) coef, *_ = np.linalg.lstsq(G, v, rcond=None) not_move = np.linalg.norm(G @ coef - v) > 1e-6 ok &= report(f'corner pi e_{c}', np.abs(B(k) @ v).max() < 1e-12 and np.abs(G.conj().T @ v).max() < 1e-12 and abs(v @ v - 31 / 16) < 1e-15 and dimker == 1 and not_move, f'B v = 0, Gamma^dagger v = 0, |v|^2 = {v @ v} = 31/16, dim ker A = {dimker}, not a vertex move') for n in (2, 3, 4): dims = set() for odd in itertools.combinations(range(4), n): k = np.zeros(4); k[list(odd)] = math.pi dims.add(int((np.abs(np.linalg.eigvalsh(A(k))) < 1e-10).sum())) ok &= report(f'corners with {n} coordinates equal to pi', dims == {0}, 'A invertible') print('Proposition 6.2: the corner cone') wK = np.linalg.eigvalsh(K) ok &= report('det K = -1013375 = -55^2 * 335, one negative eigenvalue', round(np.linalg.det(K)) == -1013375 and (wK < 0).sum() == 1, f'det K = {np.linalg.det(K):.1f}, eigenvalues {np.round(wK, 4)}') v0 = np.array([2.0 ** (d.sum() - 3) if d[0] == 0 else 0.0 for d in CLASSES]) rest = np.linalg.svd(np.eye(15) - np.outer(v0, v0) / (v0 @ v0))[0][:, :14] worst_s = worst_r = 0.0 for _ in range(6): h = rng.normal(size=4); h *= 1e-3 / np.linalg.norm(h) Ak = A(np.array([math.pi, 0, 0, 0]) + h) a11 = v0 @ Ak @ v0; a12 = v0 @ Ak @ rest; a22 = rest.T @ Ak @ rest schur = (a11 - a12 @ np.linalg.solve(a22, a12.conj())).real lead = math.sqrt(3) / 640 * (h @ K @ h) worst_s = max(worst_s, abs(schur / lead - 1)) ev = np.linalg.eigvalsh(Ak); r = ev[np.argmin(np.abs(ev))] worst_r = max(worst_r, abs(r / (math.sqrt(3) / 1240 * (h @ K @ h)) - 1)) ok &= report('Schur complement onto v_0 = (sqrt3/640) h^T K h + o(|h|^2)', worst_s < 1e-2, f'worst relative deviation {worst_s:.1e} at |h| = 1e-3') ok &= report('doubler eigenvalue = (sqrt3/1240) h^T K h + o(|h|^2)', worst_r < 1e-2, f'worst relative deviation {worst_r:.1e} at |h| = 1e-3 (16/31 of the Schur form)') print('Theorem 6.15: no common time') one = np.ones(4) q = lambda u: u @ u - u.sum() ** 2 / 3 ok &= report('record time (1,1,1,1): q = -4/3, h^T K h = 190', abs(q(one) + 4 / 3) < 1e-12 and abs(one @ K @ one - 190) < 1e-9, f'q = {q(one):.6f}, h^T K h = {one @ K @ one:.1f}') wK, VK = np.linalg.eigh(K) t = VK[:, 0] / VK[0, 0] angle = lambda x, y: math.degrees(math.acos(abs(x @ y) / np.linalg.norm(x) / np.linalg.norm(y))) print(f' [INFO] doubler time axis {np.round(t, 3)}: {angle(t, np.eye(4)[0]):.1f} deg from e_0, ' f'{angle(t, one):.1f} deg from record time (numerical, Euclidean angles in record coordinates)') worst = 0.0 for _ in range(20): k = rng.normal(size=4); S = k[1:].sum() rhs = (18 * (k[0] - S / 18) ** 2 + 85 / 3 * sum((k[i] - k[j]) ** 2 for i in range(1, 4) for j in range(i + 1, 4)) + 275 / 18 * S ** 2) worst = max(worst, abs(30 * q(k) + k @ K @ k - rhs)) ok &= report('30 q + h^T K h = 18(k_0 - S/18)^2 + (85/3) sum (k_i - k_j)^2 + (275/18) S^2', worst < 1e-9, f'max residual {worst:.1e}') print('Section 3.5 and Theorem 5.4: the static field') worst = 0.0 for _ in range(6): u = rng.normal(size=4); u -= u.mean(); u *= 1e-3 / np.linalg.norm(u) worst = max(worst, abs(newton_ratio(u) - 1)) ok &= report('chi_tt q / c_N -> 1 with G = 1/(8 pi s), every static direction', worst < 1e-5, f'max |ratio - 1| = {worst:.1e} at |k| = 1e-3 in 6 random static directions') zeta = (9 + 2 * math.sqrt(10) * 1j) / 11 ts = math.atan2(2 * math.sqrt(10), 9) W440 = {'1000': -1131 + 429 * zeta, '0100': -429 - 429 * zeta, '1100': 1760, '0010': 440, '1010': -87 + 33 * zeta, '0110': -33 - 33 * zeta, '1110': 2640, '0001': 440, '1001': -87 + 33 * zeta, '0101': -33 - 33 * zeta, '1101': 2640, '0011': 880, '1011': 1305 - 495 * zeta, '0111': 495 + 495 * zeta, '1111': 3960} W = np.zeros(15, dtype=complex) for bits, value in W440.items(): W[class_index([int(c) for c in bits])] = value / 440 u = np.array([-ts, ts, 0.0, 0.0]) resB, resG = np.abs(B(u) @ W).max(), np.abs(Gamma(u).conj().T @ W).max() ok &= report('Appendix B static resonance W at u = (-t*, t*, 0, 0), cos t* = 9/11', resB < 1e-12 and resG < 1e-12 and int(np.argmax(np.abs(W))) == TIME, f'|B W| = {resB:.1e}, |Gamma^dagger W| = {resG:.1e}, largest entry on 1111') k_res = resonance(STATIC_DIRECTIONS['(-1,1,0,0)']) ok &= report('the response along (-1,1,0,0) has its pole at |k| = sqrt2 t*', abs(k_res - math.sqrt(2) * ts) < 1e-9, f'pole at {k_res:.9f}, sqrt2 t* = {math.sqrt(2) * ts:.9f}') print('Appendix F: the cross sector on the paired direction') P01, P23 = axis_permutation([1, 0, 2, 3]), axis_permutation([0, 1, 3, 2]) cross = _range_basis((np.eye(15) - P01) @ (np.eye(15) - P23) / 4) worst = 0.0 for _ in range(20): t_, z = rng.uniform(0.02, 1.5, 2) U, V = math.sin(t_ / 2) ** 2, math.sin(z / 2) ** 2 gx = -6 * U * V + 3 * U + 4 * V * V - V val = (cross.T @ M(np.array([z + t_, z + t_, z - t_, z - t_])) @ cross)[0, 0] worst = max(worst, abs(val + 2 / math.sqrt(3) * gx)) ok &= report('cross sector one-dimensional, M = -(2/sqrt3) G_x(U, V)', cross.shape[1] == 1 and worst < 1e-12, f'max |M_cross + (2/sqrt3) G_x| = {worst:.1e} at 20 random (t, z)') print('Proposition 4.5: F1 on an axis line and a null line') f1_axis = lambda w: -8 * (w - 1) ** 4 * (w + 1) ** 2 f1_null = lambda w: ((w - 1) ** 8 * (5 * w * w - 4 * w + 5) ** 2 * (3 * w ** 6 - 6 * w ** 5 - 3 * w ** 4 + 4 * w ** 3 - 3 * w ** 2 - 6 * w + 3)) ratios = [] for f1, line in ((f1_axis, lambda s: [s, 0, 0, 0]), (f1_null, lambda s: [s, s, s, 0])): for s in (0.3, 0.7, 1.3, 2.0, 2.9): k = np.array(line(s), dtype=float) G = Gamma(k) ratios.append(np.linalg.det(A(k)) * np.exp(3j * k.sum()) / (np.linalg.det(G.conj().T @ G) ** 2 * f1(np.exp(1j * s)))) ratios = np.array(ratios) spread = np.abs(ratios / ratios[0] - 1).max() ok &= report('det A (z0 z1 z2 z3)^3 / (det(Gamma^dagger Gamma)^2 F1) is one constant on both lines', spread < 1e-9, f'constant {ratios[0].real:.6e}, relative spread {spread:.1e} at 5 points on each line') smallest = lambda k: np.abs(np.linalg.eigvalsh(A(np.array(k, dtype=float)))).min() zeros = sorted([math.acos(x.real / 2) for x in np.roots([3, -6, -12, 16]) if abs(x.imag) < 1e-12 and abs(x.real) <= 2] + [math.acos(2 / 5)]) at_zeros = max([smallest([s, s, s, 0]) for s in zeros] + [smallest([math.pi, 0, 0, 0])]) ok &= report('A singular at s = pi on the axis line and at the zeros of F1 on the null line', at_zeros < 1e-10, f'null-line zeros at s = {", ".join(f"{s:.5f}" for s in zeros)}; largest smallest |eigenvalue| there {at_zeros:.1e}') away = min(smallest([s, s, s, 0]) for s in np.linspace(0.05, math.acos(3 / 5), 12)) ok &= report('A invertible on the null line for 0 < s <= arccos(3/5)', away > 1e-8, f'smallest |eigenvalue| {away:.1e} over 12 sampled points') return ok def slice_data(path, n_a=480, n_b=240): """Sign changes of det M in the two symmetry sectors on k = a e_0 + b (1,1,1,1), a in [-pi/2, 3pi/2], b in [-pi/2, pi/2].""" perms = [[0] + list(p) for p in itertools.permutations([1, 2, 3])] sym = sum(axis_permutation(p) for p in perms) / 6 sgn = sum(np.linalg.det(np.eye(4)[p]) * axis_permutation(p) for p in perms) / 6 E = _range_basis((np.eye(15) - sym - sgn) @ (np.eye(15) + axis_permutation([0, 1, 3, 2])) / 2) T = _range_basis(sym) A_ = np.linspace(-math.pi / 2, 3 * math.pi / 2, n_a) B_ = np.linspace(-math.pi / 2, math.pi / 2, n_b) sign = {'graviton_E_sector': np.zeros((n_b, n_a)), 'doubler_symmetric_sector': np.zeros((n_b, n_a))} for i, b in enumerate(B_): for j, a in enumerate(A_): Mk = M(np.array([a + b, b, b, b])) sign['graviton_E_sector'][i, j] = np.sign(np.linalg.det(E.T @ Mk @ E)) sign['doubler_symmetric_sector'][i, j] = np.sign(np.linalg.det(T.T @ Mk @ T)) out = {'plane': 'k = a e_0 + b (1,1,1,1)', 'a': [float(A_[0]), float(A_[-1]), n_a], 'b': [float(B_[0]), float(B_[-1]), n_b], 'crossings': {}} for name, S in sign.items(): pts = [] for i in range(n_b): for j in range(n_a): if j + 1 < n_a and S[i, j] != S[i, j + 1]: pts.append([round((A_[j] + A_[j + 1]) / 2, 4), round(B_[i], 4)]) if i + 1 < n_b and S[i, j] != S[i + 1, j]: pts.append([round(A_[j], 4), round((B_[i] + B_[i + 1]) / 2, 4)]) out['crossings'][name] = pts print(f' {name}: {len(pts)} grid crossings') with open(path, 'w') as fh: json.dump(out, fh) print(f' wrote {path}') if __name__ == '__main__': parser = argparse.ArgumentParser(description='Rebuild the kappa = 1 Bloch operator from Appendix B and check it.') parser.add_argument('--slice', metavar='FILE.json', help='also write the zero-set data of the plane k = a e_0 + b(1,1,1,1)') parser.add_argument('--static', metavar='FILE.json', help='also write the static-response data along three static directions') args = parser.parse_args() good = run_checks() if args.static: print('Figure data: the static field of a resting mass') static_data(args.static) if args.slice: print('Figure data: the plane k = a e_0 + b (1,1,1,1)') slice_data(args.slice) print('ALL CHECKS PASS' if good else 'SOME CHECKS FAILED') sys.exit(0 if good else 1)