Orbits around a black hole¶
You may have seen the image of a black hole. This image comes from photons that are emitted from an accretion disk or jet around the black hole and then travel across space-time to hit our telescope detectors.
How do we understand these orbits around a black hole? This notebook will walk through Newtonian intuition, then build to GR (Schwarzschild, specifically), taking inspiration from Sean's Carroll's GR textbook (Ch. 5.4, Figs. 5.4, 5.5).
See the Schwarzschild Appendix in my overleaf.
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('lia.mplstyle')
from bh_accretion_toolkit.gr import horizons as hz
from bh_accretion_toolkit import plotting
import matplotlib.lines as mlines
from importlib import reload
Newtonian Intuition¶
Imagine a particle orbiting the sun -- it could be the Earth, or it could be our favorite "point particle".
Energy and angular momentum are conserved for this particle.
Angular momentum: $L = mr^2\dot \phi$. Energy: $E=\frac12 m(\dot r^2+r^2\dot \phi^2)-GmM/r$
Eliminating $\dot \phi$ using $L$, the last equation becomes $E=\frac12 m(\dot r^2 +L^2/(m^2r^2))-GMm/r$. Re-arranging, we have $$ \frac 12 m\dot r^2 = E - V_{\rm eff}(r),$$
where the effective potential is given by $$ V_{\rm eff}(r)=\frac 12 \frac{L^2}{mr^2}-\frac{GMm}{r}.$$
So if we know $V_{\rm eff}$, we can solve for $\dot r$ of the particle orbit for a particle with a given energy $E$.
Note that we include a constant offset $\epsilon/2$ to match the Schwarzschild energy convention used below, which we can do because it's just a constant offset and doesn't affect the dynamics.
For massive particles, the horizontal lines below show energies in the shifted convention (including the $1/2$ offset). Motion is allowed where $E \ge V_{\rm eff}$, and the marked intersections are turning points where $\dot r=0$. The potential minimum gives a circular orbit; a bound orbit has two turning points, while an unbound orbit ($E>1/2$) has one closest approach and can extend to infinity.
def effective_potential_newtonian(ang_mom, radius, epsilon=1):
"""
Calculate the effective potential energy for a particle in Schwarzschild geometry.
Follows Carroll Eq. 5.66.
Note: m = 1, M = 1, G = 1.
Parameters:
epsilon (float): 0 for massless particles, 1 for massive particles.
ang_mom (float): The angular momentum of the particle.
radius (float): The radial distance from the center of force.
Returns:
float: The effective potential energy at the given radius.
"""
return 0.5*epsilon - epsilon / radius + ang_mom**2 / (2 * radius**2)
# Plot the effective potential for both photon and massive particles
L_val = 2
r_range = np.linspace(0.5, 30, 200)
epsilon_vals = [0, 1]
particle_labels = ['Photon', 'Massive particle']
fig, ax = plt.subplots(1, 1)
particle_handles = []
for epsilon, label in zip(epsilon_vals, particle_labels):
V_eff = effective_potential_newtonian(L_val, r_range, epsilon)
l0, = ax.plot(r_range, V_eff, label=label)
particle_handles.append(l0)
# Energy levels and turning points for massive particles (epsilon = 1)
r_circular = L_val**2
E_circular = effective_potential_newtonian(L_val, r_circular, epsilon=1)
E_vals = [E_circular, 0.46, 0.55]
orbit_labels = ['Circular', 'Bound', 'Unbound']
orbit_linestyles = ['-.', '--', ':']
orbit_handles = []
for E_val, label, linestyle in zip(E_vals, orbit_labels, orbit_linestyles):
l0 = ax.axhline(E_val, color='k', ls=linestyle, lw=1, label=r"{}: $E={:.2f}$".format(label, E_val))
orbit_handles.append(l0)
# Solve E = V_eff; the circular orbit has one repeated root.
eccentricity = np.sqrt(max(0., 1. + 2.*(E_val - 0.5)*L_val**2))
turning_radii = [L_val**2 / (1. + eccentricity)]
if 0. < eccentricity < 1.:
turning_radii.append(L_val**2 / (1. - eccentricity))
ax.plot(turning_radii, [E_val]*len(turning_radii), 'o', color=l0.get_color(), markersize=5)
ax.set_xlabel(r"Radius $r/r_g$")
ax.set_ylabel(r"Effective potential $V_{\rm eff}(r)$")
ax.set_xlim(0, r_range[-1])
ax.set_ylim(0., 0.8)
particle_handles.append(mlines.Line2D([], [], color='k', marker='o', ls='none', markersize=5, label='Turning Point'))
particle_leg = ax.legend(handles=particle_handles, loc='upper right', framealpha=1)
ax.add_artist(particle_leg)
ax.legend(handles=orbit_handles[::-1], title='Orbit Type', loc='lower right', framealpha=1)
ax.grid()
plt.title(r"Newtonian potential, $L={}$".format(L_val))
plt.tight_layout()
Schwarzschild Functions¶
In the Schwarzschild metric, the radial energy equation is $$ \frac12\dot r^2 = \mathcal{E} - V_{\rm eff}(r), \qquad \mathcal{E}\equiv\frac12 E^2, $$ where $E$ is the conserved energy obtained from the timelike Killing vector. The effective potential is $$ V_{\rm eff}(r)= \frac12\epsilon -\epsilon\frac{GM}r+\frac{L^2}{2r^2}-\frac{GML^2}{r^3}, $$ where $\epsilon=1$ for massive particles and $\epsilon=0$ for photons; see Carroll Eq. 5.66.
The extra term $-GML^2/r^3$ strengthens the inward pull at small radii, overcoming the centrifugal barrier and allowing particles to plunge into the black hole.
def effective_potential_schwarzschild(ang_mom, radius, epsilon=1, gamma=1):
"""
Calculate the effective potential energy for a particle in Schwarzschild geometry.
Follows Carroll Eq. 5.66.
Parameters:
epsilon (float): 0 for massless particles, 1 for massive particles.
ang_mom (float): The angular momentum of the particle.
radius (float): The radial distance from the center of force.
gamma (float): 0 for Newtonian, 1 for GR.
Returns:
float: The effective potential energy at the given radius.
"""
return effective_potential_newtonian(ang_mom, radius, epsilon) - ang_mom**2 / (radius**3)*gamma
Larger angular momentum $L$ raises the centrifugal barrier outside the black hole, so an incoming particle needs a higher energy to cross it and plunge inward. Unlike the Newtonian case, the Schwarzschild barrier has a finite peak. At smaller radii, the negative $-GML^2/r^3$ term wins and the potential falls toward the horizon.
L_vals = [3, 4, 5]
r_range = np.linspace(hz.get_EH_location(0), 30, 200)
L_handles = []
for L_val in L_vals:
l0, = plt.plot(r_range, effective_potential_schwarzschild(L_val, r_range, epsilon=1), label=r"$L={:.0f}$".format(L_val))
L_handles.append(l0)
plt.plot(r_range, effective_potential_schwarzschild(L_val, r_range, epsilon=0), color=l0.get_color(), ls='--')
plt.axvline(hz.get_ISCO_location(0), color='k', ls=':', label='ISCO')
plt.ylabel(r"Effective potential $V_{\rm eff}(r)$")
L_leg = plt.legend(handles=L_handles, title="Angular momentum", loc='upper right')
plt.gca().add_artist(L_leg)
particle_handles = [
mlines.Line2D([], [], color='k', ls='-', label='Massive'),
mlines.Line2D([], [], color='k', ls='--', label='Massless'),
]
plt.legend(handles=particle_handles, loc='lower right')
plt.xlabel(r"$r/r_g$")
plt.title("Schwarzschild spacetime")
plt.ylim(0, 1)
(0.0, 1.0)
Circular Orbits¶
Circular orbits occur where $dV/dr=0$. Instead of doing the derivative numerically, we take the analytic expressions from Carroll Eq. 5.71.
def square_root(ang_mom):
return np.sqrt(ang_mom**4 - 12.*ang_mom**2)
def stable_circular_orbit_radius(ang_mom):
"""
Calculate the radius of a stable circular orbit for a particle with given angular momentum.
This is for massive particles in Schwarzschild geometry.
Positive root in Carroll Eq. 5.72.
Parameters:
ang_mom (float): The angular momentum of the particle.
Returns:
float: The radius of the stable circular orbit.
"""
return (ang_mom**2 + square_root(ang_mom)) / 2.
def unstable_circular_orbit_radius(ang_mom):
"""
Calculate the radius of a stable circular orbit for a particle with given angular momentum.
This is for massive particles in Schwarzschild geometry.
Negative root in Carroll Eq. 5.72.
Parameters:
ang_mom (float): The angular momentum of the particle.
Returns:
float: The radius of the stable circular orbit.
"""
return (ang_mom**2 - square_root(ang_mom)) / 2.
Reproduce Carroll Figure 5.5¶
To make sure everything is working, here I reproduce the right panel of Carroll Fig. 5.5.
# Reproduce Carroll Figure 5.5
r_range = np.linspace(0, 30, 200)
L_vals_carroll = [1, 2, 3, 4, 5]
# Value for L where stable and unstable circular orbits meet
L_critical = np.sqrt(12)
fig, ax = plt.subplots(1, 1)
for L in L_vals_carroll:
V_eff = effective_potential_schwarzschild(L, r_range, 1, 1)
l0, = ax.plot(r_range, V_eff, label=r'${:.0f}$'.format(L))
stable_r = stable_circular_orbit_radius(L)
ax.plot(stable_r, effective_potential_schwarzschild(L, stable_r, 1, 1), 's', fillstyle='full', color=l0.get_color(), markersize=8)
# print(f"Stable circular orbit for L={L} at r={stable_r:.2f}")
unstable_r = unstable_circular_orbit_radius(L)
ax.plot(unstable_r, effective_potential_schwarzschild(L, unstable_r, 1, 1), 'o', color=l0.get_color(), markersize=8)
# print(f"Unstable circular orbit for L={L} at r={unstable_r:.2f}")
# Plot stable and unstable circular orbit locations
# labelLines(axs[0].get_lines(), zorder=2.5)
L_single = ax.plot(r_range, effective_potential_schwarzschild(L_critical, r_range, 1, 1), label=r'$\sqrt{12}$', color='k', ls='--')
L_leg = ax.legend(loc='lower right', framealpha=1, title="Angular Momentum $L$", ncol=3)
ax.plot(unstable_circular_orbit_radius(L_critical), effective_potential_schwarzschild(L_critical, unstable_circular_orbit_radius(L_critical), 1, 1), 's', color='k', markersize=8)
ax.plot(stable_circular_orbit_radius(L_critical), effective_potential_schwarzschild(L_critical, stable_circular_orbit_radius(L_critical), 1, 1), 'o', color='k', markersize=8)
# radius_leg = axs[0].legend(handles=[L_single[0]], loc='lower right', framealpha=1, title="Single $L$ Value", ncol=1)
ISCO = ax.axvline(6., color='r', ls='--', label=r"Innermost {\it stable} circular orbit")
circular = ax.axvline(3., color='k', ls='-.', label=r"Innermost {\it unstable} circular orbit")
ax.set_xlabel(r"Radius $r/r_g$")
ax.grid()
ax.set_ylabel(r"Effective Potential $V_{\rm eff}(r)$")
ax.set_ylim(0., 0.8)
ax.set_xlim(0, r_range[-1])
plt.title("Massive Particles in Schwarzschild Spacetime")
plt.tight_layout()
/tmp/ipykernel_2792977/471119728.py:15: RuntimeWarning: divide by zero encountered in divide return 0.5*epsilon - epsilon / radius + ang_mom**2 / (2 * radius**2) /tmp/ipykernel_2792977/471119728.py:15: RuntimeWarning: invalid value encountered in add return 0.5*epsilon - epsilon / radius + ang_mom**2 / (2 * radius**2) /tmp/ipykernel_2792977/2932701246.py:2: RuntimeWarning: invalid value encountered in sqrt return np.sqrt(ang_mom**4 - 12.*ang_mom**2)
Exploration/Clarification¶
Just from looking at Carroll's Fig. 5.5, it was not clear to me that there were actually bound orbits, i.e. energies that particles could have and not escape from the potential well. So, the right panel here zooms in and plots values of angular momentum that show there are indeed bound orbits for $L>\sqrt{12}$.
r_range = np.linspace(0, 30, 200)
L_vals_carroll = [2, 3, 4, 5, 8]
colors = plt.rcParams['axes.prop_cycle'].by_key()['color']
# Value for L where stable and unstable circular orbits meet
L_critical = np.sqrt(12)
fig, axs = plt.subplots(1, 2, figsize=(12, 5))
for L in L_vals_carroll:
V_eff = effective_potential_schwarzschild(L, r_range, 1, 1)
if L < L_critical:
color = colors[0]
else:
color = colors[2]
l0, = axs[0].plot(r_range, V_eff, label=r'${:.0f}$'.format(L), color=color)
stable_r = stable_circular_orbit_radius(L)
axs[0].plot(stable_r, effective_potential_schwarzschild(L, stable_r, 1, 1), 's', fillstyle='full', color=l0.get_color(), markersize=8)
# print(f"Stable circular orbit for L={L} at r={stable_r:.2f}")
unstable_r = unstable_circular_orbit_radius(L)
axs[0].plot(unstable_r, effective_potential_schwarzschild(L, unstable_r, 1, 1), 'o', color=l0.get_color(), markersize=8)
# print(f"Unstable circular orbit for L={L} at r={unstable_r:.2f}")
# Plot stable and unstable circular orbit locations
# labelLines(axs[0].get_lines(), zorder=2.5)
L_single = axs[0].plot(r_range, effective_potential_schwarzschild(L_critical, r_range, 1, 1), label=r'$L_{\rm crit}=\sqrt{12}$', color=colors[1])
L_handles = [
mlines.Line2D([], [], color=colors[2], linestyle='-', label=r"$L>L_{\rm crit}$"),
L_single[0],
mlines.Line2D([], [], color=colors[0], linestyle='-', label=r"$L<L_{\rm crit}$"),
]
L_leg = axs[0].legend(handles=L_handles, loc='upper right', framealpha=1)
axs[0].plot(unstable_circular_orbit_radius(L_critical), effective_potential_schwarzschild(L_critical, unstable_circular_orbit_radius(L_critical), 1, 1), 's', color=colors[1], markersize=8)
axs[0].plot(stable_circular_orbit_radius(L_critical), effective_potential_schwarzschild(L_critical, stable_circular_orbit_radius(L_critical), 1, 1), 'o', color=colors[1], markersize=8)
axs[0].set_ylabel(r"Effective Potential $V_{\rm eff}(r)$")
axs[0].set_ylim(0., 1.4)
axs[0].set_xlim(0, r_range[-1])
# radius_leg = axs[0].legend(handles=[L_single[0]], loc='lower right', framealpha=1, title="Single $L$ Value", ncol=1)
L_vals_illustrative = [3, L_critical, 4]
for L in L_vals_illustrative:
V_eff = effective_potential_schwarzschild(L, r_range, 1, 1)
l0, = axs[1].plot(r_range, V_eff, label=r'${:.0f}$'.format(L))
stable_r = stable_circular_orbit_radius(L)
axs[1].plot(stable_r, effective_potential_schwarzschild(L, stable_r, 1, 1), 's', color=l0.get_color(), markersize=5)
# print(f"Stable circular orbit for L={L} at r={stable_r:.2f}")
unstable_r = unstable_circular_orbit_radius(L)
axs[1].plot(unstable_r, effective_potential_schwarzschild(L, unstable_r, 1, 1), 'o', color=l0.get_color(), markersize=5)
# print(f"Unstable circular orbit for L={L} at r={unstable_r:.2f}")
for ax in axs:
ISCO = ax.axvline(6., color='r', ls='--', label=r"Innermost {\it stable} circular orbit")
circular = ax.axvline(3., color='k', ls='-.', label=r"Innermost {\it unstable} circular orbit")
ax.set_xlabel(r"Radius $r/r_g$")
ax.grid()
# plt.legend(handles=[L_single[0], ISCO, circular], loc='best', framealpha=1)
axs[1].set_ylim(0.4, 0.52)
axs[1].set_xlim(2, 30)
rc_handles = [
mlines.Line2D([], [], marker='s', color='k', ls='none', label="Stable Circular Orbit"),
mlines.Line2D([], [], marker='o', color='k', ls='none', label="Unstable Circular Orbit"),
]
rc_leg = axs[1].legend(handles=rc_handles, loc='lower right', framealpha=1, ncol=1)
r_handles = [
mlines.Line2D([], [], color='k', ls='-.', label="IUCO"),
mlines.Line2D([], [], color='r', ls='--', label="ISCO"),
]
r_leg = axs[1].legend(handles=r_handles, loc='upper right', framealpha=1, ncol=1)
axs[1].add_artist(rc_leg)
plt.suptitle("Massive Particles in Schwarzschild Spacetime")
plt.tight_layout()
plotting.save_fig("schwarzschild_effective_potential", "Z_effective_potential")
/tmp/ipykernel_2792977/471119728.py:15: RuntimeWarning: divide by zero encountered in divide return 0.5*epsilon - epsilon / radius + ang_mom**2 / (2 * radius**2) /tmp/ipykernel_2792977/471119728.py:15: RuntimeWarning: invalid value encountered in add return 0.5*epsilon - epsilon / radius + ang_mom**2 / (2 * radius**2) /tmp/ipykernel_2792977/2932701246.py:2: RuntimeWarning: invalid value encountered in sqrt return np.sqrt(ang_mom**4 - 12.*ang_mom**2)
[PosixPath('figures/Z_effective_potential/schwarzschild_effective_potential.pdf'),
PosixPath('figures/Z_effective_potential/schwarzschild_effective_potential.png')]
Kerr Metric¶
See Bardeen, Press, Teukolsky 1972, Eq. 2.9a and 2.10.
For equatorial orbits ($Q=0$), $r^4\,(dr/d\tau)^2 = R(r; E, L)$, and $R$ depends on $E$ itself, so it can't be plotted as a fixed potential.
Instead, solve $R=0$ for $E$: the larger root $E_+(r; L)$ is the minimum energy allowed at radius $r$, so a particle with energy $E$ can only reach radii where $E \ge E_+(r)$.
For $a=0$, $\tfrac12 E_+^2$ is exactly Carroll's $V_{\rm eff}$ above.
See derivations/kerr_effective_potential.md for the full derivation.
Note that for retrograde orbits ($aL<0$), $E_+$ becomes negative close to the horizon: these are the negative-energy orbits inside the ergosphere that make the Penrose process possible.
from bh_accretion_toolkit.gr import orbits as cq
L = 4.
spins = [-0.9, 0., 0.5, 0.9]
fig, ax = plt.subplots(1, 1)
spin_handles = []
for spin in spins:
r_range = np.linspace(hz.get_EH_location(spin), 30, 500)
E_plus = cq.effective_potential_energy(L, r_range, spin)
l0, = ax.plot(r_range, E_plus, label=r'$a={:+.1f}$'.format(spin))
spin_handles.append(l0)
ax.axvline(hz.get_ISCO_location(spin), color=l0.get_color(), ls=':', lw=1)
ax.axhline(1., color='gray', lw=1) # E = 1: marginally bound
ax.set_xlabel(r"Radius $r/r_g$")
ax.set_ylabel(r"Minimum allowed energy $E_+(r)$")
ax.set_xlim(1., 20.)
ax.set_ylim(-0.5, 1.5)
ax.grid()
spin_leg = ax.legend(handles=spin_handles, loc='lower right', framealpha=1, title="Spin")
ax.add_artist(spin_leg)
reference_handles = [
mlines.Line2D([], [], color='k', ls=':', lw=1, label='ISCO'),
]
ax.legend(handles=reference_handles, loc='upper right', framealpha=1)
plt.title(r"Massive Particles in Kerr Spacetime, $L={:.0f}$".format(L))
plt.tight_layout()
References¶
- Carroll, Spacetime and Geometry, Sec. 5.4
- Bardeen, Press & Teukolsky 1972, ApJ 178, 347
- Misner, Thorne & Wheeler, Sec. 33.5