The script src/21-gravity/validate_2pn.py of the validation package of the paper 21-grav, rendered with an anchor per line. ▶ run the predicates on Google Colab.
1#!/usr/bin/env python32"""2PN inputs and the J0737-3039 confrontation (sub:2pn).34Exact (sympy, rational series): the isotropic-gauge deviations of the5exponential metric from Schwarzschild -- delta A = U^3/6 in -g_00,6delta B = U^2/2 in g_ij -- certified term by term.7[chart] arithmetic: the periastron excess per orbit8 delta(Domega) = pi (GM/c^2 p)^2 (2 + e^2/2)9evaluated for PSR J0737-3039 (M = 2.5871 Msun, Pb = 8834.53 s, e = 0.0878),10giving the fractional excess ~1.5e-6 relative to the 1PN advance: at current11timing precision, presently degenerate with the omega-dot -> M_tot inference.12"""13import sympy as sp14ok = True15def chk(n, c):16 global ok; ok = ok and bool(c); print(('PASS ' if c else 'FAIL ') + n)1718U = sp.symbols('U', positive=True)19exp_g00 = sp.series(sp.exp(-2*U), U, 0, 4).removeO()20schw_g00 = sp.series(((1 - U/2)/(1 + U/2))**2, U, 0, 4).removeO()21dA = sp.expand(exp_g00 - schw_g00)22chk("delta A = U^3/6 exactly (temporal, first difference at 3PN order in g_00)",23 sp.simplify(dA - U**3/6) == 0)24exp_gij = sp.series(sp.exp(2*U), U, 0, 3).removeO()25schw_gij = sp.series((1 + U/2)**4, U, 0, 3).removeO()26dB = sp.expand(exp_gij - schw_gij)27chk("delta B = U^2/2 exactly (spatial, first difference at 2PN order in g_ij)",28 sp.simplify(dB - U**2/2) == 0)29chk("1PN agreement exact: g_00 and g_ij coincide through O(U^2) and O(U) resp.",30 sp.series(exp_g00 - schw_g00, U, 0, 3).removeO() == 0 and31 sp.series(exp_gij - schw_gij, U, 0, 2).removeO() == 0)3233# ---- [chart] evaluation for J0737-3039 ----34G = 6.674e-11; c = 2.998e8; Msun = 1.989e3035M = 2.5871 * Msun; Pb = 8834.535; ecc = 0.087836import math37a = (G * M * Pb**2 / (4 * math.pi**2)) ** (1.0/3.0)38p_slr = a * (1 - ecc**2)39Up = G * M / (c**2 * p_slr)40excess_frac = Up * (2 + ecc**2/2) / 6.0 # delta(Domega)/Domega_1PN41print(f" [chart] U_p = {Up:.3e}; fractional 2PN excess = {excess_frac:.2e}")42chk("fractional excess ~1.5e-6 (within [1.0e-6, 2.0e-6])", 1.0e-6 < excess_frac < 2.0e-6)43print(" verdict: at current 1e-6 omega-dot precision, degenerate with the mass refit;")44print(" decisive once M_tot is fixed independently at 1e-6 (prediction P6).")45464748# ---- round-02: orbit-equation extraction (the analytic derivation's inputs) ----49w, E, L, M = sp.symbols('w E L M', positive=True)50BA_exp = sp.series(sp.exp(4*M*w), M, 0, 4).removeO()51B_exp = sp.series(sp.exp(2*M*w), M, 0, 4).removeO()52lnA_s = -2*M*w - (M*w)**3/653lnB_s = 2*M*w - (M*w)**2/2 + (M*w)**3/654BA_s = sp.series(sp.exp(lnB_s - lnA_s), M, 0, 4).removeO()55B_s = sp.series(sp.exp(lnB_s), M, 0, 4).removeO()56diff_eq = sp.expand((E**2*(BA_exp - BA_s) - (B_exp - B_s))/L**2)57c2 = diff_eq.coeff(w, 2).subs(E, 1 + sp.symbols('eps'))58c2 = sp.series(c2, sp.symbols('eps'), 0, 2).removeO()59c3 = diff_eq.coeff(w, 3).subs(E, 1)60chk("delta c2 leading term = eps M^2/L^2 exactly",61 sp.expand(c2).coeff(sp.symbols('eps'), 1) == M**2/L**2 and sp.expand(c2).coeff(sp.symbols('eps'), 0) == 0)62chk("delta c3 = (5/6) M^3/L^2 exactly", sp.simplify(c3 - sp.Rational(5,6)*M**3/L**2) == 0)63chk("apsidal sum identity: 5/2 - (1-e^2)/2 = 2 + e^2/2", True if sp.simplify(sp.Rational(5,2)-(1-sp.symbols('e')**2)/2-(2+sp.symbols('e')**2/2))==0 else False)6465print("\nALL PASS" if ok else "\nFAILURES PRESENT")66raise SystemExit(0 if ok else 1)67