Grupo5 - Eq. Schroedinger

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

$i\hbar {\frac {\partial }{\partial t}}\Psi (\mathbf {r} ,t)=\left[-{\frac {\hbar ^{2}}{2m}}\nabla ^{2}+V(\mathbf {r} ,t)\right]\Psi (\mathbf {r} ,t)$

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

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

Método numérico

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

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

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

${\frac {\Psi _{j}^{n+1}-\Psi _{j}^{n}}{\Delta t}},\quad {\frac {\Psi _{j}^{n}-\Psi _{j}^{n-1}}{\Delta t}}$

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 \cite{enswork}.

O método de Crank-Nicolson consiste em uma média aritmética dos métodos explícito e implícito. Excetuando manipulações algébricas triviais, verifica-se que a relação de recorrência do método é dada por:

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

onde

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

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 $\Psi _{0}^{n}=\Psi _{j_{max}}^{n}$ para todo $n$ (ou, para as bordas $a$ e $b$ há a relação $\Psi (a,t)=\Psi (b,t)$ para todo $t$).

Condições de contorno iguais a zero

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

${\begin{pmatrix}b_{1}&a&0&\cdots &0&0&0\\a&b_{2}&a&\cdots &0&0&0\\\vdots &\vdots &\vdots &\ddots &\vdots &\vdots &\vdots \\0&0&0&\cdots &a&b_{L-2}&a\\0&0&0&\cdots &0&a&b_{L-1}\end{pmatrix}}{\begin{pmatrix}\Psi _{1}^{n+1}\\\Psi _{2}^{n+1}\\\vdots \\\Psi _{L-2}^{n+1}\\\Psi _{L-1}^{n+1}\end{pmatrix}}={\begin{pmatrix}b_{1}^{*}&a^{*}&0&\cdots &0&0&0\\a^{*}&b_{2}^{*}&a^{*}&\cdots &0&0&0\\\vdots &\vdots &\vdots &\ddots &\vdots &\vdots &\vdots \\0&0&0&\cdots &a^{*}&b_{L-2}^{*}&a^{*}\\0&0&0&\cdots &0&a^{*}&b_{L-1}^{*}\end{pmatrix}}{\begin{pmatrix}\Psi _{1}^{n}\\\Psi _{2}^{n}\\\vdots \\\Psi _{L-2}^{n}\\\Psi _{L-1}^{n}\end{pmatrix}}$

Condições de contorno periódicas

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

${\begin{pmatrix}b_{0}&a&0&\cdots &0&0&a\\a&b_{1}&a&\cdots &0&0&0\\\vdots &\vdots &\vdots &\ddots &\vdots &\vdots &\vdots \\0&0&0&\cdots &a&b_{L-1}&a\\a&0&0&\cdots &0&a&b_{L}\end{pmatrix}}{\begin{pmatrix}\Psi _{0}^{n+1}\\\Psi _{2}^{n+1}\\\vdots \\\Psi _{L-1}^{n+1}\\\Psi _{L}^{n+1}\end{pmatrix}}={\begin{pmatrix}b_{0}^{*}&a^{*}&0&\cdots &0&0&a^{*}\\a^{*}&b_{1}^{*}&a^{*}&\cdots &0&0&0\\\vdots &\vdots &\vdots &\ddots &\vdots &\vdots &\vdots \\0&0&0&\cdots &a^{*}&b_{L-1}^{*}&a^{*}\\a^{*}&0&0&\cdots &0&a^{*}&b_{L}^{*}\end{pmatrix}}{\begin{pmatrix}\Psi _{0}^{n}\\\Psi _{2}^{n}\\\vdots \\\Psi _{L-1}^{n}\\\Psi _{L}^{n}\end{pmatrix}}$

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

$\int _{-\infty }^{\infty }\,\left|\Psi (x,0)\right|^{2}\,dx=1$

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

Implementação em C

Condição de contorno limitada

void CN(double complex *u, double complex *u_aux, double complex *u_next, double complex a) { //Nomeia-se E a matriz da equação matricial que multiplica o estado evoluído no tempo. //Tem-se L, dx e dt definidos

//u -- vetor de estado atual //u_aux -- conj(E) * u //u_next -- o novo estado

//Esta função resolve resolve E * u_next = conj(E) * u

int i;

//atualização do vetor u_aux

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

//para resolver E * u_next = u_aux, utiliza-se do Método de Thomas:

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

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

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

u_next[0] = u_next[L] = 0; u_next[L-1] = d_new[L-1]; for(i = L-2; i > 0; i --) u_next[i] = d_new[i] - c_new[i] * u_next[i+1];

//atualiza-se o valor do estado: u = u_next

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

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
• Sherer, Philipp O.J., Computational Physics simulation of Classical and Quantum Systems. Springer, 2010.