Una solución a la gran ecuación

La ecuación de Schrödinger 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:

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. 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 Schrödinger 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

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

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 Schrödinger es

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.

• Resolver la ecuación de Schrödinger unimensional para un potencial cuadrado.

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

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

donde

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:

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:

donde

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

que en forma explícita puede escribirse como:

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:

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:

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:

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

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 cómo obtener c y con ella F_n+1 nos queda discutir qué función de onda inicial introducimos y qué potencial usamos. La función de onda inicial que vamos a usar es una onda plana con una amplitud gausiana, esto es:

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₀ , así 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. Así 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á:

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

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 exacto en potencias de Hs, 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 Schrödinger 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).