Exibir código-fonte para Grupo - Equações de Schrödinger não-lineares acopladas

Enrique Augusto Tiran Calderoli (00276132)

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 <math> \psi(\vec{r}) </math>, e <math> \vert \psi(\vec{r})\vert^2 </math> é interpretado como a densidade de partícula. Consequentemente o número total de átomos é dado por:

<math> N = \int d\vec{r}\vert \psi(\vec{r})\vert ^2 </math>

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

<math> 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) </math>

onde <math>m</math> é a massa dos bósons, <math>V(\vec{r})</math> é o potencial externo e <math>U_o</math> representa as interações inter-partícula. A minimização desta energia no que diz respeito a variações infinitesimais de <math>\psi(\vec{r}) </math>, com um número total de átomos constante, permite obter a seguinte equação de Schrödinger não-linear

<math> 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})</math>

== Método Numérico==

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

<math> i\hbar \frac{\partial \psi}{\partial t} = H\psi    </math>

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

<math> i\hbar\frac{\psi^{n+1}_j - \psi^n_j}{\Delta t} = \sum_{k=0}^{N-1} H_{jk}\psi_{k}^n </math>

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

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

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

<math> i\hbar\frac{\psi_j^{n+1}-\psi_j^n}{\Delta t} = \sum_{k=0}^{N-1} H_{jk}\psi_{k}^{n+1}     </math>

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

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

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

<math> i\hbar \frac{\psi^{n+1}_j - \psi^{n}_j}{\Delta t} = \frac{1}{2}\sum_{k=0}^{N-1} H_{jk}(\psi_k^n + \psi_k^{n+1}) </math>

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

<math> \Psi^{n+1} = \Psi^n - \frac{i\Delta t}{2\hbar}\mathbf{H}(\Psi^n + \Psi^{n+1})                    </math>

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

<math> \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  </math>

Isolando o termo <math> \Psi^{n+1} </math> na expressão acima, encontra-se

<math> \Psi^{n+1} = \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 </math>


== Implementação e Código==

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

<math> \Psi^{n+1} = \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 </math>

pode ser reescrita na seguinte forma

<math> \Psi^{n+1} =  \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  </math>

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

<math> \Psi^{n+1} = (\mathbf{Q}^{-1} - \mathbf{I})\Psi^n = \mathbf{Q}^{-1}\Psi^n - \Psi^n                     </math>

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

<math> \mathbf{Q}\boldsymbol{\chi} = \Psi^n </math>

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

<math>  \Psi^{n+1} = \boldsymbol{\chi} - \Psi^n </math>

Uma matriz tridiagonal possui a seguinte forma 

<math> \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}
</math>

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 


<math> 
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}

</math>

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

Consequentemente, o sistema linear <math> \mathbf{Ax} = \mathbf{b}</math> 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 <math> \alpha_i </math>. Neste passo, as equações recursivas para os novos elementos da diagonal principal <math>\beta_i </math> e do vetor <math>\mathbf{b} </math> são 

<math> \beta'_i = \beta_i - \frac{\alpha_{i-1}}{\beta'_{i-1}}\gamma_{i-1}, \;\;\;\;\; i = 1, ..., N-1  </math>

e

<math>  b'_i = b_i - \frac{\alpha_{i-1}}{\beta'_{i-1}}b'_{i-1}, \;\;\;\;\; i = 1, ..., N-1                            </math>

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

<math> x_i = \frac{b'_i - \gamma_i x_{i+1}}{\beta'_i}, \;\;\;\;\; i = N-2, ..., 0                             </math>


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

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 

<math>  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_i </math>

e

<math> 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                                     </math>


e 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 <math>N=301</math> sítios, para um período total de 50 “segundos” e um passo de tempo de <math> \Delta t = 0,001 </math> “segundos”, onde as unidades de tempo possuem aspas pois são adimensionais (<math> \tau = \omega t/\hbar </math>). 

Código em Fortran para a integração numérica:

[https://github.com/Enriquecalderoli/M-todos-Computacionais-C/blob/master/nlse.f Equações de Schrödinger não-lineares acopladas]

== Aplicações==

Abaixo é possível visualizar a evolução dinâmica de uma das funções de onda (o comportamento da outra é simétrico) nos sítios centrais da rede unidimensional em questão para dois conjuntos de parâmetros físicos. No primeiro caso, tanto o campo elétrico aplicado F como o segundo potencial <math>U_2</math>, ao passo que o primeiro potencial vale <math>U_1=0.1</math>. No segundo caso o campo elétrico vale <math>F=4</math> e ambos os potenciais são iguais a <math>U_1=U_2=3</math>. O eixo vertical representa a intensidade da função de onda, que tem condição inicial unitária no sítio central. Um dos eixos horizontais representa os sítios ao redor do sítio central e o outro eixo representa a evolução temporal, de 0 a 30. Nota-se que no primeiro caso a função de onda apresenta grande dispersão, se difundindo pela rede, ao passo que no segundo caso ela permanece muito mais localizada ao redor do sítio inicial (Observação: o gráfico apresentado na segunda imagem foi girado em relação à primeira imagem pois a angulação adotada permite uma melhor visualização das propriedades de localização da função de onda nos sítios centrais). 

[[Arquivo:livre2novo.jpg|frame|200x200px|center|Gráfico 1: Evolução da função de onda na rede para <math>F=0</math>, <math>U_1 = 0.1</math> e <math>U_2=0</math>. ]]

[[Arquivo:433bnovo.jpg|frame|200x200px|center|Gráfico 2: Evolução da função de onda na rede para <math>F=4</math>, <math>U_1 = 3</math> e <math>U_2=3</math>. ]]

As oscilações de Bloch são oscilações características de uma partícula em um potencial confinante periódico e sujeita a uma força constante. Partindo da equação de movimento unidimensional para um elétron, por exemplo, submetido a um campo elétrico constante F, e usando as relações de física de estado sólido, é possível mostrar que a posição do elétron é uma função oscilatória. O período destas oscilações é dado por 

<math> T_B = \frac{2\pi\hbar}{dF} </math>

onde d representa a distância entre os poços de potencial confinante periódicos.

Uma das variáveis que podem ser estudadas com a implementação numérica mostrada acima é a função de participação de Wegner. Ela é definida por

<math> W(t) = \left[\sum_{i=1}^N \vert \psi_i(t)\vert^4 \right]^{-1} </math>

e seu valor cresce com o espalhamento espacial da probabilidade de ocupação de um sítio. Nas duas imagens abaixo, a função de Wegner é plotada como função do tempo para os casos de campo <math>F=1</math>, <math>U_1=3</math> e <math>U_2=3</math> e de campo <math>F=4</math>, <math>U_1=3</math> e <math>U_2=3</math>. É possível notar que o aumento da intensidade do campo aplicado resulta em uma frequência maior das oscilações de Bloch, em concordância com a previsão teórica. 

[[Arquivo:wegner133.png|frame|200x200px|center|Gráfico 3: Função de Wegner como função do tempo para <math>F=1</math>, <math>U_1 = 3</math> e <math>U_2=3</math>. ]]

[[Arquivo:wegner433.png|frame|200x200px|center|Gráfico 4: Função de Wegner como função do tempo para <math>F=4</math>, <math>U_1 = 3</math> e <math>U_2=3</math>. ]]

O desvio quadrático médio, mais conhecido na literatura como mean square displacement (msd), é dado por

<math> M(t) = \sum_{i=1}^{N} \vert \psi_i (t)\vert^2 (i-i_0)^2 </math>

e representa mais uma medida das propriedades espaciais das funções de onda. Nessas figuras aqui o msd é plotado como função do tempo para dois casos: campo <math>F = 4</math>, <math>U_1=3</math> e <math>U_2=0</math> e campo <math>F= 4</math>, <math>U_1=3</math> e <math>U_2=3</math>. No primeiro caso notamos dois modos de oscilação, o primeiro sendo uma oscilação de alta frequência, que é modulada por uma oscilação de menor frequência, com uma média de 10 ou 15 oscilações de alta frequência para cada uma de baixa frequência. Quando o potencial <math>U_2</math> passa de 0 para 3, nota-se que acompanhando cada oscilação de grande amplitude existem duas de amplitude menor, com duas de amplitude menor para cada oscilação de amplitude maior.

[[Arquivo:msd430.png|frame|200x200px|center|Gráfico 5: Mean Square Displacement como função do tempo para <math>F=4</math>, <math>U_1 = 3</math> e <math>U_2=0</math>. ]]

[[Arquivo:msd433.png|frame|200x200px|center|Gráfico 6: Mean Square Displacement como função do tempo para <math>F=4</math>, <math>U_1 = 3</math> e <math>U_2=3</math>. ]]

Outro parâmetro de interesse na integração destas equações acopladas é a entropia de Shannon, definida por

<math> S(t) = -\sum_{i=1}^{N} \vert \psi_i (t) \vert^2 ln\vert\psi_i(t)\vert^2 </math>

Grosso modo, ela pode ser interpretada como uma medida do grau de concentração espacial da probabilidade de encontrar a partícula na rede. Aqui são apresentadas duas imagens referentes a casos distintos: o primeiro é dado por campo <math>F = 1</math>, <math>U_1 = 3</math> e <math>U_2 = 0</math>, onde é possível notar que as oscilações de Bloch têm uma amplitude muito grande, e o segundo caso é referente a um campo <math>F = 1</math>, <math>U_1 = 3</math> e <math>U_2 = 3</math>. Para este último conjunto de valores dos parâmetros nota-se que as oscilações de Bloch possuem amplitude muito menor do que a do caso anterior, de forma que a entropia de Shannon apresenta alto grau de localização.

[[Arquivo:Shannon130.png|frame|200x200px|center|Gráfico 7: Entropia de Shannon como função do tempo para <math>F=1</math>, <math>U_1 = 3</math> e <math>U_2=0</math>. ]]

[[Arquivo:Shannon133.png|frame|200x200px|center|Gráfico 8: Entropia de Shannon como função do tempo para <math>F=1</math>, <math>U_1 = 3</math> e <math>U_2=3</math>. ]]

==Bibliografia==

[1] JUNGES, L. Dissertação de mestrado. IF-UFRGS, 2009. “A equação de Schrödinger não linear discreta com desordem de Aubry-André e com campo elétrico DC.”

[2] GARCIA, A. L. (2017) Numerical Methods for Physics

[3] COHEN-TANNOUDJI C., DIU B., LALOE F. Quantum mechanics. Volume 1. Wiley, 1991

[4] NEWMAN, M. (2013) Computational Physics