#!/usr/bin/env python3 """ Voorspelling #18: Niet-Lineaire Lichtvertraging nabij Extreme Massa ===================================================================== Spectrum Grid Model — Pulsar-Timing Analyse Hypothese: Als c een eigenschap is van grid-dichtheid, dan verdicht massa het grid lokaal, waardoor de propagatiesnelheid afneemt. Standaard ART voorspelt Shapiro-vertraging: Δt_ART = -(2GM/c³) × ln(r₁r₂/r₀²) Het grid-model voorspelt een EXTRA niet-lineaire term: Δt_grid = Δt_ART × [1 + β × (GM/rc²)^α] waar α > 1 (niet-lineaire afhankelijkheid) en β is een vrije coëfficiënt die de grid-dichteidsrespons beschrijft. In het zonnestelsel is GM/rc² ~ 10⁻⁶ (zeer zwak) → Δt_grid ≈ Δt_ART. Nabij neutronensterren is GM/rc² ~ 0.2 → de extra term wordt meetbaar. Methode: 1. Bereken standaard Shapiro-vertraging voor binaire pulsars 2. Voeg grid-correctieterm toe 3. Analyseer residuelen als functie van compactheid GM/rc² 4. Vergelijk met NANOGrav en NICER data Databron: - NANOGrav 15-year data release (nanograv.org) - NICER neutronenster-radiusmetingen - Specifiek: PSR J0737-3039 (dubbele pulsar, sterkste Shapiro-meting) - PSR J1614-2230 (2 M☉ neutronenster) Auteur: Marald Bes / Spectrum van Alles Datum: 2026-04-08 """ import numpy as np import matplotlib.pyplot as plt from scipy.optimize import curve_fit import os # ─── Natuurconstanten ─────────────────────────────────────── G = 6.674e-11 # m³/kg/s² c = 3e8 # m/s M_sun = 1.989e30 # kg R_sun = 6.96e8 # m # ─── Bekende Binaire Pulsars ──────────────────────────────── PULSARS = { 'J0737-3039A': { 'name': 'Double Pulsar A', 'M_companion': 1.249 * M_sun, # companion mass 'a_projected': 1.415, # projected semi-major axis (lt-s) 'inclination': 88.69, # graden 'compactness': 0.18, # GM/rc² van companion 'shapiro_measured': 6.21e-6, # gemeten Shapiro delay range (s) 'shapiro_error': 0.33e-6, # 1σ fout (s) 'notes': 'Sterkste Shapiro-test. Kramér et al. 2021, PRX.' }, 'J1614-2230': { 'name': 'Massive NS', 'M_companion': 0.5 * M_sun, # WD companion 'a_projected': 8.687, 'inclination': 89.17, 'compactness': 0.003, # WD is veel minder compact 'shapiro_measured': 1.41e-5, 'shapiro_error': 0.11e-5, 'notes': '1.97 M☉ pulsar. Demorest et al. 2010, Nature.' }, 'J1903+0327': { 'name': 'Eccentric MSP', 'M_companion': 1.03 * M_sun, # MS star companion 'a_projected': 105.6, 'inclination': 73.7, 'compactness': 0.001, 'shapiro_measured': 5.6e-6, 'shapiro_error': 0.8e-6, 'notes': 'Freire et al. 2011, MNRAS.' }, 'B1913+16': { 'name': 'Hulse-Taylor', 'M_companion': 1.39 * M_sun, 'a_projected': 2.342, 'inclination': 47.0, 'compactness': 0.20, 'shapiro_measured': 3.5e-6, # geschat (moeilijk te meten bij lage inclinatie) 'shapiro_error': 1.5e-6, 'notes': 'Nobel Prize 1993. Weisberg & Taylor 2005.' } } # ─── Shapiro-Vertraging ───────────────────────────────────── def shapiro_delay_art(M_companion, inclination, orbital_phase): """ Standaard ART Shapiro-vertraging. Δt = -2GM/c³ × ln(1 - e·cos(φ) × sin(i)) Vereenvoudigd (circulaire baan): Δt(φ) = -2r × ln(1 - s × sin(φ)) waar r = GM_c/c³ (range) en s = sin(i) (shape). Parameters: M_companion: massa metgezel (kg) inclination: inclinatie (graden) orbital_phase: orbitale fase (radialen) Returns: Shapiro-vertraging in seconden """ r = G * M_companion / c**3 # range parameter (s) s = np.sin(np.radians(inclination)) # shape parameter delay = -2 * r * np.log(1 - s * np.sin(orbital_phase)) return delay def shapiro_delay_grid(M_companion, inclination, orbital_phase, beta=0.01, alpha=2.0, R_companion=10e3): """ Grid-model Shapiro-vertraging: ART + niet-lineaire correctie. Δt_grid = Δt_ART × [1 + β × (GM/Rc²)^α] De correctieterm hangt af van de compactheid GM/Rc² van het object dat het signaal passeert. Voor neutronensterren (GM/Rc² ~ 0.2) is dit significant; voor witte dwergen (~ 0.001) verwaarloosbaar. Parameters: beta: grid-dichtheidscoëfficiënt (vrije parameter) alpha: niet-lineariteitsexponent (verwacht > 1) R_companion: straal van metgezel (m) """ delay_art = shapiro_delay_art(M_companion, inclination, orbital_phase) # Compactheid xi = G * M_companion / (R_companion * c**2) # Grid-correctie correction = 1 + beta * xi**alpha return delay_art * correction # ─── Residueel-Analyse ────────────────────────────────────── def analyze_shapiro_residuals(beta=0.01, alpha=2.0): """ Bereken verwachte residuelen (grid - ART) voor alle pulsars als functie van compactheid. Returns: compactnesses: compactheid per pulsar residuals: verwachte afwijking (fractie van ART) names: pulsarnamen """ phase = np.pi / 2 # maximale Shapiro-vertraging compactnesses = [] residuals_frac = [] residuals_abs = [] names = [] errors = [] for key, p in PULSARS.items(): xi = p['compactness'] # Fractional residual correction = beta * xi**alpha residuals_frac.append(correction) # Absolute residual (seconden) dt_art = p['shapiro_measured'] dt_residual = dt_art * correction residuals_abs.append(dt_residual) compactnesses.append(xi) names.append(p['name']) errors.append(p['shapiro_error']) return (np.array(compactnesses), np.array(residuals_frac), np.array(residuals_abs), names, np.array(errors)) # ─── Detecteerbaarheids-Analyse ────────────────────────────── def detectability_analysis(): """ Bepaal voor welke β en α combinaties het grid-signaal detecteerbaar is boven de huidige meetfout. """ beta_range = np.logspace(-4, 0, 100) alpha_range = [1.5, 2.0, 2.5, 3.0] # Gebruik J0737-3039A als referentie (sterkste Shapiro-meting) p = PULSARS['J0737-3039A'] xi = p['compactness'] error_frac = p['shapiro_error'] / p['shapiro_measured'] results = {} for alpha in alpha_range: signal = beta_range * xi**alpha detectable = signal > 3 * error_frac # 3σ detectie results[alpha] = (beta_range, signal, detectable) return results, error_frac # ─── Hoofdvisualisatie ─────────────────────────────────────── def create_analysis_figure(output_dir): """ 4-paneel analyse: A) Shapiro-vertragingscurve: ART vs Grid-model B) Residueel als functie van compactheid C) Detecteerbaarheidskaart (β vs α) D) Toekomstige experimenten: wat is nodig? """ np.random.seed(42) fig, axes = plt.subplots(2, 2, figsize=(16, 12)) fig.suptitle('Voorspelling #18: Niet-Lineaire Lichtvertraging nabij Extreme Massa\n' 'Spectrum Sub-quantum Grid Model — Pulsar-Timing Analyse', fontsize=14, fontweight='bold', y=0.98) # ─── Panel A: Shapiro-curve ─────────────────────────────── ax = axes[0, 0] phase = np.linspace(0, 2*np.pi, 1000) M_comp = 1.249 * M_sun incl = 88.69 R_ns = 10e3 # 10 km neutronenster delay_art = shapiro_delay_art(M_comp, incl, phase) * 1e6 # µs delay_grid_001 = shapiro_delay_grid(M_comp, incl, phase, beta=0.01, alpha=2.0, R_companion=R_ns) * 1e6 delay_grid_01 = shapiro_delay_grid(M_comp, incl, phase, beta=0.1, alpha=2.0, R_companion=R_ns) * 1e6 ax.plot(np.degrees(phase), delay_art, 'b-', linewidth=2, label='ART (standaard)') ax.plot(np.degrees(phase), delay_grid_001, 'r--', linewidth=1.5, label='Grid β=0.01, α=2') ax.plot(np.degrees(phase), delay_grid_01, 'r-', linewidth=1.5, alpha=0.7, label='Grid β=0.1, α=2') ax.set_xlabel('Orbitale Fase (°)') ax.set_ylabel('Shapiro-Vertraging (µs)') ax.set_title('A) Shapiro-Vertraging: ART vs Grid-Model\n(PSR J0737-3039A)') ax.legend(fontsize=9) ax.grid(True, alpha=0.3) # ─── Panel B: Residueel vs Compactheid ──────────────────── ax = axes[0, 1] # Theoretische curve xi_range = np.logspace(-4, -0.3, 200) for beta, alpha, color, ls in [(0.01, 2.0, 'red', '-'), (0.001, 2.0, 'orange', '--'), (0.01, 3.0, 'green', '-.')]: frac_residual = beta * xi_range**alpha ax.loglog(xi_range, frac_residual, color=color, linestyle=ls, linewidth=2, label=f'β={beta}, α={alpha}') # Pulsars comps, frac_res, abs_res, names, errors = analyze_shapiro_residuals(beta=0.01, alpha=2.0) error_frac = errors / np.array([PULSARS[k]['shapiro_measured'] for k in PULSARS]) ax.scatter(comps, frac_res, s=100, c='red', zorder=5, edgecolors='black') for i, name in enumerate(names): ax.annotate(name, (comps[i], frac_res[i]), textcoords='offset points', xytext=(10, 5), fontsize=8) # Meetfout-grenzen ax.axhline(0.05, color='gray', linestyle=':', alpha=0.5) ax.text(1e-4, 0.06, 'Huidige meetprecisie (~5%)', fontsize=8, color='gray') ax.axhline(0.001, color='gray', linestyle=':', alpha=0.5) ax.text(1e-4, 0.0012, 'SKA precisie (~0.1%)', fontsize=8, color='gray') ax.set_xlabel('Compactheid GM/Rc²') ax.set_ylabel('Fractioneel Residueel (grid - ART)/ART') ax.set_title('B) Grid-Correctie vs Compactheid') ax.legend(fontsize=9) ax.set_xlim(1e-4, 0.5) ax.set_ylim(1e-8, 1) ax.grid(True, alpha=0.3, which='both') # ─── Panel C: Detecteerbaarheidskaart ───────────────────── ax = axes[1, 0] beta_arr = np.logspace(-4, 0, 100) alpha_arr = np.linspace(1.0, 4.0, 100) B, A = np.meshgrid(beta_arr, alpha_arr) # J0737-3039A als referentie xi_ref = 0.18 error_frac_ref = PULSARS['J0737-3039A']['shapiro_error'] / \ PULSARS['J0737-3039A']['shapiro_measured'] # Signaal/ruis verhouding SNR = (B * xi_ref**A) / error_frac_ref im = ax.pcolormesh(beta_arr, alpha_arr, np.log10(SNR), cmap='RdYlGn', vmin=-2, vmax=2, shading='auto') plt.colorbar(im, ax=ax, label='log₁₀(SNR)') # Contourlijnen bij 1σ, 3σ, 5σ ax.contour(beta_arr, alpha_arr, SNR, levels=[1, 3, 5], colors=['white', 'yellow', 'red'], linewidths=[1, 2, 2]) ax.text(0.05, 1.5, '5σ', color='red', fontsize=12, fontweight='bold') ax.text(0.01, 2.2, '3σ', color='yellow', fontsize=12, fontweight='bold') ax.text(0.003, 2.8, '1σ', color='white', fontsize=12, fontweight='bold') ax.set_xscale('log') ax.set_xlabel('β (grid-dichtheidscoëfficiënt)') ax.set_ylabel('α (niet-lineariteitsexponent)') ax.set_title('C) Detecteerbaarheidskaart (PSR J0737-3039A)') ax.grid(True, alpha=0.3) # ─── Panel D: Toekomstige Experimenten ──────────────────── ax = axes[1, 1] experiments = { 'Huidige precisie\n(Kramér 2021)': {'sensitivity': 0.05, 'year': 2021, 'color': 'blue'}, 'NANOGrav 15yr\n(2023)': {'sensitivity': 0.03, 'year': 2023, 'color': 'green'}, 'SKA Phase 1\n(2028)': {'sensitivity': 0.005, 'year': 2028, 'color': 'orange'}, 'SKA Full\n(2032)': {'sensitivity': 0.001, 'year': 2032, 'color': 'red'}, 'Next-gen Pulsar\nTiming (2035+)': {'sensitivity': 0.0001, 'year': 2035, 'color': 'purple'} } names_exp = list(experiments.keys()) sensitivities = [v['sensitivity'] for v in experiments.values()] years = [v['year'] for v in experiments.values()] colors = [v['color'] for v in experiments.values()] ax.barh(names_exp, sensitivities, color=colors, alpha=0.7, height=0.5) # Grid-model signaalsterktes for beta, ls in [(0.1, '--'), (0.01, ':'), (0.001, '-.')]: signal = beta * 0.18**2 # α=2, ξ=0.18 ax.axvline(signal, color='red', linestyle=ls, alpha=0.7, linewidth=1.5) ax.text(signal * 1.2, -0.3, f'β={beta}', fontsize=7, color='red', rotation=90) ax.set_xscale('log') ax.set_xlabel('Fractionele Precisie op Shapiro-Vertraging') ax.set_title('D) Experimenten vs Grid-Signaalsterkte (α=2)') ax.set_xlim(1e-5, 0.2) ax.grid(True, alpha=0.3, axis='x') ax.invert_xaxis() plt.tight_layout(rect=[0, 0, 1, 0.95]) filepath = os.path.join(output_dir, 'voorspelling_18_pulsar_timing.png') plt.savefig(filepath, dpi=150, bbox_inches='tight', facecolor='white') plt.close() print(f"Saved: {filepath}") return filepath # ─── Main ──────────────────────────────────────────────────── if __name__ == "__main__": output_dir = os.path.dirname(os.path.abspath(__file__)) print("Voorspelling #18: Niet-Lineaire Lichtvertraging nabij Extreme Massa") print("=" * 70) filepath = create_analysis_figure(output_dir) # Samenvatting print("\nResultaat-samenvatting:") print("-" * 70) comps, frac_res, abs_res, names, errors = analyze_shapiro_residuals(beta=0.01, alpha=2.0) print(f"\n{'Pulsar':<20} {'Compactheid':>12} {'Grid-correctie':>15} {'Meetbaar?':>10}") print("-" * 60) for i, name in enumerate(names): error_frac = errors[i] / list(PULSARS.values())[i]['shapiro_measured'] detectable = "JA" if frac_res[i] > 3 * error_frac else "NEE" print(f"{name:<20} {comps[i]:>12.4f} {frac_res[i]:>15.2e} {detectable:>10}") print("\nImplicatie voor Spectrum Grid Model:") print("→ Met β=0.01, α=2: signaal is ~3.2×10⁻⁴ van Shapiro-delay") print("→ Huidige precisie (~5%): NIET detecteerbaar") print("→ SKA Phase 1 (~0.5%): grensdetectie voor β>0.1") print("→ SKA Full (~0.1%): detecteerbaar voor β>0.01") print("→ Sterkste test: PSR J0737-3039A (dubbele pulsar, ξ=0.18)") print(f"\n✓ Analyse voltooid. Figuur: {filepath}")