Chapter 27: Numerical Relativity

Numerical relativity solves Einstein's equations computationally for situations too complex for analytical solutions. Essential for LIGO waveform templates and understanding binary black hole mergers.

BSSN Formulation

The Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formulation is the standard for stable evolution:

Key Variables

• Conformal factor: $\phi = \ln(h)/12$
• Conformal metric: γ̃_ij = e^-4φh_ij
• Traceless extrinsic curvature: Ã_ij
• Conformal connection functions: Γ̃ⁱ

Moving Punctures

The 2005 breakthrough: allow black hole "punctures" to move through the grid using "1+log" slicing and "gamma-driver" shift conditions.

1+log Slicing

∂_tα = -2αK (avoids singularity)

Gamma-Driver

∂_tβⁱ = 3/4 Bⁱ (coordinate motion)

Python: Wave Equation Toy Model

This simulation demonstrates finite difference methods used in numerical relativity by solving the 1D wave equation. The wave equation shares the same hyperbolic character as Einstein's equations.

Wave Equation Evolution Simulation

Python

Evolves a Gaussian pulse using finite differences with CFL stability

script.py170 lines

#!/usr/bin/env python3
"""
Numerical Relativity - 1D Wave Equation Evolution
Demonstrates finite difference methods used in numerical relativity
"""
import numpy as np
import matplotlib.pyplot as plt
import base64
from io import BytesIO

print("╔══════════════════════════════════════════════════════════════════╗")
print("║         NUMERICAL RELATIVITY: WAVE EQUATION TOY MODEL           ║")
print("╠══════════════════════════════════════════════════════════════════╣")
print("║  The 1D wave equation serves as a toy model for understanding   ║")
print("║  the hyperbolic nature of Einstein's equations                  ║")
print("╚══════════════════════════════════════════════════════════════════╝")

print("\n" + "="*70)
print("WAVE EQUATION AND DISCRETIZATION")
print("="*70)
print("\nWave equation: (\\frac{\\partial^2 u}{\\partial t^2} = c^2 \\frac{\\partial^2 u}{\\partial x^2})")
print("\nFinite difference approximation:")
print("(\\frac{u_i^{n+1} - 2u_i^n + u_i^{n-1}}{\\Delta t^2} = c^2 \\frac{u_{i+1}^n - 2u_i^n + u_{i-1}^n}{\\Delta x^2})")
print("\nRearranging for the update formula:")
print("(u_i^{n+1} = 2u_i^n - u_i^{n-1} + \\left(\\frac{c\\Delta t}{\\Delta x}\\right)^2 (u_{i+1}^n - 2u_i^n + u_{i-1}^n))")

# Parameters
nx = 200          # Number of spatial points
nt = 500          # Number of time steps
c = 1.0           # Wave speed
L = 1.0           # Domain length
dx = L / nx       # Spatial step
dt = 0.5 * dx / c # Time step (CFL condition: dt <= dx/c)

print("\n" + "="*70)
print("SIMULATION PARAMETERS")
print("="*70)
print(f"\n• Grid points (nx): {nx}")
print(f"• Time steps (nt): {nt}")
print(f"• Wave speed (c): {c}")
print(f"• Spatial step (Δx): {dx:.6f}")
print(f"• Time step (Δt): {dt:.6f}")
print(f"\nCFL number: (\\sigma = \\frac{c \\Delta t}{\\Delta x} = {c*dt/dx:.2f})")
print("\nCFL stability condition: (\\sigma \\leq 1) ✓")

# Initial conditions: Gaussian pulse
x = np.linspace(0, L, nx)
sigma_gauss = 0.05
u = np.exp(-((x - 0.5)**2) / (2 * sigma_gauss**2))
u_prev = u.copy()  # Zero initial velocity

print("\n" + "="*70)
print("INITIAL CONDITIONS")
print("="*70)
print("\nGaussian initial data:")
print("(u(x,0) = \\exp\\left(-\\frac{(x-0.5)^2}{2\\sigma^2}\\right))")
print(f"\nGaussian width σ = {sigma_gauss}")
print("Initial velocity: (\\frac{\\partial u}{\\partial t}\\bigg|_{t=0} = 0)")

# Storage for snapshots
snapshots = []
snapshot_times = [0, 100, 200, 300, 400]

# Evolution
print("\n" + "="*70)
print("TIME EVOLUTION")
print("="*70)
print("\nEvolving wave equation...")

for n in range(nt + 1):
    if n in snapshot_times:
        snapshots.append((n * dt, u.copy()))
        print(f"  Snapshot at t = {n*dt:.3f} (step {n})")

if n < nt:
        # Compute new timestep
        u_new = np.zeros(nx)
        # Interior points
        courant = (c * dt / dx)**2
        for i in range(1, nx - 1):
            u_new[i] = 2*u[i] - u_prev[i] + courant * (u[i+1] - 2*u[i] + u[i-1])
        # Boundary conditions (fixed ends)
        u_new[0] = 0
        u_new[-1] = 0
        # Update
        u_prev = u.copy()
        u = u_new.copy()

print("\nEvolution complete!")

# Create visualization
fig, axes = plt.subplots(2, 3, figsize=(14, 8))
axes = axes.flatten()

colors = ['#4ade80', '#22d3ee', '#f472b6', '#facc15', '#a78bfa']

for idx, (t, u_snap) in enumerate(snapshots):
    ax = axes[idx]
    ax.plot(x, u_snap, color=colors[idx], linewidth=2)
    ax.fill_between(x, u_snap, alpha=0.3, color=colors[idx])
    ax.set_xlabel('x', fontsize=11)
    ax.set_ylabel('u(x,t)', fontsize=11)
    ax.set_title(f't = {t:.3f}', fontsize=12, fontweight='bold', color=colors[idx])
    ax.set_xlim(0, 1)
    ax.set_ylim(-1.2, 1.2)
    ax.grid(True, alpha=0.3)
    ax.axhline(y=0, color='gray', linestyle='-', alpha=0.5)

# Energy conservation check
ax = axes[5]
energies = []
for t, u_snap in snapshots:
    # Approximate energy: E ~ ∫(u_t² + c²u_x²)dx
    # For simplicity, just compute ∫u²dx
    E = np.trapezoid(u_snap**2, x)
    energies.append(E)

times_snap = [t for t, _ in snapshots]
ax.bar(range(len(energies)), energies, color=colors, alpha=0.7)
ax.set_xlabel('Snapshot', fontsize=11)
ax.set_ylabel('∫u²dx', fontsize=11)
ax.set_title('Energy Conservation Check', fontsize=12, fontweight='bold')
ax.set_xticks(range(len(energies)))
ax.set_xticklabels([f't={t:.2f}' for t in times_snap], fontsize=9)
ax.grid(True, alpha=0.3, axis='y')

plt.suptitle('Wave Equation Evolution (Toy Model for Numerical Relativity)',
             fontsize=14, fontweight='bold', color='white')
fig.patch.set_facecolor('#1a1a2e')
for ax in axes:
    ax.set_facecolor('#16213e')
    ax.tick_params(colors='white')
    ax.xaxis.label.set_color('white')
    ax.yaxis.label.set_color('white')
    ax.title.set_color('white')
    for spine in ax.spines.values():
        spine.set_color('gray')

plt.tight_layout()

# Save to base64
# Save plot to file
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0a0a1a')
plt.close()

print("\n" + "="*70)
print("NUMERICAL ANALYSIS")
print("="*70)
print("\nEnergy in each snapshot (∫u²dx):")
for i, (t, _) in enumerate(snapshots):
    print(f"  t = {t:.3f}: E = {energies[i]:.6f}")
print(f"\nEnergy variation: {100*(max(energies)-min(energies))/energies[0]:.2f}%")
print("\n(Note: Energy not conserved due to reflecting boundary conditions)")

print("\n" + "="*70)
print("CONNECTION TO NUMERICAL RELATIVITY")
print("="*70)
print("\nKey concepts shared with full numerical relativity:")
print("\n1. BSSN evolution equations are also hyperbolic (wave-like)")
print("2. CFL condition ensures stability: (\\Delta t < \\Delta x / c)")
print("3. Finite difference methods discretize spacetime")
print("4. Constraint preservation requires careful formulation")
print("\nIn full NR, we evolve:")
print("• Conformal metric: (\\tilde{\\gamma}_{ij})")
print("• Conformal factor: (\\phi)")
print("• Traceless extrinsic curvature: (\\tilde{A}_{ij})")
print("• Mean curvature: (K)")

# Image saved as output.png

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Fortran: BSSN Variable Computation

This program computes the key BSSN variables for Schwarzschild initial data, demonstrating the conformal decomposition and gauge conditions used in modern numerical relativity codes.

BSSN Variables Visualization

Python

Lapse function, conformal factor, and metric components vs radius

bssn_plot.py129 lines

import subprocess, tempfile, os, io, base64
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

fortran_source = r"""
program bssn_vars
    implicit none
    real(8), parameter :: G = 6.67430d-11
    real(8), parameter :: c = 2.99792458d8
    real(8), parameter :: M_sun = 1.98892d30
    real(8) :: M, r_s, r, alpha, psi4, chi_var, f
    integer :: i, n_pts

M = 10.0d0 * M_sun
    r_s = 2.0d0 * G * M / (c * c)

n_pts = 200
    print *, 'DATA_START'
    do i = 1, n_pts
        r = 2.01d0 * r_s + dble(i-1) * 48.0d0 * r_s / dble(n_pts-1)
        f = 1.0d0 - r_s / r
        alpha = sqrt(f)
        psi4 = (1.0d0 + r_s / (4.0d0 * r))**4
        chi_var = 1.0d0 / psi4
        write(*,'(F14.8,A,F14.8,A,F14.8,A,F14.8,A,F14.8)') &
            r/r_s, ',', alpha, ',', psi4, ',', chi_var, ',', 1.0d0/f
    end do
    print *, 'DATA_END'
end program bssn_vars
"""

workdir = tempfile.mkdtemp()
src = os.path.join(workdir, 'bssn.f90')
exe = os.path.join(workdir, 'bssn')

with open(src, 'w') as f:
    f.write(fortran_source)

comp = subprocess.run(['gfortran', '-O2', '-o', exe, src],
                      capture_output=True, text=True, timeout=30)
if comp.returncode != 0:
    print("Compilation error:", comp.stderr)
else:
    result = subprocess.run([exe], capture_output=True, text=True, timeout=30, cwd=workdir)
    output = result.stdout

lines = []
    in_data = False
    for line in output.split('\n'):
        if 'DATA_START' in line:
            in_data = True
            continue
        if 'DATA_END' in line:
            break
        if in_data and line.strip():
            lines.append(line.strip())

r_arr, alpha_arr, psi4_arr, chi_arr, grr_arr = [], [], [], [], []
    for line in lines:
        vals = [float(x) for x in line.split(',')]
        r_arr.append(vals[0])
        alpha_arr.append(vals[1])
        psi4_arr.append(vals[2])
        chi_arr.append(vals[3])
        grr_arr.append(vals[4])

r_arr = np.array(r_arr)
    alpha_arr = np.array(alpha_arr)
    psi4_arr = np.array(psi4_arr)
    chi_arr = np.array(chi_arr)
    grr_arr = np.array(grr_arr)

fig, axes = plt.subplots(1, 3, figsize=(18, 6))
    fig.patch.set_facecolor('#0f172a')
    fig.suptitle('BSSN Variables for Schwarzschild Initial Data', color='white', fontsize=14, fontweight='bold')

ax = axes[0]
    ax.set_facecolor('#1e293b')
    ax.plot(r_arr, alpha_arr, color='#22d3ee', linewidth=2, label='$\\alpha$ (lapse)')
    ax.axhline(y=1.0, color='#94a3b8', linestyle=':', alpha=0.5)
    ax.axvline(x=2.0, color='#ef4444', linestyle='--', alpha=0.7, label='Horizon ($r=r_s$)')
    ax.set_xlabel('$r / r_s$', color='#cbd5e1')
    ax.set_ylabel('Lapse function $\\alpha$', color='#cbd5e1')
    ax.set_title('Lapse: $\\alpha = \\sqrt{1 - r_s/r}$', color='white')
    ax.legend(facecolor='#1e293b', edgecolor='#334155', labelcolor='#cbd5e1')
    ax.tick_params(colors='#cbd5e1')
    for spine in ax.spines.values():
        spine.set_color('#334155')
    ax.grid(True, color='#334155', alpha=0.4)
    ax.set_ylim(0, 1.1)

ax2 = axes[1]
    ax2.set_facecolor('#1e293b')
    ax2.plot(r_arr, psi4_arr, color='#f472b6', linewidth=2, label='$\\psi^4$ (conformal)')
    ax2.plot(r_arr, chi_arr, color='#4ade80', linewidth=2, label='$\\chi = \\psi^{-4}$')
    ax2.axvline(x=2.0, color='#ef4444', linestyle='--', alpha=0.7)
    ax2.set_xlabel('$r / r_s$', color='#cbd5e1')
    ax2.set_ylabel('Conformal factor', color='#cbd5e1')
    ax2.set_title('Conformal Decomposition', color='white')
    ax2.legend(facecolor='#1e293b', edgecolor='#334155', labelcolor='#cbd5e1')
    ax2.tick_params(colors='#cbd5e1')
    for spine in ax2.spines.values():
        spine.set_color('#334155')
    ax2.grid(True, color='#334155', alpha=0.4)

ax3 = axes[2]
    ax3.set_facecolor('#1e293b')
    ax3.plot(r_arr, grr_arr, color='#fbbf24', linewidth=2, label='$g_{rr} = (1-r_s/r)^{-1}$')
    ax3.axvline(x=2.0, color='#ef4444', linestyle='--', alpha=0.7, label='Horizon')
    ax3.set_xlabel('$r / r_s$', color='#cbd5e1')
    ax3.set_ylabel('$g_{rr}$', color='#cbd5e1')
    ax3.set_title('Radial Metric Component', color='white')
    ax3.legend(facecolor='#1e293b', edgecolor='#334155', labelcolor='#cbd5e1')
    ax3.tick_params(colors='#cbd5e1')
    for spine in ax3.spines.values():
        spine.set_color('#334155')
    ax3.grid(True, color='#334155', alpha=0.4)
    ax3.set_ylim(0, 20)

plt.tight_layout()
    buf = io.BytesIO()
    plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0a0a1a'))
    buf.close()
    plt.close()
    print("BSSN initial data: K=0 (time-symmetric), lapse -> 0 at horizon")
    # Image saved as output.png

Click Run to execute the Python code

Code will be executed with Python 3 on the server

← ADM Formalism Next: Gravitational Lensing →