Grupo - Eq. d Schrödinger não-lineares acopladas

Bósons são partículas de spin inteiro e possuem funções de onda simétricas. Tais partículas obedecem às estatísticas de Bose-Einstein, de forma que a existência de múltiplas partículas (deste tipo) indistinguíveis em um mesmo estado quântico é possível.

Condensado de Bose Einstein e Equação de Schrödinger Não-Linear

Para temperaturas muito próximas de zero absoluto, uma parte considerável de um sistema de muitos bósons vai se encontrar no seu estado de menor energia, ou seja, no mesmo estado de menor energia.

A acumulação de bósons neste estado fundamental é chamada de condensado de Bose-Einstein. Um dos fatos mais interessantes do ponto de vista experimental deste tipo de condensado é que fenômenos quânticos se tornam visíveis macroscopicamente.

No caso de temperaturas ultrafrias (em concordância com a maioria dos experimentos envolvendo átomos alcalinos), a dinâmica do estado de partículas condensadas pode ser modelada por meio da equação de Schrödinger não-linear. Neste caso o condensado é descrito por uma única função de onda $\psi ({\vec {r}})$, e $\vert \psi ({\vec {r}})\vert ^{2}$ é interpretado como a densidade de partícula, e consequentemente o número total de átomos é dado por:

$N=\int d{\vec {r}}\vert \psi ({\vec {r}})\vert ^{2}$

A energia associada a este estado $\psi ({\vec {r}})$ do condensado para bósons no estado fundamental é, de acordo com a teoria de campo médio, igual a

$E=\int d{\vec {r}}\left({\frac {\hbar ^{2}}{2m}}\vert \nabla \psi ({\vec {r}})\vert ^{2}+V({\vec {r}})\vert \psi ({\vec {r}})\vert ^{2}+{\frac {1}{2}}U_{o}\vert \psi ({\vec {r}})\vert ^{4}\right)$

A minimização desta energia no que diz respeito a variações infinitesimais de $\psi ({\vec {r}})$, com um número total de átomos constante, obtém-se a seguinte equação de Schrödinger não-linear

$i\hbar {\frac {\partial \psi ({\vec {r}})}{\partial t}}=-{\frac {\hbar ^{2}}{2m}}\nabla ^{2}\psi ({\vec {r}})+V({\vec {r}})\psi ({\vec {r}})+U_{o}\vert \psi ({\vec {r}})\vert ^{2}\psi ({\vec {r}})$

Método Numérico

A equação de Schrödinger dependente do tempo pode ser escrita em notação de operadores como

$i\hbar {\frac {\partial \psi }{\partial t}}=H\psi$

onde H representa o operador Hamiltoniano do sistema. Como H é um operador linear, uma possível discretização da equação acima é dada por

$i\hbar {\frac {\psi _{j}^{n+1}-\psi _{j}^{n}}{\Delta t}}=\sum _{k=0}^{N-1}H_{jk}\psi _{k}^{n}$

Esta forma é conhecida como uma forma explícita no tempo, uma vez que o valor futuro de $\psi$ depende exclusivamente do valor atual de $\psi$. Em notação matricial, o desenvolvimento algébrico da expressão acima fornece

$\Psi ^{n+1}=\left(\mathbf {I} -{\frac {i\Delta t}{\hbar }}\mathbf {H} \right)\Psi ^{n}$

onde $\Psi ^{n}$ é o vetor coluna com os valores de $ \psi $ no tempo n, e $\mathbf {I}$ é a matriz identidade. Outra forma possível de discretização é a discretização implícita no tempo, dada por

$i\hbar {\frac {\psi _{j}^{n+1}-\psi _{j}^{n}}{\Delta t}}=\sum _{k=0}^{N-1}H_{jk}\psi _{k}^{n+1}$

em que o valor futuro de $\psi$ depende tanto dos valores atual e futuro de $\psi$. É possível mostrar que a forma matricial que corresponde à equação acima é dada por

$\Psi ^{n+!}=\left(\mathbf {I} +{\frac {i\Delta t}{\hbar }}\mathbf {H} \right)^{-1}\Psi ^{n}$

O método de Crank-Nicolson, por sua vez, consiste em tomar a média dos esquemas explícito e implícito, de forma que

$i\hbar {\frac {\psi _{j}^{n+1}-\psi _{j}^{n}}{\Delta t}}={\frac {1}{2}}\sum _{k=0}^{N-1}H_{jk}(\psi _{k}^{n}+\psi _{k}^{n+1})$

Tal método, além de possuir alta acurácia, é incondicionalmente estável na integração de muitas equações diferenciais parciais. Tomando a versão matricial da equação de Schrödinger discretizada por Crank-Nicolson, obtém-se

$\Psi ^{n+1}=\Psi ^{n}-{\frac {i\Delta t}{2\hbar }}\mathbf {H} (\Psi ^{n}+\Psi ^{n+!})$

ou ainda, por manipulação algébrica,

$\left(\mathbf {I} +{\frac {i\Delta t}{2\hbar }}\mathbf {H} \right)\Psi ^{n+1}=\left(\mathbf {I} -{\frac {i\Delta t}{2\hbar }}\mathbf {H} \right)\Psi ^{n}$

Isolando o termo $\Psi ^{n+1}$ na expressão acima, encontra-se

$\Psi ^{n}=\left(\mathbf {I} +{\frac {i\Delta t}{2\hbar }}\mathbf {H} \right)^{-1}\left(\mathbf {I} -{\frac {i\Delta t}{2\hbar }}\mathbf {H} \right)\Psi ^{n}$

Implementação e Código

A integração da equação de Schrödinger pelo método de Crank-Nicolson

$\Psi ^{n}=\left(\mathbf {I} +{\frac {i\Delta t}{2\hbar }}\mathbf {H} \right)^{-1}\left(\mathbf {I} -{\frac {i\Delta t}{2\hbar }}\mathbf {H} \right)\Psi ^{n}$

pode ser reescrita na seguinte forma

$\Psi ^{n+!}=\left(\mathbf {I} +{\frac {i\Delta t}{2\hbar }}\mathbf {H} \right)^{-1}\left[2\mathbf {I} -\left(\mathbf {I} +{\frac {i\Delta t}{2\hbar }}\mathbf {H} \right)\right]\Psi ^{n}=\left[2\left(\mathbf {I} +{\frac {i\Delta t}{2\hbar }}\mathbf {H} \right)^{-1}-\mathbf {I} \right]\Psi ^{n}$

Definindo $\mathbf {Q} ={\frac {1}{2}}\left[\mathbf {I} +{\frac {i\Delta t}{2\hbar }}\mathbf {H} \right]$, a equação acima pode ser expressa como

$\Psi ^{n}=(\mathbf {Q} ^{-1}-\mathbf {I} )\Psi ^{n}=\mathbf {Q} ^{-1}\Psi ^{n}-\Psi ^{n}$

Portanto, a evolução dinâmica do sistema pode ser avaliada separando o processo em duas partes. Primeiramente, o sistema linear

$\mathbf {Q} {\boldsymbol {\chi }}=\Psi ^{n}$

é resolvido para o vetor ${\boldsymbol {\chi }}$ e, em sequência, os valores da função de onda são atualizados de acordo com

$\Psi ^{n+1}={\boldsymbol {\chi }}-\Psi ^{n}$

Uma matriz tridiagonal possui a seguinte forma

$\mathbf {A} ={\begin{bmatrix}\beta _{0}&\gamma _{0}&0&\cdots &0\\\alpha _{0}&\beta _{1}&\gamma _{1}&\cdots &0\\0&\alpha _{1}&\beta _{2}&\cdots &0\\\vdots &\vdots &\vdots &\ddots &\vdots \\0&0&0&\cdots &\beta _{N-1}\end{bmatrix}}$

ou seja, apenas os termos da diagonal principal e os termos imediatamente acima ou abaixo são não-nulos. Ela pode ser armazenada de forma mais compacta no formato

$A_{c}={\begin{bmatrix}*&\beta _{0}&\gamma _{0}\\\alpha _{0}&\beta _{1}&\gamma _{1}\\\alpha _{1}&\beta _{2}&\gamma _{2}\\\vdots &\vdots &\vdots \\\alpha _{N-2}&\beta _{N-1}&*\end{bmatrix}}$

onde os elementos denotados por um asterisco não são utilizados.

Consequentemente, o sistema linear $ \mathbf{Ax} = \mathbf{b}$ pode ser resolvido pelo método de eliminação Gaussiana, também conhecido como algoritmo de Thomas neste caso. O procedimento é dividido em duas partes: por primeiro, a eliminação progressiva é realizada, que consiste em determinar equações lineares equivalentes que independam de $\alpha _{i}$. Neste passo, as equações recursivas para os novos elementos da diagonal principal $\beta _{i}$ e do vetor $\mathbf {b}$ são

$\beta '_{i}=\beta _{i}-{\frac {\alpha _{i-1}}{\beta '_{i-1}}}\gamma _{i-1},\;\;\;\;\;i=1,...,N-1$

e

$b'_{i}=b_{i}-{\frac {\alpha _{i-1}}{\beta '_{i-1}}}b'_{i-1},\;\;\;\;\;i=1,...,N-1$

com $\beta '_{0}=\beta _{0}$ e $b'_{0}=b_{0}$. Na etapa seguinte, em que ocorre a chamada substituição regressiva, a última equação é resolvida para obter $x_{N-1}=b'_{N-1}/\beta '_{N-1}$, e este resultado é inserido na equação anterior, gerando uma relação recursiva da forma

$x_{i}={\frac {b'_{i}-\gamma _{i}x_{i+1}}{\beta '_{i}}},\;\;\;\;\;i=N-2,...,0$

Como a matriz $\mathbf {Q}$ é tridiagonal, o sistema linear $\mathbf {Q} {\boldsymbol {\chi }}=\Psi ^{n}$ pode ser resolvido com o uso do algoritmo de Thomas para obter ${\boldsymbol {\chi }}$ e, em sequência, ${\boldsymbol {\chi }}$ é utilizado para obter $\Psi ^{n+1}$.

O código abaixo foi implementado para integrar numericamente duas equações de Schrödinger não-lineares acopladas, que foram discretizadas sobre uma rede unidimensional da seguinte forma

$i\hbar {\frac {\partial \psi _{i}}{\partial t}}=-\omega (\psi _{i-1}+\psi _{i+1})+V_{i}\psi _{i}+U_{1}\vert \psi _{i}\vert ^{2}\psi _{i}+U_{2}\vert \phi _{i}\vert ^{2}\psi$

e

$i\hbar {\frac {\partial \phi _{i}}{\partial t}}=-\omega (\phi _{i-1}+\phi _{i+1})+V_{i}\phi _{i}+U_{1}\vert \phi _{i}\vert ^{2}\phi _{i}+U_{2}\vert \psi _{i}\vert ^{2}\phi _{i}$

que representam a evolução dinâmica de dois condensados de Bose-Einstein nesta mesma rede. O último termo de ambas as equações denota o acoplamento do sistema.

Para esta integração numérica, foram utilizados $N=301$ sítios, para um período total de 50 “segundos” e um passo de tempo de $\Delta t=0,001$ “segundos”, onde as unidades de tempo possuem aspas pois são adimensionais $\tau =\omega t/\hbar$.