#!/usr/bin/env python3 """ verify_P138.py — Verification script for Addendum P138 Pion FSI and rho-omega Interference Corrections to Delta_alpha_had Verified claims: 1. FSI Omnes correction: +2.182e-4 2. rho-omega cross-term: -2.316e-4 3. Net correction: -1.34e-5 4. P138 total Delta_alpha_had: 3.718e-3 5. Residual gap: -9.55% 6. W-boson mass: 80.383 GeV Author: verification script for León Fernando Vlegels, P138 Copyright: León Fernando Vlegels. License: MIT. """ import math import sys # --------------------------------------------------------------------------- # TOE constants # --------------------------------------------------------------------------- ALPHA_INV = 4 * math.pi**3 + math.pi**2 + math.pi ALPHA = 1.0 / ALPHA_INV PHI = (1.0 + math.sqrt(5.0)) / 2.0 # --------------------------------------------------------------------------- # Physical constants (PDG 2022) # --------------------------------------------------------------------------- M_Z = 91.1876 # GeV M_PI = 0.13957 # GeV, charged pion mass M_RHO = 0.77526 # GeV, rho(770) mass G_RHO = 0.1478 # GeV, rho total width M_OME = 0.78265 # GeV, omega mass G_OME = 0.00849 # GeV, omega total width G_F = 1.1663787e-5 # GeV^-2, Fermi constant ALPHA_S_MZ = 0.1186 # strong coupling at M_Z (P106) # Derived S_THRESH = 4.0 * M_PI**2 # 4m_pi^2 (GeV^2) M_RHO2 = M_RHO**2 M_OME2 = M_OME**2 M_Z2 = M_Z**2 # --------------------------------------------------------------------------- # P118 rho-omega mixing amplitude # --------------------------------------------------------------------------- EPS = -ALPHA_S_MZ / 3.0 # = -0.03953... # --------------------------------------------------------------------------- # P137 baseline and targets (from P138 text) # --------------------------------------------------------------------------- DA_HAD_P137 = 3.731e-3 # GS baseline DA_HAD_FERMI = 4.110e-3 # Fermi target (P135) DA_LEP = 0.031422 # leptonic running (P134) PDG_MW = 80.377 # GeV # --------------------------------------------------------------------------- # Helper functions # --------------------------------------------------------------------------- def beta(s): """Pion velocity beta(s) = sqrt(1 - 4m_pi^2/s).""" arg = 1.0 - S_THRESH / s if arg <= 0.0: return 0.0 return math.sqrt(arg) def bw_ff_sq(s): """|F_pi|^2_BW = M_rho^4 / [(s-M_rho^2)^2 + M_rho^2 Gamma_rho^2].""" denom = (s - M_RHO2)**2 + M_RHO2 * G_RHO**2 return M_RHO2**2 / denom def omnes_factor(s): """ Omnes ratio |F_pi|^2_Om / |F_pi|^2_BW = s^(-s/M_rho^2) using linear phase-shift approximation delta_1^1(s) ~ (s/M_rho^2)*(pi/2) and Lambda = 1 GeV cutoff. = exp[(s/M_rho^2) * ln(1/s)] with s in GeV^2. """ exp_arg = (s / M_RHO2) * math.log(1.0 / s) return math.exp(exp_arg) def re_b_omega(s): """Re[B_omega(s)] = M_omega^2 (M_omega^2 - s) / [(M_omega^2-s)^2 + (M_omega Gamma_omega)^2].""" num = M_OME2 * (M_OME2 - s) denom = (M_OME2 - s)**2 + (M_OME * G_OME)**2 return num / denom def dispersive_kernel(s): """Common dispersive weight: M_Z^2 / (M_Z^2 - s).""" return M_Z2 / (M_Z2 - s) # --------------------------------------------------------------------------- # Gaussian-Legendre quadrature (pure Python, no scipy dependency) # --------------------------------------------------------------------------- def gauss_legendre_nodes_weights(n): """ Compute n-point Gauss-Legendre nodes and weights on [-1, 1]. Uses Newton iteration on Legendre polynomial recurrence. """ nodes = [0.0] * n weights = [0.0] * n for i in range((n + 1) // 2): # Initial guess xi = math.cos(math.pi * (i + 0.75) / (n + 0.5)) for _ in range(100): p0, p1 = 1.0, xi for k in range(2, n + 1): p0, p1 = p1, ((2*k - 1)*xi*p1 - (k-1)*p0) / k dp = n * (p0 - xi * p1) / (1.0 - xi**2) dx = p1 / dp xi -= dx if abs(dx) < 1e-15: break w = 2.0 / ((1.0 - xi**2) * dp**2) nodes[i] = -xi nodes[n - 1 - i] = xi weights[i] = w weights[n - 1 - i] = w return nodes, weights def integrate_gl(f, a, b, n): """Integrate f on [a,b] using n-point Gauss-Legendre.""" nodes, weights = gauss_legendre_nodes_weights(n) mid = 0.5 * (a + b) half = 0.5 * (b - a) return half * sum(w * f(mid + half * x) for x, w in zip(nodes, weights)) def integrate_subintervals(f, breakpoints, n_per_sub): """Integrate f over subintervals defined by a list of breakpoints.""" total = 0.0 for k in range(len(breakpoints) - 1): total += integrate_gl(f, breakpoints[k], breakpoints[k+1], n_per_sub) return total # --------------------------------------------------------------------------- # Claim 1: FSI Omnes correction # --------------------------------------------------------------------------- def da_fsi_integrand(s): b = beta(s) ff2 = bw_ff_sq(s) om = omnes_factor(s) dR = 0.25 * b**3 * ff2 * (om - 1.0) return dR / s * dispersive_kernel(s) # Integration range: 4m_pi^2 to (0.6 GeV)^2 = 0.36 GeV^2 S_HI_FSI = 0.36 # GeV^2 # Five sub-intervals matching the paper's Table 1 (sqrt(s) boundaries in GeV) SQRT_S_BREAKS_FSI = [math.sqrt(S_THRESH), 0.346, 0.400, 0.458, 0.520, 0.600] S_BREAKS_FSI = [x**2 for x in SQRT_S_BREAKS_FSI] DA_FSI = (ALPHA / (3.0 * math.pi)) * integrate_subintervals( da_fsi_integrand, S_BREAKS_FSI, n_per_sub=400 ) # --------------------------------------------------------------------------- # Claim 2: rho-omega cross-term # --------------------------------------------------------------------------- def da_rw_integrand(s): b = beta(s) ff2 = bw_ff_sq(s) reBw = re_b_omega(s) dR = 0.25 * b**3 * ff2 * 2.0 * EPS * reBw return dR / s * dispersive_kernel(s) # Integration range: 4m_pi^2 to 1 GeV^2, with fine grid around omega pole S_HI_RW = 1.0 # GeV^2 DELTA_OMEGA = 5.0 * M_OME * G_OME # half-window around omega pole S_BREAKS_RW = [ S_THRESH, max(S_THRESH + 1e-6, M_OME2 - DELTA_OMEGA), M_OME2 - 0.5 * M_OME * G_OME, M_OME2, M_OME2 + 0.5 * M_OME * G_OME, M_OME2 + DELTA_OMEGA, S_HI_RW, ] # Remove any breakpoints below threshold or above upper limit, keep sorted unique S_BREAKS_RW = sorted(set( max(S_THRESH + 1e-9, min(S_HI_RW - 1e-9, bp)) for bp in S_BREAKS_RW if S_THRESH < bp < S_HI_RW )) S_BREAKS_RW = [S_THRESH] + S_BREAKS_RW + [S_HI_RW] DA_RW = (ALPHA / (3.0 * math.pi)) * integrate_subintervals( da_rw_integrand, S_BREAKS_RW, n_per_sub=600 ) # --------------------------------------------------------------------------- # Derived claims # --------------------------------------------------------------------------- DA_NET = DA_FSI + DA_RW DA_P138 = DA_HAD_P137 + DA_NET GAP_ABS = DA_P138 - DA_HAD_FERMI GAP_PCT = GAP_ABS / DA_HAD_FERMI * 100.0 # W-boson mass from Fermi on-shell constraint DA_TOTAL_P138 = DA_LEP + DA_P138 ALPHA_MZ = ALPHA / (1.0 - DA_TOTAL_P138) C = math.pi * ALPHA_MZ / (math.sqrt(2.0) * G_F * M_Z2) x = 0.5 * (1.0 + math.sqrt(1.0 - 4.0 * C)) # larger root gives M_W > M_Z*sin(theta_W) MW_P138 = M_Z * math.sqrt(x) MW_RESIDUAL_PCT = (MW_P138 - PDG_MW) / PDG_MW * 100.0 # --------------------------------------------------------------------------- # Report # --------------------------------------------------------------------------- PASS = FAIL = 0 def check(n, desc, computed, claimed, tol_rel, unit=""): global PASS, FAIL rel_err = abs(computed - claimed) / abs(claimed) if claimed != 0 else abs(computed) ok = rel_err <= tol_rel PASS += ok; FAIL += (not ok) print(f" [{'PASS' if ok else 'FAIL'}] {n:>2}. {desc}") print(f" computed={computed:>12.5g} claimed={claimed:>12.5g} rel err={rel_err * 100:>7.3f}%") return ok print(" P138 Verification: Pion FSI + rho-omega Corrections to Delta_alpha_had") print() # Core constants print(f" ALPHA_INV = {ALPHA_INV:.10f}") print(f" ALPHA = {ALPHA:.10e}") print(f" epsilon = {EPS:.6f} (claimed -0.03953)") print(f" S_thresh = {S_THRESH:.6f} GeV^2 (4m_pi^2)") print() # Sub-interval FSI breakdown print(" FSI sub-interval contributions (x1e-4):") fsi_subs = [] for k in range(len(S_BREAKS_FSI) - 1): sub = (ALPHA / (3.0 * math.pi)) * integrate_gl( da_fsi_integrand, S_BREAKS_FSI[k], S_BREAKS_FSI[k+1], 400 ) claimed_str = [0.052, 0.151, 0.304, 0.520, 1.155][k] fsi_subs.append(sub) lo = math.sqrt(S_BREAKS_FSI[k]) hi = math.sqrt(S_BREAKS_FSI[k+1]) print(f" [{lo:.3f}, {hi:.3f}] GeV: computed={sub*1e4:.3f} claimed={claimed_str:.3f} (x1e-4)") print() check(1, "epsilon (rho-omega mixing)", EPS, -0.03953, 5e-3) check(2, "FSI Omnes correction (Da_FSI)", DA_FSI, +2.182e-4, 2e-2) check(3, "rho-omega cross-term (Da_rw)", DA_RW, -2.316e-4, 2e-2) check(4, "Net correction (Da_FSI+Da_rw)", DA_NET, -1.34e-5, 0.20) # looser: near-cancellation check(5, "P138 total (Da_had_P138)", DA_P138, 3.718e-3, 2e-3) check(6, "Residual gap fraction (%)", GAP_PCT, -9.55, 2e-2) check(7, "W-boson mass MW (GeV)", MW_P138, 80.383, 5e-4) print() print(f" DA_FSI = {DA_FSI:.4e} (claimed +2.182e-4)") print(f" DA_RW = {DA_RW:.4e} (claimed -2.316e-4)") print(f" Net = {DA_NET:.4e} (claimed -1.34e-5)") print(f" DA_had(P138) = {DA_P138:.6e} (claimed 3.718e-3)") print(f" Gap = {GAP_PCT:.4f}% (claimed -9.55%)") print(f" MW = {MW_P138:.6f} GeV (claimed 80.383 GeV)") print(f"\n{'='*60}\nRESULT: {PASS} PASS / {FAIL} FAIL") sys.exit(0 if FAIL == 0 else 1)