05-ilang-free-fall / 14-pair-level-next-order
14pair level next orderverified
Determines the bound level of a pair of two bodies one order beyond leading order: binding cost, curvature, rest-term weights and curvature inertia, and whether the inertia follows the binding cost.
check_pair_level.py17 KB · raw
"""Check of 14-pair-level-next-order: numerical diagonalization of the fibre contract H(P) of (11.1)
for a diagonal on-site binding, compared with the leading forms (11.14) and the next-order forms
(14.6)-(14.11). Also checks the effective operator (14.4) on a ring and the small-u forms (14.12)-(14.13)
(exact series with sympy; leading terms on single and double cancellation sets in 50-digit arithmetic).
Output: output/check_pair_level.log (all numbers cited in derivation.v2.md)."""
from pathlib import Path
import numpy as np
import sympy as sp
import mpmath as mp
OUT = Path(__file__).parent / "output"
OUT.mkdir(exist_ok=True)
LOG = []
def log(s=""):
print(s)
LOG.append(s)
# hatted data sets: (tau1, tau2, M1, M2, w11, w12, w21, w22); all w scale as eps^2, tau as eps
DATA = {
"D1": (1.0, 1.6, 1.0, 2.5, -0.40, -0.50, 0.30, 0.70),
"D2": (1.3, 0.8, 1.5, 0.7, +0.60, 0.40, -0.20, -0.50),
"D3": (1.2, 0.9, 0.8, 1.3, -0.12, 0.20, 0.10, 0.30),
}
EPS = [0.2, 0.1, 0.05]
ETA = 0.5
H_FD = 0.05 # finite-difference step in P (7-point stencil); 2*H_FD as control
P_LIST = [0.0, H_FD, 2 * H_FD, 3 * H_FD, 4 * H_FD, 6 * H_FD]
def fibre_sector(P, L, t1, t2, M1, M2, W, eta=ETA):
"""H(P) of (11.1) on r in [-L, L], restricted to the sector Theta = (-1)^r s3(1) s3(2) = +1."""
sig3 = (1, -1) # branch index 0 = branch 1 (upper), 1 = branch 2
states = [(r, s1, s2) for r in range(-L, L + 1) for s1 in (0, 1) for s2 in (0, 1)
if (-1) ** (r % 2) * sig3[s1] * sig3[s2] == 1]
idx = {st: k for k, st in enumerate(states)}
n = len(states)
H = np.zeros((n, n), complex)
e1p, e1m = np.exp(1j * eta * P), np.exp(-1j * eta * P)
e2p, e2m = np.exp(1j * (1 - eta) * P), np.exp(-1j * (1 - eta) * P)
for (r, s1, s2), k in idx.items():
H[k, k] += M1 * sig3[s1] + M2 * sig3[s2] + (W[s1][s2] if r == 0 else 0.0)
# A1 = (i t1/2) sigma1^(1) (e^{i eta P} S^dag - e^{-i eta P} S); (S^dag psi)(r) = psi(r+1)
for rr, c in ((r + 1, 0.5j * t1 * e1p), (r - 1, -0.5j * t1 * e1m)):
k2 = idx.get((rr, 1 - s1, s2))
if k2 is not None:
H[k, k2] += c
# A2 = (i t2/2) sigma1^(2) (e^{i(1-eta)P} S - e^{-i(1-eta)P} S^dag); (S psi)(r) = psi(r-1)
for rr, c in ((r - 1, 0.5j * t2 * e2p), (r + 1, -0.5j * t2 * e2m)):
k2 = idx.get((rr, s1, 1 - s2))
if k2 is not None:
H[k, k2] += c
s3a = np.array([sig3[s1] for (_, s1, _) in states], float)
s3b = np.array([sig3[s2] for (_, _, s2) in states], float)
return H, s3a, s3b
def bound_level(P, L, t1, t2, M1, M2, W, vec=False):
H, s3a, s3b = fibre_sector(P, L, t1, t2, M1, M2, W)
E0 = M1 + M2
ev, U = np.linalg.eigh(H)
top = np.where(ev > E0 - min(M1, M2))[0] # cluster of the band B_++ and its bound level
k = top[-1] if W[0][0] > 0 else top[0]
if not vec:
return ev[k]
v = U[:, k]
return ev[k], float(np.vdot(v, s3a * v).real), float(np.vdot(v, s3b * v).real)
def coef(a1, a2, M1, M2, u, sg, sqrt):
"""eps^4 coefficients of (14.7)-(14.10) at given u and sign sg of w11 (numpy, mpmath or sympy)."""
a, b1, b2 = a1 + a2, a1 / M1, a2 / M2
beta, sig, nu = b1 + b2, sqrt(1 + u * u), a1 * a2 / (a1 + a2) ** 2
gs, gc, gm = a1 ** 2 / M1 + a2 ** 2 / M2, a2 ** 2 / M1 + a1 ** 2 / M2, a1 * a2 * (1 / M1 + 1 / M2)
sps = u * u / (sig + 1) if sg < 0 else sig + 1 # sig + sg, without cancellation for small u
phA = 2 * sps / sig - u ** 2 + u ** 3 / sig
phB = u ** 2 * (sg + sig - u) / sig
dphA = -2 * u + 3 * u ** 2 / sig - u * (2 * sg + u ** 3) / sig ** 3
dphB = 2 * u - 3 * u ** 2 / sig + (sg * u * (2 + u ** 2) + u ** 4) / sig ** 3
E4 = -gs / 2 * phA - (gs + gm) / 2 * phB
C = (2 * u ** 2 / sig ** 3 * (gs * (sg + sig ** 3 - u - u ** 3) - gm * (sg + sig ** 3 - 2 * u - u ** 3))
+ 4 * gm * (1 + sg / sig - 2 * u ** 2 + 2 * u ** 3 / sig) + 2 * u ** 2 * (sig - u) / sig * (gc + gm))
Xi = [3 * bi * phA / 2 - gs / (2 * a) * u * dphA + (beta + 2 * a / Mi) * phB / 2 - beta / 2 * u * dphB
for bi, Mi in ((b1, M1), (b2, M2))]
return dict(a=a, b=(b1, b2), sig=sig, nu=nu, E4=E4, C=C, Xi=Xi, rel=sg * sig * C / (4 * a))
def predictions(t1, t2, M1, M2, w):
"""Leading (11.14) and next-order (14.7)-(14.11) forms; unhatted inputs."""
a1, a2 = t1 ** 2 / (4 * M1), t2 ** 2 / (4 * M2)
sg, u = np.sign(w), abs(w) / (a1 + a2)
c = coef(a1, a2, M1, M2, u, sg, np.sqrt)
a, (b1, b2), sig, nu, E4, C, Xi = c["a"], c["b"], c["sig"], c["nu"], c["E4"], c["C"], c["Xi"]
Lam = 1 + sg / sig
Eb0, Epp0 = a * (1 + sg * sig), -4 * sg * nu * a / sig
return dict(u=u, Eb0=Eb0, Eb1=Eb0 + E4, Epp0=Epp0, Epp1=Epp0 + nu * C,
s30=[1 - b1 * Lam, 1 - b2 * Lam], s31=[1 - b1 * Lam + b1 * Xi[0], 1 - b2 * Lam + b2 * Xi[1]],
mb0=1 / Epp0, mb1=(1 / Epp0) * (1 + c["rel"]), D=-(sig ** 2 * C - 4 * E4) / (16 * a1 * a2))
def Eb_closed(P, t1, t2, M1, M2, w):
"""E_b(P) - (M1+M2) of (14.6): leading part and order-eps^4 part."""
a1, a2 = t1 ** 2 / (4 * M1), t2 ** 2 / (4 * M2)
b1, b2 = a1 / M1, a2 / M2
c2, c4 = np.cos(2 * P), np.cos(4 * P)
s = np.sqrt(w * w + a1 ** 2 + a2 ** 2 + 2 * a1 * a2 * c2)
aw, sg = abs(w), np.sign(w)
k1, k2 = a1 + a2 * c2, a2 + a1 * c2
l1, l2 = a1 ** 2 + 2 * a1 * a2 * c2 + a2 ** 2 * c4, a2 ** 2 + 2 * a1 * a2 * c2 + a1 ** 2 * c4
g = (s + 2 * aw) / (s * (s + aw) ** 2)
lead = a1 + a2 + sg * s
corr = (-(a1 * b1 / 8) * (6 + 8 * sg * k1 / s + 2 * g * l1) - (a2 * b2 / 8) * (6 + 8 * sg * k2 / s + 2 * g * l2)
- w * w / (2 * s) * (sg * (b1 + b2) + (b1 * k1 + b2 * k2) / (s + aw)))
return lead, corr
def d2_7pt(f, h):
"""second derivative at 0 of an even function from f(0), f(h), f(2h), f(3h)"""
return (4 * f[3] - 54 * f[2] + 540 * f[1] - 490 * f[0]) / (180 * h * h)
def slope(dev, eps):
return [np.log(abs(dev[k] / dev[k + 1])) / np.log(eps[k] / eps[k + 1]) for k in range(len(dev) - 1)]
# ---------------------------------------------------------------- Part A: bound level of H(P)
log("Part A: bound level of H(P) (11.1), sector Theta=+1, open chain r in [-L, L], eta = %.2f" % ETA)
log("FD: 7-point stencil for E_b''(0), step h = %.3f; control step 2h = %.3f" % (H_FD, 2 * H_FD))
summary = {}
for name, (T1, T2, M1, M2, W11, W12, W21, W22) in DATA.items():
a_hat = T1 ** 2 / (4 * M1) + T2 ** 2 / (4 * M2)
u = abs(W11) / a_hat
rho = np.sqrt(1 + u * u) - u
L = 2 * int(np.ceil(40 / np.log(1 / rho))) # |psi(L)| ~ rho^(L/2) <= e^-40
log("")
log("%s: tau^=(%.2f,%.2f) M=(%.2f,%.2f) w^=(w11 %.2f, w12 %.2f, w21 %.2f, w22 %.2f) u=%.4f L=%d"
% (name, T1, T2, M1, M2, W11, W12, W21, W22, u, L))
rows = []
for eps in EPS:
t1, t2, w = eps * T1, eps * T2, eps ** 2 * W11
W = [[eps ** 2 * W11, eps ** 2 * W12], [eps ** 2 * W21, eps ** 2 * W22]]
res = {}
for LL in (L, int(1.5 * L)):
E = []
for P in P_LIST:
if P == 0.0:
e, s1, s2 = bound_level(P, LL, t1, t2, M1, M2, W, vec=True)
else:
e = bound_level(P, LL, t1, t2, M1, M2, W)
E.append(e - (M1 + M2))
Epp_h = d2_7pt([E[0], E[1], E[2], E[3]], H_FD)
Epp_2h = d2_7pt([E[0], E[2], E[4], E[5]], 2 * H_FD)
res[LL] = dict(Eb=E[0], Epp=Epp_h, Epp2h=Epp_2h, s3=(s1, s2), E03=E[5])
r0, r1 = res[L], res[int(1.5 * L)]
conv = max(abs(r0["Eb"] - r1["Eb"]), abs(r0["Epp"] - r1["Epp"]),
abs(r0["s3"][0] - r1["s3"][0]), abs(r0["s3"][1] - r1["s3"][1]))
pr = predictions(t1, t2, M1, M2, w)
lead03, corr03 = Eb_closed(0.3, t1, t2, M1, M2, w)
row = dict(eps=eps, Eb=r0["Eb"], Epp=r0["Epp"], s3=r0["s3"], pr=pr, conv=conv,
fd=abs(r0["Epp"] - r0["Epp2h"]), E03=r0["E03"], E03p=(lead03, lead03 + corr03))
rows.append(row)
log(" eps=%.3f E_bind=%.12e E''(0)=%.12e <s3_1>-1=%.6e <s3_2>-1=%.6e" %
(eps, row["Eb"], row["Epp"], row["s3"][0] - 1, row["s3"][1] - 1))
log(" |change L->1.5L| max=%.1e |E''(h)-E''(2h)|=%.1e" % (conv, row["fd"]))
summary[name] = rows
# deviations
log(" deviation from leading form (11.14) | from next-order form (14.6)-(14.9); slopes in eps:")
for lab, num, p0, p1 in (
("E_bind", lambda r: r["Eb"], lambda r: r["pr"]["Eb0"], lambda r: r["pr"]["Eb1"]),
("E''(0)", lambda r: r["Epp"], lambda r: r["pr"]["Epp0"], lambda r: r["pr"]["Epp1"]),
("<s3_1>", lambda r: r["s3"][0], lambda r: r["pr"]["s30"][0], lambda r: r["pr"]["s31"][0]),
("<s3_2>", lambda r: r["s3"][1], lambda r: r["pr"]["s30"][1], lambda r: r["pr"]["s31"][1]),
("E_b(0.3)", lambda r: r["E03"], lambda r: r["E03p"][0], lambda r: r["E03p"][1])):
d0 = [num(r) - p0(r) for r in rows]
d1 = [num(r) - p1(r) for r in rows]
log(" %-8s lead: %s slopes %s" % (lab, " ".join("%+.3e" % x for x in d0),
" ".join("%.2f" % s for s in slope(d0, EPS))))
log(" %-8s next: %s slopes %s" % ("", " ".join("%+.3e" % x for x in d1),
" ".join("%.2f" % s for s in slope(d1, EPS))))
# item 3: m_b - (m1 + m2 - 4 m1 m2 E_bind)
log(" item 3: D = m_b - (m1+m2-4 m1 m2 E_bind), numerical vs (14.11), both divided by m_b^(0):")
for r in rows:
eps = r["eps"]
m1, m2 = M1 / (eps * T1) ** 2, M2 / (eps * T2) ** 2
Dn = 1 / r["Epp"] - (m1 + m2 - 4 * m1 * m2 * r["Eb"])
log(" eps=%.3f D_num/m_b0=%+.6e D_pred/m_b0=%+.6e difference=%+.2e" %
(eps, Dn / r["pr"]["mb0"], r["pr"]["D"] / r["pr"]["mb0"], (Dn - r["pr"]["D"]) / r["pr"]["mb0"]))
log(" relative correction of m_b: numerical m_b/m_b^(0)-1 vs (14.10):")
for r in rows:
log(" eps=%.3f num=%+.6e pred=%+.6e" % (r["eps"], (1 / r["Epp"]) / r["pr"]["mb0"] - 1,
r["pr"]["mb1"] / r["pr"]["mb0"] - 1))
# ---------------------------------------------------------------- Part B: effective operator on a ring
log("")
log("Part B: spectrum of the B_++ cluster of H(P) on a ring of N_r = 40 sites vs the effective operator")
log(" (14.4) built from its explicit coefficients; P = 0.37, eta = 0.3, data D1")
NR, PB, EB = 40, 0.37, 0.3
def ring_full(P, t1, t2, M1, M2, W, eta):
S = np.roll(np.eye(NR), 1, axis=0) # (S psi)(r) = psi(r-1)
Sd = S.T
T1 = 0.5j * t1 * (np.exp(1j * eta * P) * Sd - np.exp(-1j * eta * P) * S)
T2 = 0.5j * t2 * (np.exp(1j * (1 - eta) * P) * S - np.exp(-1j * (1 - eta) * P) * Sd)
s1, s3, I2 = np.array([[0, 1], [1, 0]]), np.diag([1.0, -1.0]), np.eye(2)
P0 = np.zeros((NR, NR)); P0[0, 0] = 1
H = (np.kron(T1, np.kron(s1, I2)) + np.kron(T2, np.kron(I2, s1)) + M1 * np.kron(np.eye(NR), np.kron(s3, I2))
+ M2 * np.kron(np.eye(NR), np.kron(I2, s3)) + np.kron(P0, np.diag([W[0][0], W[0][1], W[1][0], W[1][1]])))
return np.linalg.eigvalsh(H)
def ring_eff(P, t1, t2, M1, M2, W, eta, order):
a1, a2 = t1 ** 2 / (4 * M1), t2 ** 2 / (4 * M2)
b1, b2 = a1 / M1, a2 / M2
f1, f2 = np.exp(2j * eta * P), np.exp(-2j * (1 - eta) * P)
w = W[0][0]
H = np.zeros((NR, NR), complex)
if order == 2: # (11.11)
diag, J2, J4, Jb, Om, Lc = a1 + a2, a1 * f1 + a2 * f2, 0, 0, 0, 0
else: # (14.4); order 3: without the w12, w21 terms
diag = a1 * (1 - 0.75 * b1) + a2 * (1 - 0.75 * b2)
J2 = a1 * (1 - b1) * f1 + a2 * (1 - b2) * f2
J4 = a1 * b1 * f1 ** 2 + a2 * b2 * f2 ** 2
Jb = b1 * f1 + b2 * f2
Om = W[1][0] * b1 + W[0][1] * b2 if order == 4 else 0
Lc = W[1][0] * b1 * f1 + W[0][1] * b2 * f2 if order == 4 else 0
for r in range(NR):
H[r, r] += M1 + M2 + diag
H[r, (r + 2) % NR] += -0.5 * J2
H[r, (r - 2) % NR] += -0.5 * np.conj(J2)
H[r, (r + 4) % NR] += -0.125 * J4
H[r, (r - 4) % NR] += -0.125 * np.conj(J4)
beta = b1 + b2 if order > 2 else 0
H[0, 0] += w * (1 - beta / 2)
H[0, 2] += w / 8 * Jb; H[0, NR - 2] += w / 8 * np.conj(Jb)
H[NR - 2, 0] += w / 8 * Jb; H[2, 0] += w / 8 * np.conj(Jb)
H[1, 1] += Om / 4; H[NR - 1, NR - 1] += Om / 4
H[NR - 1, 1] += -0.25 * Lc; H[1, NR - 1] += -0.25 * np.conj(Lc)
assert np.allclose(H, H.conj().T)
return np.linalg.eigvalsh(H)
T1, T2, M1, M2, W11, W12, W21, W22 = DATA["D1"]
devs = {2: [], 3: [], 4: []}
for eps in EPS:
W = [[eps ** 2 * W11, eps ** 2 * W12], [eps ** 2 * W21, eps ** 2 * W22]]
ev = ring_full(PB, eps * T1, eps * T2, M1, M2, W, EB)
clus = np.sort(ev[ev > M1 + M2 - min(M1, M2)])
for o in (2, 3, 4):
devs[o].append(np.max(np.abs(clus - ring_eff(PB, eps * T1, eps * T2, M1, M2, W, EB, o))))
for o, lab in ((2, "(11.11) "), (3, "(14.4) without w12,w21"),
(4, "(14.4) complete ")):
log(" max |eigenvalue difference|, %s: %s slopes %s" %
(lab, " ".join("%.3e" % x for x in devs[o]), " ".join("%.2f" % s for s in slope(devs[o], EPS))))
# ---------------------------------------------------------------- Part C: small-u forms (14.12)-(14.13)
log("")
log("Part C: ordered limit u -> 0 (w11 < 0) of the eps^4 coefficients of (14.7)-(14.10)")
def expansion(a1, a2, x, y):
"""coefficients of u^0..u^n displayed in (14.12), sigma = -1; x = 1/M1, y = 1/M2"""
a, nu = a1 + a2, a1 * a2 / (a1 + a2) ** 2
gs, gc, gm = a1 ** 2 * x + a2 ** 2 * y, a2 ** 2 * x + a1 ** 2 * y, a1 * a2 * (x + y)
def X(ai, aj, xi, xj):
return [0, 0, 0, aj / (2 * a) * ((ai - 2 * aj) * xi + (2 * ai - aj) * xj),
(ai ** 2 - 7 * ai * aj + 4 * aj ** 2) * xi / (8 * a) + 3 * aj * (aj - ai) * xj / (4 * a),
aj / (4 * a) * ((2 * aj - 3 * ai) * xi + (aj - 4 * ai) * xj)]
return dict(E4=[0, 0, 0, gm / 2],
nuC=[0, 0, 2 * nu * (gc - 2 * gm), 2 * nu * (5 * gm - gs - gc), nu * (3 * gs - 9 * gm / 2)],
rel=[0, 0, -(gc - 2 * gm) / (2 * a), -(5 * gm - gs - gc) / (2 * a),
-(6 * gs + 2 * gc - 13 * gm) / (8 * a)],
Xi1=X(a1, a2, x, y), Xi2=X(a2, a1, y, x))
def lines(c):
return dict(E4=c["E4"], nuC=c["nu"] * c["C"], rel=c["rel"], Xi1=c["Xi"][0], Xi2=c["Xi"][1])
# C1: exact series of the closed forms against (14.12), symbolic a1, a2 > 0, x = 1/M1 > 0, y = 1/M2 > 0
log(" C1: exact Taylor coefficients in u of the closed forms (sympy) vs the coefficients of (14.12):")
a1s, a2s, xs, ys, us = sp.symbols("a1 a2 x y u", positive=True)
exact, disp = lines(coef(a1s, a2s, 1 / xs, 1 / ys, us, -1, sp.sqrt)), expansion(a1s, a2s, xs, ys)
for k in exact:
n = len(disp[k]) - 1
ser = sp.expand(sp.series(exact[k], us, 0, n + 1).removeO())
ok = all(sp.simplify(ser.coeff(us, m) - disp[k][m]) == 0 for m in range(n + 1))
log(" %-4s through u^%d: identical = %s" % (k, n, ok))
# C2: leading terms (14.13) on data sets with no, single and double cancellations, 50 digits
mp.mp.dps = 50
ts = mp.findroot(lambda t: 2 * t ** 3 - 9 * t ** 2 - 3 * t + 8, 0.86)
r21 = mp.sqrt(21)
log(" C2: t_* (root of 2t^3-9t^2-3t+8 in (1/2,2)) = %.10f, (2-t_*)/(2t_*-1) = %.10f; 50-digit arithmetic"
% (float(ts), float((2 - ts) / (2 * ts - 1))))
log(" leading power = first displayed coefficient with |c| > 1e-30; ratio = closed form / leading term")
def hat(T1, T2, M1, M2):
T1, T2, M1, M2 = (mp.mpf(str(v)) for v in (T1, T2, M1, M2))
return T1 ** 2 / (4 * M1), T2 ** 2 / (4 * M2), M1, M2
SETS = {"D3": hat(*DATA["D3"][:4]), "equal": hat(1, 1, 1, 1),
"S2": (mp.mpf(1) / 4, mp.mpf(1), mp.mpf(8) / 7, mp.mpf(1)), # gc = 2gm, gs != 3gm
"S3": ((5 - r21) / 2, mp.mpf(1), (11 + r21) / 10, mp.mpf(1)), # gc = 2gm, gs = 3gm
"X3": (ts, mp.mpf(1), (2 - ts) / (2 * ts - 1), mp.mpf(1))} # X3 = X4 = 0 for Xi_1
for lab, (a1, a2, M1, M2) in SETS.items():
disp = expansion(a1, a2, 1 / M1, 1 / M2)
lead, zmax = {}, mp.mpf(0)
for k, v in disp.items():
n = next(m for m, c in enumerate(v) if abs(c) > mp.mpf(10) ** -30)
lead[k] = n
zmax = max([zmax] + [abs(c) for c in v[:n]])
gs, gc, gm = a1 ** 2 / M1 + a2 ** 2 / M2, a2 ** 2 / M1 + a1 ** 2 / M2, a1 * a2 * (1 / M1 + 1 / M2)
log(" %-5s a^=(%.6f,%.6f) M=(%.6f,%.6f) (gc-2gm)/gm=%+.3e (gs-3gm)/gm=%+.3e max|skipped c|=%.1e"
% (lab, float(a1), float(a2), float(M1), float(M2), float((gc - 2 * gm) / gm),
float((gs - 3 * gm) / gm), float(zmax)))
log(" leading powers: %s" % " ".join("%s u^%d (c=%+.6e)" % (k, n, float(disp[k][n]))
for k, n in lead.items()))
for u in (mp.mpf("1e-2"), mp.mpf("1e-3")):
val = lines(coef(a1, a2, M1, M2, u, -1, mp.sqrt))
log(" u=%.0e ratios: %s" % (float(u), " ".join(
"%s %.6f" % (k, float(val[k] / (disp[k][n] * u ** n))) for k, n in lead.items())))
(OUT / "check_pair_level.log").write_text("\n".join(LOG) + "\n", encoding="utf-8")
output/check_pair_level.log10 KB · raw
Part A: bound level of H(P) (11.1), sector Theta=+1, open chain r in [-L, L], eta = 0.50
FD: 7-point stencil for E_b''(0), step h = 0.050; control step 2h = 0.100
D1: tau^=(1.00,1.60) M=(1.00,2.50) w^=(w11 -0.40, w12 -0.50, w21 0.30, w22 0.70) u=0.7905 L=112
eps=0.200 E_bind=-5.537956235082e-03 E''(0)=1.588965060821e-02 <s3_1>-1=-2.168838e-03 <s3_2>-1=-8.828018e-04
|change L->1.5L| max=1.5e-12 |E''(h)-E''(2h)|=8.4e-10
eps=0.100 E_bind=-1.388688933774e-03 E''(0)=3.969821660361e-03 <s3_1>-1=-5.396400e-04 <s3_2>-1=-2.206895e-04
|change L->1.5L| max=5.9e-12 |E''(h)-E''(2h)|=2.0e-10
eps=0.050 E_bind=-3.474340889640e-04 E''(0)=9.922902739865e-04 <s3_1>-1=-1.347498e-04 <s3_2>-1=-5.517182e-05
|change L->1.5L| max=2.8e-12 |E''(h)-E''(2h)|=5.1e-11
deviation from leading form (11.14) | from next-order form (14.6)-(14.9); slopes in eps:
E_bind lead: +2.238e-05 +1.396e-06 +8.723e-08 slopes 4.00 4.00
next: +5.791e-08 +9.031e-10 +1.410e-11 slopes 6.00 6.00
E''(0) lead: +1.389e-05 +8.819e-07 +5.533e-08 slopes 3.98 3.99
next: -2.914e-07 -4.547e-09 -7.229e-11 slopes 6.00 5.98
<s3_1> lead: -1.370e-05 -8.544e-07 -5.338e-08 slopes 4.00 4.00
next: -3.410e-08 -5.349e-10 -8.365e-12 slopes 5.99 6.00
<s3_2> lead: -5.545e-08 -2.872e-09 -1.703e-10 slopes 4.27 4.08
next: -1.263e-08 -1.953e-10 -3.043e-12 slopes 6.01 6.00
E_b(0.3) lead: +2.289e-05 +1.428e-06 +8.924e-08 slopes 4.00 4.00
next: +4.492e-08 +7.008e-10 +1.095e-11 slopes 6.00 6.00
item 3: D = m_b - (m1+m2-4 m1 m2 E_bind), numerical vs (14.11), both divided by m_b^(0):
eps=0.200 D_num/m_b0=-6.629865e-06 D_pred/m_b0=-2.799318e-05 difference=+2.14e-05
eps=0.100 D_num/m_b0=-5.663204e-06 D_pred/m_b0=-6.998295e-06 difference=+1.34e-06
eps=0.050 D_num/m_b0=-1.664859e-06 D_pred/m_b0=-1.749574e-06 difference=+8.47e-08
relative correction of m_b: numerical m_b/m_b^(0)-1 vs (14.10):
eps=0.200 num=-8.742482e-04 pred=-8.933671e-04
eps=0.100 num=-2.221467e-04 pred=-2.233418e-04
eps=0.050 num=-5.575948e-05 pred=-5.583544e-05
D2: tau^=(1.30,0.80) M=(1.50,0.70) w^=(w11 0.60, w12 0.40, w21 -0.20, w22 -0.50) u=1.1759 L=80
eps=0.200 E_bind=5.136989357678e-02 E''(0)=-1.275968735072e-02 <s3_1>-1=-1.199641e-02 <s3_2>-1=-2.033606e-02
|change L->1.5L| max=6.5e-13 |E''(h)-E''(2h)|=8.0e-10
eps=0.100 E_bind=1.294382462252e-02 E''(0)=-3.249295987270e-03 <s3_1>-1=-3.069918e-03 <s3_2>-1=-5.303397e-03
|change L->1.5L| max=4.2e-13 |E''(h)-E''(2h)|=2.0e-10
eps=0.050 E_bind=3.242458562345e-03 E''(0)=-8.161338742853e-04 <s3_1>-1=-7.720324e-04 <s3_2>-1=-1.340283e-03
|change L->1.5L| max=7.5e-13 |E''(h)-E''(2h)|=5.0e-11
deviation from leading form (11.14) | from next-order form (14.6)-(14.9); slopes in eps:
E_bind lead: -5.444e-04 -3.474e-05 -2.183e-06 slopes 3.97 3.99
next: +1.546e-05 +2.482e-07 +3.905e-09 slopes 5.96 5.99
E''(0) lead: +3.189e-04 +2.036e-05 +1.279e-06 slopes 3.97 3.99
next: -9.094e-06 -1.453e-07 -2.283e-09 slopes 5.97 5.99
<s3_1> lead: +3.806e-04 +2.433e-05 +1.529e-06 slopes 3.97 3.99
next: -1.168e-05 -1.869e-07 -2.938e-09 slopes 5.97 5.99
<s3_2> lead: +1.187e-03 +7.725e-05 +4.879e-06 slopes 3.94 3.98
next: -6.696e-05 -1.091e-06 -1.724e-08 slopes 5.94 5.98
E_b(0.3) lead: -5.304e-04 -3.385e-05 -2.127e-06 slopes 3.97 3.99
next: +1.507e-05 +2.419e-07 +3.807e-09 slopes 5.96 5.99
item 3: D = m_b - (m1+m2-4 m1 m2 E_bind), numerical vs (14.11), both divided by m_b^(0):
eps=0.200 D_num/m_b0=+4.227327e-02 D_pred/m_b0=+4.284997e-02 difference=-5.77e-04
eps=0.100 D_num/m_b0=+1.067556e-02 D_pred/m_b0=+1.071249e-02 difference=-3.69e-05
eps=0.050 D_num/m_b0=+2.675800e-03 D_pred/m_b0=+2.678123e-03 difference=-2.32e-06
relative correction of m_b: numerical m_b/m_b^(0)-1 vs (14.10):
eps=0.200 num=+2.499426e-02 pred=+2.508015e-02
eps=0.100 num=+6.264615e-03 pred=+6.270038e-03
eps=0.050 num=+1.567169e-03 pred=+1.567510e-03
D3: tau^=(1.20,0.90) M=(0.80,1.30) w^=(w11 -0.12, w12 0.20, w21 0.10, w22 0.30) u=0.1981 L=408
eps=0.200 E_bind=-4.699931522243e-04 E''(0)=1.816024850353e-02 <s3_1>-1=-4.283860e-04 <s3_2>-1=-9.169974e-05
|change L->1.5L| max=1.6e-12 |E''(h)-E''(2h)|=2.0e-09
eps=0.100 E_bind=-1.176597100687e-04 E''(0)=4.540273684594e-03 <s3_1>-1=-1.071905e-04 <s3_2>-1=-2.286127e-05
|change L->1.5L| max=5.3e-13 |E''(h)-E''(2h)|=4.8e-10
eps=0.050 E_bind=-2.942500201675e-05 E''(0)=1.135081467153e-03 <s3_1>-1=-2.680349e-05 <s3_2>-1=-5.711356e-06
|change L->1.5L| max=4.9e-13 |E''(h)-E''(2h)|=1.2e-10
deviation from leading form (11.14) | from next-order form (14.6)-(14.9); slopes in eps:
E_bind lead: +8.606e-07 +5.373e-08 +3.357e-09 slopes 4.00 4.00
next: +1.308e-09 +2.039e-11 +3.179e-13 slopes 6.00 6.00
E''(0) lead: -1.124e-06 -6.953e-08 -4.335e-09 slopes 4.02 4.00
next: -1.606e-08 -2.582e-10 -6.230e-12 slopes 5.96 5.37
<s3_1> lead: +5.012e-07 +3.132e-08 +1.957e-09 slopes 4.00 4.00
next: +9.991e-11 +1.581e-12 +2.442e-14 slopes 5.98 6.02
<s3_2> lead: -3.392e-07 -2.113e-08 -1.319e-09 slopes 4.00 4.00
next: -1.538e-09 -2.397e-11 -3.744e-13 slopes 6.00 6.00
E_b(0.3) lead: +6.392e-07 +3.991e-08 +2.494e-09 slopes 4.00 4.00
next: +7.306e-10 +1.139e-11 +1.786e-13 slopes 6.00 5.99
item 3: D = m_b - (m1+m2-4 m1 m2 E_bind), numerical vs (14.11), both divided by m_b^(0):
eps=0.200 D_num/m_b0=+9.675148e-05 D_pred/m_b0=+9.581065e-05 difference=+9.41e-07
eps=0.100 D_num/m_b0=+2.401307e-05 D_pred/m_b0=+2.395266e-05 difference=+6.04e-08
eps=0.050 D_num/m_b0=+5.993875e-06 D_pred/m_b0=+5.988166e-06 difference=+5.71e-09
relative correction of m_b: numerical m_b/m_b^(0)-1 vs (14.10):
eps=0.200 num=+6.191205e-05 pred=+6.102416e-05
eps=0.100 num=+1.531314e-05 pred=+1.525604e-05
eps=0.050 num=+3.819513e-06 pred=+3.814010e-06
Part B: spectrum of the B_++ cluster of H(P) on a ring of N_r = 40 sites vs the effective operator
(14.4) built from its explicit coefficients; P = 0.37, eta = 0.3, data D1
max |eigenvalue difference|, (11.11) : 2.618e-04 1.657e-05 1.039e-06 slopes 3.98 4.00
max |eigenvalue difference|, (14.4) without w12,w21: 6.825e-06 2.367e-07 1.231e-08 slopes 4.85 4.26
max |eigenvalue difference|, (14.4) complete : 4.233e-06 6.754e-08 1.060e-09 slopes 5.97 5.99
Part C: ordered limit u -> 0 (w11 < 0) of the eps^4 coefficients of (14.7)-(14.10)
C1: exact Taylor coefficients in u of the closed forms (sympy) vs the coefficients of (14.12):
E4 through u^3: identical = True
nuC through u^4: identical = True
rel through u^4: identical = True
Xi1 through u^5: identical = True
Xi2 through u^5: identical = True
C2: t_* (root of 2t^3-9t^2-3t+8 in (1/2,2)) = 0.8625413912, (2-t_*)/(2t_*-1) = 1.5687293044; 50-digit arithmetic
leading power = first displayed coefficient with |c| > 1e-30; ratio = closed form / leading term
D3 a^=(0.450000,0.155769) M=(0.800000,1.300000) (gc-2gm)/gm=-6.852e-01 (gs-3gm)/gm=-1.080e+00 max|skipped c|=0.0e+00
leading powers: E4 u^3 (c=+7.077016e-02) nuC u^2 (c=-3.705082e-02) rel u^2 (c=+8.004808e-02) Xi1 u^3 (c=+9.585799e-02) Xi2 u^3 (c=-2.769231e-01)
u=1e-02 ratios: E4 0.999750 nuC 0.974152 rel 0.974201 Xi1 0.990132 Xi2 0.990624
u=1e-03 ratios: E4 0.999980 nuC 0.997423 rel 0.997424 Xi1 0.999029 Xi2 0.999063
equal a^=(0.250000,0.250000) M=(1.000000,1.000000) (gc-2gm)/gm=-1.000e+00 (gs-3gm)/gm=-2.000e+00 max|skipped c|=0.0e+00
leading powers: E4 u^3 (c=+6.250000e-02) nuC u^2 (c=-6.250000e-02) rel u^2 (c=+1.250000e-01) Xi1 u^4 (c=-3.125000e-02) Xi2 u^4 (c=-3.125000e-02)
u=1e-02 ratios: E4 0.997450 nuC 0.970077 rel 0.970126 Xi1 1.039844 Xi2 1.039844
u=1e-03 ratios: E4 0.999750 nuC 0.997001 rel 0.997001 Xi1 1.003998 Xi2 1.003998
S2 a^=(0.250000,1.000000) M=(1.142857,1.000000) (gc-2gm)/gm=+0.000e+00 (gs-3gm)/gm=-7.500e-01 max|skipped c|=0.0e+00
leading powers: E4 u^3 (c=+2.343750e-01) nuC u^3 (c=+1.125000e-01) rel u^3 (c=-1.406250e-01) Xi1 u^3 (c=-8.125000e-01) Xi2 u^3 (c=+2.031250e-01)
u=1e-02 ratios: E4 1.000575 nuC 1.014814 rel 1.014865 Xi1 0.991945 Xi2 0.992550
u=1e-03 ratios: E4 1.000062 nuC 1.001498 rel 1.001499 Xi1 0.999197 Xi2 0.999268
S3 a^=(0.208712,1.000000) M=(1.558258,1.000000) (gc-2gm)/gm=+7.800e-51 (gs-3gm)/gm=+1.560e-50 max|skipped c|=2.2e-51
leading powers: E4 u^3 (c=+1.713258e-01) nuC u^4 (c=+2.202760e-01) rel u^4 (c=-3.189205e-01) Xi1 u^3 (c=-7.165151e-01) Xi2 u^3 (c=+1.495454e-01)
u=1e-02 ratios: E4 1.002450 nuC 0.995417 rel 0.995467 Xi1 0.990726 Xi2 0.993474
u=1e-03 ratios: E4 1.000249 nuC 0.999554 rel 0.999555 Xi1 0.999075 Xi2 0.999361
X3 a^=(0.862541,1.000000) M=(1.568729,1.000000) (gc-2gm)/gm=-1.022e+00 (gs-3gm)/gm=-1.956e+00 max|skipped c|=8.4e-53
leading powers: E4 u^3 (c=+7.061879e-01) nuC u^2 (c=-7.177271e-01) rel u^2 (c=+3.874586e-01) Xi1 u^5 (c=-3.791529e-01) Xi2 u^4 (c=-1.688121e-01)
u=1e-02 ratios: E4 0.997560 nuC 0.970927 rel 0.970975 Xi1 0.998819 Xi2 1.022334
u=1e-03 ratios: E4 0.999760 nuC 0.997086 rel 0.997087 Xi1 0.999895 Xi2 1.002245