El problema restringido de tres cuerpos y los puntos de Lagrange

1. Planteamiento

Dos cuerpos masivos, los primarios de masas $m_1 \ge m_2$, describen órbitas circulares en torno a su centro de masas. Un tercer cuerpo de masa despreciable (la partícula de prueba) se mueve en su campo gravitatorio sin perturbarlos. Restringimos el movimiento al plano orbital de los primarios.

Elegimos unidades adimensionales (Szebehely 1967, cap. 1):

• unidad de masa: $m_1 + m_2 = 1$, y definimos el cociente de masas $\mu = m_2/(m_1+m_2) \in (0, \tfrac12]$;
• unidad de longitud: la separación constante de los primarios;
• unidad de tiempo: tal que $G = 1$. Por la tercera ley de Kepler, $n^2 \cdot 1^3 = G(m_1+m_2) = 1$, así que el movimiento medio de los primarios es $n = 1$ y su periodo es $2\pi$.

En el marco sinódico (que gira con velocidad angular $n=1$ alrededor del centro de masas) los primarios están en reposo en $(-\mu, 0)$ y $(1-\mu, 0)$. Las distancias de la partícula $(x,y)$ a ellos son

Es el problema que Euler (1767) y Lagrange (1772) estudiaron buscando soluciones particulares del problema de tres cuerpos; Lagrange encontró las cinco configuraciones de equilibrio que llevan su nombre, y un siglo después se descubrieron los primeros troyanos de Júpiter exactamente donde él había predicho.

2. Ecuaciones de movimiento en el marco rotante

Sea $\mathbf r$ la posición de la partícula en el marco rotante y $\boldsymbol\omega = \hat{\mathbf z}$ su velocidad angular constante. La aceleración inercial se descompone como

donde $2\boldsymbol\omega\times\dot{\mathbf r}$ es el término de Coriolis y $\boldsymbol\omega\times(\boldsymbol\omega\times\mathbf r) = -(x, y, 0)$ el centrífugo. Igualando a la fuerza gravitatoria por unidad de masa, $-\nabla\big(-\tfrac{1-\mu}{r_1} - \tfrac{\mu}{r_2}\big)$, y agrupando el término centrífugo en un potencial, se obtiene

con el potencial efectivo (o pseudo-potencial)

Explícitamente,

Convención: Murray & Dermott (1999, §3.3) usan la misma $\Omega$ (allí $U$); algunos autores añaden la constante $\tfrac12\mu(1-\mu)$, lo que desplaza todos los valores de $C_J$ en esa cantidad.

3. La integral de Jacobi

Proposición. La cantidad $C_J = 2\Omega(x,y) - (\dot x^2 + \dot y^2)$ se conserva a lo largo de cualquier solución.

Demostración. Multiplicamos la primera ecuación por $\dot x$, la segunda por $\dot y$ y sumamos. Los términos de Coriolis se cancelan, $-2\dot x\dot y + 2\dot y\dot x = 0$ (la fuerza de Coriolis es perpendicular a la velocidad y no realiza trabajo):

porque $\Omega$ no depende explícitamente de $t$. Integrando, $\tfrac12(\dot x^2+\dot y^2) = \Omega - \tfrac12 C_J$. $\blacksquare$

$C_J$ es la única integral primera conocida del problema; Poincaré demostró que no existen otras integrales analíticas uniformes, lo que es la raíz de la no integrabilidad del problema de tres cuerpos. Nótese que en el marco inercial ni la energía ni el momento angular de la partícula se conservan (el campo depende del tiempo); $C_J$ es la combinación $-2(E - \omega L_z)$ de ambas, que sí lo hace.

4. Curvas de velocidad cero (Hill)

Como $v^2 = 2\Omega - C_J \ge 0$, una partícula con constante $C_J$ solo puede ocupar la región de Hill $\{(x,y) : 2\Omega(x,y) \ge C_J\}$. Su frontera, $2\Omega(x,y) = C_J$, es la curva de velocidad cero.

• $C_J$ grande: la región permitida son tres piezas disjuntas — dos «burbujas» alrededor de cada primario y el exterior lejano. Una partícula atrapada cerca de $m_2$ no puede escapar: es la base del criterio de estabilidad de Hill.
• Al disminuir $C_J$, las burbujas se tocan primero en $L_1$ (cuando $C_J = C_{L_1}$), luego la región interior se abre al exterior por $L_2$, después por $L_3$ y finalmente la región prohibida desaparece en $L_4$ y $L_5$.

Para $\mu \le \tfrac12$ se cumple $C_{L_1} > C_{L_2} \ge C_{L_3} > C_{L_4} = C_{L_5} = 3 - \mu(1-\mu)$. Con $\mu = 0.01215$ (Tierra–Luna) los valores calculados con src/lib/math/crtbp.ts son

En el simulador, la zona roja es la región prohibida para la $C_J$ elegida. Prueba a lanzar una partícula en reposo cerca de la Luna con $C_{L_2} < C_J < C_{L_1}$: puede pasar por el «cuello» de $L_1$ hacia la Tierra, pero no escapar al exterior.

5. Los puntos de equilibrio

Los equilibrios en el marco rotante satisfacen $\dot x = \dot y = 0$ y $\nabla\Omega = 0$.

5.1 Puntos triangulares $L_4$, $L_5$

Si $y \ne 0$, $\Omega_y = 0$ exige $1 = \frac{1-\mu}{r_1^3} + \frac{\mu}{r_2^3}$. Escribiendo $x$ e $y^2$ en función de $r_1$ y $r_2$ (se tiene $r_1^2 - r_2^2 = 2x + 2\mu - 1$, e $x^2+y^2 = (1-\mu)r_1^2 + \mu r_2^2 - \mu(1-\mu)$), es inmediato comprobar que

En las variables $(r_1, r_2)$, $\partial\Omega/\partial r_1 = (1-\mu)(r_1 - r_1^{-2}) = 0$ y $\partial\Omega/\partial r_2 = \mu(r_2 - r_2^{-2}) = 0$, de donde $r_1 = r_2 = 1$: la partícula forma un triángulo equilátero con los primarios,

(El cambio $(x,y)\mapsto(r_1,r_2)$ es un difeomorfismo local fuera del eje $y=0$, así que los puntos críticos coinciden.)

5.2 Puntos colineales $L_1$, $L_2$, $L_3$

Sobre $y = 0$, $\Omega_y \equiv 0$ y los equilibrios son las raíces de

Proposición. $f$ tiene exactamente una raíz en cada intervalo $(-\infty,-\mu)$, $(-\mu, 1-\mu)$, $(1-\mu,\infty)$.

Demostración. En cada intervalo, $f'(x) = 1 + \dfrac{2(1-\mu)}{|x+\mu|^3} + \dfrac{2\mu}{|x-1+\mu|^3} > 0$, luego $f$ es estrictamente creciente. En los extremos: $f \to +\infty$ al acercarse a un primario por la izquierda y $f \to -\infty$ por la derecha, mientras que $f(x) \sim x \to \pm\infty$ cuando $x \to \pm\infty$. Por el teorema de Bolzano hay una raíz en cada intervalo y, por monotonía, solo una. $\blacksquare$

Esta demostración justifica el algoritmo de la librería: bisección en cada intervalo, que converge con garantía (el error se divide por dos en cada iteración; ~50 iteraciones dan precisión de máquina). Llamando $\rho = 1 - \mu - x_{L_1}$ la distancia de $L_1$ a $m_2$, el desarrollo para $\mu \ll 1$ da

donde $r_H$ es el radio de Hill; $L_2$ está a $\approx r_H(1 + \tfrac13 r_H)$ al otro lado de $m_2$.

6. Estabilidad lineal

Sea $(x_0,y_0)$ un equilibrio y $x = x_0 + \xi$, $y = y_0 + \eta$. A primer orden,

con las derivadas segundas evaluadas en el equilibrio. Buscando soluciones $\propto e^{\lambda t}$, el determinante

da la ecuación característica

Es bicuadrática: con $\Lambda = \lambda^2$, $\Lambda^2 + b\Lambda + c = 0$. La estabilidad lineal (todas las $\lambda$ imaginarias puras y distintas) exige que ambas raíces $\Lambda$ sean reales, negativas y distintas: $b > 0$, $c > 0$, $b^2 - 4c > 0$.

Puntos colineales. En $y=0$, $\Omega_{xy} = 0$ y, con $A = \frac{1-\mu}{r_1^3} + \frac{\mu}{r_2^3}$, $\Omega_{xx} = 1 + 2A > 0$ y $\Omega_{yy} = 1 - A$. La condición de equilibrio $f(x)=0$ implica $A > 1$ (Problema 2), luego $\Omega_{yy} < 0$ y $c = \Omega_{xx}\Omega_{yy} < 0$: una raíz $\Lambda$ es positiva y hay un $\lambda$ real positivo. $L_1$, $L_2$, $L_3$ son inestables para todo $\mu$ (son puntos silla × centro).

Puntos triangulares. En $L_4$, con $r_1 = r_2 = 1$:

de modo que $b = 1$, $c = \tfrac{27}{4}\mu(1-\mu) > 0$ y

Criterio de Routh (1875). $L_4$ y $L_5$ son linealmente estables si y solo si $27\mu(1-\mu) < 1$, es decir,

Para $\mu < \mu_R$ hay dos frecuencias, $\omega_\pm = \sqrt{-\Lambda_\mp}$. Para $\mu \ll 1$: la frecuencia de libración lenta $\omega_s \approx \sqrt{27\mu/4}$ y la rápida $\omega_f \approx 1$ (epiciclo). Que $L_4$ sea un máximo de $\Omega$ y aun así estable es un efecto genuino de la fuerza de Coriolis (estabilidad giroscópica). La estabilidad lineal no implica estabilidad no lineal; el análisis completo (Deprit & Deprit-Bartholomé 1967, vía KAM) excluye dos valores resonantes de $\mu$ (aquellos en que $\omega_f/\omega_s = 2$ y $3$) y prueba la estabilidad para el resto (Murray & Dermott 1999, §3.8).

7. Aplicaciones

Troyanos de Júpiter. Para Sol–Júpiter, $\mu = GM_J/(GM_\odot + GM_J) = 9.54\times10^{-4} \ll \mu_R$ (GM de Horizons): $L_4$ y $L_5$ son estables y albergan miles de asteroides troyanos que libran en órbitas de «renacuajo» alrededor de $60°$ por delante y por detrás de Júpiter. El periodo de libración lineal es $T_s = 2\pi/\omega_s \approx 2\pi/\sqrt{27\mu/4}$ en unidades adimensionales; dividido entre el periodo de Júpiter ($2\pi$) da $1/\sqrt{27\mu/4} = 12.5$ periodos jovianos, unos 148 años. Pulsa «Órbita troyana (L4)» en el simulador con el preajuste Sol–Júpiter.

Tierra–Luna. Con $\mu = 0.01215$, la Luna está a $384\,400$ km de la Tierra, y $L_1$, $L_2$ quedan a $58\,000$ km y $64\,500$ km de la Luna (valores obtenidos con las abscisas de la tabla anterior: $384\,400\times(1 - \mu - 0.83692)$ y $384\,400\times(1.15568 - 1 + \mu)$). Como $\mu > \mu_R$ no se cumple ($0.01215 < 0.0385$), $L_4$ y $L_5$ Tierra–Luna son linealmente estables, pero la perturbación solar impide que acumulen material de forma duradera (las «nubes de Kordylewski» siguen siendo objeto de debate).

JWST en Sol–Tierra L2. Para el sistema Sol–(Tierra+Luna), $\mu \approx 3.04\times10^{-6}$ y $L_2$ está a $\approx 0.01004\,(1 + 0.00335)$ UA $= 0.01008$ UA $\approx 1.5\times10^{6}$ km de la Tierra, en dirección antisolar ($r_H = 0.01004$ UA), de acuerdo con los «1.5 millones de kilómetros» que da la NASA. Como $L_2$ es inestable (un $\lambda$ real positivo con tiempo de $e$-plegado de unas semanas), el JWST no está en $L_2$ sino en una órbita de halo a su alrededor que, según la NASA, «tarda unos 6 meses en completarse» y es «de tamaño similar a la órbita de la Luna alrededor de la Tierra»; la deriva se corrige con pequeñas maniobras periódicas de mantenimiento de estación. Así el Sol, la Tierra y la Luna quedan siempre del mismo lado del parasol y el telescopio nunca entra en la sombra de la Tierra.

8. Nota numérica

El simulador integra con RK4 (rk4Step) con un paso adaptado a la distancia al primario más cercano, $h = \min(h_{\max},\ \eta\, r_{\min}^{3/2})$, con $h_{\max} = 5\times10^{-3}$ y $\eta = 0.02$. RK4 no es simpléctico: el error en $C_J$ crece lentamente (típicamente de forma lineal en $t$, $\propto h^4$). La lectura «ΔC/C» del simulador es la medida natural del error: con estos parámetros queda por debajo de $10^{-9}$ en cientos de unidades de tiempo para órbitas que no se acercan a los primarios. Dos observaciones:

• El paso adaptado a $r_{\min}^{3/2}$ sigue la escala de tiempo de caída libre local, $\sqrt{r^3/\mu}$: cerca de un primario el campo varía rápido y hace falta resolverlo. Sin esa adaptación, un paso fijo que funciona lejos de los primarios produce errores de $C_J$ de orden unidad en un encuentro cercano.
• La constante de Jacobi es un diagnóstico necesario pero no suficiente: una integración puede conservar $C_J$ a $10^{-9}$ y aun así tener errores de fase grandes en órbitas caóticas (lección 5). Para medir el orden del método hay que comparar con una solución de referencia (paso mucho menor) y verificar que el error escala como $h^4$ (Problema 6).

Problemas

1. (Coriolis no trabaja) Demuestra que si la velocidad angular del marco $n$ fuera una función de $t$, $C_J$ dejaría de conservarse, e identifica el término responsable.

2. (Inestabilidad de los colineales) Demuestra que en $L_1$, $L_2$, $L_3$ se cumple $A = \frac{1-\mu}{r_1^3} + \frac{\mu}{r_2^3} > 1$.

3. (Routh) Comprueba que $\Omega_{xx}+\Omega_{yy} = 3$ en $L_4$ para todo $\mu$ y deduce $b = 1$. Demuestra que $c$ es máximo en $\mu = \tfrac12$ y que $27\mu(1-\mu) = 1$ tiene una única raíz en $(0, \tfrac12)$.

4. (Radio de Hill) A partir de $f(x) = 0$ con $x = 1-\mu-\rho$, $\rho \ll 1$, $\mu \ll 1$, demuestra que $\rho^3 \approx \mu/3$ al orden más bajo, y obtén la corrección $\rho \approx r_H(1 - r_H/3)$.

5. (Simulador) Con el preajuste Tierra–Luna, lanza una partícula en reposo cerca de la Luna con $C_J$ entre $C_{L_1}$ y $C_{L_2}$. Mide cuánto tarda en cruzar a la región terrestre. Explica por qué, aunque el cruce por $L_1$ esté permitido energéticamente, la partícula puede tardar muchas órbitas en encontrarlo.

6. (Computacional: orden del integrador) Integra con RK4 y paso fijo $h$ una órbita que no se acerque a los primarios (por ejemplo, una órbita troyana con el preajuste Sol–Júpiter) durante $t = 50$ y mide $|\Delta C_J|$ para $h = 0.02, 0.01, 0.005, 0.0025$. Representa $\log|\Delta C_J|$ frente a $\log h$ y estima la pendiente. Repite con leapfrog y comenta.

import numpy as np
def f(s, mu):
    x, y, vx, vy = s; r1 = np.hypot(x+mu, y); r2 = np.hypot(x-1+mu, y)
    Ox = x - (1-mu)*(x+mu)/r1**3 - mu*(x-1+mu)/r2**3; Oy = y - (1-mu)*y/r1**3 - mu*y/r2**3
    return np.array([vx, vy, 2*vy + Ox, -2*vx + Oy])
def CJ(s, mu):
    x, y, vx, vy = s; r1 = np.hypot(x+mu, y); r2 = np.hypot(x-1+mu, y)
    return x*x + y*y + 2*(1-mu)/r1 + 2*mu/r2 - vx*vx - vy*vy
mu = 9.54e-4; s0 = np.array([0.5 - mu + 0.02, np.sqrt(3)/2, 0.0, 0.0])  # cerca de L4
for h in [0.02, 0.01, 0.005, 0.0025]:
    s = s0.copy(); C0 = CJ(s, mu)
    for _ in range(int(50/h)):
        k1 = f(s, mu); k2 = f(s + h/2*k1, mu); k3 = f(s + h/2*k2, mu); k4 = f(s + h*k3, mu)
        s = s + h/6*(k1 + 2*k2 + 2*k3 + k4)
    print(h, abs(CJ(s, mu) - C0))

Autoevaluación

Lecturas

• Murray & Dermott (1999), cap. 3: deducción completa de las ecuaciones, la integral de Jacobi, los puntos de Lagrange, la estabilidad (§3.6–3.8) y las órbitas de renacuajo y herradura (§3.9–3.12).
• Szebehely (1967), Theory of Orbits: el tratado de referencia del problema restringido, con los puntos de equilibrio (cap. 4) y su estabilidad (cap. 5).
• Danby (1988), cap. 8: presentación clásica con énfasis en las órbitas periódicas.
• Routh (1875), DOI 10.1112/plms/s1-6.1.86: el criterio original; Deprit & Deprit-Bartholomé (1967), DOI 10.1086/110213: estabilidad no lineal de los puntos triangulares.
• NASA Science, «Webb’s Orbit»: descripción oficial de la órbita del JWST en torno a L₂.