The Eddington eclipse expedition: starlight bent by the Sun

General relativity predicts that a light ray passing a mass \(M\) at impact parameter \(b\) is deflected by \(\delta = 4GM/(c^2b)\), twice what Newtonian gravity gives for a particle moving at \(c\). For a ray grazing the Sun that is 1.75 arcseconds. On 29 May 1919 expeditions to Sobral and Príncipe photographed the stars around the eclipsed Sun and compared their positions with night-time plates of the same field. The stars had shifted outward by the Einstein amount.

This example computes the deflection with light_deflection_angle() in solar units, draws the displacement pattern of a field of stars around the eclipsed Sun, and uses exact_deflection_angle() to show that the weak-field formula is essentially exact at the Sun but fails near a black hole’s photon sphere.

import matplotlib.pyplot as plt
import numpy as np

from physicskit.relativity.chapters.lensing import exact_deflection_angle
from physicskit.relativity.chapters.schwarzschild import SchwarzschildBlackHole
from physicskit.relativity.utils import constants as const

ARCSEC = np.pi / 180 / 3600

The grazing deflection

In geometrized units the Sun’s mass is \(GM_\odot/c^2 = 1.477\) km.

M_sun = const.solar_masses_to_geometrized(1.0)
R_sun = 6.957e8
sun = SchwarzschildBlackHole(M=M_sun)
delta_gr = sun.light_deflection_angle(R_sun) / ARCSEC
print(f"GM_sun / c^2 = {M_sun / 1e3:.3f} km")
print(f"deflection at the solar limb: Einstein {delta_gr:.3f} arcsec, Newtonian {delta_gr / 2:.3f} arcsec")
print("1919 results: Sobral 1.98 +/- 0.12 arcsec, Principe 1.61 +/- 0.30 arcsec")

GM_sun / c^2 = 1.477 km
deflection at the solar limb: Einstein 1.751 arcsec, Newtonian 0.876 arcsec
1919 results: Sobral 1.98 +/- 0.12 arcsec, Principe 1.61 +/- 0.30 arcsec

What the plates showed

Every star is pushed radially away from the Sun by \(\delta(b) = 1.75'' \times R_\odot/b\). The shifts, hugely exaggerated here, fall off as \(1/b\): only stars near the limb move measurably, which is why an eclipse was needed.

rng = np.random.default_rng(1919)
r_star = R_sun * rng.uniform(1.5, 6.0, 40)
phi = rng.uniform(0, 2 * np.pi, 40)
shift = np.array([sun.light_deflection_angle(b) for b in r_star]) / ARCSEC
arrow_per_arcsec = 0.6  # drawn arrow length, in solar radii, per arcsecond of shift
exaggeration = arrow_per_arcsec * 959.6  # the Sun's radius subtends 959.6 arcsec
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
ax1.add_patch(plt.Circle((0, 0), 1.0, color="k"))
ax1.add_patch(plt.Circle((0, 0), 1.6, color="gold", alpha=0.25, lw=0))
x0, y0 = r_star / R_sun * np.cos(phi), r_star / R_sun * np.sin(phi)
ax1.plot(x0, y0, "o", color="0.6", ms=4, label="night-time position")
ax1.quiver(
    x0,
    y0,
    np.cos(phi) * shift,
    np.sin(phi) * shift,
    angles="xy",
    scale_units="xy",
    scale=1 / arrow_per_arcsec,
    color="firebrick",
    width=0.004,
    label=f"eclipse shift (x{exaggeration:.0f})",
)
ax1.set_xlim(-6.5, 6.5)
ax1.set_ylim(-6.5, 6.5)
ax1.set_aspect("equal")
ax1.set_xlabel(r"x [$R_\odot$]")
ax1.set_ylabel(r"y [$R_\odot$]")
ax1.set_title("Star field around the eclipsed Sun")
ax1.legend(fontsize=8, loc="upper right")

b_over_R = np.linspace(1, 6, 100)
ax2.plot(b_over_R, 1.75 / b_over_R, color="firebrick", label=r"Einstein $4GM/c^2b$")
ax2.plot(b_over_R, 0.875 / b_over_R, "--", color="steelblue", label="Newtonian (half)")
ax2.errorbar([1.0, 1.0], [1.98, 1.61], yerr=[0.12, 0.30], fmt="s", color="k", capsize=3, label="1919, extrapolated to the limb")
ax2.set_xlabel(r"impact parameter $b / R_\odot$")
ax2.set_ylabel("deflection [arcsec]")
ax2.set_title("Deflection falls off as 1/b")
ax2.legend(fontsize=8)
fig.tight_layout()

Where the weak-field formula holds

Integrating the exact null geodesic around a unit-mass black hole shows \(4M/b\) is accurate for \(b \gg M\). The Sun’s limb is at \(b \approx 470{,}000\,M\); near the photon sphere, \(b_c=3\sqrt3M\approx5.2M\), the exact deflection grows without bound.

bh = SchwarzschildBlackHole(M=1.0)
b_vals = np.array([6.0, 8.0, 12.0, 20.0, 50.0])
exact = np.array([exact_deflection_angle(bh, b) for b in b_vals])
for b, ex in zip(b_vals, exact):
    print(f"b = {b:6.1f} M: exact {ex:.5f} rad, weak-field {4 / b:.5f} rad, ratio {ex / (4 / b):.3f}")
print(f"the Sun's limb: b = {R_sun / M_sun:.3e} M, so the next-order correction is of relative size ~ M/b = {M_sun / R_sun:.1e}")

plt.show()

b =    6.0 M: exact 1.72003 rad, weak-field 0.66667 rad, ratio 2.580
b =    8.0 M: exact 0.85870 rad, weak-field 0.50000 rad, ratio 1.717
b =   12.0 M: exact 0.45260 rad, weak-field 0.33333 rad, ratio 1.358
b =   20.0 M: exact 0.23614 rad, weak-field 0.20000 rad, ratio 1.181
b =   50.0 M: exact 0.08508 rad, weak-field 0.08000 rad, ratio 1.064
the Sun's limb: b = 4.711e+05 M, so the next-order correction is of relative size ~ M/b = 2.1e-06

Check

4 G M_sun / c^2 R_sun = 1.75 arcsec (Newton: half); the exact deflection approaches 4M/b from above as b grows.

assert abs(delta_gr - 1.751) < 0.002
assert np.all(exact / (4 / b_vals) > 1) and np.all(np.diff(exact / (4 / b_vals)) < 0)
series = 4 / b_vals + 15 * np.pi / (4 * b_vals**2) + 128 / (3 * b_vals**3) + 3465 * np.pi / (64 * b_vals**4)
assert abs(exact[-1] / series[-1] - 1) < 1e-4