Geodesics

Geodesics are the straightest possible curves on a Riemannian manifold—the paths of freely falling particles in general relativity and the shortest curves connecting nearby points. The exponential map converts tangent vectors into geodesic endpoints, while Jacobi fields measure how nearby geodesics spread apart or converge.

Historical Context

The study of geodesics dates to Johann Bernoulli's brachistochrone problem (1696) and Euler's work on shortest paths on surfaces (1728). Gauss studied geodesics on surfaces systematically in the Disquisitiones (1827), establishing that geodesic curvature is an intrinsic quantity.

Carl Gustav Jacob Jacobi introduced what we now call Jacobi fields in his study of conjugate points and the second variation of arc length (1837). These ideas were crucial for understanding when geodesics minimize length and when they fail to do so.

In general relativity, geodesics describe the motion of freely falling particles and the propagation of light. Karl Schwarzschild's 1916 solution gave the first exact geodesics in a curved spacetime, explaining Mercury's perihelion precession and light bending.

Derivation 1: The Geodesic Equation

A geodesic is a curve $\gamma(t)$ whose tangent vector is parallel-transported along itself:

$\nabla_{\dot{\gamma}}\dot{\gamma} = 0$

In local coordinates $x^i(t)$, writing $\dot{\gamma}^i = dx^i/dt$:

$\boxed{\frac{d^2 x^i}{dt^2} + \Gamma^i_{jk}\frac{dx^j}{dt}\frac{dx^k}{dt} = 0}$

Variational Derivation

Alternatively, geodesics arise as critical points of the length functional:

$L[\gamma] = \int_a^b \sqrt{g_{ij}\frac{dx^i}{dt}\frac{dx^j}{dt}}\,dt$

It is more convenient to work with the energy functional (which gives the same curves with affine parameterization):

$E[\gamma] = \frac{1}{2}\int_a^b g_{ij}\frac{dx^i}{dt}\frac{dx^j}{dt}\,dt$

The Euler-Lagrange equation for $E$ reproduces the geodesic equation with Christoffel symbols arising from the metric derivatives.

GR interpretation: In general relativity, freely falling particles follow timelike geodesics, and light rays follow null geodesics ($g_{\mu\nu}\dot{x}^\mu\dot{x}^\nu = 0$). The equivalence principle states that locally, geodesic motion is indistinguishable from inertial motion.

Derivation 2: The Exponential Map

For each point $p \in M$ and tangent vector $v \in T_pM$, let $\gamma_v(t)$ be the unique geodesic with $\gamma_v(0) = p$ and $\dot{\gamma}_v(0) = v$. The exponential map is:

$\exp_p: T_pM \supset U \to M, \quad \exp_p(v) = \gamma_v(1)$

Normal Coordinates

Choosing an orthonormal basis $\{e_i\}$ for $T_pM$, the map $(x^1,\ldots,x^n) \mapsto \exp_p(x^i e_i)$defines normal (Riemann) coordinates centered at $p$. In these coordinates:

$g_{ij}(0) = \delta_{ij}, \quad \Gamma^k_{ij}(0) = 0$

The metric is Euclidean to first order at the center

$g_{ij}(x) = \delta_{ij} - \frac{1}{3}R_{ikjl}(0)\,x^k x^l + O(|x|^3)$

Curvature appears at second order

Injectivity Radius

The injectivity radius $\text{inj}(p)$ is the largest radius for which $\exp_p$ is a diffeomorphism onto its image. For the unit sphere,$\text{inj}(p) = \pi$ (distance to the antipodal point). For compact manifolds, the global injectivity radius $\text{inj}(M) = \inf_p \text{inj}(p) > 0$.

Geodesic completeness: A Riemannian manifold is geodesically complete if $\exp_p$ is defined on all of $T_pM$. By the Hopf-Rinow theorem (1931), this is equivalent to completeness as a metric space, and also to the property that every pair of points can be joined by a length-minimizing geodesic.

Derivation 3: Jacobi Fields and Geodesic Deviation

Consider a one-parameter family of geodesics $\gamma_s(t)$. The variation vector field $J(t) = \frac{\partial\gamma_s}{\partial s}\big|_{s=0}$measures the infinitesimal separation between neighboring geodesics. It satisfies the Jacobi equation:

$\boxed{\frac{D^2 J}{dt^2} + R(J, \dot{\gamma})\dot{\gamma} = 0}$

Derivation from the Variation

Starting from the geodesic equation $\nabla_{\dot{\gamma}_s}\dot{\gamma}_s = 0$ for each $s$, differentiate with respect to $s$:

$\frac{D}{ds}\frac{D}{dt}\dot{\gamma}_s = 0$

Swapping the order of covariant differentiation introduces the Riemann tensor:

$\frac{D}{ds}\frac{D}{dt}V = \frac{D}{dt}\frac{D}{ds}V + R\left(\frac{\partial\gamma}{\partial s}, \frac{\partial\gamma}{\partial t}\right)V$

This yields the Jacobi equation. The curvature controls whether geodesics focus or defocus.

Constant Curvature Solutions

For a space of constant sectional curvature $K$, with initial conditions $J(0) = 0$, $J'(0) = e$:

$K > 0: \quad J(t) = \frac{\sin(\sqrt{K}\,t)}{\sqrt{K}}\,e \quad \text{(converging)}$

$K = 0: \quad J(t) = t\,e \quad \text{(constant spreading)}$

$K < 0: \quad J(t) = \frac{\sinh(\sqrt{|K|}\,t)}{\sqrt{|K|}}\,e \quad \text{(diverging)}$

Geodesic deviation in GR: The Jacobi equation in spacetime gives the equation of geodesic deviation $\ddot{\xi}^\mu = -R^\mu{}_{\nu\alpha\beta}u^\nu\xi^\alpha u^\beta$, where $u^\mu$ is the 4-velocity and $\xi^\mu$ the separation vector. This describes tidal forces—the physical effect detected by LIGO when gravitational waves pass through.

Derivation 4: Conjugate Points and Minimality

A point $q = \gamma(t_0)$ is conjugate to $p = \gamma(0)$ along a geodesic $\gamma$ if there exists a non-zero Jacobi field $J$ with $J(0) = J(t_0) = 0$. Conjugate points are where the exponential map becomes singular:

$q \text{ conjugate to } p \iff d(\exp_p)_{t_0 v} \text{ is singular}$

Second Variation Formula

The second variation of the energy functional for a geodesic $\gamma$ with variation field $W$:

$\frac{d^2E}{ds^2}\bigg|_{s=0} = \int_0^L \left[\left|\frac{DW}{dt}\right|^2 - R(W,\dot{\gamma},\dot{\gamma},W)\right]dt$

A geodesic fails to minimize length past its first conjugate point. On the unit sphere, every geodesic has its first conjugate point at distance $\pi$ (the antipodal point), so great circle arcs longer than $\pi$ are not length-minimizing.

Morse index theorem: The number of conjugate points (counted with multiplicity) along a geodesic from $p$ to $q$ equals the index (number of negative eigenvalues) of the second variation operator. This connects geodesic geometry to Morse theory.

Derivation 5: Geodesics in Schwarzschild Spacetime

The Schwarzschild metric for a mass $M$ (in units $G = c = 1$):

$ds^2 = -\left(1-\frac{2M}{r}\right)dt^2 + \left(1-\frac{2M}{r}\right)^{-1}dr^2 + r^2 d\Omega^2$

The symmetries yield two conserved quantities along geodesics: energy $E = (1-2M/r)\dot{t}$and angular momentum $L = r^2\dot{\phi}$. The radial equation becomes:

$\frac{1}{2}\dot{r}^2 + V_{\text{eff}}(r) = \frac{1}{2}E^2$

where the effective potential is:

$\boxed{V_{\text{eff}}(r) = \frac{1}{2}\left(1-\frac{2M}{r}\right)\left(\epsilon + \frac{L^2}{r^2}\right)}$

Here $\epsilon = 1$ for massive particles and $\epsilon = 0$ for photons. The extra$-ML^2/r^3$ term (absent in Newtonian gravity) causes perihelion precession and the innermost stable circular orbit (ISCO) at $r = 6M$.

Observational tests: Mercury's perihelion precesses by 43 arcseconds per century due to the geodesic deviation from Newtonian orbits. Light bending near the Sun (1.75 arcseconds) was confirmed by Eddington's 1919 eclipse expedition, making Einstein famous.

Interactive Simulation

This simulation numerically solves the geodesic equation on the 2-sphere, visualizes the exponential map metric factor for the three model geometries, computes Jacobi fields showing geodesic focusing and defocusing, and plots the Schwarzschild effective potential for massive particle orbits.

Geodesics: Geodesic Equation, Exponential Map & Jacobi Fields

Python

script.py222 lines

#!/usr/bin/env python3
"""
geodesics.py - Geodesic Equation, Exponential Map, and Jacobi Fields
Solves geodesic equations on sphere and Schwarzschild spacetime
"""
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
import base64
from io import BytesIO

# ============================================================
# 1. Geodesics on the 2-sphere
# ============================================================
print("=" * 72)
print("  GEODESICS: SHORTEST PATHS ON CURVED SPACES")
print("=" * 72)
print()

def sphere_geodesic(t, y):
    """Geodesic equation on unit S^2: (theta, phi, dtheta/dt, dphi/dt)"""
    th, ph, dth, dph = y
    ddth = np.sin(th) * np.cos(th) * dph**2
    ddph = -2.0 * np.cos(th) / np.sin(th) * dth * dph
    return [dth, dph, ddth, ddph]

print("  GEODESICS ON S^2 (GREAT CIRCLES)")
print("  " + "-" * 50)
print()
print("  Geodesic equations on the unit sphere:")
print("    d^2(theta)/dt^2 - sin(theta)cos(theta)(dphi/dt)^2 = 0")
print("    d^2(phi)/dt^2 + 2 cot(theta) (dtheta/dt)(dphi/dt) = 0")
print()

# Solve several geodesics with different initial conditions
t_span = (0, 2 * np.pi)
t_eval = np.linspace(0, 2 * np.pi, 500)

initial_conditions = [
    [np.pi/2, 0, 0, 1],        # equator
    [np.pi/4, 0, 0.5, 0.8],    # tilted
    [np.pi/3, 0, 0.7, 0.5],    # another tilt
    [np.pi/6, 0, 0.3, 1.0],    # near pole
]

geodesic_solutions = []
for ic in initial_conditions:
    sol = solve_ivp(sphere_geodesic, t_span, ic, t_eval=t_eval,
                    method='RK45', max_step=0.01)
    geodesic_solutions.append(sol)

print("  Computed {} geodesics on S^2".format(len(geodesic_solutions)))
for i, sol in enumerate(geodesic_solutions):
    th0, ph0 = sol.y[0][0], sol.y[1][0]
    print("    Geodesic {}: theta_0 = {:.4f}, phi_0 = {:.4f}".format(
        i+1, th0, ph0))
print()

# ============================================================
# 2. Exponential map
# ============================================================
print("  EXPONENTIAL MAP AT NORTH POLE")
print("  " + "-" * 50)
print()
print("  exp_p: T_pM -> M maps tangent vectors to geodesic endpoints")
print("  exp_p(v) = gamma_v(1) where gamma_v is the geodesic with")
print("    gamma_v(0) = p and gamma_v'(0) = v")
print()

# Exponential map from north pole of S^2
# At theta=epsilon (near north pole), the exponential map is approximately
# the identity map in polar coordinates
radii = np.linspace(0.01, np.pi, 100)
# On S^2, exp maps a disk of radius pi to the entire sphere
# The metric in exponential coordinates: ds^2 = dr^2 + sin^2(r) dphi^2
metric_factor = np.sin(radii) / radii  # compared to flat: r

print("  In normal coordinates at p, the metric expands as:")
print("    g_ij(x) = delta_ij - (1/3) R_{ikjl} x^k x^l + O(|x|^3)")
print()
print("  For S^2 (K=1), in polar normal coords:")
print("    ds^2 = dr^2 + sin^2(r) dphi^2")
print()

for r in [0.5, 1.0, np.pi/2, 2.0, np.pi]:
    ratio = np.sin(r) / r if r > 0 else 1.0
    print("    r = {:.4f}: sin(r)/r = {:.6f}  (flat = 1.000)".format(r, ratio))
print()

# ============================================================
# 3. Jacobi fields
# ============================================================
print("  JACOBI FIELDS ON CONSTANT-CURVATURE SPACES")
print("  " + "-" * 50)
print()
print("  Jacobi equation: D^2 J/dt^2 + R(J, gamma')gamma' = 0")
print()

t_jac = np.linspace(0, np.pi, 300)

# K > 0 (sphere): J(t) = sin(sqrt(K) t) / sqrt(K)
K_pos = 1.0
J_sphere = np.sin(np.sqrt(K_pos) * t_jac) / np.sqrt(K_pos)

# K = 0 (flat): J(t) = t
J_flat = t_jac

# K < 0 (hyperbolic): J(t) = sinh(sqrt|K| t) / sqrt|K|
K_neg = 1.0
J_hyp = np.sinh(np.sqrt(K_neg) * t_jac) / np.sqrt(K_neg)

print("  For constant curvature K, initially parallel geodesics:")
print("    K > 0 (sphere):     J(t) = sin(sqrt(K) t) / sqrt(K)")
print("    K = 0 (flat):       J(t) = t")
print("    K < 0 (hyperbolic): J(t) = sinh(sqrt|K| t) / sqrt|K|")
print()
print("  Conjugate points (where J=0) on the sphere: t = pi/sqrt(K)")
print("    For unit sphere: first conjugate point at t = pi (antipodal)")
print()

# ============================================================
# 4. Schwarzschild geodesics
# ============================================================
print("  SCHWARZSCHILD GEODESICS (ORBITS)")
print("  " + "-" * 50)
print()

# Effective potential for massive particle
r_sch = np.linspace(2.5, 30.0, 500)
M = 1.0

L_values = [3.5, 4.0, 4.5]
print("  Effective potential: V_eff(r) = (1 - 2M/r)(1 + L^2/r^2) / 2")
print()

for L in L_values:
    V_eff = 0.5 * (1 - 2*M/r_sch) * (1 + L**2/r_sch**2)
    r_circ = L**2 / (2*M) * (1 - np.sqrt(1 - 12*M**2/L**2)) if L**2 > 12*M**2 else None
    if r_circ:
        print("    L = {:.1f}: Circular orbit at r = {:.4f}".format(L, r_circ))
    else:
        print("    L = {:.1f}: No stable circular orbit".format(L))
print()

# ============================================================
# Visualization
# ============================================================
fig, axes = plt.subplots(2, 2, figsize=(14, 11))
fig.patch.set_facecolor('#0a0a0a')

for ax in axes.flat:
    ax.set_facecolor('#111111')
    for spine in ax.spines.values():
        spine.set_color('#2dd4bf')
    ax.tick_params(colors='#99f6e4')
    ax.grid(True, color='#134e4a', alpha=0.4)

# Plot 1: Geodesics on S^2 (projected to 2D)
ax1 = axes[0, 0]
colors = ['#2dd4bf', '#f472b6', '#a78bfa', '#f59e0b']
for i, sol in enumerate(geodesic_solutions):
    th, ph = sol.y[0], sol.y[1]
    x = np.sin(th) * np.cos(ph)
    y = np.sin(th) * np.sin(ph)
    z = np.cos(th)
    # Mollweide-like projection
    ax1.plot(ph % (2*np.pi), th, color=colors[i], linewidth=1.5,
             alpha=0.8, label='Geodesic {}'.format(i+1))
ax1.set_xlabel('phi', color='#5eead4')
ax1.set_ylabel('theta', color='#5eead4')
ax1.set_title('Geodesics on S^2 (coord view)', color='#99f6e4', fontsize=12, fontweight='bold')
ax1.legend(facecolor='#111111', edgecolor='#134e4a', labelcolor='#99f6e4', fontsize=8)
ax1.invert_yaxis()

# Plot 2: Exponential map metric factor
ax2 = axes[0, 1]
ax2.plot(radii, metric_factor, color='#2dd4bf', linewidth=2.5, label='S^2: sin(r)/r')
ax2.plot(radii, np.ones_like(radii), color='#6b7280', linewidth=1.5, linestyle='--', label='Flat: 1')
ax2.plot(radii, np.sinh(radii)/radii, color='#f472b6', linewidth=2.5, label='H^2: sinh(r)/r')
ax2.set_xlabel('Geodesic distance r', color='#5eead4')
ax2.set_ylabel('Metric factor', color='#5eead4')
ax2.set_title('Exponential Map Metric Factor', color='#99f6e4', fontsize=12, fontweight='bold')
ax2.legend(facecolor='#111111', edgecolor='#134e4a', labelcolor='#99f6e4')
ax2.set_ylim(0, 3)

# Plot 3: Jacobi fields
ax3 = axes[1, 0]
ax3.plot(t_jac, J_sphere, color='#34d399', linewidth=2.5, label='K>0 (sphere)')
ax3.plot(t_jac, J_flat, color='#f59e0b', linewidth=2.5, label='K=0 (flat)')
ax3.plot(t_jac, np.clip(J_hyp, 0, 5), color='#f472b6', linewidth=2.5, label='K<0 (hyperbolic)')
ax3.axhline(y=0, color='#6b7280', linewidth=0.8)
ax3.set_xlabel('Parameter t', color='#5eead4')
ax3.set_ylabel('Jacobi field J(t)', color='#5eead4')
ax3.set_title('Jacobi Fields: Geodesic Deviation', color='#99f6e4', fontsize=12, fontweight='bold')
ax3.legend(facecolor='#111111', edgecolor='#134e4a', labelcolor='#99f6e4')
ax3.set_ylim(-0.5, 5)

# Plot 4: Schwarzschild effective potential
ax4 = axes[1, 1]
for i, L in enumerate(L_values):
    V_eff = 0.5 * (1 - 2*M/r_sch) * (1 + L**2/r_sch**2)
    ax4.plot(r_sch, V_eff, color=colors[i], linewidth=2.5,
             label='L = {:.1f}'.format(L))
ax4.axhline(y=0.5, color='#6b7280', linewidth=1, linestyle='--', alpha=0.5)
ax4.set_xlabel('r / M', color='#5eead4')
ax4.set_ylabel('V_eff', color='#5eead4')
ax4.set_title('Schwarzschild Effective Potential', color='#99f6e4', fontsize=12, fontweight='bold')
ax4.legend(facecolor='#111111', edgecolor='#134e4a', labelcolor='#99f6e4')
ax4.set_ylim(0.3, 0.7)

plt.tight_layout(pad=2.0)
buf = BytesIO()
plt.savefig(buf, format='png', dpi=150, bbox_inches='tight', facecolor='#0a0a0a')
buf.seek(0)
img_b64 = base64.b64encode(buf.read()).decode()
print('<img src="data:image/png;base64,' + img_b64 + '" alt="Geodesics Visualization" />')
buf.close()
print()
print("Visualization complete.")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Summary

Geodesic Equation

Geodesics satisfy $\ddot{x}^i + \Gamma^i_{jk}\dot{x}^j\dot{x}^k = 0$, arising from parallel transport of the velocity vector or extremizing the energy functional.

Exponential Map

Maps tangent vectors to geodesic endpoints, providing normal coordinates where the metric is Euclidean to first order and curvature appears at second order.

Jacobi Fields

Measure geodesic deviation: positive curvature focuses geodesics (sphere), negative curvature defocuses them (hyperbolic space). Conjugate points mark where geodesics stop minimizing.

Schwarzschild Geodesics

The effective potential approach reveals perihelion precession, light bending, and the ISCO—all consequences of the relativistic $1/r^3$ correction.

← Connections & Curvature Next: Gauss-Bonnet Theorem →