Geodésicas de Schwarzschild: órbitas y luz alrededor de un agujero negro

1. La métrica de Schwarzschild

Trabajamos en unidades geométricas $G = c = 1$: una masa $M$ equivale a la longitud $GM/c^2$ (para el Sol, $GM_\odot/c^2 \approx 1.477\ \text{km}$) y al tiempo $GM/c^3$. La solución de vacío, estática y esféricamente simétrica de las ecuaciones de Einstein es, en coordenadas $(t, r, \theta, \varphi)$,

donde $M$ es la masa del cuerpo central y $r$ es la coordenada areal (la esfera $r = \text{cte}$ tiene área $4\pi r^2$). En $r = 2M$ se anula $g_{tt}$: es el horizonte de sucesos (radio de Schwarzschild $r_s = 2GM/c^2$). La singularidad de $g_{rr}$ allí es de coordenadas; la de curvatura está en $r = 0$.

Una partícula de prueba sigue una geodésica $x^\alpha(\tau)$, parametrizada por el tiempo propio $\tau$ si es masiva (con $u^\alpha = dx^\alpha/d\tau$ y $u\cdot u = -1$) o por un parámetro afín $\lambda$ si es luz ($u\cdot u = 0$). Por simetría esférica podemos restringirnos al plano ecuatorial $\theta = \pi/2$ sin pérdida de generalidad.

2. Cantidades conservadas a partir de los vectores de Killing

La métrica no depende de $t$ ni de $\varphi$, así que $\xi = \partial_t$ y $\eta = \partial_\varphi$ son vectores de Killing. Para todo vector de Killing $K$ y toda geodésica, $K\cdot u$ es constante:

Definimos entonces la energía y el momento angular por unidad de masa:

Para $r\to\infty$, $E \to dt/d\tau = \gamma$: es la energía total (incluida la de reposo) por unidad de masa, y $E < 1$ caracteriza las órbitas ligadas.

3. Ecuación radial y potencial efectivo

Sustituyendo $dt/d\tau$ y $d\varphi/d\tau$ en la normalización $g_{\alpha\beta}u^\alpha u^\beta = -1$:

Los tres primeros términos reproducen el problema newtoniano (con $\tfrac12(V-1)$ como energía potencial efectiva por unidad de masa: $-M/r + L^2/2r^2$). El término $-2ML^2/r^3$ es nuevo: domina a radios pequeños y hace que, en lugar de crecer sin límite como en Newton, $V$ alcance un máximo y luego caiga hasta $V(2M) = 0$: la barrera centrífuga es finita y una partícula con $E^2$ por encima de ella cae al agujero.

Derivando respecto a $\tau$ se obtiene la forma de segundo orden, que es la que integra el simulador:

Órbitas circulares y la ISCO

Los extremos de $V$ cumplen $dV/dr = 0$, es decir

$r_+$ es un mínimo (órbita circular estable) y $r_-$ un máximo (circular inestable). Despejando $L$ en la ecuación cuadrática obtenemos, para una órbita circular de radio $r$,

que solo existen para $r > 3M$. Las dos raíces se funden cuando $L^2 = 12M^2$, en

la órbita circular estable más interna (ISCO). Allí $E = \sqrt{8/9}$, de modo que el gas de un disco de acreción que espiralea lentamente hasta la ISCO y cae libera una fracción $1-\sqrt{8/9} \approx 5.7\,\%$ de su energía en reposo, frente al $0.7\,\%$ de la fusión del hidrógeno. Para $r$ entre $3M$ y $6M$ las órbitas circulares existen pero son inestables: cualquier perturbación las lleva al horizonte o al infinito.

Frecuencia y periodo. La frecuencia angular de una órbita circular medida en tiempo coordenado $t$ (el de un observador lejano) es, con $L = r^2\,d\varphi/d\tau$ y $E = (1-2M/r)\,dt/d\tau$,

exactamente la tercera ley de Kepler, un hecho notable de la coordenada areal $r$ de Schwarzschild. El periodo de la ISCO es $T_{\rm ISCO} = 2\pi\sqrt{216}\,M = 92.3\,GM/c^3$: con $GM_\odot/c^3 = 4.925\,\mu$s,

La variabilidad más rápida posible de un disco de acreción está fijada por $T_{\rm ISCO}$: las oscilaciones cuasiperiódicas de cientos de Hz en binarias de rayos X y las llamaradas de ~30 min de Sgr A* (GRAVITY) encajan con estas cifras. La velocidad orbital medida por un observador estático local en la ISCO es $v = \sqrt{M/(r-2M)} = c/2$. Por último, la ISCO depende del espín: en la métrica de Kerr se acerca a $M$ (prógrada, espín máximo) y la eficiencia de acreción sube hasta el 42 % (Bardeen, Press & Teukolsky 1972), lo que explica la luminosidad de los cuásares.

4. La ecuación de Binet relativista

Para la forma de la órbita usamos $u = 1/r$ y $\varphi$ como variable independiente. Como $dr/d\tau = (dr/d\varphi)(L/r^2) = -L\,du/d\varphi$, la ecuación radial queda

Derivando respecto a $\varphi$ y dividiendo por $2\,du/d\varphi$ (para órbitas no circulares):

Sin el término $3Mu^2$ es la ecuación de Binet newtoniana, cuya solución es la cónica $u = (M/L^2)(1+e\cos\varphi)$. En unidades SI el cociente entre el término nuevo y $M/L^2$ es $3L^2u^2/c^2 = 3(v_\perp/c)^2$: diminuto en el Sistema Solar, de orden unidad cerca de un agujero negro.

5. Precesión del periastro

Escribimos $u = u_0 + u_1$ con $u_0 = (M/L^2)(1+e\cos\varphi)$ y tratamos $3Mu^2 \approx 3Mu_0^2$ como perturbación. Usando $\cos^2\varphi = \tfrac12(1+\cos2\varphi)$:

El término constante y el de $\cos2\varphi$ dan correcciones acotadas y periódicas. El término en $\cos\varphi$ está en resonancia con el oscilador $u'' + u$: la solución particular de $u'' + u = A\cos\varphi$ es $\tfrac{A}{2}\varphi\sin\varphi$ (compruébalo), un término secular que crece con $\varphi$. Con $A = 6M^3e/L^4$:

donde hemos usado $\cos(\varphi - \epsilon\varphi) \approx \cos\varphi + \epsilon\varphi\sin\varphi$ para $\epsilon\varphi \ll 1$. El periastro (máximo de $u$) se repite cuando $(1-3M^2/L^2)\varphi = 2\pi n$, de modo que avanza, por órbita,

donde usamos $L^2 \simeq Ma(1-e^2)$, la relación newtoniana entre el momento angular y el semilatus rectum.

Mercurio

Con $GM_\odot/c^2 = 1476.6\ \text{m}$, $a = 5.791\times10^{10}\ \text{m}$ y $e = 0.2056$:

En un siglo juliano hay $36\,525/87.969 = 415.2$ órbitas, así que la precesión relativista es $\approx 42.98''$ por siglo. Este es exactamente el exceso que Le Verrier (1859) no pudo atribuir a las perturbaciones de los demás planetas, y su explicación por Einstein (1915) fue la primera prueba de la relatividad general. (El avance total observado, unas $5600''$ por siglo, se debe sobre todo a la precesión del equinoccio del sistema de referencia y a los demás planetas.)

Comprobación numérica

En src/lib/math/schwarzschild.ts, precesionPorOrbita integra la ecuación de Binet con RK4 en $\varphi$ desde el periastro ($u' = 0$) hasta el siguiente cambio de signo $u' : + \to -$, localizado por interpolación lineal. Como $u''' = -u'(1-6Mu)$ se anula en el periastro, el error de la interpolación es $O(h^3)$ y el del RK4, $O(h^4)$ global. Con $h = 10^{-3}$ rad, para Mercurio la integración y la fórmula de primer orden difieren en $3\times10^{-5}$ en términos relativos: la discrepancia la dominan los términos $O(M/p)$, con $p = a(1-e^2)$, que la fórmula ignora. Esos términos crecen al acercarse al agujero: para $a = 10^3M$, $e = 0.3$ ($p = 910M$) la diferencia ya es del $0.5\,\%$; para $a = 100M$, $e = 0.5$, del $6.5\,\%$; y cerca de la ISCO la precesión por órbita supera los $100^\circ$ y la órbita dibuja una roseta. Compruébalo en el simulador.

En el modo partícula la partícula parte de $r_0$ con $dr/d\tau = 0$, de modo que $E^2 = V(r_0)$. Fíjate en el panel inferior: la órbita está ligada si la línea $E^2$ corta a $V(r)$ en dos puntos fuera de la barrera; si $L < \sqrt{12}\,M$ no hay barrera y toda partícula termina cruzando el horizonte. La órbita newtoniana (mismos $r_0$ y $L$) se cierra sobre sí misma; la de Schwarzschild no.

6. La luz: esfera de fotones y parámetro de impacto crítico

Para la luz $u\cdot u = 0$ y $E$, $L$ se definen igual con un parámetro afín $\lambda$. Repitiendo el paso 3:

donde $b$ es el parámetro de impacto (distancia del rayo incidente a la recta paralela que pasa por el centro). $W'(r) = 0$ da un único máximo en $r = 3M$, la esfera de fotones, con $W(3M) = 1/(27M^2)$. Es una órbita circular inestable de la luz. Un fotón llegando desde el infinito:

• con $b > b_c$ rebota en la barrera en el radio $r_0$, la raíz mayor de $r^3 - b^2r + 2Mb^2 = 0$, $\;r_0 = \tfrac{2b}{\sqrt3}\cos\!\left[\tfrac13\arccos\!\left(-3\sqrt3 M/b\right)\right]$;
• con $b < b_c$ supera la barrera y cae al agujero, donde

Propiedades de la esfera de fotones. La frecuencia angular coordenada es $\Omega_c = \sqrt{M/27M^3} = 1/(3\sqrt3M)$ (la misma ley de Kepler), con periodo $2\pi\cdot3\sqrt3\,M = 32.6\,GM/c^3$: 0.6 ms para $10\,M_\odot$, 11 min para Sgr A*, 12 días para M87*. Un observador estático en $r = 3M$ ve el tiempo lejano dilatado por $\sqrt{1-2/3} = 0.577$. La órbita es inestable con exponente de Lyapunov $\lambda = \Omega_c$ en tiempo coordenado (Cardoso et al. 2009): una desviación radial inicial $\delta r_0$ crece como $\delta r_0e^{\lambda t}$, de modo que en cada vuelta se multiplica por $e^{2\pi}\approx535$. Esa misma inestabilidad gobierna cómo «suena» un agujero negro perturbado: las frecuencias de sus modos cuasinormales son $\omega\approx\Omega_c(\ell+\tfrac12) - i\lambda(n+\tfrac12)$ en el límite eikonal, lo que conecta la esfera de fotones con el ringdown de GW150914 (lección 4).

La ecuación de Binet de la luz se obtiene igual que la masiva, sin el término $M/L^2$ (pues $1/L^2 \to 0$ al reescalar con $\lambda$):

7. Desviación de la luz

A orden cero, $u_0 = \sin\varphi/b$: una recta a distancia $b$ del origen. A primer orden en $M/b$,

El rayo llega del infinito ($u = 0$) en $\varphi = -\delta/2$ y sale en $\varphi = \pi + \delta/2$. Imponiendo $u(-\delta/2) = 0$ con $\delta \ll 1$: $-\delta/(2b) + 2M/b^2 = 0$, es decir

el doble del valor «newtoniano» que se obtiene tratando la luz como un proyectil. Para un rayo rasante al Sol, $b = R_\odot = 6.957\times10^8\ \text{m}$:

Las expediciones de Dyson, Eddington y Davidson a Sobral y Príncipe durante el eclipse total del 29 de mayo de 1919 midieron desplazamientos estelares compatibles con $1.75''$ y no con $0.87''$; hoy la interferometría de radio (VLBI) confirma la predicción con precisión mejor que $10^{-4}$. La función deflexionLuz integra la ecuación de Binet desde el periastro $r_0$ hasta $u = 0$ y usa la simetría $\delta = 2\varphi_{1/2} - \pi$; en el limbo solar coincide con $4M/b$ con error relativo menor que $10^{-3}$, mientras que para $b \to b_c$ la deflexión diverge logarítmicamente: el fotón da vueltas cerca de la esfera de fotones antes de escapar.

8. La sombra de un agujero negro

Un observador lejano a distancia $D$ que mira hacia el agujero ve oscura toda dirección cuyo rayo, trazado hacia atrás, tiene $b < b_c$: esos rayos terminan en el horizonte. La sombra es un disco de radio angular

notablemente mayor que el horizonte ($2M$) e incluso que la esfera de fotones ($3M$), por el efecto de lente del propio agujero. Para M87* ($M \approx 6.5\times10^9 M_\odot$, $D \approx 16.8\ \text{Mpc}$) el diámetro $2\theta_{\text{sh}}$ resulta $\approx 40\ \mu\text{as}$, comparable con el anillo de $42 \pm 3\ \mu\text{as}$ fotografiado por el Event Horizon Telescope en 2019; para Sgr A* ($M \approx 4\times10^6 M_\odot$, $D \approx 8.15\ \text{kpc}$) se predicen $\approx 50\ \mu\text{as}$, frente a $51.8 \pm 2.3\ \mu\text{as}$ medidos en 2022. (El anillo observado es emisión del plasma lensada cerca de $b_c$, no la sombra misma; la rotación —métrica de Kerr— modifica $b_c$ en pocos por ciento.)

En el modo luz del simulador, el haz de rayos paralelos con $|b| = 1, 2, \dots, 12\,M$ muestra la frontera: los rayos con $|b| \le 5$ son capturados y los de $|b| \ge 6$ escapan. Lleva el rayo destacado a $b = 5.20\,M$ para ver cómo orbita varias veces la esfera de fotones.

9. Problemas

1. Órbitas circulares (demostración). Demuestra que $L^2 = Mr^2/(r-3M)$ y $E = (1-2M/r)/\sqrt{1-3M/r}$, y que no hay órbitas circulares (de ningún tipo) para $r$ menor que $3M$.

2. Estabilidad y ISCO (demostración). Demuestra que la órbita circular de radio $r$ es estable si y solo si $r$ es mayor que $6M$.

3. Precesión y energía. Usando el resultado de la sección 5, demuestra que en unidades SI la precesión por órbita es $\Delta\varphi = 24\pi^3a^2/[c^2T^2(1-e^2)]$, con $T$ el periodo kepleriano. Evalúala para Mercurio.

4. Máxima aproximación de un fotón (demostración). Demuestra que para $b$ mayor que $b_c$ el periastro $r_0$ satisface $b^2 = r_0^3/(r_0-2M)$ y que $r_0 \to 3M$ cuando $b \to b_c^+$.

5. Periodos de la ISCO. Demuestra que $\Omega = \sqrt{M/r^3}$ para toda órbita circular de Schwarzschild y calcula $T_{\rm ISCO}$ y la frecuencia de la onda gravitacional ($f_{\rm GW} = 2/T$) para un sistema de masa total $65\,M_\odot$ como GW150914.

6. Deflexión a segundo orden (computacional). Con deflexionLuz calcula $\delta(b)$ para $b \in [6, 1000]\,M$ y ajusta $\delta - 4M/b \approx C\,M^2/b^2$. Compara con el valor analítico $C = 15\pi/4$.

import { deflexionLuz } from '@/lib/math/schwarzschild';
for (const b of [6, 10, 20, 50, 100, 1000]) {
  const d = deflexionLuz(b, 1, 1e-4);
  console.log(b, (d - 4 / b) * b * b);
}

7. Análisis de error (computacional). Estudia el error de precesionPorOrbita frente al paso $h$ para $a = 100M$, $e = 0.5$ usando como referencia $h = 10^{-4}$. Comprueba en una gráfica log-log que la pendiente es $\approx 4$ (RK4) mientras el error de interpolación no domine, y explica por qué este último es $O(h^3)$.

10. Autoevaluación

11. Lecturas

• Hartle (2003), Gravity, cap. 9 (geodésicas de Schwarzschild, potencial efectivo, ISCO) y cap. 15 (agujeros negros). Schutz (2009), cap. 11. Carroll (2004), §5.4–5.5. Weinberg (1972), §8.4–8.6.
• Bardeen, Press & Teukolsky (1972), ApJ 178, 347: órbitas circulares en Kerr y eficiencia de acreción.
• Cardoso et al. (2009), Phys. Rev. D 79, 064016: esfera de fotones, exponentes de Lyapunov y modos cuasinormales.
• Event Horizon Telescope Collaboration (2019), ApJL 875, L1; (2022), ApJL 930, L12.
• Dyson, Eddington & Davidson (1920), Phil. Trans. R. Soc. A 220, 291: el eclipse de 1919 (lección 3).