8.12. El potencial modificado

En la sección anterior comprobamos teorica y numéricamente que las superficies de cero velocidad actúan como espejos contra los que se reflejan las trayectorias de partículas de prueba en el CRTBP. La razón de ello es que, por definición, cuando la partícula alcanza estas superficies imaginarias se detiene. Esto es consistente con el hecho además de que más allá de ellas se encuentran regiones que no pueden ser atravesadas por la partícula. Pero surge la pregunta de ¿por qué una vez la partícula llega a la superficie de cero velocidad y se detiene, vuelve después a ganar velocidad para desprenderse de ella? ¿en qué dirección se acelera después de alcanzar la superficie de cero velocidad?.

Para resolver estas preguntas podemos apelar a la ecuación de movimiento de la partícula que escribimos en la Sección El problema circular restringido de los tres cuerpos:

\begin{equation} \label{eq:crtbp_edm} \ddot{\vec{r}} = -(1-\alpha)\frac{\vec{r}_{1}}{r^3_{1}} - \alpha\frac{\vec{r}_{2}}{r^3_{2}} - \hat{e}_z\times(\hat{e}_z\times \vec{r}) -2\hat{e}_z\times \dot{\vec{r}} \end{equation}

Para el asunto particular que nos ocupa, nos interesa la forma que adopta esta ecuación de movimiento unicamente cuando la partícula esta en reposo, es decir se encuentra sobre una superficie de cero velocidad. En estas condiciones el término de la aceleración de Coriolís desaparece y la ecuación se convierte en:

\begin{equation} \label{eq:edm_crtbp_cero_velocidad} \ddot{\vec{r}} = -(1-\alpha)\frac{\vec{r}_{1}}{r^3_{1}} - \alpha\frac{\vec{r}_{2}}{r^3_{2}} - \hat{e}_z\times(\hat{e}_z\times \vec{r}) \end{equation}

donde el subíndice 0 en $\ddot{\vec{r}}_0$ lo hemos usado para indicar que esta ecuación solo aplica para partículas que se encuentran en reposo en el sistema rotante.

Nota: Reposo no implica equilibrio. Aunque para este punto el lector debe seguramente haber reflexionado ampliamente sobre el malentendido común de que el reposo es lo mismo a la ausencia total de movimiento, no sobra, dada la importancia en este punto del libro, recordar esta distinción conceptual fundamental.

Por reposo entendemos el estado de movimiento caracterizado exclusivamente por la condición $\vec v=\vec o$. Ahora bien, el estado completo de movimiento de un cuerpo no solo esta determinado por su velocidad. Otras propiedades cinemáticas también son importantes, siendo la más relevante para esta discusión la aceleración $\vec a$.

Si una partícula esta en reposo, $\vec v=\vec o$ pero también tiene una aceleración nula $\vec a=\vec o$, entonces, en virtud del teorema de inercia la partícula se mantendrá en reposo hasta que una fuerza la saque de esa situación. A esta condición la llamamos equilibrio.

Sin embargo hay estados de movimiento en los que $\vec v=\vec o$ pero $\vec a\neq\vec o$ (reposo pero no equilibrio). Este es el caso por ejemplo de una partícula lanzada hacia arriba en un campo gravitacional casi uniforme (el campo gravitacional muy cerca a la superficie de la Tierra). En la parte alta de su trayectoria la partícula alcanza el reposo (se detiene) pero su aceleración (la aceleración de la gravedad) es distinta de cero (no hay equilibrio); un instante de tiempo después, la partícula estará nuevamente en movimiento hacia abajo.

Esto es justamente lo que pasa cuando en el CRTBP la partícula de prueba alcanza una superficie de cero velocidad: al hacerlo, se detiene, pero si la aceleración en aquel lugar es distinta de cero (dependiendo de lo que dicte justamente la Ec. edm_crtbp_cero_velocidad) entonces un momento después estará otra vez en movimiento.

Es fácil mostrar (ver Problemas al final del capítulo) que la ecuación de movimiento en (edm_crtbp_cero_velocidad) puede escribirse en la forma:

\begin{equation} \label{eq:edm_crtbp_reposo} \ddot{\vec{r_0}} = - \vec\nabla \left[-\frac{(1-\alpha)}{r_{1}}-\frac{\alpha}{r_{2}}-\frac{1}{2}(x^2 + y^2)\right] \end{equation}

es decir, en el sistema de referencia rotante, una partícula en reposo experimentará una aceleración producto de un potencial escalar diferente del simple potencial gravitacional newtoniano. A este potencial lo llamaremos aquí el potencial modificado del CRTBP y tiene la forma:

\begin{equation} \label{eq:Vmod_crtbp} V_\mathrm{mod} = -\frac{(1-\alpha)}{r_{1}} - \frac{\alpha}{r_{2}} - \frac{1}{2}(x^2 + y^2) \end{equation}

En términos del potencial modificado las ecuación de movimiento del sistema (Ec. crtbp_edm) incluso cuando la velocidad de la partícula de prueba no es 0, se pueden escribir como:

\begin{eqnarray} \label{eq:crtbp_edm_Vmod} \ddot{x}-2\dot y & = & -\frac{\partial V_\mathrm{mod}}{\partial x} \\ \nonumber \ddot{y}+2\dot x & = & -\frac{\partial V_\mathrm{mod}}{\partial y} \\ \nonumber \ddot{z} & = & -\frac{\partial V_\mathrm{mod}}{\partial z} \\ \end{eqnarray}

¿Cómo luce este potencial modificado?. Naturalmente el potencial $V_\mathrm{mod}$ depende de todas $(x,y,z)$ y visualizarlo sería muy complejo. Sin embargo podemos, como hicimos en secciones anteriores, concentrarnos en el valor del potencial sobre el plano $x-y$ (el plano de la órbita de las partículas masivas). Visualizarlo allí puede darnos alguna intuición de su estructura 3d.

En el algoritmo a continuación se calcula el valor del potencial modificado en una malla coordenada como la que hemos usado a lo largo de este capítulo y se visualiza en el espacio de tres dimensiones:

In [1]:
#Parámetros del sistema
alfa=0.3

#Malla coordenada
from numpy import meshgrid,zeros_like,linspace
rango=1.5
NG=100
X,Y=meshgrid(linspace(-rango,rango,NG),
             linspace(-rango,rango,NG))
Z=zeros_like(X)
#Factor de suavizado (ver nota)
eps=0.8*alfa

#Distancia relativa
from numpy import sqrt
r1=sqrt((X+alfa)**2+Y**2+Z**2+eps**2)
r2=sqrt((X-1+alfa)**2+Y**2+Z**2+eps**2)

#Calcula el potencial
Vmod=-(1-alfa)/r1-alfa/r2-0.5*(X**2+Y**2)
Vmod = 
[[-2.72 -2.68 -2.64 ... -2.64 -2.68 -2.72]
 [-2.68 -2.64 -2.6  ... -2.6  -2.64 -2.68]
 [-2.65 -2.6  -2.57 ... -2.56 -2.6  -2.64]
 ...
 [-2.65 -2.6  -2.57 ... -2.56 -2.6  -2.64]
 [-2.68 -2.64 -2.6  ... -2.6  -2.64 -2.68]
 [-2.72 -2.68 -2.64 ... -2.64 -2.68 -2.72]]

Nótese que el valor absoluto del potencial es del orden de la mitad de la constante de Jacobi. Esto era de esperarse del hecho de que las expresiones que definen ambas cantidades son casi idénticas (Ecs. Vmod_crtbp y constante_Jacobi).

Nota: el factor de suavizado. En el cálculo del potencial hemos usado un artefacto matemático, que no procede de la teoría del CRTBP. Puesto que el valor del potencial es teóricamente infinito cuando $r_1=0$ o $r_2=0$ hemos introducido en las expresiones para el potencial una cantidad ficticia eps proporcional al parámetro $\alpha$ del sistema (relacionado con la masa de las partículas), y que hace que el potencial mantenga un valor finito cerca al centro de las partículas. La introducción de esta cantidad artificial implica también que los valores calculados del potencial no son los correctos a una cierta distancia de las partículas (típicamente similar al valor de eps).

Este artefacto numérico facilita la representación gráfica del potencial, pero también es importante en otras aplicaciones de la mecánica celeste. Así por ejemplo cuando se simula numéricamente la interacción gravitacional de partículas en sistemas de N cuerpos tales como galaxias o nubes de gas, el facto eps recibe el nombre de factorde suavizado y se usa para simular de alguna manera el hecho de que las *partículas de materia** tienen en realidad una dimensión finita (no son realmente puntuales).

Una gráfica del potencial modificado sobre el plano $x-y$ en el CRTBP se puede crear con el algoritmo a continuación:

In [3]:
%matplotlib inline

In [4]:
#Grafico
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

fig=plt.figure()
ax=fig.add_subplot(111, projection='3d')
ax.plot_surface(X,Y,Vmod,cmap="jet",rstride=2,cstride=2)
ax.set_zlim((-3,0))
ax.view_init(30,-15)

ax.set_xlabel("$x$")
ax.set_ylabel("$y$")
ax.set_zlabel("$V_\mathrm{mod}$")
ax.set_title(f"Potencial modificado en el CRTBP con $\\alpha={alfa}$")
fig.tight_layout()

Figura 8.122.

Otra manera de visualizar el mismo potencial es usando un gráfico de contornos:

In [5]:
fig=plt.figure()
ax=fig.gca()

#Contornos
ax.contour(X,Y,Vmod,levels=15,colors='k',linestyles='solid',alpha=0.3)
c=ax.contourf(X,Y,Vmod,levels=100,cmap="jet")
fig.colorbar(c)

#Decoración
ax.set_xlabel("$x$")
ax.set_ylabel("$y$")
ax.set_title(f"Potencial modificado en el CRTBP con $\\alpha={alfa}$")
fig.tight_layout()

#Puntos de interés
ax.plot([0],[0],'ko')
ax.text(0,0,'A',fontsize=18)
ax.plot([1],[1],'ko')
ax.text(1,1,'B',fontsize=18)
ax.plot([0.8*alfa],[-1+0.7*alfa],'ko')
ax.text(0.8*alfa,-1+0.7*alfa,'C',fontsize=18)
ax.plot([1+0.5*alfa],[0],'ko')
ax.text(1+0.5*alfa,0,'D',fontsize=18);

Figura 8.123. Contornos del potencial modificado en el CRTBP (equipotenciales) y algunos puntos de interés para analizar.

¿Cómo podemos interpretar la estructura del potencial modificado?. Como es corriente en mecánica, el potencial asociado a una fuerza, dicta la manera como la partícula se acelera. En particular, la dirección de la aceleración es perpendicular a las curvas de equipotencial (dirección del gradiente) y va en el sentido en el que el potencial disminuye.

En la gráfica de equipotenciales de la Figura (code:Vmod_contornos) hemos resaltado algunos puntos que merecen atención y que nos permiten interpretar mejor el significado de $V_\mathrm{mod}$:

  • Punto A. Un punto cercano a la partícula más masiva. Al examinar la posición de este punto entre los equipotenciales nos damos cuenta que la aceleración instantánea de una partícula de prueba puesta en reposo allí, apuntara en el sentido de la partícula (hacia donde disminuye el potencial). En otros términos, las partículas que en su movimiento en el CRTBP alcanzan superficies de cero velocidad cercanas a las partículas, tenderán a reflejarse en dirección a ellas. Eso es precisamente lo que vimos al estudiar numéricamente la dinámica en los ejemplos de la Sección Las regiones de exclusión.
  • Punto B. Un punto lejano a ambas partículas. Allí la situación también es bastante siemple. Dado que en puntos lejanos a las partículas el potencial disminute hacia afuera (como vemos por las faldas del potencial modificado en el gráfico en 3d de la Figura (code:crtbp_Vmod_3d)), una partícula dejada en reposo en el punto B se acelerará inicialmente hacia afuera, como si fuera repelida por el sistema. La razón aparente de este contrasentido es que debemos recordar que en el sistema rotante, incluso si la partícula está en reposo, actúa la aceleración centrífuga que se hace mayor en tanto más lejos esté la partícula del origen de coordenadas.
  • Punto C. Un punto no colineal con las partículas masivas pero aproximadamente equidistante de ellas. Como vemos en Figura (code:crtbp_Vmod_3d)) en este punto se produce un máximo del potencial modificado. Esto implica que allí cerca el gradiente del potencial es cero y una partícula en reposo tenderá a quedarse en reposo. Llamamos a este un punto de equilibrio en el sistema (ver próxima sección).
  • Punto D. Un punto situado cerca al punto colineal $L_2$, que definimos en la Sección Las regiones de exclusión. Este punto, junto con aquellos cercanos también a $L_1$ y $L_3$, tiene una característica peculiar. Como vemos en Figura (code:crtbp_Vmod_3d)) en este punto se produce un máximo local del potencial modificado. Esto implica, como sucedio en el caso de C, que allí cerca el gradiente del potencial también es cero y por lo tanto se trata también de un punto de equilibrio. Pero hay una diferencia fundamental entre el punto D y el C. No importa la dirección en la que se desplace la partícula después de salir de su posición en C, la aceleración se dirigirá en cualquier dirección menos hacia C. C actúa como una especie de punto repulsivo. En su lugar en el punto D ocurre algo curioso. Si la partícula se desplaza ligeramente en dirección perpendicular a la línea que une las partículas masivas ($+y$ o $-y$) la aceleración tendrá a restituirla allí de donde salió. En su lugar si se desplaza acercándose o alejándose de las partículas masivas el potencial tenderá a acercar o alejar rápidamente a la partícula de los cuerpos masivos. Decimos que D está en un punto de silla del potencial.