The Schwarzschild solution and orbits around a mass

In December 1915, a month after Einstein published his field equations, Karl Schwarzschild sent him an exact solution: the spacetime outside a spherical, non-rotating mass. Schwarzschild was serving as an artillery officer on the Russian front and computing in his spare time; he died of an illness contracted there in May 1916. Einstein, who had found only approximate solutions, was astonished that the problem could be solved exactly. The Schwarzschild metric describes the spacetime around the Sun, the Earth to good accuracy, neutron stars and non-rotating black holes. This lesson derives its main features, studies its geometry, and solves for the orbits of planets and of light around it. Lesson 7 uses the orbit equation to compute the classic tests of general relativity; lesson 8 follows the metric inside its horizon.

The metric

For a static, spherically symmetric spacetime, coordinates can be chosen so that the metric takes the form \(ds^2 = -A(r)\,c^2dt^2 + B(r)\,dr^2 + r^2(d\theta^2 + \sin^2\theta\,d\phi^2)\), where \(r\) is defined so that a sphere of radius \(r\) has area \(4\pi r^2\). Solving the vacuum equations \(R_{\mu\nu} = 0\) gives \(A = 1/B = 1 - r_s/r\), with a constant \(r_s\) that the Newtonian limit fixes as \(2GM/c^2\). The Schwarzschild metric is therefore

Far from the mass, where \(r \gg r_s\), the time part is \(-(1 + 2\Phi/c^2)c^2dt^2\) with \(\Phi = -GM/r\), exactly the weak-field metric of lesson 4, so Newtonian gravity is recovered. A remarkable fact, Birkhoff's theorem (1923), states that the Schwarzschild metric is the unique spherically symmetric vacuum solution, even if the source is not static. A spherical star that pulsates or collapses, keeping its spherical symmetry, produces no change in the field outside it. In particular, a spherical star cannot emit gravitational waves, just as a spherically pulsating charge cannot emit electromagnetic radiation.

The coordinates have a physical meaning that is easy to misread. The coordinate \(t\) is the time measured by a clock at rest very far away. A clock at rest at radius \(r\) ticks proper time \(d\tau = \sqrt{1 - r_s/r}\,dt\), slower than the distant clock, as lesson 3 predicted; light from it reaches distant observers with redshift \(1 + z = 1/\sqrt{1 - r_s/r}\). The coordinate \(r\) is not the distance from the middle; the proper radial distance between two spheres is \(\int dr/\sqrt{1 - r_s/r}\), larger than the difference of their \(r\) values. Space near the mass is stretched radially.

The shape of space

The spatial geometry at a fixed time is easier to picture through an embedding diagram. Take the equatorial plane, \(\theta = \pi/2\), at constant \(t\). Its metric \(dr^2/(1 - r_s/r) + r^2d\phi^2\) is the metric of a curved two-dimensional surface, and that surface can be drawn in ordinary three-dimensional space as a surface of revolution \(z(r)\). Matching the metrics gives \(z = 2\sqrt{r_s(r - r_s)}\), a parabola rotated about the vertical axis, called Flamm's paraboloid after Ludwig Flamm, who found it in 1916. MATH-430's capstone computed its Gaussian curvature, \(K = -r_s/(2r^3)\), in three different ways. Distances measured within the surface are the true radial distances; the circumference at radius \(r\) is still \(2\pi r\). The funnel shape shows how the proper distance between neighboring circles grows as \(r_s\) is approached.

The embedding diagram is a picture of space only; it says nothing about time, and it is time that dominates the motion of slow bodies. The diagram should not be read as a rubber sheet down which planets roll.

Conserved quantities

The Schwarzschild metric does not depend on \(t\) or on \(\phi\). By the same argument that gave Clairaut's relation on the sphere in lesson 4, each of these symmetries gives a conserved quantity along every geodesic. With \(\dot{x} = dx/d\tau\), they are the energy per unit rest energy and the angular momentum per unit mass:

The first is the relativistic generalization of energy conservation; for a particle at rest at infinity, \(\tilde{E} = 1\). The second is Kepler's second law in proper time. The normalization of the four-velocity, \(g_{\mu\nu}\dot{x}^\mu\dot{x}^\nu = -c^2\), confined to the equatorial plane, then gives a single equation for the radial motion:

This looks exactly like the one-dimensional problem of PHYS-244, a particle moving in an effective potential, with two familiar terms and one new one. The first is the Newtonian gravitational potential and the second the centrifugal barrier. The third, which falls off as \(1/r^3\), is the general-relativistic correction. It is attractive and grows faster than the centrifugal barrier at small \(r\), so close enough to the mass it always wins. In Newtonian gravity, any particle with angular momentum is kept away from the middle by the centrifugal barrier. In general relativity the barrier has a maximum, and a particle with enough energy goes over the top and falls in.

Circular orbits and the innermost stable orbit

Circular orbits sit where \(dV_\text{eff}/dr = 0\). Working in units with \(G = c = 1\), where \(r_s = 2M\), the condition is \(Mr^2 - \tilde{L}^2 r + 3M\tilde{L}^2 = 0\), with two roots for each \(\tilde{L}\): a stable orbit at the bottom of the well and an unstable one at the top of the peak. The roots merge when \(\tilde{L}^2 = 12M^2\), at

Inside this innermost stable circular orbit no stable circular orbits exist: any small disturbance sends a particle spiralling in. For a black hole of ten solar masses the ISCO is at 89 km. A particle orbiting there has \(\tilde{E} = \sqrt{8/9} = 0.943\), so a particle that spirals down slowly through a disc to the ISCO and then plunges releases \(1 - 0.943\), about 5.7 parts in a hundred, of its rest energy as heat and light. That is roughly eight times the efficiency of nuclear fusion, which turns 0.7 parts in a hundred of hydrogen's mass into energy, and it is why accreting black holes power quasars, the brightest steady sources in the universe. For a rapidly spinning black hole, lesson 8 shows, the efficiency rises to as much as 42 parts in a hundred.

The unstable circular orbits extend inward to \(r = 3M\), where \(\tilde{L} \to \infty\) and \(\tilde{E} \to \infty\): only massless particles can orbit there. This is the photon sphere, at one and a half Schwarzschild radii. Light passing a black hole close to the photon sphere can loop around it several times before escaping, which is why the Event Horizon Telescope images of M87* and Sgr A* show a bright ring of radius about \(\sqrt{27}M = 5.2M\), the apparent size of the photon sphere as seen from far away, rather than the horizon itself.

The orbit equation

For the shape of orbits, use \(u = 1/r\) and \(\phi\) as the independent variable, as in the classical Kepler problem. Differentiating the energy equation with respect to \(\phi\) and using \(\tilde{L} = r^2 d\phi/d\tau\) gives

Without the last term this is Binet's equation of PHYS-244, whose solutions are the conic sections of Kepler. The extra term, smaller than the others by the factor \(3GMu/c^2 = 3r_s/(2r)\), is the whole relativistic correction for planetary orbits. Lesson 7 solves it by perturbation theory to find the perihelion advance. Here we integrate it numerically, together with the full geodesic equations, and measure the angle through which the orbit turns between successive perihelia.

import math
G, c, Msun = 6.674e-11, 2.998e8, 1.989e30
def precession(M, a, e, relativistic=True, steps_per_rad=4000):
    """Integrate u'' + u = GM/L^2 + (3GM/c^2) u^2 from perihelion to the next perihelion.
    Returns the advance per orbit in radians."""
    p = a * (1 - e * e)                  # semi-latus rectum
    k = G * M / c**2                     # GM/c^2 (m)
    A = 1 / p                            # Newtonian GM/L^2 so that u = (1 + e cos phi)/p
    eps = 3 * k if relativistic else 0.0
    u, du = (1 + e) / p, 0.0             # start at perihelion
    h = 1 / steps_per_rad
    f = lambda u: A + eps * u * u - u
    phi, prev = 0.0, du
    while True:
        # velocity Verlet for u'' = f(u)
        du_half = du + 0.5 * h * f(u)
        u += h * du_half
        du = du_half + 0.5 * h * f(u)
        phi += h
        if phi > 3.0 and prev > 0 and du <= 0 and u > (1 + e) / p * 0.5:   # u at a maximum: perihelion
            # refine the crossing linearly
            frac = prev / (prev - du)
            return phi - h + frac * h - 2 * math.pi
        prev = du
# Mercury: subtract the Newtonian run so that the small step-size error cancels
adv = precession(Msun, 5.791e10, 0.2056) - precession(Msun, 5.791e10, 0.2056, relativistic=False)
orbits_per_century = 36525 / 87.969
arcsec = adv * 180 / math.pi * 3600 * orbits_per_century
print(f"Mercury: {adv:.3e} rad per orbit = {arcsec:.1f} arcsec per century")
formula = 6 * math.pi * G * Msun / (c**2 * 5.791e10 * (1 - 0.2056**2))
print(f"formula 6 pi GM / (c^2 a (1 - e^2)) = {formula:.3e} rad")
# the star S2 around Sgr A*: a = 1031 AU, e = 0.885, M = 4.3 million suns
AU = 1.496e11
s2 = precession(4.3e6 * Msun, 1031 * AU, 0.885, steps_per_rad=20000) - precession(4.3e6 * Msun, 1031 * AU, 0.885, False, 20000)
print(f"S2: {s2*180/math.pi*60:.1f} arcmin per orbit")

Mercury: 5.020e-07 rad per orbit = 43.0 arcsec per century
formula 6 pi GM / (c^2 a (1 - e^2)) = 5.020e-07 rad
S2: 12.3 arcmin per orbit

The numerical integration recovers Mercury's 43 arcseconds per century, matching the perturbative formula derived in lesson 7. The code subtracts an identical run without the correction term, which would close exactly; this cancels the small error that the finite step size introduces into the timing of each perihelion. For the star S2, which passes within 120 AU of the black hole Sgr A* every 16 years, the advance is about 12 arcminutes per orbit, about seven thousand times larger per orbit than Mercury's 0.1 arcseconds. The GRAVITY collaboration at the European Southern Observatory measured exactly this Schwarzschild precession for S2 in 2020, after following the star through its 2018 closest approach, and Reinhard Genzel shared the 2020 Nobel Prize in Physics with Andrea Ghez for this program of measurements.

Falling in

A body dropped from rest at infinity has \(\tilde{E} = 1\) and \(\tilde{L} = 0\), and the energy equation becomes \((dr/d\tau)^2 = c^2 r_s/r\), formally the Newtonian free-fall equation. Integrating, the proper time to fall from radius \(r_0\) to the middle is \(\tau = \frac{2}{3}\sqrt{r_0^3/(c^2 r_s)}\), finite. A falling observer reaches \(r = r_s\) and then \(r = 0\) in a finite time on their own watch. The coordinate time \(t\) tells a different story: since \(dt/d\tau = 1/(1 - r_s/r)\) diverges at \(r_s\), a distant observer computes that the infall takes infinite coordinate time, and sees the falling object slow down, redden and fade as it approaches the Schwarzschild radius, never quite crossing. Both descriptions are correct. They reflect the fact that the coordinates \(t\) and \(r\) misbehave at \(r_s\), a coordinate singularity, not a physical one. Lesson 8 introduces coordinates that pass smoothly through it.

import math
G, c, Msun = 6.674e-11, 2.998e8, 1.989e30
M = 10 * Msun
rs = 2 * G * M / c**2
r0 = 10 * rs
# proper time from r0 to r for radial infall from rest at infinity: tau = (2/3)(r0^1.5 - r^1.5)/(c sqrt(rs))
tau = lambda r: (2/3) * (r0**1.5 - r**1.5) / (c * math.sqrt(rs))
# coordinate time by numerical integration of dt/dr = -(r/rs)^0.5 / (c (1 - rs/r))
def coord_time(r_end, n=200000):
    t, r, dr = 0.0, r0, (r0 - r_end) / n
    for _ in range(n):
        rm = r - dr / 2
        t += math.sqrt(rm / rs) / (c * (1 - rs / rm)) * dr
        r -= dr
    return t
print(f"rs = {rs/1e3:.2f} km; start at r0 = {r0/1e3:.1f} km")
for frac in (5, 2, 1.1, 1.01, 1.001):
    r = frac * rs
    print(f"r = {frac:6.3f} rs: proper time {tau(r)*1e3:8.4f} ms, far-away coordinate time {coord_time(r)*1e3:9.4f} ms")
print(f"proper time to reach the horizon {tau(rs)*1e3:.4f} ms, and r = 0 {tau(0)*1e3:.4f} ms")
# ISCO and photon sphere for this hole
print(f"ISCO {3*rs/1e3:.1f} km, photon sphere {1.5*rs/1e3:.1f} km, binding at ISCO {1-math.sqrt(8/9):.4f} of rest energy")

rs = 29.54 km; start at r0 = 295.4 km
r =  5.000 rs: proper time   1.3428 ms, far-away coordinate time    1.5556 ms
r =  2.000 rs: proper time   1.8914 ms, far-away coordinate time    2.3450 ms
r =  1.100 rs: proper time   2.0014 ms, far-away coordinate time    2.7215 ms
r =  1.010 rs: proper time   2.0105 ms, far-away coordinate time    2.9619 ms
r =  1.001 rs: proper time   2.0114 ms, far-away coordinate time    3.1900 ms
proper time to reach the horizon 2.0115 ms, and r = 0 2.0771 ms
ISCO 88.6 km, photon sphere 44.3 km, binding at ISCO 0.0572 of rest energy