validate_strongfield.py — src/21-gravity

1"""Strong field: exact profile, critical radii, photon sphere, shadow, ringdown, entropy, echoes.23Validates: Prop (exact static profile)         existence boundary r_* = sqrt(Gm);4           Prop (horizonless, operationally black)  r_f = 2Gm/ln(Omega);5           Prop (photon sphere and shadow)     r_ph = 2Gm, b_c = 2e Gm, delta = +4.63%;6           ringdown shift -4.4%; S/S_BH = ln^-2(Omega); echo delay ~ Omega (no echoes);7           Sgr A* / M87* angular numbers;8           Prop (ISCO and accretion efficiency)  r_ISCO = (3+sqrt5) Gm,9               efficiency 5.48% vs Schwarzschild 5.72%, Gm*Omega ratio 0.931.10"""11import math12import numpy as np1314ok = True15Gm = 1.0e6                       # Planck units; corrections (Gm/r^2)^2 negligible at r ~ Gm16lnOm = math.log(1e122)1718# exact potential u(r) = int_r^inf arcsin(Gm/s^2) ds  (numeric, log grid to 1e12)19def u_exact(r):20    s = np.logspace(np.log10(r), 12, 400000)21    return np.trapezoid(np.arcsin(np.clip(Gm/s**2, 0, 1)), s)2223rstar = math.sqrt(Gm)24print(f"  slip core r_* = sqrt(Gm) = {rstar:.0f}  (arcsin domain boundary)")25print(f"  u(2Gm) = {u_exact(2*Gm):.6f}  vs Gm/r = 0.5  (exact profile ~ Newtonian at photon sphere)")26if abs(u_exact(2*Gm) - 0.5) > 1e-3: ok = False2728# photon sphere: extremise b(r) = r exp(2 u(r)) with the exact potential29rr = np.linspace(1.6*Gm, 2.4*Gm, 81)30bb = [r*math.exp(2*u_exact(r)) for r in rr]31r_ph = rr[int(np.argmin(bb))]32b_c = min(bb)33print(f"  photon sphere r_ph = {r_ph/Gm:.3f} Gm (target 2),  b_c = {b_c/Gm:.4f} Gm (target 2e = {2*math.e:.4f})")34if abs(r_ph/Gm - 2) > 0.05 or abs(b_c/Gm - 2*math.e) > 0.02: ok = False3536delta = 2*math.e/(3*math.sqrt(3)) - 137ring = 3*math.sqrt(3)/(2*math.e) - 138print(f"  shadow deviation delta = {delta*100:.2f}%  (target +4.63), ringdown {ring*100:.2f}% (target -4.42)")39if abs(delta - 0.04627) > 1e-4: ok = False4041# operational horizon, entropy ratio, echo delay42r_f = 2*Gm/lnOm43print(f"  r_f = 2Gm/ln(Omega) = r_s/{2*Gm/r_f:.0f},  S/S_BH = {1/lnOm**2:.2e},  "44      f"echo delay factor Omega/ln^2(Omega) = {1e122/lnOm**2:.1e} (no observable echoes)")4546# Sgr A* and M87* angular sizes47Msun_km = 1.4766248for name, M, D_kpc in (("Sgr A*", 4.297e6, 8.277), ("M87*", 6.5e9, 16.8e3)):49    theta_g = M*Msun_km*1e3/(D_kpc*3.0857e19)*206264.806e650    gr, frc = 2*3*math.sqrt(3)*theta_g, 2*2*math.e*theta_g51    print(f"  {name}: shadow GR {gr:.1f} muas, FRC {frc:.1f} muas")52    if name == "Sgr A*" and not (52 < gr < 55 and 54 < frc < 57): ok = False5354# --- ISCO and accretion efficiency (exponential metric, dimensionless M = Gm/c^2 = 1) ---55# circular orbits: rdot^2 = E^2 - W,  W = e^{-2/r} + L^2 e^{-4/r}/r^2 ; ISCO: W' = W'' = 056def _W(r, L2):  return math.exp(-2/r) + L2*math.exp(-4/r)/r**257def _L2(r):58    h = 1e-459    dA = (math.exp(-2/(r+h)) - math.exp(-2/(r-h)))/(2*h)60    dB = (math.exp(-4/(r+h))/(r+h)**2 - math.exp(-4/(r-h))/(r-h)**2)/(2*h)61    return -dA/dB62def _Wrr(r, L2):63    h = 1e-464    return (_W(r+h, L2) - 2*_W(r, L2) + _W(r-h, L2))/h**265r_isco = 3 + math.sqrt(5)                       # exact closed form to verify66L2 = _L2(r_isco); E = math.sqrt(_W(r_isco, L2)); L = math.sqrt(L2)67MOm = L*math.exp(-4/r_isco)/(E*r_isco**2)       # Gm*Omega (M=1)68eff = (1 - E)*100; ratio = MOm*6*math.sqrt(6)69print(f"  ISCO r = (3+sqrt5) Gm = {r_isco:.4f} Gm,  |W''| = {abs(_Wrr(r_isco, L2)):.1e} (marginally stable)")70print(f"  E/m = {E:.5f},  efficiency = {eff:.3f}% (Schwarzschild {(1-math.sqrt(8/9))*100:.3f}%),  "71      f"Gm*Omega = {MOm:.5f}, ratio to Schwarzschild = {ratio:.4f}")72if abs(eff - 5.479) > 0.02 or abs(ratio - 0.9308) > 0.005 or abs(_Wrr(r_isco, L2)) > 1e-3: ok = False7374print("PASS" if ok else "FAIL", "- strong field")75