Problema del átomo de hidrógeno

De luz-wiki

Resumen

En este texto propongo el método de diferencias finitas para obtener una solución numérica del problema del átomo de hidrógeno, es decir, resolver la ecuación de Schrödinger independiente del tiempo. Finalmente implementé el método en python para obtener una representación gráfica de la probabilidad de densidad de la posición del electrón en función de su distancia al núcleo.

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 $m_p$. 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,

$R=\frac{m_e r_e + m_p r_p}{m_e + m_p}$

donde $r_e$ y $r_p$ 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 $r=r_e=R_p$, donde

$r_e = R + \frac{m_p}{m_e + m_p} r$

$r_p = R - \frac{m_e}{m_e + m_p} r$


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

$μ = \frac{m_e m_p}{m_e + m_p}$

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?

Solución cuántica

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

$(\frac{ħ}{2μ}ᐁ^{2} - \frac{e^{2}}{4 π ε_0 r}) Ψ(r,t) = iħ \frac{∂}{∂t} Ψ(r,t)$

donde $ħ$ es la constante de Planck reducida, $ᐁ$ es el operador diferencial nabla, el término $\frac{e^{2}}{4 π ε_0 r}$ corresponde al potencial de Coulomb y $ Ψ(r,t)$ es la función de onda de la partícula en una dimensión.

Ya que el potencial de Coulomb no depende explícitamente del tiempo, propondré una solución "separada", un producto de funciones,

$Ψ(r,t) = Ω(r) e^{\frac{i E t}{ħ}}$


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, Ω(r) y es por lo tanto la ecuación que me interesa resolver, enunciada cómo:

$(\frac{ħ}{2μ}ᐁ^{2} - \frac{e^{2}}{4 π ε_0 r}) Ω(r) = EΩ(r)$

Ahora es conveniente expresar la ecuación en términos de coordenadas esféricas $θ,Φ,r$.

$[\frac{ħ^{2}}{2μ r^2} ( r^{2} \frac{∂}{∂r}) + \frac{1}{\mathrm{sen}(θ)} \frac{∂}{∂θ} (\mathrm{sen}(θ)\frac{∂}{∂θ} + \frac{1}{\mathrm{sen}^{2}(θ)} \frac{∂^{2}}{∂Φ^{2}} - \frac{e^{2}}{4 π ε_0 r}]Ψ = EΨ $

La solución de esta ecuación tiene 3 partes, la solución azimutal $Ψ(Φ)$, la solución polar $Ψ(θ)$ y la solución radial $Ψ(r)$. Por el alcance de este texto nos interesa sólo la parte radial de la solución.

Solución radial

La solución radial de la ecuación de Schrödinger introduce una variable nueva, el número cuántico orbital $l$, que adopta valores enteros positivos. Así, la ecuación radial queda,


$\frac{1}{R} \frac{∂}{∂r} ( r^{2} \frac{∂R}{∂r}) - \frac{2μ r^{2}}{ħ^{2}} (\frac{e^{2}}{4 π ε_0 r}) = -l(l+1) $


Con algo de álgebra [1] podemos reescribir la ecuación anterior cómo, \begin{equation} \frac{∂^{2} R}{∂ r^{2}} + \frac{2}{r} \frac{∂R}{∂r} - \frac{2μ}{ħ^{2}} (\frac{e^{2}}{4 π ε_0 r} - \frac{-l(l+1) ħ^{2}}{2μ r^{2}} - E) R = 0 \label{radial} \end{equation}

Ahora me tomaré un momento para explorar cómo se comporta ésta solución para valores muy grande de $r$, donde hago despreciables los términos $r^{-1}$ y $r^{-2}$, es decir

$r → ∞: \frac{∂^{2} R_∞}{∂r^{2}} + \frac{2μE}{ħ^{2}} R_∞ = 0$

la solución de la ecuación anterior es una combinación lineal de exponenciales reales pero para evitar divergencia para valores valores muy grandes de $r$, me quedo sólo con una exponencial decreciente,

$R_∞ (r) = C e^{-\frac{\sqrt{2μ|E|}}{ħ}r} = \frac{1Ry}{n^{2}}$

La forma final de la solución radial se obtiene de una serie de potencias. No es mi intención abundar en cómo se soluciona la ecuación analíticamente así que bastará mostrar que (\ref{radial}) queda cómo

$R = C e^{-\frac{r}{n a*_0}} (\frac{2r}{n a*_0})^{l} L_{n-l-1}(\frac{2r}{n a*_0})$

donde $a*_0$ es el radio de Bohr reducido,

$a*_0 = \frac{4π ε_0 ħ^{2} }{μ e^{2}}$

$n$ es el número cuántico principal $nЄ(1,2,3,...)$ y la energía queda cómo

$E = \frac{μ e^{4}}{32 π^{2} ε_0^{2} ħ^{2} n^{2}} = \frac{1Ry}{n^{2}} $

donde Ry se refiere a 1 Rydberg, $1 Ry = 13.605 eV$ y es claro que la energía está cuantizada.

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:

$\frac{df(x)}{dx} = lim_{h→0} \frac{f(x+h)-f(x)}{h} $


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,

\begin{equation} \frac{d^{2}f(x)}{dx^{2}} ≈ \frac{\frac{f(x+h) - f(x)}{h} - \frac{f(x) - f(x-h)}{h}}{h} = \frac{f(x-h) - 2f(x) + f(x+h)}{h^{2}} \label{segunda} \end{equation}


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$ y $x+h$. Puedo elegir más puntos para la aproximación pero ésto requiere expandir el operador en serie de Taylor. Por ahora bastará mostrar a manera de ejemplo cómo se ve la solución numérica que obtengo con el método. Para ello, considero la ecuación de onda de una partícula en una caja,

$\frac{∂^{2}}{∂x^{2}} ψ(x) = -\frac{2m}{ħ} Eψ(x) $

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 con índices,

$\frac{∂^{2}}{∂x^{2}} ψ_i = -\frac{2m}{ħ} Eψ_i$ → $ᗄ_i Є (1,...,N)$


Donde $N+2$ es el número total de puntos en el dominio discretizado. Utilizando la ecuación (\ref{segunda}) para reescribir la ecuación de onda,

$\frac{ψ_{i-1} - 2ψ_i + ψ_i+1}{h^{2}} = -\frac{2m}{ħ^{2}} Eψ_i$

Puedo expresar este sistema de ecuaciones de forma matricial,

$\begin{bmatrix} -2 & 1 & & & \\ 1 & -2 & 1 & & \\ & 1 & -2 & \ddots & \\ & & \ddots & \ddots & 1\\ & & & 1 & -2 \\ \end{bmatrix} * \begin{bmatrix} ψ_1 \\ ψ_2 \\ \vdots \\ ψ_n \end{bmatrix} = \frac{-2mE}{ħ} * \begin{bmatrix} ψ_1 \\ ψ_2 \\ \vdots \\ ψ_n \end{bmatrix}$


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 cuantizados de energía. Con respecto a las condiciones de frontera es pertinente un comentario adicional: la función de onda de una partícula en una caja debe cumplir las condiciones de frontera de Dirilecht, $ ψ_0 = ψ_N+1 = 0$. Ésto implica dos soluciones triviales pero separadas para los valores $ ψ_0$ y $ψ_N+1 $. Cómo el valor es 0, el operador diferencial no se ve afectado por su presencia y se pueden implementar condiciones de frontera igual a 0 ignorando los valores $ ψ_0 y ψ_N+1 $, y me quedo con las N ecuaciones del dominio discretizado.

Discretización de la coordenada radial por diferencias finitas

Retomando la ecuación (\ref{radial}), es conveniente sustituir $ρ = rR → R=\frac{ρ}{r}$, por lo que reescribo el operador cómo,

$\frac{∂}{∂r} (r^{2} \frac{∂}{∂r} \frac{ρ}{r}) = \frac{∂}{∂r} [r^{2} \frac{1}{r} \frac{∂ρ}{∂r} + r^{2} ρ (-\frac{1}{r^{2}})] = \frac{∂ρ}{∂r} + r \frac{∂^{2}ρ}{∂r^{2}}-\frac{∂ρ}{∂r} = r \frac{∂^{2}ρ}{∂r^{2}} $

Por lo tanto,

$r \frac{∂^{2}ρ}{∂r^{2}} - \frac{2μr^{2}}{ħ^{2}}(\frac{e^{2}}{4 π ε_0 r} - \frac{l(l+1)ħ^{2}}{2μr^{2}} - E ) \frac{ρ}{r} = 0 $

dividido por r,

$ \frac{∂^{2}ρ}{∂r^{2}} - \frac{2μ}{ħ^{2}}\frac{e^{2}}{4 π ε_0 r} - \frac{l(l+1)ħ^{2}}{2μr^{2}} - E ) ρ = 0 $

Divido por el término constante y pasó el término que involucra a la energía a la derecha,

$- \frac{ħ}{2μ} \frac{∂^{2}ρ}{∂r^{2}} + \frac{e^{2}}{4 π ε_0 r}ρ - \frac{l(l+1)ħ^{2}}{2μr^{2}}ρ = Eρ $

Discretizo 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 corresponde al operador de Laplace,

$-\frac{ħ^{2}}{2μ}*\begin{bmatrix} -2 & 1 & & & \\ 1 & -2 & 1 & & \\ & 1 & -2 & \ddots & \\ & & \ddots & \ddots & 1\\ & & & 1 & -2 \\ \end{bmatrix}$


junto con otras dos matrices diagonales,

$\frac{e^{2}}{4πε_0}\begin{bmatrix} \frac{1}{r_1} & & & \\ & \frac{1}{r_2} & & \\ & & \ddots & \\ & & &\frac{1}{r_N} \\ \end{bmatrix}$

por último,

$\frac{ħ^{2}}{2μ}\begin{bmatrix} \frac{1}{r_1^{2}} & & & \\ & \frac{1}{r_2^{2}} & & \\ & & \ddots & \\ & & &\frac{1}{r_N^{2}} \\ \end{bmatrix}$

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 electrón en el átomo de hidrógeno gracias a la sustitución que adopté antes, $ρ= rR$.

Implementación en Python

Para la implementación definimos funciones que calculan las matrices diagonales señaladas en la sección anterior. Estas funciones dependen de la coordenada discretizada $r$. Con las funciones construimos un hamiltoniano con forma matricial. El siguiente paso es establecer el número de elementos $N$ de nuestro espacio discretizado, el número cuántico orbital $l$ y una región radial de largo 20 $nm$ donde excluimos el origen coordenado para evitar una indeterminación. A grandes rasgos, el código consiste en un algoritmo iterativo que diagonaliza sólo una fracción de todos los eigenestados. En presencia del potencial de Coulomb obtendremos eigenvalores negativos (bound states). Por la finitud del espacio discretizado y las condiciones de frontera, el espectro de energías positivas es análogo al de una partícula en una caja. Despreciamos estos estados porque nuestro enfoque está en el potencial negativo, es decir que sólo me quedo con eigenestados de menor magnitud. Con éstos últimos el programa calcula la probabilidad de densidad correspondiente a cada estado energético. La probabilidad de densidad se refiere a la probabilidad de encontrar el electrón a una distancia $x$ del núcleo.

Probabilidad de densidad para los primeros 3 niveles energéticos del átomo de hidrógeno

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 (línea azul) , 2s (línea verde) y 3s (línea roja). Se sabe que para estos estados energéticos 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$.

Conclusión

A pesar de que se obtuvieron resultados esperados por la magnitud de la energía para cada nivel, valdría la pena explorar la incertidumbre asociada a éste método de solución y compararlo con otros métodos numéricos como el de elementos finitos. Por ahora hemos discutido cómo el método de diferencias finitas se utiliza para discretizar un operador diferencial y lo he utilizado para resolver el problema del átomo de hidrógeno. La implementación en código nos permite tomar la solución radial de la ecuación de Schrödinger para visualizar los niveles de energía y su probabilidad de densidad asociada.

Referencias

CODATA Value: Rydberg constant times hc in eV. (2018). The NIST Reference on Constants, Units, and Uncertainty. https://physics.nist.gov/cgi-bin/cuu/Value?rydhcev

French, A. P., & Taylor, E. F. (1978). Introduction to Quantum Physics (M.I.T. Introductory Physics Series) (1st ed.). W. W. Norton & Company.

Kong, Q., Siauw, T., & Bayen, A. (2020). Python Programming and Numerical Methods: A Guide for Engineers and Scientists (1st ed.). Academic Press.

Schroeder, D. V. (2015). The radial equation. In The radial equation (pp. 1–2). https://physics.weber.edu/schroeder/quantum/RadialEquation.pdf

  1. Schroeder, D. V. (2015). The radial equation. In The radial equation (pp. 1–2).https://physics.weber.edu/schroeder/quantum/RadialEquation.pdf