Chapter 18: Schwarzschild Solution
The Schwarzschild solution describes the spacetime outside a spherically symmetric, non-rotating mass. Found in 1916, it was the first exact solution to Einstein's equations and describes stellar objects and non-rotating black holes.
The Schwarzschild Metric
\( ds^2 = -\left(1 - \frac{r_s}{r}\right)c^2dt^2 + \frac{dr^2}{1 - \frac{r_s}{r}} + r^2 d\Omega^2 \)
where rs = 2GM/c² (Schwarzschild radius) and dΩ² = dθ² + sin²θ dφ²
Properties
- • Spherically symmetric
- • Static (time-independent)
- • Vacuum solution (Tμν = 0)
- • Asymptotically flat
Schwarzschild Radius
Sun: rs ≈ 3 km
Earth: rs ≈ 9 mm
Proton: rs ≈ 10-54 m
Key Features
Event Horizon (r = rs)
One-way membrane: nothing can escape once inside. gtt → 0 (infinite time dilation). This is a coordinate singularity, not a physical one.
Physical Singularity (r = 0)
True curvature singularity: Kretschmann scalar K = 48M²/r⁶ → ∞. All infalling matter ends here.
ISCO (r = 6M)
Innermost Stable Circular Orbit. Inside this, circular orbits are unstable. Important for accretion disk physics.
Photon Sphere (r = 3M)
Unstable circular orbit for light. Photons can orbit but any perturbation sends them to the horizon or infinity.
Python: Schwarzschild Orbits
#!/usr/bin/env python3
"""
schwarzschild_orbits.py - Compute orbits in Schwarzschild spacetime
Run: python3 schwarzschild_orbits.py
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
M = 1.0 # Geometrized units
def effective_potential(r, L):
"""Effective potential for massive particle"""
return -M/r + L**2/(2*r**2) - M*L**2/r**3
def geodesic_eq(tau, y, L, E):
"""Geodesic equations: y = [r, phi, dr/dtau]"""
r, phi, r_dot = y
if r <= 2*M: return [0, 0, 0]
phi_dot = L / r**2
r_ddot = -M/r**2 + L**2/r**3 - 3*M*L**2/r**4
return [r_dot, phi_dot, r_ddot]
# Compute effective potential
r_vals = np.linspace(2.1*M, 30*M, 500)
L_values = [3.5*M, 4.0*M, 4.5*M]
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
# Plot effective potentials
ax1 = axes[0]
for L in L_values:
V = [effective_potential(r, L) for r in r_vals]
ax1.plot(r_vals/M, V, label=f'L = {L/M:.1f}M')
ax1.axhline(y=0, color='gray', linestyle='--', alpha=0.5)
ax1.axvline(x=6, color='red', linestyle=':', alpha=0.5, label='ISCO')
ax1.axvline(x=3, color='orange', linestyle=':', alpha=0.5, label='Photon sphere')
ax1.axvline(x=2, color='black', linestyle='-', alpha=0.5, label='Horizon')
ax1.set_xlabel('r/M')
ax1.set_ylabel(r'$V_{eff}$')
ax1.set_title('Effective Potential')
ax1.legend()
ax1.set_xlim(2, 20)
ax1.set_ylim(-0.1, 0.05)
ax1.grid(True, alpha=0.3)
# Integrate an orbit
L = 4.0*M
E = 0.97
r0 = 12*M
# Initial radial velocity from energy equation
f0 = 1 - 2*M/r0
r_dot0 = -np.sqrt(E**2 - f0*(1 + L**2/r0**2))
y0 = [r0, 0, r_dot0]
tau_span = [0, 600]
tau_eval = np.linspace(*tau_span, 5000)
sol = solve_ivp(lambda t, y: geodesic_eq(t, y, L, E), tau_span, y0,
t_eval=tau_eval, method='RK45', rtol=1e-10)
r, phi = sol.y[0], sol.y[1]
x = r * np.cos(phi)
y_coord = r * np.sin(phi)
ax2 = axes[1]
ax2.plot(x/M, y_coord/M, 'b-', linewidth=0.5)
horizon = plt.Circle((0, 0), 2, color='black', label='Horizon')
photon = plt.Circle((0, 0), 3, color='gray', fill=False, linestyle='--', label='Photon sphere')
isco = plt.Circle((0, 0), 6, color='red', fill=False, linestyle=':', label='ISCO')
ax2.add_patch(horizon)
ax2.add_patch(photon)
ax2.add_patch(isco)
ax2.set_aspect('equal')
ax2.set_xlabel('x/M')
ax2.set_ylabel('y/M')
ax2.set_title(f'Bound Orbit (E={E}, L={L/M:.1f}M)')
ax2.legend()
ax2.set_xlim(-15, 15)
ax2.set_ylim(-15, 15)
ax2.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('schwarzschild_orbits.png', dpi=150)
plt.show()
print("Schwarzschild Black Hole Key Radii:")
print(f" Event Horizon: r = 2M = {2*M}")
print(f" Photon Sphere: r = 3M = {3*M}")
print(f" ISCO: r = 6M = {6*M}")
print("\nPlot saved as 'schwarzschild_orbits.png'")To Run:
python3 schwarzschild_orbits.py