Pasado y Futuro del Sistema Solar: Introducción a la Dinámica Molecular

• Teoría: Mecánica clásica, ecuaciones de Hamilton, campos centrales, colisiones inelásticas.
• Ténica numérica: Dinámica Molecular: resolución de sistemas de ecuaciones diferenciales por el método de Verlet.
• Programación: Fortran, C, Introducción a X-windows.

Introducción Teórica:

Hace 200 años, Pierre Simon, Marqués de Laplace comenzó a pensar sobre el origen del sistema solar. Su hipótesis fué presentada en 1796 en Exposition du système du monde. De acuerdo con la teoría de Laplace, el sistema solar evolucionó a partir de un gas incasdencente girando alrededor de un eje. Según se iba enfriando, el gas se iba contrayendo lo que originaba un rotación más rápida debido a la conservación del momento angular. En algun momento la parte exterior del gas se separó formando anillos que posteriormente condensaron para formar planetas. La parte interior se condensó en un sólo cuerpo: el sol. La idea detrás de esta hipótesis no fue del todo de Laplace, ya fue propuesta antes en una forma menos elaborada por Thomas Wright e Immanuel Kant. En cualquier caso Laplace elaboró formalmente la idea cuyos detalles técnicos se encuentran en su gran obra Mécanique céleste 1799-1825 (5 volúmenes). Esta obra se piensa que es la culminación del punto de vista newtoniano sobre la gravitación. De hecho, Laplace mostró que el sistema solar era estable bajo la acción de pequeñas perturbaciones. Su trabajo mostraba la no necesidad de intervención divina en la formación y estabilidad del sistema solar. Se dice que Napoleón le comentó este hecho y Laplace le respondió: No necesito introducir dicha hipótesis.

Esta imágen genérica de Laplace es la que actualmente aceptamos como la teoría más viable que explica la formación del sistema solar. Sin embargo, Laplace no consiguió explicar con precisión como a partir de discos de gas se podían formar planetas. Hoy sabemos lo que pudo ocurrir para la formación de los planetas: las partículas del gas comenzaron a colisionar entre si formando partículas cada vez mayores. En cuanto alcanzaron tamaños suficientemente grandes, su atracción gravitatoria favoreció su acumulación. Esos objetos se denominan planetesimales. Los cometas y los asteroides se cree que son planetesimales sobrantes durante la formación del sistema solar. En una última etapa, los planetesimales se fueron uniendo para formar los planetas.

Podemos finalmente resumir la teoría actual sobre la formación del sol y de los planetas:

• 1. Inicialmente existe una nube de Hidrógeno y Helio girando alrededor de un eje.
• 2. El gas comienza a contraerse debido a la fuerza de la gravedad.
• 3. La nube forma un disco. Las partículas tienden a dirigirse hacia el centro formando un proto-sol.
• 4. Debido a la rotación, algo de polvo y gas orbita alrededor del centro del protosol como un disco plano.
• 5. El proto-sol colapsa, la temperatura se incrementa hasta que los átomos de Hidrógeno comienzan a fusionarse.
• 6. Los proto-planetas comienzan a formarse a partir del polvo y gas que se encuentra orbitando el proto-sol.
• 7. El sistema solar se va vaciando de escombros.
• 8. Se separan los elementos sólidos de los gaseosos: los elementos gaseosos fueron expulsados de los planetas interiores por el viento solar.
• 9. Diferenciación: se separan los elementos pesados de los ligeros. Los elementos pesados (Fe, Ni) se depositan en los centros de los planetas, los elementos ligeros (Al, Si, Mg, Na, K, Ca) se sitúan el la corteza y manto.

Hay evidencias observacionales del anterior esquema evolutivo:

El futuro del sistema solar puede ser estudiado a partir de simulaciones numéricas exhaustivas. Como ejemplo ver el artí culo de J. Sackmann, A.I. Boothroyd y K. Kraemer Our Sun. III. Present and Future in The Astrophysical Journal 418 457-468 (1993). Allí se nos realata como el sol se convertirá en gigante roja, y después de varios episodios de inestabilidad y pérdida de masa, acabará su vida como enana blanca.

Problemas:

• Obligatorio: Diseñar y escribir el programa que simula el comportamiento del sistema solar. Supondremos sólo la existencia del Sol y los nueve planetas. Suponer que el movimiento planetario se realiza en un plano. Usar datos reales de excentricidad de las órbitas, masas y distancias.

• Voluntario (1): Formación de los planetas a partir de planetesimales. Partimos de la condición inicial de un sol masivo y un gran número de pequeños cuerpos orbitando alrededor de éste con diferentes energías y excentricidades. Los pequeños planetesimales interaccionan entre si gravitatoriamente. La formación de planetas se va realizando mediante choques totalmente inelásticos entre planetesimales. Dada una condición inicial de planetesimales, estudiar la distribución de tamaños de los planetas finales. Puesto que la colisión entre planetesimales es inelástica, parte de su energía cinética antes de su colisión se convierte en calor. Mientras que en ese mismo proceso las energías caloríficas que tenían se suman. Seguir la pista de esas energías y obtener las energías internas de la distribución de planetas finales. Comparar las propiedades obtenidas con las del sistema solar en la actualidad.

• Voluntario (2): Futuro del sistemas solar. A partir de las condiciones iniciales actuales (posición, velocidad y excentricidad de la órbita) de los planetas, estudiar su comportamiento futuro sabiendo que el Sol irá sufriendo cambios estructurales, asi realizar la simulación siguiendo el siguiente proceso:

  • Secuencia Principal: Dejar evolucionar la condición inicial un tiempo (para que se ajusten las órbitas). La masa del Sol es M = 1 y su radio R = 1 (en unidades apropiadas.
  • Se va consumiendo el Hidrógeno: Cambiar la masa y el radio del Sol a M = 1, R = 1.575. Dejar evolucionar para que se ajusten las órbitas.
  • Se va consumiendo el Hidrógeno: Cambiar la masa y el radio del Sol a M = 0.9998, R = 2.3. Dejar evolucionar para que se ajusten las órbitas.
  • Se va consumiendo el Hidrógeno: Cambiar la masa y el radio del Sol a M = 0.9935, R = 9.5. Dejar evolucionar para que se ajusten las órbitas.
  • Ignición del Helio: Cambiar la masa y el radio del Sol a M = 0.7249, R = 166. Dejar evolucionar para que se ajusten las órbitas.
  • Se acaba el Helio del núcleo: Cambiar la masa y el radio del Sol a M = 0.708, R = 20. Dejar evolucionar para que se ajusten las órbitas.
  • Pulsos inestables  de la corteza de Helio: Cambiar la masa y el radio del Sol a M = 0.591, R = 100. Dejar evolucionar para que se ajusten las órbitas.
  • Enana Blanca: Cambiar la masa y el radio del Sol a M = 0.54137, R = 0.058. Dejar evolucionar para que se ajusten las órbitas.

• Voluntario (3): Sol binario. Suponer que dos soles, equidistantes y de masa igual a la masa solar, giran respecto a su centro de masa con una velocidad angular w. Colocar 9 planetas con condiciones iniciales idénticas a los correspondientes al sistema solar y estudiar su evolución. Estudiar la estabilidad del sistema. Tener en cuenta que un planeta desaparece cuando colisiona con uno de los soles. Estudiar el efecto de w. 

Modelo y método numérico:

Nuestro sistema se compone de N esferas de diferentes masas m_i i = 1, N. Los cuerpos interaccionan entre si mediante fuerzas gravitatorias newtonianas. La ecuación de evolución del cuerpo i es:

A diferencia del problema de la nave espacial aqui hemos de resolver (genéricamente) un gran número de ecuaciones diferenciales (sobre todo para el problema Voluntario (2)). En este caso, el uso del algoritmo de Runge-Kutta es excesivamente costoso desde el punto de vista computacional. En estos casos, uno puede renunciar a la precisión del cálculo de cada una de las trayectorias de las partículas por la rápidez del cálculo. Este hecho, per se, justifica el uso de dichos algoritmos simplificados. Pero además, desde un punto de vista más formal este cambio de estrategia está permitido siempre que las propiedades estadísticas del sistema lo permita. Esto es, sea un punto inicial del espacio de las fases x. Supongamos que su evolución mediante sus ecuaciones del movimiento lo llevan a x_t cuando transcurre un tiempo t. Si en lugar de las ecuaciones del movimiento, utilizamos una aproximación discreta de las mismas, es natural pensar que partiendo de x no llegaremos exáctamente a x_t sino a un x_t¢ que será cercano a x_t si los incrementos del algoritmo son lo suficientemente pequeños. Es obvio que al usar el algoritmo discreto estamos cometiendo un error en la evaluación de las trayectorias, pero estaríamos salvados si x_t¢ perteneciese a una trayectoria exácta correspondiente a una condición inicial x¢. Esta propiedad ocurre si el problema es suficientemente insensible a las condiciones iniciales, suficientemente caótico y tiene muchos grados de libertad. Para estos casos podemos concluir que dado un punto x y su evolución algorítmica x_t¢, existe una trayectoria típica exácta que es la sombra de la trayectoria algorítmica en el espacio de las fases (sombra significa matemáticamente que para todo tiempo la distancia entre la trayectoria algorítmica y la sombra está acotada y es pequeña). Obviamente la demostración rigurosa de estas propiedades y las condiciones en que ocurren es extraordinariamente compleja. En cualquier caso, sea debido a la complejidad del problema, o porque la existencia de trayectorias sombra nos lo permita, el uso de estos algoritmos simplificados sólo tendrá sentido cuando nos interesen propiedades globales o macroscópicas del sistema puesto que en ese caso lo que importa no es el detalle del camino seguido (como ocurre en el problema del cohete donde saber si orbita la Luna o se estrella en ella es importante) sino sus propiedades promedio. En cualquier caso, el algoritmo utilizado debe de tener ciertas propiedades de estabilidad. Por ejemplo, si la dinámica conserva la energía y el momento angular, es importante y necesario que el algoritmo sea diseñado para que los conserve.

Desde el punto de vista computacional los algoritmos más eficaces son los llamados algoritmos de un paso, esto es, dadas las posiciones y momentos de la partículas en un instante dado, podemos a partir de ellas calcular sus respectivas posiciones y momentos en un paso posterior (observar que el algoritmo de Runge-Kutta utiliza dos pasos). Veamos como construir uno de estos algoritmos, el llamado algoritmo de Verlet. Este algoritmo se basa en la discretización más obvia de las primeras y segundas derivadas, esto es:

En resumen, el algoritmo de Verlet es:

• (0) Dar una h inicial, t = 0, y las posiciones y velocidades de las partículas inicialmente.
• (1) Evaluar F(t).
• (2) Evaluar r(t+h) y w = v(t)+hF(t)/2.
• (3) Evaluar F(t+h).
• (4) Evaluar v(t+h) = w+hF(t+h)/2.
• (5) t = t+h. Ir a (2).

Para aplicar el algoritmo de Verlet en nuestro caso, hemos de tener en cuenta que vamos a trabajar con números muy grandes y muy pequeños simultáneamente (G = 6.67×10^-11 y la Masa solar es M_s = 1.99×10³0). Esto nos hace pensar que va a ser conveniente reescalar las ecuaciones del movimiento. En particular, se puede elegir la siguiente reescala:

donde elegimos c = 1.496×10¹¹ que no es más que la distancia entre la tierra y el Sol. En estas unidades, Neptuno se encuentra aproximadamente a 30 unidades c de distancia al Sol. Observamos que una unidad t¢ corresponden a unos 58.1 dias. En estas unidades, las ecuaciones del movimiento dadas por (1) se convierten en:

Para el caso del problema de la formación de los planetas a partir de los planetesimales suponemos que la agregación de materia se produce mediante colisiones totalmente inelásticas. Asi, el algoritmo de colisión entre las partículas i y j lo elegimos:

• (0) Dados r_i(t), v_i(t), r_j(t) y v_j(t).
• (1) t = t+h
• (2) Evaluar r_i(t), v_i(t), r_j(t) y v_j(t).
• (3) Si |r_i(t)-r_j(t)| > R_i+R_j entonces ir a (1).
• (4) Hay colisión: v_i = (m_iv_i+m_jv_j)/(m_i+m_j). Desaparece la partícula j. R¢_i = R_i[(m_i+m_j)/m_i]^1/3 y m_i = m_i+m_j.
• (5) Ir a (1).

Notar que al introducir el choque inelástico, la energía cinética no se conserva generándose calor en el proceso. Esto hay que tenerlo en cuenta si queremos utilizar la energía total del sistema para controlar los errores de redondeo.

El último punto que hemos de tener en cuenta es la elección de la condición inicial. Para ello notemos que el momento angular de sistema se conserva aún con la presencia de la interacción inelástica. Asi pues, para garantizar que el estado final van a ser planetas girando en el mismo sentido y con momento angular definido, vamos a tomar como condición inicial que todos los planetesimales tienen el mismo momento angular inicial. Además les supondremos que tienen igual masa e igual excentricidad en la órbita, igual foco y el ángulo del eje mayor distribuido de forma aleatoria uniformemente. Asi, dado el momento angular l, la excentricidad e < 1 (para tener órbitas cerradas, y masa m del planetesimal, utilizaremos el siguiente algoritmo para darle una órbita.

• (0) Dados l, e y m.
• (1) Evaluar la energía total de la órbita: E = -m³(1-e²)/2l².
• (2) Elegir una v en el intervalo [m(1-e)/l,m(1+e)/l] mediante un número aleatorio.
• (3) Elegir un f en el intervalo [0,2p] mediante un número aleatorio. Construir v_x = vcos(f) y v_y = vsin(f).
• (4) Evaluar r = 1/(v²/2-E/m).
• (5) Elegir con probabilidad 1/2 uno de los dos valores: y = -lv_x/mv²±v_y[(mrv+l)(mrv-l)]^1/2/mv².
• (6) Evaluar x = (r²-y²)^1/2.
• (7) Evaluar w = xv_y-yv_x. Si w < 0 entonces x = -x.

Pulsar aquí para los datos astronómicos de masa, posición, ....de los planetas del sistema solar.