Are the courses on CoursesHub.World free?

Yes, all courses on CoursesHub.World are completely free and open access. We believe in democratizing education and making university-level science courses available to everyone worldwide.

What subjects are covered on CoursesHub.World?

CoursesHub.World offers courses in physics (Quantum Mechanics, General Relativity, QFT, Plasma Physics, Cosmology), biology (Molecular Biology, Cell Physiology, Pharmacology), and earth sciences (Oceanography, Atmospheric Science, Climatology).

What level are the courses?

Our courses are designed at the graduate and advanced undergraduate level. They include rigorous mathematical derivations and are suitable for physics students, researchers, and serious self-learners.

Where do the video lectures come from?

Our 400+ video lectures come from world-renowned sources including MIT OpenCourseWare, Stanford, lectures by Nobel laureates, and other leading universities and educators.

← Back to Part I: Differential Geometry

Chapter 6: Geodesics

Geodesics are the straightest possible paths in curved spacetime. Free-falling particles follow timelike geodesics; light rays follow null geodesics. This is the essence of Einstein's equivalence principle: gravity is geometry.

The Geodesic Equation

A geodesic is a curve whose tangent vector is parallel transported along itself:

$ \frac{d^2x^\mu}{d\lambda^2} + \Gamma^\mu_{\nu\rho} \frac{dx^\nu}{d\lambda} \frac{dx^\rho}{d\lambda} = 0 $

Second-order ODE for the worldline x^μ(λ)

Timelike Geodesics

$ ds^2 < 0 $: massive particles, λ = proper time τ

Null Geodesics

$ ds^2 = 0 $: light rays, λ is affine parameter

Spacelike Geodesics

$ ds^2 > 0 $: tachyons (hypothetical), λ = proper length

Extremal Property

Geodesics extremize the action $ S = \int ds $

Constants of Motion

Killing vectors generate symmetries that yield conserved quantities along geodesics:

Energy Conservation

Time translation symmetry ∂_t gives $ E = -g_{t\mu} \frac{dx^\mu}{d\tau} = \left(1 - \frac{2M}{r}\right)\frac{dt}{d\tau} $

Angular Momentum Conservation

Axial symmetry ∂_φ gives $ L = g_{\phi\mu} \frac{dx^\mu}{d\tau} = r^2 \frac{d\phi}{d\tau} $

Normalization

$ g_{\mu\nu} \frac{dx^\mu}{d\tau} \frac{dx^\nu}{d\tau} = -1 $ for timelike geodesics

Effective Potential Analysis

Using conserved quantities, radial motion reduces to a 1D problem with effective potential:

$ \frac{1}{2}\left(\frac{dr}{d\tau}\right)^2 + V_{eff}(r) = \frac{1}{2}(E^2 - 1) $

where $ V_{eff} = -\frac{M}{r} + \frac{L^2}{2r^2} - \frac{ML^2}{r^3} $

The last term (-ML²/r³) is the relativistic correction, responsible for:

Perihelion precession of Mercury
Unstable circular orbits at r = 3M (photon sphere)
Innermost stable circular orbit (ISCO) at r = 6M

Python: Geodesic Integration in Schwarzschild

Python

geodesics_schwarzschild.py162 lines

#!/usr/bin/env python3
"""
╔══════════════════════════════════════════════════════════════════════════════╗
║           Geodesic Integration in Schwarzschild Spacetime                    ║
║                                                                              ║
║  The geodesic equation:                                                      ║
║    $\frac{d^2x^\mu}{d\tau^2} + \Gamma^\mu_{\nu\rho} \frac{dx^\nu}{d\tau} \frac{dx^\rho}{d\tau} = 0$              ║
║                                                                              ║
║  Effective potential for radial motion:                                      ║
║    $V_{eff}(r) = -\frac{M}{r} + \frac{L^2}{2r^2} - \frac{ML^2}{r^3}$                    ║
║                                                                              ║
║  The last term is the GR correction responsible for perihelion precession.  ║
╚══════════════════════════════════════════════════════════════════════════════╝
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
import base64
from io import BytesIO

# Constants (geometrized units: G = c = 1)
M = 1.0

def geodesic_equations(tau, y):
    """
    Geodesic equations for Schwarzschild metric
    y = [t, r, theta, phi, dt/dtau, dr/dtau, dtheta/dtau, dphi/dtau]
    """
    t, r, theta, phi, t_dot, r_dot, theta_dot, phi_dot = y

if r <= 2*M:
        return [0]*8  # Inside horizon

f = 1 - 2*M/r

# Christoffel symbols
    Gamma_t_tr = M / (r * (r - 2*M))
    Gamma_r_tt = M * (r - 2*M) / r**3
    Gamma_r_rr = -M / (r * (r - 2*M))
    Gamma_r_thth = -(r - 2*M)
    Gamma_r_phiphi = -(r - 2*M) * np.sin(theta)**2
    Gamma_th_rth = 1/r
    Gamma_th_phiphi = -np.sin(theta) * np.cos(theta)
    Gamma_phi_rphi = 1/r
    Gamma_phi_thphi = np.cos(theta) / np.sin(theta) if np.sin(theta) != 0 else 0

# Geodesic equations: d²x^μ/dτ² + Γ^μ_νρ (dx^ν/dτ)(dx^ρ/dτ) = 0
    t_ddot = -2 * Gamma_t_tr * t_dot * r_dot
    r_ddot = -Gamma_r_tt * t_dot**2 - Gamma_r_rr * r_dot**2 - Gamma_r_thth * theta_dot**2 - Gamma_r_phiphi * phi_dot**2
    theta_ddot = -2 * Gamma_th_rth * r_dot * theta_dot - Gamma_th_phiphi * phi_dot**2
    phi_ddot = -2 * Gamma_phi_rphi * r_dot * phi_dot - 2 * Gamma_phi_thphi * theta_dot * phi_dot

return [t_dot, r_dot, theta_dot, phi_dot, t_ddot, r_ddot, theta_ddot, phi_ddot]

def effective_potential(r, L, M=1):
    """Effective potential for massive particle"""
    return -M/r + L**2/(2*r**2) - M*L**2/r**3

# Initial conditions for bound orbit
r0 = 10.0 * M
theta0 = np.pi / 2
phi0 = 0

# Set up for slightly eccentric orbit
L = 4.0 * M  # Angular momentum
E = 0.97     # Energy < 1 for bound orbit

# Calculate initial velocities from conserved quantities
f0 = 1 - 2*M/r0
t_dot0 = E / f0
phi_dot0 = L / r0**2

# dr/dtau from normalization: g_μν u^μ u^ν = -1
# -f*t_dot^2 + (1/f)*r_dot^2 + r^2*theta_dot^2 + r^2*sin^2(theta)*phi_dot^2 = -1
r_dot0 = np.sqrt(f0 * (E**2/f0 - 1 - L**2/r0**2))  # Positive for outward motion

y0 = [0, r0, theta0, phi0, t_dot0, r_dot0, 0, phi_dot0]

# Integrate
tau_span = [0, 500]
tau_eval = np.linspace(*tau_span, 5000)

sol = solve_ivp(geodesic_equations, tau_span, y0, t_eval=tau_eval, method='RK45', rtol=1e-10)

t, r, theta, phi = sol.y[0], sol.y[1], sol.y[2], sol.y[3]

print("╔══════════════════════════════════════════════════════════════════════╗")
print("║     Geodesic Integration Results                                     ║")
print("╚══════════════════════════════════════════════════════════════════════╝")
print()
print(f"Initial conditions:")
print(f"  $r_0 = {r0:.2f}M$, $E = {E}$, $L = {L}M$")
print()
print(f"Orbital extrema:")
print(f"  $r_{{min}} = {r.min():.4f}M$")
print(f"  $r_{{max}} = {r.max():.4f}M$")
print()

# Calculate perihelion precession
# Find perihelion passages (local minima of r)
perihelion_idx = []
for i in range(1, len(r)-1):
    if r[i] < r[i-1] and r[i] < r[i+1]:
        perihelion_idx.append(i)

if len(perihelion_idx) >= 2:
    delta_phi = phi[perihelion_idx[-1]] - phi[perihelion_idx[0]]
    n_orbits = len(perihelion_idx) - 1
    precession_per_orbit = delta_phi / n_orbits - 2*np.pi
    print(f"Perihelion precession:")
    print(f"  Measured: $\Delta\phi = {precession_per_orbit:.6f}$ rad/orbit")
    print(f"  Expected (weak field): $\Delta\phi \\approx \\frac{{6\pi M}}{{r_0}} = {6*np.pi*M/r0:.6f}$ rad/orbit")
    print()

# Plot orbit
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# r vs phi (orbital shape)
ax1 = axes[0]
x = r * np.cos(phi)
y_coord = r * np.sin(phi)
ax1.plot(x, y_coord, 'b-', linewidth=0.8, label='Geodesic orbit')
circle = plt.Circle((0, 0), 2*M, color='black', label='Event horizon (r=2M)')
ax1.add_patch(circle)
circle2 = plt.Circle((0, 0), 3*M, color='gray', fill=False, linestyle='--', linewidth=2, label='Photon sphere (r=3M)')
ax1.add_patch(circle2)
circle3 = plt.Circle((0, 0), 6*M, color='orange', fill=False, linestyle=':', linewidth=2, label='ISCO (r=6M)')
ax1.add_patch(circle3)
ax1.set_aspect('equal')
ax1.set_xlabel('x/M', fontsize=12)
ax1.set_ylabel('y/M', fontsize=12)
ax1.set_title('Geodesic Orbit in x-y Plane', fontsize=14, fontweight='bold')
ax1.legend(loc='upper right')
ax1.grid(True, alpha=0.3)

# Effective potential
ax2 = axes[1]
r_pot = np.linspace(3*M, 20*M, 500)
V_eff = [effective_potential(ri, L) for ri in r_pot]
ax2.plot(r_pot/M, V_eff, 'b-', linewidth=2, label=f'$V_{{eff}}$ (L={L}M)')
ax2.axhline(y=0.5*(E**2 - 1), color='r', linestyle='--', linewidth=2, label=f'E = {E}')
ax2.axvline(x=r.min()/M, color='gray', linestyle=':', alpha=0.7, linewidth=2)
ax2.axvline(x=r.max()/M, color='gray', linestyle=':', alpha=0.7, linewidth=2)
ax2.fill_between([r.min()/M, r.max()/M], ax2.get_ylim()[0], ax2.get_ylim()[1], alpha=0.2, color='green', label='Allowed region')
ax2.set_xlabel('r/M', fontsize=12)
ax2.set_ylabel('$V_{eff}$', fontsize=12)
ax2.set_title('Effective Potential Analysis', fontsize=14, fontweight='bold')
ax2.legend(loc='upper right')
ax2.grid(True, alpha=0.3)

plt.tight_layout()

# Convert plot to base64
buf = BytesIO()
plt.savefig(buf, format='png', dpi=150, bbox_inches='tight')
buf.seek(0)
img_base64 = base64.b64encode(buf.read()).decode('utf-8')
buf.close()
plt.close()

print(f'<img src="data:image/png;base64,{img_base64}" alt="Geodesic Orbit" />')

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Fortran: Null Geodesics (Light Rays)

Fortran

null_geodesics.f90139 lines

! ╔══════════════════════════════════════════════════════════════════════════════╗
! ║                   Null Geodesics in Schwarzschild Spacetime                  ║
! ║                                                                              ║
! ║  For light rays in the equatorial plane (θ = π/2):                          ║
! ║    Impact parameter: $b = \frac{L}{E}$                                        ║
! ║    Critical value: $b_{crit} = 3\sqrt{3}M ≈ 5.196M$                            ║
! ║                                                                              ║
! ║  Deflection angle in weak field limit:                                      ║
! ║    $\Delta\phi \approx \frac{4M}{b}$ (for $b \gg M$)                                  ║
! ╚══════════════════════════════════════════════════════════════════════════════╝

program null_geodesics
    implicit none
    integer, parameter :: dp = selected_real_kind(15)
    real(dp), parameter :: M = 1.0_dp
    real(dp), parameter :: pi = 4.0_dp * atan(1.0_dp)

real(dp) :: r, phi, r_dot, phi_dot
    real(dp) :: b, lambda, dlambda
    real(dp) :: r_min, deflection
    integer :: i, n_steps
    real(dp) :: y(4), dy(4)

! Impact parameter b = L/E (in units of M)
    b = 5.0_dp * M  ! > 3√3 M ≈ 5.196 M so ray escapes

print '(A)', '╔══════════════════════════════════════════════════════════════════╗'
    print '(A)', '║     Null Geodesics (Light Rays) in Schwarzschild Spacetime      ║'
    print '(A)', '╚══════════════════════════════════════════════════════════════════╝'
    print '(A)', ''
    print '(A,F8.4,A)', 'Impact parameter: $b = ', b, 'M$'
    print '(A,F8.4,A)', 'Critical value: $b_{crit} = ', 3.0_dp*sqrt(3.0_dp)*M, 'M$'
    print '(A)', ''
    print '(A)', 'For $b > b_{crit}$, light ray deflects but escapes to infinity.'
    print '(A)', 'For $b < b_{crit}$, light ray spirals into the black hole.'
    print '(A)', ''

! Initial conditions: light coming from infinity
    r = 100.0_dp * M
    phi = pi  ! Coming from left

! From b = r²(dφ/dλ)/(dt/dλ) and null condition
    ! For large r: dr/dλ ≈ -1 (incoming), dφ/dλ = b/r²
    r_dot = -1.0_dp
    phi_dot = b / r**2

y = [r, phi, r_dot, phi_dot]

n_steps = 100000
    dlambda = 0.01_dp
    r_min = r

! Integrate until r increases again (ray escaping)
    do i = 1, n_steps
        call rk4_step(y, dlambda, dy)
        y = y + dy

r = y(1)
        phi = y(2)

if (r < r_min) r_min = r

! Stop if ray escapes to large r
        if (r > 100.0_dp * M .and. y(3) > 0) exit

! Stop if ray hits horizon
        if (r < 2.1_dp * M) then
            print '(A)', 'Ray captured by black hole!'
            stop
        end if
    end do

print '(A)', '──────────────────────────────────────────────────────────────────'
    print '(A)', 'Results:'
    print '(A)', '──────────────────────────────────────────────────────────────────'
    print '(A,F12.6,A)', 'Minimum approach: $r_{min} = ', r_min, 'M$'
    print '(A)', ''

! Deflection angle
    deflection = abs(phi - pi) - pi

print '(A,F12.6,A)', 'Final angle: $\phi = ', phi, '$ rad'
    print '(A,F12.6,A)', 'Deflection: $\Delta\phi = ', deflection, '$ rad'
    print '(A,F12.6,A)', '           = ', deflection * 180.0_dp / pi, '°'
    print '(A)', ''
    print '(A)', 'Weak field approximation:'
    print '(A,F12.6,A)', '  $\Delta\phi \approx \frac{4M}{b} = ', 4.0_dp*M/b, '$ rad'
    print '(A,F12.6,A)', '  Error: ', abs(deflection - 4.0_dp*M/b)/deflection * 100, '%'

contains

subroutine geodesic_deriv(y, dydt)
        real(dp), intent(in) :: y(4)
        real(dp), intent(out) :: dydt(4)
        real(dp) :: r, phi, rdot, phidot, f

r = y(1)
        phi = y(2)
        rdot = y(3)
        phidot = y(4)

f = 1.0_dp - 2.0_dp*M/r

! dr/dlambda = rdot
        dydt(1) = rdot

! dphi/dlambda = phidot
        dydt(2) = phidot

! d(rdot)/dlambda = -Gamma^r terms
        ! For null geodesic in equatorial plane:
        dydt(3) = -M*(r-2*M)/r**3 + (r-2*M)*phidot**2

! d(phidot)/dlambda = -2*Gamma^phi_r_phi * rdot * phidot
        dydt(4) = -2.0_dp/r * rdot * phidot

end subroutine

subroutine rk4_step(y, h, dy)
        real(dp), intent(in) :: y(4), h
        real(dp), intent(out) :: dy(4)
        real(dp) :: k1(4), k2(4), k3(4), k4(4), ytmp(4)

call geodesic_deriv(y, k1)

ytmp = y + 0.5_dp*h*k1
        call geodesic_deriv(ytmp, k2)

ytmp = y + 0.5_dp*h*k2
        call geodesic_deriv(ytmp, k3)

ytmp = y + h*k3
        call geodesic_deriv(ytmp, k4)

dy = h * (k1 + 2*k2 + 2*k3 + k4) / 6.0_dp

end subroutine

end program null_geodesics

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

← Parallel Transport Next: Part II - Curvature →