Resolución de la ecuación de Schrodinger en una dimensión.

Teoría: Física Cuántica
Técnica numérica: Esquemas implícitos en al resolución de ecuaciones en derivadas parciales parabólicas.
Programación: Fortran, C, Dynamiclattice.

Introducción teórica:

La ecuación de Schrodinger describe la evolución (determinista) de un campo F llamado función de ondas. Esta ecuación en una dimensión es de la forma:

_
h

¶F(x,t)

¶t

é
ê
ê
ê
ë

_
h

¶²

¶x²

+V(x)

ù
ú
ú
ú
û

F(x,t) = HF

(1)

donde `h es la constante de Planck, F(x,t) es el campo que representa a una partícula en el instante t, m es la masa de dicha partícula, V(x) es el potencial al que se ve sometida y H es el operador Hamiltoniano (que es Hermítico). La función de ondas es la magnitud que contiene toda la información accesible a un observador externo al sistema (notar que es una variable compleja). Para evitar el uso continuado de m y `h, reescalamos el tiempo t® t`h y el espacio x® x`h/Ö[2m] lo que nos da una ecuación de Schrodinger equivalente a (1) pero con `h = 1 y m = 1/2. Esta ecuación reescalada será la que utilicemos a partir de ahora.

El carácter físico de la función de ondas fue dado por Max Born al interpretar que

P(x,t) = |F(x,t)|²

(2)

es proporcional a la probabilidad de encontrar la partícula en el punto x en el instante t. Para que sea una probabilidad, debemos de garantizar que

ó
õ

¥

-¥

dx P(x,t) = 1

(3)

esto es, las funciones de onda han de ser de cuadrado integrable luego han de pertenecer a un espacio de Hilbert. Notar que la integral de probabilidad es una constante durante la evolución.

Finalmente decir que la solución formal de la ecuación de Schrodinger es

F(x,t) = e^-itH F(x,0)

(4)

donde el operador e^itH es unitario al ser H hermítico. Esta última propiedad nos será útil a la hora de hallar la discretización más adecuada.

Problemas:

Obligatorio: Resolver la ecuación de Schrodinger unimensional para un potencial cuadrado.
Voluntario: Estudiar el coeficiente de transmisión (la probabilidad de encontrar a la partícula al otro lado del obstáculo para tiempos largos pero antes de que se refleje en la pared derecha) en función del parámetro del potencial l, para valores positivos y negativos. Estudiar lo mismo para un potencial escalón y uno triangular. Comparar con la teoría.

Método numérico.

Discretizaremos el espacio y el tiempo de forma que la función de ondas es:

F(x,t)®F(jh,ns) = f_j,n

(5)

donde n = 0,1,2,...¥, j = 0,1,...,N, h es el espaciado del retículo y s es el espaciado temporal (ambos supuestamente pequeños. Supondremos que las condiciones de contorno para la función de ondas en j = 0 y j = N son las correspondientes a la existencia de un potencial infinito en esos puntos, esto es, la probabilidad de encontrar la partícula en dichos puntos es cero. Asi, F_0,n = F_N,n = 0, para todo tiempo. En este contexto, podemos definir de forma inmediata un primer algoritmo discreto a partir de la expresión (4). Según esa expresión, si t es muy pequeña, podemos aproximar el operador de evolución exp(-itH) por su desarrollo de Taylor quedándonos a primer orden en t. Si además observamos que el operador Hamiltoniano no es más que la aplicación de una derivada segunda, su discretización espacial es obvia y en total obtenemos el siguiente algoritmo

F(j,n+1) = (1-isH_D+O((sH_D)²)F(j,n)

(6)

donde

H_D f_j = (-

h²

d_j²+V_j)f_j = -

h²

(f_j+1-2f_j+f_j-1)+V_j

(7)

y V_j = V(jh). Este algoritmo tan simple y natural tiene un problema grave: el operador (1-isH) no es unitario. Esto implica que durante la iteración del algoritmo para encontrar la evolución de la función de ondas, su normalización, å_j h |F_j,n|², irá variando con el tiempo lo que es incompatible con el carácter de probabilidad que se le tiene que dar. Asi pues, nuestro objetivo será encontrar un operador evolución similar al anterior pero que sea unitario. Esto lo conseguimos con el siguiente algoritmo:

F_j,n+1 =

1-isH_D/2

1+isH_D/2

F_j,n

(8)

Notar que además de ser unitario, este operador es exácto hasta orden (sH_D)³. El único problema que nos resta con este algoritmo es el de diseñar la estrategia utilizar eficazmente el operador 1/(1+isH_D). Para ello reescribimos la anterior ecuación como:

F_j,n+1 =

é
ê
ë

1+isH_D/2

-1

ù
ú
û

F_j,n = c_j,n-F_j,n

(9)

donde

1+isH_D/2

F_j,n = c_j,n

(10)

o lo que es lo mismo, dada F_j,n, j = 0,..,N, c_j,n es la solución de la ecuación

[ 1+isH_D/2 ]c_j,n = 2F_j,n

(11)

que en forma explícita puede escribirse como:

c_j+1,n+

é
ê
ê
ê
ë

-2+

~
s

~
V

ù
ú
ú
ú
û

c_j,n+c_j-1,n =

~
s

F_j,n

(12)

donde [s\tilde] = s/h² y [V\tilde]_j = h²V_j. De esta forma el algoritmo básico de evolución queda bastante claro: dada una función de ondas en el instante n para toda posición j se resuelve el conjunto de ecuaciones (12) y se obtiene c_j,n j = 0,..,N. Con esta solución vamos a la ecuación (9) y obtenemos las nuevas F_j,n+1 j = 0,...,N, volviendo a iterar el proceso. Lo único que nos queda por conocer es como resolver ecuaciones del tipo (12), osea, como invertir matrices tridiagonales.

En nuestro caso hemos de resolver un conjunto de ecuaciones del tipo:

A_j^-c_j-1+A_j⁰c_j+A_j⁺c_j+1 = b_j j = 1,¼,N-1

(13)

donde A_j^- = 1, A_j⁰ = -2+2i/[s\tilde]-[V\tilde]_j, A_j⁺ = 1 y b_j = 4iF_j,n/[s\tilde]. Las condiciones de contorno c₀ = c_N = 0 (notar que estas condiciones de contorno implican que en (9) F₀ y F_N son cero para todo tiempo.

Para resolver la recurrencia (13) suponemos que su solución es del tipo:

c_j+1 = a_jc_j+b_j j = 0,¼,N-1

(14)

donde para garantizar que se cumple c_N = 0 tomaremos que a_N-1 = b_N-1 = 0. Sustituyendo esta expresión en (13) obtenemos:

c_j = -

A_j^-

A_j⁰+A_j⁺a_j

c_j-1+

b_j-A_j⁺b_j

A_j⁰+A_j⁺a_j

(15)

Si identificamos las ecuaciones (14) y (15) podemos definir las recurrencias para los coeficientes a y b, asi

a_j-1 = -A_j ^-g_j

(16)

b_j-1 = g_j(b_j-A_j⁺b_j )

(17)

donde g_j^-1 = A_j⁰+A_j⁺a_j. Estas ecuaciones nos dan la forma de obtener todas las a y b partiendo de j = N-1 y obteniendo en orden decreciente a_j y b_j j = N-2, ...,1. Una vez que tenemos todas las a y b, utilizamos la recurrencia (14) para hallar las c_j en orden creciente en j.

Una vez ya sabemos como obtener c y con ella F_n+1 nos queda discutir que función de onda inicial introducimos y que potencial usamos. La función de onda inicial que vamos a usar es una onda plana con una amplitud gausiana, esto es:

F(x,0) = e^ik₀xe^{-(x-x₀)²/2s²}

(18)

notar que con esta elección, la probabilidad de encontrar inicialmente a la partícula en un punto x es una gausiana centrada en x₀ y de anchura s. El número de oscilaciones completas que la función de ondas tiene sobre la red depende de k₀ asi k₀Nh = 2pn_ciclos. En lugar de dar como parámetro inicial k₀, darémos n_ciclos. Obviamente n_ciclos = 0,1,....,N pero físicamente no tendría mucho sentido que una oscilación completa de la función de ondas se produjese entre dos puntos del retículo pues querría decir que no estamos haciendo lo suficientemente fina nuestra discretización. Asi pues el parámetro n_ciclos lo restringiremos a valores 1,...,N/4. De esta forma, un ciclo, como mínimo tendrá cuatro puntos. La posición media inicial y la anchura de la gausiana las tomamos: x₀ = Nh/4 y s = Nh/8 respectivamente. Por útimo, el potencial que usaremos será de anchura N/5, centrado en N/2 y altura proporcional a la energía de la función de ondas incidente: lk₀² donde l = 0.3 por ejemplo.

En resumen, utilizando las constantes reescaladas, la función de ondas en el retículo se tomará:

F_j,0 = e^{i[k]₀ j} e^{-2(4j-N)²/N²}

(19)

donde [k]₀ = k₀ h = 2pn_ciclos/N donde n_ciclos = 1,...,N/4. El potencial

~
V

= V_j h² = 0 si j not in [2N/5,3N/5]

(20)

= l

[k]

2
0

si j Î [2N/5,3N/5]

(21)

Lo único que nos queda por fijar es el parámetro [s] = s/h². Puesto que la energía es proporcional a k₀² y nuestro operador dinámico discreto tiende a ser exácto en potencias de H s, lo óptimo es elegir Hs < 1, esto es, k₀² s < 1, asi deducimos que [s] < 1/[k\tilde]₀². En particular tomaremos [s] = 1/4[k]₀². Notar que los parámetro que hemos de dar al sistema inicialmente son N, n_ciclos y l pues todos los demás se fijan a partir de ellos.

En resumen, el algoritmo para resolver la ecuación de Schrodinger unidimensional es:

(0) Dar los parámetros iniciales: N, n_ciclos y l. Generar [s], [k]₀, [V]_j y F_j,0 (incluyendo las condiciones de contorno F_0,0 = F_N,0 = 0.
(1) Calcular a's y b's utilizando la recurrencias (16) y (17) respectivamente.
(2) Calcular c's de (14).
(3) Calcular F_j,n+1 de (9).
(4) n=n+1, Ir a (2).

Números complejos en Fortran.

Declaración: implicit double complex (h)
Asignación de valores constantes: h = (1.0d+0,0.3d+0) (=1+i0.3).
Asignación de variables: h = complex(a,b) (=a+ib).