Problema del átomo de hidrógeno

Planteamiento del problema

Considere un átomo de hidrógeno con un electrón de carga -e y masa m_e y un protón con carga +e y masa mp. Dadas las propiedades del sistema, en particular, que mp>>me, tomaremos al protón cómo una partícula cargada con un comportamiento cómo se describiría clásicamente. Así, obtenemos un problema clásico de mecánica de dos cuerpos. Ahora definimos el centro de masa de nuestro sistema,

Archivo:Ecuacioncentrodemasa.gif

donde r^e y rp se refieren a la posición del electrón y el protón respectivamente. Para un sistema aislado, sin fuerzas externas actuando sobre él, el centro de masa se mueve a velocidad constante por lo tanto si expresamos la posición de ambas partículas en función de la posición centro de masa R y sus posiciones relativas,

Archivo:Posición electrón.gif

Archivo:Posicionproton.gif

es claro que hemos reducido el problema a una variable, la posición relativa r, por lo que podemos resolver el sistema clásicamente introduciendo la variable de masa reducida

Archivo:CodeCogsEqn (3).gif

Aunque la complejidad del problema se ha reducido es necesario hacer algunas aclaraciones. Aplicando métodos de solución clásicos encontraríamos una respuesta más bien trivial en la medida que el electrón y el protón se atraen hasta colisionar. Por otro lado, si le otorgamos al electrón suficiente momento para establecer una órbita alrededor del núcleo, ésta no duraría mucho tiempo ya que las partículas cargadas aceleradas emiten radiación por lo que un protón en órbita circular bajo los efecto de la aceleración centrípeta, rápidamente se precipitaría al núcleo por la pérdida de energía cinética y el resultante decaimiento de su órbita.

Así, hemos planteado uno de los problemas fundacionales de lo que hoy se conoce cómo mecánica cuántica: si la teoría clásica nos dice que el átomo de hidrógeno no debería ser estable, ¿por qué lo es?

Mi objetivo es mostrar un método de solución para la función de onda del átomo de hidrógeno e implementarlo en Python para obtener una representación visual de la probabilidad de densidad del electrón en la órbita del átomo de hidrógeno.

Solución cuántica

A partir de aquí es necesario introducir la ecuación de Schrödinger dependiente del tiempo,

Archivo:CodeCogsEqn (4).gif

sin embargo el potencial de Coulomb no depende explícitamente del tiempo así que propondré una solución de la forma

Archivo:CodeCogsEqn (5).gif

La parte temporal de la solución es una función periódica exponencial y compleja mientras que la ecuación de Schrödinger independiente del tiempo se asocia a la parte espacial y es por lo tanto la ecuación que me interesa resolver, enunciada cómo:

Archivo:CodeCogsEqn (6).gif

Método de diferencias finitas

El método de diferencias finitas se utiliza para resolver numéricamente ecuaciones diferenciales. La definición de la derivada en forma de límite es:

Archivo:CodeCogsEqn.gif

Lo que busco es aproximar el operador diferencial con pequeñas diferencias. Bajo éste esquema, el dominio de la función está discretizado en escalones finitos de h. Haciendo el desplazamiento h cada vez menor se obtiene una aproximación arbitrariamente precisa. Así, para aumentar la precisión de éste método podemos considerar una cantidad mayor de valores cercanos. Para la segunda derivada,

lo que obtengo es una diferencia de diferencias, por lo tanto es conveniente considerar al menos 3 puntos en mi aproximación: x-h, x, x+h.

Ahora puedo, a manera de ejemplo, utilizar éste método para tratar la ecuación de onda de una partícula en una caja,

Archivo:Part en caja.svg

Voy a discretizar el dominio utilizando escalones de valor h. La solución es la lista de valores Ψ(0), Ψ(h), Ψ(2h),...,Ψ((N+1)h). Para estos valores, puedo expresar una ecuación de la forma,

Archivo:Part en caja 2.svg

Donde N+2 es el número total de puntos en el dominio discretizado. Utilizando la forma del método de diferencias finitas para segundas derivadas se obtiene,

Archivo:Disc ej.svg

Notemos que el sistema de ecuaciones puede expresarse como la ecuación matricial,

Archivo:Eqn matr.svg

Se ha obtenido un problema de eigenvalores; al resolverlo se determinan los eigenvectores Ψ que representan la solución así como los eigenvalores que representan los niveles escalonados de energía. Con respecto a las condiciones de frontera, la función de onda de una partícula en una caja debe cumplir las condiciones de frontera de Dirilecht Ψ₀=Ψ_N+1=0. Ésto implica dos soluciones triviales pero separadas para los valores Ψ₀ y Ψ_N+1. Cómo el valor es 0 por diseño, el operador no se ve afectado por su presencia y se pueden implementar condiciones de frontera igual a 0 ignorando los valores Ψ₀ y Ψ_N+1, reteniendo las N ecuaciones del dominio discretizado.

Solución radial de la ecuación de Schrödinger independiente del tiempo

Ahora voy a aplicar el método para la ecuación radial que se deriva resolviendo la ecuación de Schrödinger,

Archivo:CodeCogsEqn (1).gif

Utilizaré la sustitución,

Archivo:Sust rad shr.gif

Así, el operador diferencial queda cómo

Archivo:Svg dif op.svg

Archivo:Svg 2.svg

Dividiendo entre r,

Archivo:Svgentrer.svg

y ahora dividiendo por la constante del segundo término,

Archivo:Svg 3 final.svg

Discretización de la coordenada radial por diferencias finitas

Discretizaré la coordenada radial usando una red equidistante r_i con N elementos y un desplazamiento igual a h=r_i+1-r_i. El Hamiltoniano discretizado consiste de 3 términos donde el primero se corresponde con el operador de Laplace

mat

junto con otras dos matrices diagonales,

Archivo:Matc.svg

mat disc

La suma de estas 3 matrices representa la matriz hamiltoniana que se quiere diagonalizar para resolver el problema de eigenvalores. Los eigenvalores son los estados energéticos y sus eigenvectores ρ están relacionados a la parte radial de la ecuación de onda del hidrógeno gracias a la sustitución que adoptamos antes, ρ= rR.

Implementación en Python

Para esta simulación elegí el número de momento angular cuántico l = 0. Por lo tanto en la gráfica se muestran resultados para los primeros 3 niveles energéticos del átomo de hidrógeno (1s,2s,3s,...). Se sabe que para estos estados corresponden energías de 1 Ry, 1/4 Ry, 1/9 Ry respectivamente lo cual es exactamente lo que se obtiene en la simulación.

Modificar el valor de l no influye en el espectro de energía obtenido, con la excepción de que elimina ciertos estados energéticos. Esto se debe a que para un número cuántico principal n, los valores permitidos de l son 0,1,...,n-1. Es decir que si elegimos l=1 la solución de menor energía corresponde a n=2.