Equação de Dirac

Grupo: André Luis Della Valentina, Lucas dos Santos Assmann, Vinícius Bayne Müller

Introdução

A equação de Dirac descreve uma partícula relativística de spin ${\frac {1}{2}}$, como o elétron, com estrutura análoga a da equação de Schrödinger. Ela surgiu inicialmente como tentativa de explicar discrepâncias entre experimentos e teoria para o espectro do átomo de hidrogênio, possibilitando correções para o cálculo da energia do elétron em diferentes níveis (as chamadas correções de estrutura fina). Além disso, por meio dela foi possível prever a existência de antimatéria: descrevendo o elétron, ela também descreve o pósitron.

A fim de compatibilizar a Mecânica Quântica com a Relatividade Especial, a equação diferencial parcial é de primeira ordem tanto no tempo quanto no espaço (diferentemente da equação de Schrödinger, que é de segunda ordem no espaço). A equação de Dirac pode ser escrita de diversas formas; aqui, apresentamo-la explicitamente como um sistema de EDPs acopladas, mais conveniente para os propósitos do trabalho.

Assim como a equação de Schrödinger, a construção da equação de Dirac vem a partir do operador Hamiltoniano, que descreve a evolução temporal do estado quântico em questão:

$i\hbar {\frac {\partial }{\partial t}}\Psi ({\boldsymbol {x}},t)=H\Psi ({\boldsymbol {x}},t)$

onde, como anteriormente, os autovalores de $H$ correspondem aos valores possíveis de energia que o sistema pode assumir.

A mudança com relação à Mecânica Quântica não relativística acontece quando consideramos a energia relativística da partícula:

$E^{2}=p^{2}c^{2}+m^{2}c^{4}$

Assim, o Hamiltoniano é modificado para representar a medida da energia relativística total.

Diferentemente da equação de Schrödinger, aqui $\Psi ({\boldsymbol {x}},t)$ não representa apenas uma função de onda, mas sim um conjunto de quatro delas. Usando a notação

$\Psi ={\begin{bmatrix}\Phi _{1}\\\Phi _{2}\\\Phi _{3}\\\Phi _{4}\end{bmatrix}}$,

as componentes de $\Psi$ representam as funções de onda associadas ao elétron e ao pósitron: $\Phi _{1}$ ($\Phi _{2}$) representa a função de onda do elétron com spin up (down), e $\Phi _{3}$ ($\Phi _{4}$) representa a função de onda do pósitron com spin up (down). O objeto $\Psi ({\boldsymbol {x}},t)$ é chamado de spinor.

Dedução da equação de Dirac em uma e duas dimensões

Consideraremos neste trabalho a equação de Dirac em uma e duas dimensões, $x$ e $y$. A escolha dessas coordenadas se dá pela conveniência do acoplamento das EDPs: nesse caso, as quatro equações acopladas passam a ser acopladas duas a duas, facilitando o estudo do sistema.

Construção do Hamiltoniano completo

Consideremos uma partícula sob ação de um potencial $V({\boldsymbol {x}};t)$ (onde ${\boldsymbol {x}}=(x,y,z)^{T}$), que afeta a energia potencial da partícula, e de um potencial "escalar" $V_{sc}({\boldsymbol {x}};t)$, que afeta a massa de repouso da mesma. Seguindo uma das propostas possíveis para o Hamiltoniano, temos

$H=c{\boldsymbol {\alpha }}\cdot {\boldsymbol {p}}+\beta (mc^{2}+V_{sc})+VI_{4}$

onde ${\boldsymbol {\alpha }}=\alpha _{x}{\hat {i}}+\alpha _{y}{\hat {j}}+\alpha _{z}{\hat {k}}$; $\alpha _{i}$ e $\beta$ são matrizes 4x4 adimensionais e ${\boldsymbol {p}}$ é o vetor momento linear da partícula.

Pode-se mostrar que ${\boldsymbol {\alpha }}$ e $\beta$ devem satisfazer certos vínculos, limitando as escolhas possíveis para essas matrizes. A escolha mais simples e usualmente adotada consiste em tomar

$\beta ={\begin{pmatrix}I_{2}&0\\0&-I_{2}\end{pmatrix}}={\begin{pmatrix}1&0&0&0\\0&1&0&0\\0&0&-1&0\\0&0&0&-1\end{pmatrix}}$

$\alpha _{x}={\begin{pmatrix}0&\sigma _{x}\\\sigma _{x}&0\end{pmatrix}}={\begin{pmatrix}0&0&0&1\\0&0&1&0\\0&1&0&0\\1&0&0&0\end{pmatrix}}$

$\alpha _{y}={\begin{pmatrix}0&\sigma _{y}\\\sigma _{y}&0\end{pmatrix}}={\begin{pmatrix}0&0&0&-i\\0&0&i&0\\0&-i&0&0\\i&0&0&0\end{pmatrix}}$

$\alpha _{z}={\begin{pmatrix}0&\sigma _{z}\\\sigma _{z}&0\end{pmatrix}}={\begin{pmatrix}0&0&1&0\\0&0&0&-1\\1&0&0&0\\0&-1&0&0\end{pmatrix}}$

Sendo ${\boldsymbol {p}}=-i\hbar \nabla$, podemos escrever o produto escalar ${\boldsymbol {\alpha }}\cdot {\boldsymbol {p}}$ como

${\boldsymbol {\alpha }}\cdot {\boldsymbol {p}}=-i\hbar \left(\alpha _{x}{\frac {\partial }{\partial x}}+\alpha _{y}{\frac {\partial }{\partial y}}+\alpha _{z}{\frac {\partial }{\partial z}}\right)$

Portanto, em notação matricial o Hamiltoniano $H$ pode ser escrito como

$H=-i\hbar c{\begin{pmatrix}0&0&{\frac {\partial }{\partial z}}&{\frac {\partial }{\partial x}}-i{\frac {\partial }{\partial y}}\\0&0&{\frac {\partial }{\partial x}}+i{\frac {\partial }{\partial y}}&-{\frac {\partial }{\partial z}}\\{\frac {\partial }{\partial z}}&{\frac {\partial }{\partial x}}-i{\frac {\partial }{\partial y}}&0&0\\{\frac {\partial }{\partial x}}+i{\frac {\partial }{\partial y}}&-{\frac {\partial }{\partial z}}&0&0\\\end{pmatrix}}+{\begin{pmatrix}V+mc^{2}+V_{sc}&0&0&0\\0&V+mc^{2}+V_{sc}&0&0\\0&0&V-mc^{2}-V_{sc}&0\\0&0&0&V-mc^{2}-V_{sc}\\\end{pmatrix}}$ $H={\begin{pmatrix}V+mc^{2}+V_{sc}&0&-i\hbar c{\frac {\partial }{\partial z}}&-i\hbar c{\frac {\partial }{\partial x}}-\hbar c{\frac {\partial }{\partial y}}\\0&V+mc^{2}+V_{sc}&-i\hbar c{\frac {\partial }{\partial x}}+\hbar c{\frac {\partial }{\partial y}}&i\hbar c{\frac {\partial }{\partial z}}\\-i\hbar c{\frac {\partial }{\partial z}}&-i\hbar c{\frac {\partial }{\partial x}}-\hbar c{\frac {\partial }{\partial y}}&V-mc^{2}-V_{sc}&0\\-i\hbar c{\frac {\partial }{\partial x}}+\hbar c{\frac {\partial }{\partial y}}&i\hbar c{\frac {\partial }{\partial z}}&0&V-mc^{2}-V_{sc}\\\end{pmatrix}}$

Unidades naturais e redução para duas dimensões

A fim de simplificar o formalismo, adotamos as chamadas "unidades naturais", onde $\hbar =c=m=1$. Note que isso equivale a reescalar as quantidades físicas do problema por um fator adequado. Ao fazer $c=1$, também assumimos que a partícula está no limite relativístico.

Além disso, queremos estudar o problema em duas dimensões. Observamos que $\Psi (x,y,z)=\Psi (x,y)$; logo, ${\frac {\partial \Psi }{\partial z}}=0$. Portanto, temos o Hamiltoniano simplificado

$H={\begin{pmatrix}V+1+V_{sc}&0&0&-i{\frac {\partial }{\partial x}}-{\frac {\partial }{\partial y}}\\0&V+1+V_{sc}&-i{\frac {\partial }{\partial x}}+{\frac {\partial }{\partial y}}&0\\0&-i{\frac {\partial }{\partial x}}-{\frac {\partial }{\partial y}}&V-1-V_{sc}&0\\-i{\frac {\partial }{\partial x}}+{\frac {\partial }{\partial y}}&0&0&V-1-V_{sc}\\\end{pmatrix}}$

Forma explícita final

Retornando ao problema original, queremos resolver

$i\hbar {\frac {\partial }{\partial t}}\Psi =H\Psi \to \left[iI_{4}{\frac {\partial }{\partial t}}-H\right]\Psi =0$

Novamente utilizando a notação matricial, obtemos

${\begin{pmatrix}i{\frac {\partial }{\partial t}}-V-V_{sc}-1&0&0&i{\frac {\partial }{\partial x}}+{\frac {\partial }{\partial y}}\\0&i{\frac {\partial }{\partial t}}-V-V_{sc}-1&i{\frac {\partial }{\partial x}}-{\frac {\partial }{\partial y}}&0\\0&i{\frac {\partial }{\partial x}}+{\frac {\partial }{\partial y}}&i{\frac {\partial }{\partial t}}-V+V_{sc}+1&0\\i{\frac {\partial }{\partial x}}-{\frac {\partial }{\partial y}}&0&0&i{\frac {\partial }{\partial t}}-V+V_{sc}+1\end{pmatrix}}{\begin{pmatrix}\Phi _{1}\\\Phi _{2}\\\Phi _{3}\\\Phi _{4}\end{pmatrix}}=0$

Realizando a multiplicação matricial, pode-se ver que se obtém dois sistemas acoplados: $\Phi _{1}$ com $\Phi _{4}$ e $\Phi _{2}$ com $\Phi _{3}$. Escolhendo o sistema de $\Phi _{1}$ com $\Phi _{4}$:

${\begin{cases}\left(i{\dfrac {\partial }{\partial t}}-V-V_{sc}-1\right)\Phi _{1}+\left(i{\dfrac {\partial }{\partial x}}+{\dfrac {\partial }{\partial y}}\right)\Phi _{4}=0\\\left(i{\dfrac {\partial }{\partial x}}-{\dfrac {\partial }{\partial y}}\right)\Phi _{1}+\left(i{\dfrac {\partial }{\partial t}}-V+V_{sc}+1\right)\Phi _{4}=0\end{cases}}$

Simplificando e isolando a derivada temporal, obtemos por fim

${\begin{cases}{\dfrac {\partial \Phi _{1}}{\partial t}}=-i(V+V_{sc}+1)\Phi _{1}-{\dfrac {\partial \Phi _{4}}{\partial x}}+i{\dfrac {\partial \Phi _{4}}{\partial y}}\\{\dfrac {\partial \Phi _{4}}{\partial t}}=-i(V-V_{sc}-1)\Phi _{4}-{\dfrac {\partial \Phi _{1}}{\partial x}}-i{\dfrac {\partial \Phi _{1}}{\partial y}}\end{cases}}$

Por fim, a equação em uma dimensão ($x$) é facilmente obtida: basta fazer ${\dfrac {\partial \Phi _{1}}{\partial y}}={\dfrac {\partial \Phi _{4}}{\partial y}}=0$

${\begin{cases}{\dfrac {\partial \Phi _{1}}{\partial t}}=-i(V+V_{sc}+1)\Phi _{1}-{\dfrac {\partial \Phi _{4}}{\partial x}}\\{\dfrac {\partial \Phi _{4}}{\partial t}}=-i(V-V_{sc}-1)\Phi _{4}-{\dfrac {\partial \Phi _{1}}{\partial x}}\end{cases}}$

Neste trabalho, serão considerados principalmente problemas apenas com o potencial $V$; assim, caso indicado o contrário, será assumido $V_{sc}=0$.

Discretização

A equação de Dirac 1D pode ser escrita, na forma matricial, como

$i\partial _{t}\mathbf {\Phi } =[-i\sigma _{1}\partial _{x}+\sigma _{3}]\mathbf {\Phi } +[V(t,x)I_{2}]\mathbf {\Phi }$

onde $\mathbf {\Phi } =(\phi _{1},\phi _{4})^{T}$ e $I_{2}$ é matriz identidade de dimensão 2. As matrizes de Pauli ${\vec {\sigma }}$ são escritas, aqui, como $\sigma _{1}=\sigma _{x}$ e $\sigma _{3}=\sigma _{z}$.

Para discretizar uma equação diferencial parcial como essa, é necessário discretizar o espaço e o tempo. Convenciona-se $\Delta t$ como um passo finito de tempo e $h$ como um passo finito no espaço, de tal forma que $x_{j}=x_{0}+jh,t_{n}=t_{0}+n\Delta t$, onde $j,n$ são números inteiros. Define-se a notação $\mathbf {\Phi } (t_{n},x_{j})=\mathbf {\Phi } _{j}^{n}$ e também $V(t_{n},x_{n})=V_{j}^{n}$. Discretiza-se as derivadas parciais explicitamente com uma expansão em série de taylor da própria função:

$\mathbf {\Phi _{j}^{n+1}} =\mathbf {\Phi _{j}^{n}} +\partial _{t}\mathbf {\Phi _{j}^{n}} \Delta t+{\mathcal {O}}(\Delta t^{2})$

Considerando uma derivada discretizada $\delta _{t}\approx \partial _{t}$ e truncando na primeira ordem:

$\delta _{t}\mathbf {\Phi _{j}^{n}} ={\frac {\mathbf {\Phi _{j}^{n+1}} -\mathbf {\Phi _{j}^{n}} }{\Delta t}}$

O processo é completamente análogo para a derivada espacial, porém para facilitar a aplicação do método mantém-se o espaço centrado em $x_{j}$, em outras palavras faz-se uma expansão em torno de $x_{j-1}$, obtendo:

$\delta _{x}\mathbf {\Phi _{j}^{n}} ={\frac {\mathbf {\Phi _{j+1}^{n}} -\mathbf {\Phi _{j-1}^{n}} }{2h}}$

Com isso, obtém-se uma equação para um método explícito no tempo da equação de Dirac 1D.

$i\delta _{t}\mathbf {\Phi _{j}^{n}} =[-i\sigma _{1}\delta _{x}+\sigma _{3}]\mathbf {\Phi _{j}^{n}} +[V_{j}^{n}I_{2}]\mathbf {\Phi _{j}^{n}}$

$i{\frac {\mathbf {\Phi _{j}^{n+1}} -\mathbf {\Phi _{j}^{n}} }{\Delta t}}=[-i\sigma _{1}\delta _{x}+\sigma _{3}]\mathbf {\Phi _{j}^{n}} +[V_{j}^{n}I_{2}]\mathbf {\Phi _{j}^{n}}$

$i\mathbf {\Phi _{j}^{n+1}} =\mathbf {\Phi _{j}^{n}} +\Delta t[-i\sigma _{1}\delta _{x}+\sigma _{3}]\mathbf {\Phi _{j}^{n}} +\Delta t[V_{j}^{n}I_{2}]\mathbf {\Phi _{j}^{n}}$

Pode-se também desenvolver um método implícito no tempo fazendo a expansão de $\mathbf {\Phi _{j}^{n-1}}$ em torno de $t_{n}$, obtendo:

$\delta _{t}\mathbf {\Phi _{j}^{n}} ={\frac {\mathbf {\Phi _{j}^{n}} -\mathbf {\Phi _{j}^{n-1}} }{\Delta t}}$

Ao aplicar esta aproximação na equação discretizada basta dar um passo a frente em todos os elementos, obtendo um método implícito no tempo, já que há dependência com $t_{n+1}$.

$i\mathbf {\Phi _{j}^{n+1}} =\mathbf {\Phi _{j}^{n}} +\Delta t[-i\sigma _{1}\delta _{x}+\sigma _{3}]\mathbf {\Phi _{j}^{n+1}} +\Delta t[V_{j}^{n+1}I_{2}]\mathbf {\Phi _{j}^{n+1}}$

Neste trabalho, foram utilizados dois métodos de integração numérica diferentes: o de Crank-Nicolson, para a equação de Dirac em uma dimensão, e o de Leap-Frog, para a equação em duas dimensões. Deixamos a referência 3 para estes e outros métodos numéricos aplicados ao problema.

Método de Crank-Nicolson

O método de Crank-Nicolson (CN) consiste em uma média entre um método explícito e outro implícito no espaço. Utilizar-se-á a notação $\mathbf {\Phi _{j}^{n+1/2}}$ para representar justamente a média entre ambos os métodos, ou seja:

$\mathbf {\Phi _{j}^{n+1/2}} ={\frac {\mathbf {\Phi _{j}^{n+1}} +\mathbf {\Phi _{j}^{n}} }{2}}$

Define-se a notação:

$V_{j}^{n+1/2}\mathbf {\Phi _{j}^{n+1/2}} ={\frac {V_{j}^{n+1}\mathbf {\Phi _{j}^{n+1}} +V_{j}^{n}\mathbf {\Phi _{j}^{n}} }{2}}$

Dessa maneira, enuncia-se o método CN para a equação de Dirac 1D como:

$i\delta _{t}\mathbf {\Phi _{j}^{n}} =[-i\sigma _{1}\delta _{x}+\sigma _{3}]\mathbf {\Phi _{j}^{n+1/2}} +[V_{j}^{n+1/2}I_{2}]\mathbf {\Phi _{j}^{n+1/2}}$,

onde $\delta _{t},\delta _{x}$ são as discretizações explícitas das derivadas.

Para que seja possível aplicar e estudar o método, é necessário passar da notação matricial para escalar:

$i\delta _{t}\mathbf {\Phi _{j}^{n}} =-i\sigma _{1}\delta _{x}\left({\frac {\mathbf {\Phi _{j}^{n+1}} +\mathbf {\Phi _{j}^{n}} }{2}}\right)+\sigma _{3}\left({\frac {\mathbf {\Phi _{j}^{n+1}} +\mathbf {\Phi _{j}^{n}} }{2}}\right)+I_{2}\left({\frac {V_{j}^{n+1}\mathbf {\Phi _{j}^{n+1}} +V_{j}^{n}\mathbf {\Phi _{j}^{n}} }{2}}\right)$

$i{\frac {\mathbf {\Phi _{j}^{n+1}} -\mathbf {\Phi _{j}^{n}} }{\Delta t}}=-{\frac {i}{2}}\sigma _{1}\left[{\frac {\mathbf {\Phi _{j+1}^{n+1}} -\mathbf {\Phi _{j-1}^{n+1}} }{2h}}+{\frac {\mathbf {\Phi _{j+1}^{n}} -\mathbf {\Phi _{j-1}^{n}} }{2h}}\right]+\sigma _{3}\left({\frac {\mathbf {\Phi _{j}^{n+1}} +\mathbf {\Phi _{j}^{n}} }{2}}\right)+I_{2}\left({\frac {V_{j}^{n+1}\mathbf {\Phi _{j}^{n+1}} +V_{j}^{n}\mathbf {\Phi _{j}^{n}} }{2}}\right)$

Isolando cada tempo em um lado da igualdade:

$\left[I_{2}+{\frac {i}{2}}\Delta t\sigma _{3}+{\frac {i}{2}}\Delta tV_{j}^{n+1}I_{2}\right]\mathbf {\Phi _{j}^{n+1}} +{\frac {\Delta t}{4h}}\sigma _{1}\left[\mathbf {\Phi _{j+1}^{n+1}} -\mathbf {\Phi _{j-1}^{n+1}} \right]=$ $\left[I_{2}-{\frac {i}{2}}\Delta t\sigma _{3}-{\frac {i}{2}}\Delta tV_{j}^{n}I_{2}\right]\mathbf {\Phi _{j}^{n}} -{\frac {\Delta t}{4h}}\sigma _{1}\left[\mathbf {\Phi _{j+1}^{n}} -\mathbf {\Phi _{j-1}^{n}} \right]$

Abrindo as matrizes $\sigma _{1},\sigma _{3}$ e $I_{2}$ e operando-as sobre o vetor $\mathbf {\Phi }$ na equação, tem-se:

$\left({\begin{bmatrix}1&0\\0&1\\\end{bmatrix}}+{\frac {i\Delta t}{2}}{\begin{bmatrix}1&0\\0&-1\\\end{bmatrix}}+{\frac {i}{2}}\Delta tV_{j}^{n+1}{\begin{bmatrix}1&0\\0&1\\\end{bmatrix}}\right){\begin{bmatrix}\psi _{1,j}^{n+1}\\\psi _{4,j}^{n+1}\\\end{bmatrix}}+{\frac {\Delta t}{4h}}{\begin{bmatrix}0&1\\1&0\\\end{bmatrix}}{\begin{bmatrix}\psi _{1,j+1}^{n+1}-\psi _{1,j-1}^{n+1}\\\psi _{4,j+1}^{n+1}-\psi _{4,j-1}^{n+1}\\\end{bmatrix}}=\left({\begin{bmatrix}1&0\\0&1\\\end{bmatrix}}-{\frac {i\Delta t}{2}}{\begin{bmatrix}1&0\\0&-1\\\end{bmatrix}}-{\frac {i}{2}}\Delta tV_{j}^{n}{\begin{bmatrix}1&0\\0&1\\\end{bmatrix}}\right){\begin{bmatrix}\psi _{1,j}^{n}\\\psi _{4,j}^{n}\\\end{bmatrix}}-{\frac {\Delta t}{4h}}{\begin{bmatrix}0&1\\1&0\\\end{bmatrix}}{\begin{bmatrix}\psi _{1,j+1}^{n}-\psi _{1,j-1}^{n}\\\psi _{4,j+1}^{n}-\psi _{4,j-1}^{n}\\\end{bmatrix}}$

Pode-se realizar as operações matriciais e escrever duas equações escalares. Para facilitar a notação, utiliza-se $f_{j}^{n}=\psi _{1,j}^{n}$ e $g_{j}^{n}=\psi _{4,j}^{n}$:

${\begin{cases}\left[1+{\frac {i\Delta t}{2}}(V_{j}^{n+1}+1)\right]f_{j}^{n+1}+{\frac {\Delta t}{4h}}(g_{j+1}^{n+1}-g_{j-1}^{n+1})=\left[1-{\frac {i\Delta t}{2}}(V_{j}^{n}+1)\right]f_{j}^{n}-{\frac {\Delta t}{4h}}(g_{j+1}^{n}-g_{j-1}^{n})\\\left[1+{\frac {i\Delta t}{2}}(V_{j}^{n+1}-1)\right]g_{j}^{n+1}+{\frac {\Delta t}{4h}}(f_{j+1}^{n+1}-f_{j-1}^{n+1})=\left[1-{\frac {i\Delta t}{2}}(V_{j}^{n}-1)\right]g_{j}^{n}-{\frac {\Delta t}{4h}}(f_{j+1}^{n}-f_{j-1}^{n})\end{cases}}$

Tem-se então um número $n$ de equações onde $n$ é o número de elementos do espaço discretizado. Portanto o primeiro termo das duas equações gera uma matriz diagonal, pois multiplica os termos espaciais dependentes de $x_{j}$; já o segundo termo gera uma matriz tridiagonal com diagonal principal nula. Nota-se que os primeiros termos dos dois lados da igualdade são o conjugado um do outro: define-se, portanto, $\alpha ^{n}=1+{\frac {i\Delta t}{2}}(V_{j}^{n+1}+1)$ e $\beta =1+{\frac {i\Delta t}{2}}(V_{j}^{n+1}-1)$.

${\begin{cases}\alpha ^{n+1}f_{j}^{n+1}+{\dfrac {\Delta t}{4h}}(g_{j+1}^{n+1}-g_{j-1}^{n+1})=\alpha ^{n^{*}}f_{j}^{n}-{\dfrac {\Delta t}{4h}}(g_{j+1}^{n}-g_{j-1}^{n})\\\beta ^{n+1}g_{j}^{n+1}+{\dfrac {\Delta t}{4h}}(f_{j+1}^{n+1}-f_{j-1}^{n+1})=\beta ^{n^{*}}g_{j}^{n}-{\dfrac {\Delta t}{4h}}(f_{j+1}^{n}-f_{j-1}^{n})\end{cases}}$

Considerando que o potencial $V$ é só função da posição, escreve-se o método como:

${\begin{cases}Af^{n+1}+Bg^{n+1}=A^{*}f^{n}-Bg^{n}\\Cg^{n+1}+Bf^{n+1}=C^{*}g^{n}-Bf^{n}\\\end{cases}}$,

onde

$A={\begin{bmatrix}\alpha &0&0&\cdots &0\\0&\alpha &0&\cdots &0\\0&0&\alpha &\cdots &0\\\vdots &\vdots &\vdots &\ddots &\vdots \\0&\cdots &\cdots &\cdots &\alpha \\\end{bmatrix}};\quad B={\begin{bmatrix}0&{\frac {\Delta t}{4h}}&0&\cdots &0\\-{\frac {\Delta t}{4h}}&0&{\frac {\Delta t}{4h}}&\cdots &0\\0&-{\frac {\Delta t}{4h}}&0&\cdots &0\\\vdots &\vdots &\ddots &\ddots &\vdots \\0&\cdots &\cdots &-{\frac {\Delta t}{4h}}&0\\\end{bmatrix}};\quad C={\begin{bmatrix}\beta &0&0&\cdots &0\\0&\beta &0&\cdots &0\\0&0&\beta &\cdots &0\\\vdots &\vdots &\vdots &\ddots &\vdots \\0&\cdots &\cdots &\cdots &\beta \\\end{bmatrix}}$.

Por fim, pode-se escrever o método resolvendo o sistema

${\begin{cases}f^{n+1}=F^{-1}Df^{n}-F^{-1}Eg^{n}\\g^{n+1}=J^{-1}Gf^{n}-J^{-1}Hg^{n}\\\end{cases}}$,

onde $D=(B^{-1}A^{*}+C^{-1}B)$, $E=(I+C^{-1}C^{*})$, $F=(B^{-1}A-C^{-1}B)$, $G=(A^{-1}A^{*}+I)$, $H=(A^{-1}B+B^{-1}C^{*})$ e $J=(A^{-1}B-B^{-1}C)$.

Com isso, e com condições iniciais e de contorno bem estabelecidas, já é possível aplicar o método, dado que todas essas matrizes dependem somente dos parâmetros do sistema.

Estabilidade Crank-Nicolson

Utilizar-se-á o método de von Neumann para analisar a estabilidade do método de Crank-Nicolson para a equação de Dirac unidimensional. Para tanto, supõe-se que a função $\mathbf {\Phi _{j}^{n}}$ pode ser dada pela série de Fourier

$\mathbf {\Phi _{j}^{n}} =\sum _{k=0}^{\inf }\mathbf {A} ^{n}e^{ikqjh}$

Devido à independência linear de cada termo do somatório, ao substituir na equação do método haverá uma equação para cada ente do somatório. Se $\left|{\frac {\mathbf {A} ^{n+1}}{\mathbf {A} ^{n}}}\right|\leq 1$, então pode-se dizer que o método estável, já que dessa forma garante-se uma não divergência.

Aplica-se um termo geral da série de índice $k$ no método CN para a equação de Dirac 1D:

$\left[I_{2}+{\frac {i}{2}}\Delta t\sigma _{3}+{\frac {i}{2}}\Delta tV_{j}^{n+1}I_{2}\right]\mathbf {A} ^{n+1}e^{ikqjh}+{\frac {\Delta t}{4h}}\sigma _{1}\mathbf {A} ^{n+1}\left[e^{ikq(j+1)h}-e^{ikq(j-1)h}\right]=\left[I_{2}-{\frac {i}{2}}\Delta t\sigma _{3}-{\frac {i}{2}}\Delta tV_{j}^{n}I_{2}\right]\mathbf {A} ^{n}e^{ikqjh}-{\frac {\Delta t}{4h}}\sigma _{1}\mathbf {A} ^{n}\left[e^{ikq(j+1)h}-e^{ikq(j-1)h}\right]$

Divide-se tudo por $e^{ikqjh}$ e isola-se $\mathbf {A}$:

$\left[I_{2}+{\frac {i}{2}}\Delta t\sigma _{3}+{\frac {i}{2}}\Delta tV_{j}^{n+1}I_{2}+{\frac {\Delta t}{4h}}\sigma _{1}\left(e^{ikqh}-e^{-ikqh}\right)\right]\mathbf {A} ^{n+1}=\left[I_{2}-{\frac {i}{2}}\Delta t\sigma _{3}-{\frac {i}{2}}\Delta tV_{j}^{n}I_{2}-{\frac {\Delta t}{4h}}\sigma _{1}\left(e^{ikqh}-e^{-ikqh}\right)\right]\mathbf {A} ^{n}$

$\left[I_{2}+{\frac {i}{2}}\Delta t\sigma _{3}+{\frac {i}{2}}\Delta tV_{j}^{n+1}I_{2}+{\frac {\Delta t}{4h}}\sigma _{1}\left({\frac {\sin(kqh)}{2i}}\right)\right]\mathbf {A} ^{n+1}=\left[I_{2}-{\frac {i}{2}}\Delta t\sigma _{3}-{\frac {i}{2}}\Delta tV_{j}^{n}I_{2}-{\frac {\Delta t}{4h}}\sigma _{1}\left({\frac {\sin(kqh)}{2i}}\right)\right]\mathbf {A} ^{n}$

$\left[I_{2}+{\frac {i}{2}}\Delta t\sigma _{3}+{\frac {i}{2}}\Delta tV_{j}^{n+1}I_{2}-{\frac {i\Delta t}{4h}}\sigma _{1}\left({\frac {\sin(kqh)}{2}}\right)\right]\mathbf {A} ^{n+1}=\left[I_{2}-{\frac {i}{2}}\Delta t\sigma _{3}-{\frac {i}{2}}\Delta tV_{j}^{n}I_{2}+{\frac {i\Delta t}{4h}}\sigma _{1}\left({\frac {\sin(kqh)}{2}}\right)\right]\mathbf {A} ^{n}$

Nota-se que os termos que multiplicam o fator $\mathbf {A}$ são o conjugado um do outro. Define-se $z_{l}$ como a componente l escalarda multiplicação da matriz $I_{2}+{\frac {i}{2}}\Delta t\sigma _{3}+{\frac {i}{2}}\Delta tV_{j}^{n+1}I_{2}-{\frac {i\Delta t}{4h}}\sigma _{1}\left({\frac {\sin(kqh)}{2}}\right)$ pelo vetor $A$:

${\frac {{A_{l}}^{n+1}}{{A_{l}}^{n}}}={\frac {z^{*}}{z}}$

$\left|{\frac {{A_{l}}^{n+1}}{{A_{l}}^{n}}}\right|=\left|{\frac {z^{*}}{z}}\right|={\frac {|z^{*}|}{|z|}}=1$,

onde $z$ é sempre diferente de zero, visto que a parte real é dada por uma matriz identidade constante.

Mostra-se, portanto, que a razão entre os coeficientes da série de Fourier nunca diverge, ou seja, o método é incondicionalmente estável.

Simulações em Julia

Equação 1D

Fez-se simulações do método de Crank-Nicolson para a equação de Dirac unidimensional em três configurações de potenciais diferentes: Nulo (Partícula Livre), Poço Infinito e Oscilador Harmônico Simples. Para todos os casos utilizou-se uma condição inicial de um pacote gaussiano de desvio padrão $\sigma$ para uma das componentes de $\mathbf {\Phi }$, sendo a outra componente nula. O módulo quadrado do pacote gaussiano deve ter área unitária dentro da malha utilizada (de $x=0$ até $x=50$), por isso a constante de normalização deve ser ${\frac {1}{\sqrt {\sigma {\sqrt {\pi }}}}}$.

Segue trecho do código comum a todos as simulações realizadas; a única diferença é que em um caso do OHS o pacote gaussiano é deslocado do meio da malha para a posição $x=18$.

using Plots
using LaTeXStrings
using LinearAlgebra
using Integrals

function init(x, sigma)
	arg1 = @. ((-(x-25)^2)/(2*sigma^2))
	arg2 = 1/sqrt(sqrt(pi)*sigma)
	return arg2*exp.(arg1) ##Inicializo com um pacote gaussiano
end

function OHS(V,x, k)
	V = @. 0.5*k*(x-25)^2
	return V

end

function poco(V, h)
	
	#Quero colocar nas posições 18 e 32
	V[round(Int64, 18/h)] = 100000
	V[round(Int64, 32/h)] = 100000
	
	return V
end

##Matrizes do método de Crank-Nicholson

function matriz(dt,h,L, V)

	A = Diagonal(ones(L).*(1 .+ im*dt*0.5.*(V.+1)))  
	B = LinearAlgebra.Tridiagonal(ones(L-1).*(-dt/(4*h)), zeros(L), ones(L-1).*(dt/(4*h)))
	C = Diagonal(ones(L).*(1 .+ im*dt*0.5.*(V.-1)))
	
	A = Array(A)
	B = Array(B)
	C = Array(C)
	
	A_i = inv(A) 
	B_i = inv(B)
	C_i = inv(C)
	AiB = A_i*B
	CiB = C_i*B
	
	
	D = B_i*conj(A) + CiB
	E = I(L) + C_i*conj(C)
	F = B_i*A -CiB
	G = A_i*conj(A) + I(L)
	H = AiB + B_i*conj(C)
	J = AiB - B_i*conj(C)
	
	F_i = inv(F)
	J_i = inv(J)
	
	K = F_i*D
	L = F_i*E
	
	M = J_i*G
	N = J_i*H
	
	return K, L, M, N

end

function CN(size, h, tsize, V, dt, ci)

	f = zeros(Complex, tsize, size) #Cada linha representa vetor posição em um tempo diferente. ## phi 1
	g = zeros(Complex, tsize, size) ##phi 4
	
	A, B, C, D = matriz(dt, h, size, V) #Matrizes do método Crank-Nicholson
	
	f[1, :] = ci ##Condição Inicial
	g[1, :] = zeros(size)
	
	##Condições de Contorno
	f[:, 1] = zeros(tsize)
	f[:, end] = zeros(tsize)
	
	g[:, 1] = zeros(tsize)
	g[:, end] = zeros(tsize)
	
	i=2 #Não interfiro nos contornos
	while i<tsize
		
		f[i, :] = A*f[i-1, :] - B*g[i-1, :]
		g[i, :] = C*f[i-1, :] - D*g[i-1, :]
		i+=1
		
	end
	return f, g
end

Partícula Livre

Aplicou-se o código demonstrado com potencial nulo para simular o caso da partícula livre. Segue abaixo código utilizado para gerar a animação:

Nota-se que mesmo que $\phi _{4}$ seja nula no início, a existência da partícula $\phi _{1}$ (neste caso, o elétron com spin up) gera a outra (pósitron com spin down). O comportamento é de dispersão, ou seja, a tendência é que a posição fique cada vez menos definida: a partícula livre, por estar fora da ação de um potencial, apresenta momento bem definido; pelo Princípio da Incerteza, a posição torna-se incerta. Nota-se também que as bordas estão longe o suficiente na escala de tempo utilizada, de maneira que as condições de contorno não afetam a simulação.

function main()

	
	L = 50 #Dimensão espacial matriz
	k=0.1 #"Constante Elástica" do oscilador harmônico
	sigma = 1 #Largura do pacote gaussiano
	dt = 0.02
	size= 10*L
	h = L/size ##Passo espacial
	tmax = 15
	tsize = round(Int64, tmax/dt) #Tamanho do vetor tempo
	
	t = LinRange(0, tmax, tsize)
	x = LinRange(0, L, size)
	V = zeros(size)

	xi = init(LinRange(0, L, size), sigma) #Condição inicial pacote gaussiano


	f,g =CN(size,h,tsize, V, dt,xi)

	
	i=1
	anim = @animate while i<=tsize
		f_real = real(f[i, :])
		f_imag = imag(f[i, :])
		g_real = real(g[i, :])
		g_imag = imag(g[i, :])
		probf = abs2.(f[i, :])
		probg = abs2.(g[i, :])
		prob = probf + probg
		

		p=plot(x, probf, label=L"$|\phi_1|^2$", legendfont=font(15))
		plot!(x, probg, label=L"$|\phi_4|^2$", legendfont=font(15))
		titulo = "Partícula Livre"*", "*L"$\sigma=$"*string(sigma)*", "* L"$dt=$"*string(dt)* ", " *L"$t=$" * string(round(t[i], digits=3))
		plot!(title=titulo, 
		xlabel = "x", 
		xlim=(0,50),
		xticksfont = font(13),
		ticksfontsize = 10,
		ylim=(0,1),
		yticksfont=font(13),
		)

		i+=1
	end
	
	arquivo = "diracfree.gif"
	gif(anim, arquivo, fps=30)

end

main()

Poço Infinito

O poço infinito foi simulado colocando dois potenciais muito grandes em $x=18$ e $x=32$. Abaixo está reproduzida a animação.

Nota-se que a função simulada é a soma dos quadrados dos módulos de cada função de onda (elétron e pósitron); é a integral de $|\Phi |^{2}=|\phi _{1}|^{2}+|\phi _{4}|^{2}$ que deve sempre ser unitária, o que concorda com o obtido. Na animação é possível perceber um comportamento semelhante ao de uma onda estacionária, onde tem-se vários "harmônicos" associados ao pacote de ondas: como esperado, a função de onda é uma combinação linear dos autoestados do poço infinito.

Segue trecho do código utilizado para gerar a animação:

function main()

	
	L = 50 #Dimensão espacial matriz
	sigma = 1 #Largura do pacote gaussiano
	dt = 0.05 #Passo temporal
	size= 30*L #Utilizado para definir o passo espacial
	h = L/size #Passo espacial
	tmax = 30 #
	tsize = round(Int64, tmax/dt) tamanho do vetor tempo.
	
	t = LinRange(0, tmax, tsize)
	x = LinRange(0, L, size)
	V = zeros(size) 

	xi = init(LinRange(0, L, size), sigma) #Condição inicial de pacote gaussiano centralizado em x=25
	#V = OHS(V,x,k)
	V = poco(V, h) #Inicializa o potencial



	f,g =CN(size,h,tsize, V, dt,xi) #Obtém as funções de onda e sua respectiva evolução temporal

		

	
	
	##Produz-se a animação
	i=1
	anim = @animate while i<=tsize
		f_real = real(f[i, :])
		f_imag = imag(f[i, :])
		g_real = real(g[i, :])
		g_imag = imag(g[i, :])
		probf = abs2.(f[i, :])
		probg = abs2.(g[i, :])
		prob = probf + probg
		
		
		problem = SampledIntegralProblem(prob, x) #Resolve a integral em cada tempo, calculando a área e mostrando na animação
		method = TrapezoidalRule()
		integral = solve(problem, method)
		
		integral = round(integral[1]; digits=3) #Arredonda.
		integral = string(integral)
		
		p = plot(x, prob,label=L"$|\Phi|^2$", legendfont=font(15), fill=true, fillalpha=0.2)
		titulo = "Poço Infinito"*", "*", "*L"$\sigma=$"*string(sigma)*", "* L"$dt=$"*string(dt)* ", " *L"$t=$" * string(round(t[i], digits=3))*", Área="*integral
		plot!(x,V, label="Potencial", legendfont=font(15))
		plot!(title=titulo, 
		xlabel = "x", 
		xlim=(15,35),
		xticksfont = font(13),
		ticksfontsize = 10,
		ylim=(0,0.6),
		yticksfont=font(13),
		)

		i+=2 ##O passo é dado de 2 em 2 para deixar gif mais leve.
	end
	
	arquivo = "diracpoco.gif"
	gif(anim, arquivo, fps=30)

end

main()

Oscilador Harmônico Simples

Inicialmente, realizou-se a simulação deslocando o pacote gaussiano para a posição inicial $x=18$, com a finalidade de observar a oscilação. Testou-se a simulação com diferentes $k$'s; escolheu-se o que está mostrado abaixo pois facilita a visualização do comportamento oscilatório.

Nota-se um comportamento análogo ao caso clássico; porém, em alguns momentos é possível observar certos harmônicos do pacote de onda. Além disso, nota-se nas extremidades um valor máximo para $|\mathbf {\Phi } |^{2}$, também condizente com o caso clássico.

Durante testes nas simulações, notou-se que a área sob a curva possui variações quando se utiliza um passo $\Delta t$ muito alto. Nesse caso, a convergência do método de Crank-Nicolson para os valores teóricos depende da malha utilizada, embora a estabilidade dele seja incondicional. A diferença pode ser vista comparando-se a animação abaixo com a anterior:

Em seguida, passou-se para uma condição inicial com um pacote gaussiano centralizado:

Mesmo que a posição inicial esteja no centro, devido ao fato da posição não estar bem definida ocorrem pequenas oscilações em torno do ponto de equilíbrio.

Por fim, segue o código utilizado para o oscilador harmônico simples (a diferença entre cada caso é só as constantes e a condição inicial).

function main()

	
	L = 50 #Dimensão espacial matriz
	k=0.1 #"Constante Elástica" do oscilador harmônico
	sigma = 1 #Largura do pacote gaussiano
	dt = 0.02
	size= 10*L
	h = L/size ##Passo espacial
	tmax = 50
	tsize = round(Int64, tmax/dt) #Tamanho do vetor tempo
	
	t = LinRange(0, tmax, tsize)
	x = LinRange(0, L, size)
	V = zeros(size)

	xi = init(LinRange(0, L, size), sigma) #Condição inicial pacote gaussiano
	V = OHS(V,x,k) ##inicializa o potencial

	f,g =CN(size,h,tsize, V, dt,xi)

	i=1
	#Faz o gif
	anim = @animate while i<=tsize
	
		f_real = real(f[i, :])
		f_imag = imag(f[i, :])
		g_real = real(g[i, :])
		g_imag = imag(g[i, :])
		probf = abs2.(f[i, :])
		probg = abs2.(g[i, :])
		prob = probf + probg
		
		
		problem = SampledIntegralProblem(prob, x) ##Calcula a área sob a curva
		method = TrapezoidalRule()
		integral = solve(problem, method)
		integral = round(integral[1]; digits=3)
		integral = string(integral)
		
		p = plot(x, prob,label=L"$|\Phi|^2$", legendfont=font(15), fill=true, fillalpha=0.2)
		
		#=
		p=plot(x, probf, label=L"$|\phi_1|^2$", legendfont=font(15))
		plot!(x, probg, label=L"$|\phi_4|^2$", legendfont=font(15))
		=#
		
		titulo = "OHS"*", "*L"$k=$"*string(k)*", "*L"$\sigma=$"*string(sigma)*", "* L"$dt=$"*string(dt)* ", " *L"$t=$" * string(round(t[i], digits=3))*", Área="*integral
		plot!(x,V, label="Potencial", legendfont=font(15))
		plot!(title=titulo, 
		xlabel = "x", 
		xlim=(0,50),
		xticksfont = font(13),
		ticksfontsize = 10,
		ylim=(0,1),
		yticksfont=font(13),
		)
		i+=5 ##Feito de 5 em 5 para deixar .gif mais leve
	end
	
	arquivo = "diracOHS.gif"
	gif(anim, arquivo, fps=30)

end

main()

Equação 1D com potencial "escalar"

Por fim, consideraremos o problema com o potencial "escalar" $V_{sc}$. Dentro do formalismo da relatividade, os potenciais podem ser classificados de acordo com o seu comportamento frente a uma transformação de Poincaré: o potencial $V$ se transforma como um vetor (por isso chamado de potencial "vetor"), e o potencial $V_{sc}$, como um escalar (potencial "escalar"). Como exemplo de potencial do tipo vetor, temos os potenciais eletromagnéticos; e de potencial do tipo escalar, pode-se citar modelos de confinamento. Na prática, pode-se dizer que o potencial $V_{sc}$ altera a massa de repouso "efetiva" da partícula.

Por ser um tópico bastante especializado, será considerado aqui apenas o caso do poço infinito. O problema é apenas uma extensão do que foi exposto nessa seção; a discretização e a estabilidade do método de Crank-Nicolson são análogas, bastando apenas fazer $\alpha ^{n}=1+{\frac {i\Delta t}{2}}(V_{j}^{n}+V_{{sc}_{j}}^{n}+1)$ e $\beta ^{n}=1+{\frac {i\Delta t}{2}}(V_{j}^{n}-V_{{sc}_{j}}^{n}-1)$.

Poço infinito

Fez-se uma simulação completamente análoga à anterior, mas usando $V=0$ e definindo o poço infinito em $V_{sc}$. A animação está reproduzida abaixo.

Não há diferença notável com relação ao poço infinito em $V$: o confinamento da partícula não muda, mesmo que a sua massa de repouso "efetiva" sim. No entanto, ressalta-se que esse é o comportamento conjunto das duas componentes $\phi _{1}$ e $\phi _{4}$; o comportamento individual de cada uma pode ser diferente.

Em seguida, foi feita uma simulação colocando o poço infinito nos dois potenciais, ou seja, fazendo $V=V_{sc}$, conforme animação abaixo.

Dessa vez, o comportamento é bastante diferente: a partícula fica apenas parcialmente confinada no poço, sendo parte da função de onda transmitida. Assim, observa-se um fenômeno que não ocorre na formulação não relativística da Mecânica Quântica: o de tunelamento em uma barreira infinita.

Segue abaixo o código utilizado para produzir as animações, adaptado do código utilizado aqui. A diferença entre as duas animações é de apenas uma linha, definindo $V=0$ ou $V=V_{sc}$.

using Plots
using LaTeXStrings
using LinearAlgebra
using Integrals

function init(x, sigma)
    arg1 = @. ((-(x-25)^2)/(2*sigma^2))
    arg2 = 1/sqrt(sqrt(pi)*sigma)
    return arg2*exp.(arg1) ##Inicializo com um pacote gaussiano
end


function poco(V, h)
    
    #Quero colocar nas posições 18 e 32
    V[round(Int64, 18/h)] = 100000
    V[round(Int64, 32/h)] = 100000
    
    return V
end

##Matrizes do método de Crank-Nicholson

function matriz(dt,h,L, V, Sc)

    A = Diagonal(ones(L).*(1 .+ im*dt*0.5.*(V+Sc.+1)))  
    B = LinearAlgebra.Tridiagonal(ones(L-1).*(-dt/(4*h)), zeros(L), ones(L-1).*(dt/(4*h)))
    C = Diagonal(ones(L).*(1 .+ im*dt*0.5.*(V-Sc.-1)))
    
    A = Array(A)
    B = Array(B)
    C = Array(C)
    
    A_i = inv(A) 
    B_i = inv(B)
    C_i = inv(C)
    AiB = A_i*B
    CiB = C_i*B
    
    
    D = B_i*conj(A) + CiB
    E = I(L) + C_i*conj(C)
    F = B_i*A -CiB
    G = A_i*conj(A) + I(L)
    H = AiB + B_i*conj(C)
    J = AiB - B_i*conj(C)
    
    F_i = inv(F)
    J_i = inv(J)
    
    K = F_i*D
    L = F_i*E
    
    M = J_i*G
    N = J_i*H
    
    return K, L, M, N

end

function CN(size, h, tsize, V, Sc, dt, ci)

    f = zeros(Complex, tsize, size) #Cada linha representa vetor posição em um tempo diferente. ## phi 1
    g = zeros(Complex, tsize, size) ##phi 4
    
    A, B, C, D = matriz(dt, h, size, V, Sc) #Matrizes do método Crank-Nicholson
    
    f[1, :] = ci ##Condição Inicial
    g[1, :] = zeros(size)
    
    ##Condições de Contorno
    f[:, 1] = zeros(tsize)
    f[:, end] = zeros(tsize)
    
    g[:, 1] = zeros(tsize)
    g[:, end] = zeros(tsize)
    
    i=2 #Não interfiro nos contornos
    while i<tsize
        
        f[i, :] = A*f[i-1, :] - B*g[i-1, :]
        g[i, :] = C*f[i-1, :] - D*g[i-1, :]
        i+=1
        
    end
    return f, g
end


function main()

    
    L = 50 #Dimensão espacial matriz
    sigma = 1 #Largura do pacote gaussiano
    dt = 0.02
    size= 10*L
    h = L/size ##Passo espacial
    tmax = 50
    tsize = round(Int64, tmax/dt) #Tamanho do vetor tempo
    
    t = LinRange(0, tmax, tsize)
    x = LinRange(0, L, size)
    V = zeros(size)
    Sc = copy(V)

    xi = init(LinRange(0, L, size), sigma) #Condição inicial pacote gaussiano
    Sc = poco(Sc,h) ##potencial "escalar"
    V = copy(Sc) ##potencial "vetor"

    f,g = CN(size,h,tsize,V,Sc,dt,xi)

    i=1
    #Faz o gif
    anim = @animate while i<=tsize
    
        f_real = real(f[i, :])
        f_imag = imag(f[i, :])
        g_real = real(g[i, :])
        g_imag = imag(g[i, :])
        probf = abs2.(f[i, :])
        probg = abs2.(g[i, :])
        prob = probf + probg
        
        
        problem = SampledIntegralProblem(prob, x) ##Calcula a área sob a curva
        method = TrapezoidalRule()
        integral = solve(problem, method)
        integral = round(integral[1]; digits=3)
        integral = string(integral)
        
        p = plot(x, prob,label=L"$|\Phi|^2$", legendfont=font(15), fill=true, fillalpha=0.2)
        
        titulo = "Poço"*", "*L"$\sigma=$"*string(sigma)*", "* L"$dt=$"*string(dt)* ", " *L"$t=$" * string(round(t[i], digits=3))*", Área="*integral
        plot!(x,Sc, label=L"$V=V_{sc}$", legendfont=font(15))
        plot!(title=titulo, 
        xlabel = "x", 
        xlim=(0,50),
        xticksfont = font(13),
        ticksfontsize = 10,
        ylim=(0,1),
        yticksfont=font(13),
        )
        i+=5 ##Feito de 5 em 5 para deixar .gif mais leve
    end
    
    arquivo = "diracVscVPoco.gif"
    gif(anim, arquivo, fps=30)

end

main()

Referências

1. The quantum theory of the electron. Proceedings of the Royal Society of London A, v. 117, n. 778, p. 610–624, fev. 1928.
2. SAKURAI, J. J. Mecânica quântica moderna. 2. ed. Porto Alegre: Bookman, 2012.
3. BAO, W. et al. Numerical Methods and Comparison for the Dirac Equation in the Nonrelativistic Limit Regime. Journal of Scientific Computing, v. 71, n. 3, p. 1094–1134, jun. 2017.
4. SOFF, G. et al. Solution of the Dirac Equation for Scalar Potentials and its Implications in Atomic Physics. Zeitschrift für Naturforschung A, v. 28, n. 9, p. 1389–1396, 1 set. 1973.
5. THALLER, B. The Dirac equation. Berlin: Springer, 2010.