Equação de Schroedinger

A evolução temporal do estado quântico $\mathrm{\Psi }(\mathbf{𝐫},t)$ é dada pela equação de Schrödinger, a qual é postulada como:

$i\hslash \frac{\partial }{\partial t}\mathrm{\Psi }(\mathbf{𝐫},t)=[-\frac{{\hslash }^2}{2m}{\mathrm{\nabla }}^2+V(\mathbf{𝐫},t)]\mathrm{\Psi }(\mathbf{𝐫},t)$

Posto em unidades atômicas (onde $m_e$ e $\hslash$ são unitários), o caso unidimensional de um elétron num potencial independente do tempo reduz-se a:

$\frac{\partial }{\partial t}\mathrm{\Psi }(x,t)=[\frac{i}{2}\frac{{\partial }^2}{\partial x^2}-iV(x)]\mathrm{\Psi }(x,t)$

Método numérico

Buscando resolver a equação numericamente, tem-se a discretização de $\frac{{\partial }^2\mathrm{\Psi }}{\partial x^2}$ :

$\frac{{\mathrm{\Psi }}_{j-1}^n-2{\mathrm{\Psi }}_j^n+{\mathrm{\Psi }}_{j+1}^n}{{\left(\mathrm{\Delta }x\right)}^2}$

e as discretizações de $\frac{\partial \mathrm{\Psi }}{\partial t}$ (explícita e implícita, respectivamente):

$\frac{\partial \mathrm{\Psi }}{\partial t}=\frac{{\mathrm{\Psi }}_j^{n+1}-{\mathrm{\Psi }}_j^n}{\mathrm{\Delta }t}$ (explícita)

$\frac{\partial \mathrm{\Psi }}{\partial t}=\frac{{\mathrm{\Psi }}_j^n-{\mathrm{\Psi }}_j^{n-1}}{\mathrm{\Delta }t}$ (implícita)

No método implícito as diferenças são tomadas no tempo $n+1$ ao invés de tomá-las no tempo $n$ como no método explícito.

$\frac{{\mathrm{\Psi }}_j^{n+1}-{\mathrm{\Psi }}_j^n}{\mathrm{\Delta }t}=\frac{{\mathrm{\Psi }}_{j+1}^{n+1}-2{\mathrm{\Psi }}_j^{n+1}+{\mathrm{\Psi }}_{j-1}^{n+1}}{{\left(\mathrm{\Delta }x\right)}^2}$

fazendo $\lambda =\frac{\mathrm{\Delta }t}{(\mathrm{\Delta }x{)}^2}$, temos:

${\mathrm{\Psi }}_j^n=-\lambda {\mathrm{\Psi }}_{j+1}^{n+1}+(1+2\lambda ){\mathrm{\Psi }}_j^{n+1}-\lambda {\mathrm{\Psi }}_{j-1}^{n+1}$

Considerando $0\le j\le M$ temos $M-1$ equações simultâneas. O método implícito converge à solução da EDP desde que $\mathrm{\Delta }x\to 0$ e $\mathrm{\Delta }t\to 0$ independente do valor de $\lambda$.

Para obter os valores no dado tempo $n+1$ se resolve o conjunto de equações simultaneamente dado pela equação acima que pode ser escrita na forma matricial,

$\left(\begin{array}{ccccccc}-\lambda & (1+2\lambda )& -\lambda & \cdots & 0& 0& 0\\ 0& -\lambda & (1+2\lambda )& -\lambda & \cdots & 0& 0\\ ⋮& ⋮& \ddots & \ddots & \ddots & ⋮& ⋮\\ 0& 0& 0& \cdots & -\lambda & (1+2\lambda )& -\lambda \end{array}\right)\left(\begin{array}{c}{\mathrm{\Psi }}_0^{n+1}\\ {\mathrm{\Psi }}_1^{n+1}\\ ⋮\\ {\mathrm{\Psi }}_{j-1}^{n+1}\end{array}\right)=\left(\begin{array}{c}{\mathrm{\Psi }}_0^n\\ {\mathrm{\Psi }}_1^n\\ ⋮\\ {\mathrm{\Psi }}_{j-1}^n\end{array}\right)$

Tanto no método explícito quanto no método implícito não é conservada a norma do estado (o que é estritamente necessário, já que o estado pode ser interpretado como uma onda de probabilidade). Por esse motivo, utiliza-se o método de Crank-Nicolson, o qual tem essa propriedade.

O método de Crank-Nicolson consiste em uma média aritmética dos métodos explícito e implícito. Essa combinação dá a estabilidade do método implícito e a precisão do método explícito, apesar de causar oscilações numéricas (podendo ser contornadas usando passos de tempo menores). Excetuando manipulações algébricas triviais, verifica-se que a relação de recorrência do método é dada por:

$a{\mathrm{\Psi }}_{j-1}^{n+1}+b_j{\mathrm{\Psi }}_j^{n+1}+a{\mathrm{\Psi }}_{j+1}^{n+1}=a^{*}{\mathrm{\Psi }}_{j-1}^n+b_j^{*}{\mathrm{\Psi }}_j^n+a^{*}{\mathrm{\Psi }}_{j+1}^n,$

onde

$a\equiv -\frac{i\mathrm{\Delta }t}{4(\mathrm{\Delta }x{)}^2}$ e $b_j\equiv 1+\frac{i\mathrm{\Delta }t}{2}[\frac{1}{(\mathrm{\Delta }x{)}^2}+V(j\mathrm{\Delta }x)]$.

A integração numérica depende, portanto, do potencial em que o elétron está sujeito, bem como da sua condição inicial e suas das condições de contorno.

Que condições podemos impor para a fronteira? Quando se trata do problema analiticamente, costuma-se considerar que a função de onda tende a zero no infinito. Numericamente, pode-se fazer uma transposição disso, criando uma condição para bordas em pontos suficientemente distantes do centro da distribuição da função de onda, igualando-as a zero. Outra forma de tratar o problema numericamente é criando condições de contorno periódicas, em que para as bordas vale ${\mathrm{\Psi }}_0^n={\mathrm{\Psi }}_{j_{max}}^n$ para todo $n$ (ou, para as bordas $a$ e $b$ há a relação $\mathrm{\Psi }(a,t)=\mathrm{\Psi }(b,t)$ para todo $t$).

Condições de contorno iguais a zero

Para as condições de contorno ${\mathrm{\Psi }}_0^n={\mathrm{\Psi }}_L^n=0$, a iteração reduz-se à equação matricial:

$\left(\begin{array}{ccccccc}{b}_1& a& 0& \cdots & 0& 0& 0\\ a& b_2& a& \cdots & 0& 0& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ 0& 0& 0& \cdots & a& b_{L-2}& a\\ 0& 0& 0& \cdots & 0& a& b_{L-1}\end{array}\right)\left(\begin{array}{c}{\mathrm{\Psi }}_1^{n+1}\\ {\mathrm{\Psi }}_2^{n+1}\\ ⋮\\ {\mathrm{\Psi }}_{L-2}^{n+1}\\ {\mathrm{\Psi }}_{L-1}^{n+1}\end{array}\right)=\left(\begin{array}{ccccccc}{b}_1^{*}& a^{*}& 0& \cdots & 0& 0& 0\\ a^{*}& b_2^{*}& a^{*}& \cdots & 0& 0& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ 0& 0& 0& \cdots & a^{*}& b_{L-2}^{*}& a^{*}\\ 0& 0& 0& \cdots & 0& a^{*}& b_{L-1}^{*}\end{array}\right)\left(\begin{array}{c}{\mathrm{\Psi }}_1^n\\ {\mathrm{\Psi }}_2^n\\ ⋮\\ {\mathrm{\Psi }}_{L-2}^n\\ {\mathrm{\Psi }}_{L-1}^n\end{array}\right)$

Condições de contorno periódicas

De maneira semelhante, a iteração do caso das condições de contorno periódicas - ${\mathrm{\Psi }}_0^n={\mathrm{\Psi }}_L^n$ - reduz-se à equação matricial:

$\left(\begin{array}{ccccccc}{b}_0& a& 0& \cdots & 0& 0& a\\ a& b_1& a& \cdots & 0& 0& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ 0& 0& 0& \cdots & a& b_{L-1}& a\\ a& 0& 0& \cdots & 0& a& b_L\end{array}\right)\left(\begin{array}{c}{\mathrm{\Psi }}_0^{n+1}\\ {\mathrm{\Psi }}_2^{n+1}\\ ⋮\\ {\mathrm{\Psi }}_{L-1}^{n+1}\\ {\mathrm{\Psi }}_L^{n+1}\end{array}\right)=\left(\begin{array}{ccccccc}{b}_0^{*}& a^{*}& 0& \cdots & 0& 0& a^{*}\\ a^{*}& b_1^{*}& a^{*}& \cdots & 0& 0& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ 0& 0& 0& \cdots & a^{*}& b_{L-1}^{*}& a^{*}\\ a^{*}& 0& 0& \cdots & 0& a^{*}& b_L^{*}\end{array}\right)\left(\begin{array}{c}{\mathrm{\Psi }}_0^n\\ {\mathrm{\Psi }}_2^n\\ ⋮\\ {\mathrm{\Psi }}_{L-1}^n\\ {\mathrm{\Psi }}_L^n\end{array}\right)$

Condição inicial

Já a condição inicial é arbitrária, pois define o estado inicial do sistema que queremos tratar. Fazendo uma referência ao tratamento de sistemas clássicos, seria como definir posição e momento iniciais. É claro que, para ter o sentido físico de uma função de onda, deve-se ter o cuidado de criar uma condição inicial normalizada, satisfazendo

$\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{|\mathrm{\Psi }(x,0)|}^{2}dx=1$

bastando, então, inseri-la no programa.

Implementação em C

#include <stdio.h>
#include <math.h>
#include <complex.h> 

#define N_steps 100000
#define L 1000
#define dt 1
#define dx 4.0
#define w 0.002

double V(int x_rel)
{
	return - pow(w,2) * pow(x_rel - 500,2) / 2.0; //SHO

	//return 0; 
}

double complex b(int i)
{
	return 0.5 * dt * I * (1.0 / pow(dx,2.0) + V(i)) + 1.0;
}

double u_0(int x_rel)
{
	//return sqrt(2 / L) * sin(5*x_rel * M_PI / L); //eigenstate of infinite square well

	//return (1.0/ sqrt(5 * M_PI)) * exp(I * 0.5 * x_rel) * exp(-pow(x_rel - 500, 2) / (2 * 25)); //gaussian packet

	return  (w / M_PI) * exp(-w * pow(x_rel - 500, 2) / 2.0);
}

void CN(double complex *u, double complex *u_aux, double complex *u_next, double complex a, int bi)
{
	//bi = 1 (bounded); bi = 0 (periodic)

	int i;

	for(i = bi; i <= L - bi; i++) u_aux[i] = conj(a) * (u[(i-1+L)%L] + u[(i+1+L)%L]) + conj(b(i)) * u[i];

	//thomas algorithm

	double complex c_new[L+1], d_new[L+1];

	c_new[bi] = a / b(bi);
	for(i = 1 + bi; i <= L - bi; i++) c_new[i] = a / (b(i) - c_new[i-1] * a);

	d_new[bi] = u_aux[bi] / b(bi);
	for(i = 1 + bi; i <= L - bi; i++) d_new[i] = (u_aux[i] - d_new[i-1] * a) / (b(i) - c_new[i-1] * a);

	if (bi == 1) u_next[0] = u_next[L] = 0;
	u_next[L-bi] = d_new[L-bi];
	for(i = L-1-bi; i >= bi; i --) u_next[i] = d_new[i] - c_new[i] * u_next[i+1];
	
	//u = u_next

	for(i = 0; i <= L; i++) u[i] = u_next[i];
} 

int main(void)
{
	int i, j, n = 0;

	double complex u[L+1], u_next[L+1], u_aux[L+1], a = - 0.25 * I * dt / pow(dx,2.0);

	//Initial Contition
	for(i = 0; i <= L; i ++) u[i] = u_0(i);

	while(n < N_steps)
	{
		CN(u, u_aux, u_next, a, 0);
		printf("set title 'Time = %d'\nplot \'-' w lp pt 7 ps 0.1", n);
		for(i = 0; i <= L; i ++) printf("%d\t%.10lf\n",i,pow(creal(u[i]),2) + pow(cimag(u[i]),2));
		//for(i = 0; i <= L; i ++) printf("%d\t%lf\n",i,cimag(u[i]));
		//for(i = 0; i <= L; i ++) printf("%d\t%lf\n",i,creal(u[i]));
		printf("e\npause 0.1\n"); 
		n ++;
	}

	return 0;
	
}

Aplicações

Linearidade da função de onda

Pacote gaussiano estacionário num oscilador harmônico

Pacote gaussiano tunelando uma berreira de potencial

Referências

• Cohen-Tannoudji C., Diu B., Laloe F. Quantum mechanics. Volume 1. Wiley, 1991.
• Numerical Resolution Of The Schrödinger Equation. Jorgensen L., Lopes Cardozo D., Thivierge E. http://web.pa.msu.edu/people/duxbury/courses/phy480/SchrodingerDynamics.pdf
• Crank, J.; Nicolson, P. (1947). "A practical method for numerical evaluation of solutions of partial differential equations of the heat conduction type". Proc. Camb. Phil. Soc. 43 (1): 50–67. doi:10.1007/BF02127704.
• Sherer, Philipp O.J., Computational Physics simulation of Classical and Quantum Systems. Springer, 2010.
• Born M., Nobel lecture: The statistical interpretation of quantum mechanics. 11 de Dezembro de 1954. https://www.nobelprize.org/nobel_prizes/physics/laureates/1954/born-lecture.pdf