9.8.8. El problema de los dos cuerpos con n arbitrario¶
Como vimos experimentando con la solución numérica al problema de los dos cuerpos sometidos a fuerzas centrales, esto es con una función de energía potencial $U(r)$ que solo depende de la distancia relativa, si $n$ no es 2 las trayectorias resultantes del vector relativo tienen una estructura peculiar: no se cierran sobre sí mismas y parecen elipses que precesan alrededor del origen de coordenadas. Aunque solo en el caso de $n=1$ parece existir una solución matemáticamente simple y regular, las variacións observadas de la variable radial $r$ y la variable angular $\theta$, incluso para valores arbitrarios reales de $0<n<2$, siguen siendo relativamente regulares.
¿Qué podemos decir, desde el punto de vista analítico, sobre las trayectorias de partículas sometidas a fuerzas centrales con $n\neq 1$?. En partícular podemos preguntarnos por aquellos sistemas que tienen un $n$ valores de $n$ muy cercanos a 1.
Consideremos por ejemplo el movimiento de una partícula sometida a una fuerza específica del tipo:
\begin{equation} \label{eq:f_modificada} f(r)=-\frac{\mu}{r^2}-\frac{\sigma}{r^4} \end{equation}donde $\sigma\ll\mu r^2$ para todos los valores de $r$ relevantes para el sistema. Si bien esta elección de la fuerza parece completamente arbitraria, y bien podríamos comenzar analizando el caso de una fuerza del tipo $f(r)=-\mu/r^{2+\delta}$ con $\delta$ pequeño, hay dos razones para que la sencilla exploración que realizaremos en los siguientes párrafos utilice esta forma específica. La primera es que esta forma de la fuerza es una perturbación de la forma newtoniana que aparece explícitamente en el primer término. La segunda razón la descubriremos a lo largo del desarrollo: esta fuerza permite una descripción analítica adecuada para los propósitos del libro. Y la tercera y quizas la más importante, es el hecho de esta es la dependencia funcional que tiene la fuerza específica efectiva en la teoría general de la relatividad. Es decir, no estamos tratando aquí con una fórmula arbitraria, sino con una que puede llegar a tener un verdadero interés físico.
Si reemplazamos $f(r)$ en la ecuación de la forma orbital (Ec. forma_orbital), la ecuación radial para esta forma particular de la fuerza adopta la forma:
\begin{equation} \label{eq:forma_orbital_precesion} \frac{d^2 u}{d\theta^2}+u = \frac{\mu}{h^2}+\frac{\sigma}{h^2} u^2 \end{equation}Podemos comprobar que se trata de una ecuación diferencial no homogenea muy similar a la ecuación análoga para el caso newtoniano (Ec. forma_orbital_kepler). La única diferencia con aquella es el término adicional al lado derecho $\sigma u^2/h^2$. Como es común en la teoría de ecuaciones diferenciales, podemos proponer una solución a esta ecuación que sea la suma de la solución encontrada para la Ec. (forma_orbital_kepler) y una función adicional, desconocida que podemos llamar $g(\theta)$:
$$ u(\theta)=\frac{1}{p}(1+e\cos\theta)+g(\theta) $$donde $p=h^2/\mu$. Asumiremos que $g(\theta)$ contiene términos que dependen de $\sigma$ y que hacen que en general $g(\theta)\ll u(\theta)$.
Si reemplazamos la solución propuesta en la Ec. (forma_orbital_precesion) obtenemos:
$$ \frac{d^2 g}{dt^2}+g=q\left[\frac{1}{p}(1+e\cos\theta)+g(\theta)\right]^2 $$donde hemos introducido un nuevo parámetro $q=\sigma/h^2$
Si ahora desarrollamos el lado derecho de la ecuación y eliminamos todos los términos cuadráticos o de orden superior en $q$ o $g(\theta)$, la ecuación diferencial para $g$ nos queda:
$$ \frac{d^2 g}{dt^2}+g\approx\frac{q}{2p^2}(2+e^2+4e\cos\theta+e^2\cos 2\theta) $$Podemos absorber el factor constante $(2p^2/q)$ en la función $g$ y escribir la función resultante $g'(\theta)= 2p^2 g(\theta)/q$ como la superposición de tres funciones
$$ g'=g'_0+g'_1(\theta)+g'_2(\theta) $$donde la primera es una constante que no depende de $\theta$. La segunda depende de funciones trigonométricas con argumento $\theta$ ($\cos\theta$ o $\sin\theta$) y la tercera es una función que depende de funciones trigonométricas con argumento $2\theta$.
No es difícil mostrar (ver Problemas al final del capítulo) que las funciones resultantes son:
\begin{eqnarray} \nonumber g'_0 & = & 2+e^2\\ \nonumber g'_1 & = & 2e\theta\sin\theta\\ \nonumber g'_2 & = & -\frac{1}{3} e^2\cos 2\theta\\ \end{eqnarray}De este modo una solución general a la ecuación diferencial en Ec. (forma_orbital_precesion) se puede escribir en la forma:
\begin{equation} \label{eq:conica_precesion} u(\theta)\approx\frac{1}{p}\left[1+e\cos\theta+\frac{q}{p}\left(1+\frac{e^2}{2}\right)+\frac{q}{p}e\theta\sin\theta-\frac{q}{p}\frac{e}{6}\cos 2\theta\right] \end{equation}Como vemos $u(\theta)$ en el caso de la fuerza postulada, es casi una cónica pero esta perturbada por 3 términos pequeños (proporsionales a $q$) que producen contribuciones que podríamos clasificar como:
- Un término constante $(q/p)(1+e^2/2)$. Este término aumenta ligeramente el valor de $u(\theta)$ (disminuye el valor de $r$), pero lo hace de forma independiente de la posición angular.
- Un término cuyo valor aumenta monótonamente, $(qe/p)\theta\sin\theta$. Llamamos a este el término secular.
- Y un término pequeño, que oscila alrededor de cero con amplitud proporcional $-qe/p$ y con un período cercano a la mitad del período kepleriano (el período que tendría la órbita si la fuerza fuera la gravedad newtoniana). Llamamos a este el término oscilatorio.
El siguiente algoritmo permite comparar, visualmente, la contribución que cada uno de estos términos hace en la solución final obtenida en la Ec. (conica_precesion).
%matplotlib inline
#Parametros del sistema
q=0.05
p=1.0
e=0.1
#Valores de teta
from numpy import linspace,sin,cos,pi
teta_min=0.0
teta_max=5*pi
tetas=linspace(teta_min,teta_max,1000)
#Términos de u
u_conica=(1/p)*(1+e*cos(tetas))
u_constante=(q/p)*(1+e**2/2)
u_secular=(q/p)*e*tetas*sin(tetas)
u_oscilatorio=-(q/p)*(e/6)*cos(2*tetas)
#Gráfico
import matplotlib.pyplot as plt
fig=plt.figure()
ax=fig.gca()
ax.plot(tetas,u_conica+u_constante,'k--',label="cónica+constante")
ax.plot(tetas,u_conica+u_constante+10*u_oscilatorio,'r:',
label="cónica+constante+oscilatorio")
ax.plot(tetas,u_conica+u_constante+u_oscilatorio+u_secular,'b:',
label="cónica+constante+oscilatorio+secular")
#Decoración
ax.legend(loc='lower right')
ax.set_xlabel("$\\theta$ [rad]")
ax.set_ylabel("$u(\\theta)$ [rad]")
ax.set_xlim((teta_min,teta_max))
ax.set_ylim((0.8,1.2))
ax.grid();
En la solución representada en la Figura (code:conica_precesion) vemos que el término oscilatorio (que esta multiplicado allí artificialmente por un factor de 10) no tiene casi ningún efecto en modificar la cónica original. Como mencionamos antes el término constante simplemente modifica la distancia al periapsis y al apoapsis (extremos de la curva punteada), sin moficiar sustancialmente la forma de la misma. Solo el término secular introduce modificaciones significativas en la trayectoria y es responsable del efecto que observabamos en secciones anteriores: la precesión del periapsis.
Teniendo en cuenta estos resultados, podemos simplificar aún más la Ec. (conica_precesion) escribiéndola como:
$$ u'(\theta)\approx\frac{1}{p'}\left[1+e\cos\theta+\frac{qe}{p}\theta\sin\theta\right] $$donde $p'$ es un semilatus rectum efectivo (debido al cambio en la distancia al periapsis y al apoapsis de la órbita por efecto del término constante).
El segundo término del lado derecho lo podemos escribir también en la forma:
$$ \cos\theta+\frac{q}{p}\theta\sin\theta\approx\cos\left(\theta-\frac{q\theta}{p}\right) $$donde nos hemos valido del hecho que como $q/p\ll 1$, entonces $\sin(q\theta/p)\approx q\theta/p$ y $\cos(q\theta/p)\approx 1$
De allí, finalmente, la solución para la ecuación orbital del sistema relativo de dos cuerpos sometidos a una fuerza del tipo dado por la Ec. (f_modificada) se puede aproximar como:
$$ u'(\theta)\approx\frac{1}{p'}[1+e\cos(\theta-q\theta/p)] $$al menos para $\theta$ no muy grande (las primeras vueltas de la partícula).
La característica más notable de esta trayectoria es que los puntos de máximo (o de mínimo) no se repiten, como en el caso de la cónica kepleriana, exactamente cada $2\pi$ radianes. Después de $\theta=0$ el siguiente máximo de $u'(\theta)$, es decir, el siguiente periapsis, se produce cuando:
$$ \theta_p-\frac{q\theta_p}{p}=2\pi $$Despejando $\theta_p$ obtenemos:
$$ \theta_p=\frac{2\pi}{1-q/p}\approx2\pi(1+q/p) $$Es decir, por cada vuelta se produce un desplazamiento en la posición del perihelio correspondiente a:
\begin{equation} \label{eq:avance_perihelio} \Delta\theta_p\approx\frac{2\pi q}{p} \end{equation}Un poco de historia: el avance del perihelio de Mercurio. En el año 1915, el físico Suizo-Americano Albert Einstein ("ainshtain") estaba en una carrera contra el tiempo buscando las que se convertirían a la larga en las "verdaderas" leyes de la gravitación universal. Desde el año 1907 Einstein ya sospechaba que el fenómeno gravitacional era algo más que lo que Newton había intuído que era y que nunca había podido precisar. A la larga, Einstein comprendió que la "fuerza de gravedad" no era en realidad una fuerza aplicada, sino una fuerza resultante, una manifestación del movimiento de los cuerpos en el espacio-tiempo "distorsionado" cerca a cuerpos muy masivos. Una de las primeras pruebas contundentes de la validez de su teoría, la obtuvo el 18 de noviembre de 1915 cuando resolvió uno de los problemas abiertos mejor conocidos de la mecánica celeste: la precesión anómala del perihelio de Mercurio.
$$ f(r)=-\frac{\mu}{r^2}-\frac{3\mu h^2}{c^2 r^4} $$Desde hacía más de 50 años las observaciones precisas de la órbita del planeta Mercurio, habían revelado que el perihelio del planeta no ocurría en cada órbital en el mismo lugar que en la anterior. El perihelio se desplazaba un poco (1.38 segundos de arco o arcsec por cada vuelta, o equivalentemente 575 arcsec por siglo), alargando la llegada de Mercurio a este punto en su órbita. Usando la mejor mecánica celeste de la época, los astrónomos habían determinado que las perturbaciones gravitacionales de los planetas podían explicar la mayor parte del efecto: si la teoría de la gravitación de Newton era correcta el Perihelio de Mercurio se debía desplazar 1.28 arcsec por vuelta o lo que es lo mismo 532 arcsec por siglo. Por mucho que intentaron corregir la cifra (que era inferior a las mejores estimaciones astronómicas hechas en aquella época por 43 segundos de arco por siglo) no lo consiguieron. Intentaron poniendo un planeta perturbador entre el Sol y Mercurio (Vulcano lo llamaron) pero nunca lo encontraron. Intentaron poner un cinturon de asteroides, pero no consiguieron ver nada como eso cerca al Sol.
Lo que la teoría de Einstein decía en 1915 era, basicamente, que la teoría Newtoniana de la gravedad solo era una aproximación de primer orden, a la verdadera ley de gravitación universal. En particular, después de utilizar la nueva teoría de Einstein, puede probarse que la ecuación radial de movimiento para un partícula que se mueve cerca al Sol, tiene una forma idéntica a la ecuación de la forma orbital que dedujimos en esta sección. Sin embargo el término asociado con la fuerza específica en la mecánica celeste clásica, es diferente. Esta diferencia puede explicarse si se asume que la fuerza específica que mueve a las partículas alrededor de nuestra estrella (y en realidad alrededor de cualquier cuerpo esférico), es, en lugar de la newtoniana $f(r)=-\mu/r^2$, una fuerza de la forma:
$$ \sigma=\frac{3\mu h^2}{c^2} $$donde $h$ es el momento angular reducido específico y $c$ es la velocidad de la luz. Esta es justamente la forma funcional que tiene la fuerza central no Newtoniana que estudiamos en este capítulo donde reconocemos el valor de la constante
$$ q=\frac{3\mu}{c^2} $$De allí podemos calcular a su vez el valor del parámetro $q$:
$$ \Delta\theta_p=\frac{6\pi\mu}{c^2 a(1-e^2)} $$Con esto podemos y teniendo en cuenta que $p=a(1-e^2)$, el avance del perihelio por vuelta, predicho por la teoría que desarrollamos en esta sección y que calculamos con la Ec. (avance_perihelio) será:
$$ \Delta\theta_p=0.1030\;\mathrm{arcsec} $$Reemplazando los valores conocidos para la masa del Sol ($\mu$), la velocidad de la luz y el semijememayor y excentricidad promedio de mercurio se obtiene:
$$ \frac{36.525 \Delta\theta_p}{87,97}=42.66\;\mathrm{arcsec} $$por vuelta. Teniendo en cuenta que el período de revolución de Mercurio es 87.97 días y un siglo contiene 36.525 días, el avance por siglo, calculado según la ecuación Ec. (avance_perihelio) y usando la fuerza específica equivalente de la relatividad general, será:
que corresponde a la anomalía observada.
