Quatro abordagens numéricas para resolver a equação de Schrödinger independente do tempo

Four numerical approaches to solve the time-independent Schrödinger equation

Resumo

A equação de Schrödinger foi um marco no desenvolvimento da Mecânica Quântica, e continua sendo até hoje imprescindível no estudo de sistemas não-relativísticos. Livros-texto de Mecânica Quântica geralmente abordam problemas em que é possível encontrar expressões analíticas para as autoenergias e autoestados do sistema. No entanto, existem inúmeros sistemas físicos de interesse para os quais não conhecemos soluções analíticas. Assim, métodos numéricos mostram-se essenciais para complementar as abordagens analíticas. Nesse trabalho, apresentamos métodos numéricos para a solução de problemas de Mecânica Quântica. Começamos discutindo aspectos da teoria, os quais usamos para motivar as diferentes abordagens numéricas. Em seguida, introduzimos a discretização e a adimensionalização dos problemas, visando empregá-las nas simulações numéricas. Apresentamos os métodos da tentativa, do encontro, matricial e variacional de Monte Carlo. Aplicamos os métodos em diferentes sistemas físicos e discutimos as vantagens e desvantagens de cada um deles. Esperamos que os alunos sejam capazes de reproduzir os resultados apresentados e estendê-los para as suas aplicações de interesse.

Palavras-chave:
Equação de Schrödinger independente do tempo; Método variacional de Monte Carlo; Método da tentativa; Método do encontro

Abstract

Schrödinger’s equation was a milestone in the development of Quantum Mechanics and continues to be essential in the study of non-relativistic systems to this day. Quantum Mechanics textbooks generally address problems in which it is possible to find analytical expressions for the eigenenergies and eigenstates of the system. However, we do not know analytical solutions for numerous physical systems of interest. Thus, numerical methods are essential to complement analytical approaches. In this work, we present numerical methods for solving Quantum Mechanics problems. We begin by discussing aspects of the theory we use to motivate the different numerical approaches. Next, we introduce discretization and dimensionless approaches, aiming to use them in numerical simulations. We present the shooting, matching, matrix and variational Monte Carlo methods. We apply the techniques to different physical systems and discuss the advantages and disadvantages of each one. We hope students can reproduce the results presented and extend them to their applications of interest.

Keywords
Time-independent Schrödinger’s equation; Variational Monte Carlo; Shooting method; Matching method

1. Introdução

A teoria quântica de Schrödinger é baseada na equação que leva o seu nome, e suas soluções fornecem toda informação que desejamos sobre o sistema. A priori, essa é uma equação diferencial parcial dependente do tempo. No entanto, podemos aplicar a técnica de separação de variáveis, obtendo uma parte que varia no tempo (importante para considerações dinâmicas) e outra estacionária, regida pela equação de Schrödinger independente do tempo [¹]. As suas soluções contém características da Mecânica Quântica (MQ): a quantização da energia, a superposição de estados quânticos, e a natureza probabilística da teoria.

Neste artigo iremos focar apenas na equação de Schrödinger independente do tempo por vários motivos. Primeiramente, para um conjunto muito grande de aplicações, estamos interessados apenas nas soluções estacionárias. Em segundo lugar, em posse das autofunções e autovalores, podemos obter a função de onda total dependente do tempo de maneira direta. Por último, os métodos numéricos são diferentes para abordar as duas versões da equação de Schrödinger. Isso vem do fato da equação dependente do tempo conter uma derivada temporal de primeira ordem além da derivada espacial de segunda ordem presente em ambas. Na nomenclatura de equações diferenciais parciais de segunda ordem, dizemos que a equação dependente do tempo é parabólica, enquanto a independente do tempo é elíptica [²]. O leitor interessado em um algoritmo relativamente simples para resolver a equação dependente do tempo pode consultar a Ref. [³].

Soluções analíticas da equação de Schrödinger independente do tempo são extremamente raras. De fato, mesmo para potenciais relativamente simples (como por exemplo o oscilador harmônico e Coulomb) as soluções não são trivialmente obtidas. Além disso, nos casos em que tais soluções existem, elas podem se tornar tão complexas que é difícil compreender toda a informação codificada. Devido a isso, os métodos numéricos desempenham um papel fundamental em MQ. Vários métodos numéricos estão disponíveis, porém esses possuem diferenças entre si, e muitas vezes não são aplicáveis na mesma situação. Um exemplo é o método de Euler, capaz de fornecer uma discretização para equações diferenciais ordinárias de primeira ordem, mas falha em problemas envolvendo movimento harmônico simples [⁴].

Um dos grandes papéis da Física Computacional é traduzir os problemas físicos de maneira que eles possam ser resolvidos com algoritmos computacionais [⁴]. Adotar uma abordagem computacional não significa que podemos nos descuidar com o embasamento teórico da MQ. Temos que que garantir que os cálculos feitos pela máquina correspondam ao sistema físico do qual queremos a solução. Portanto, a escolha do algoritmo adequado e a imposição das condições físicas à simulação são de vital importância para a obtenção da resposta correta.

Não há um método numérico universal para abordar os sistemas quânticos: as peculiaridades de cada sistema devem ser exploradas para determinarmos a melhor maneira de proceder. Com isso, nosso objetivo neste artigo é demonstrar a aplicação de diferentes abordagens numéricas em sistemas ilustrativos. Escolhemos alguns sistemas que possuem solução analítica, pois assim temos uma maneira direta de conferir e testar os resultados numéricos provenientes dos códigos desenvolvidos.

Nesse trabalho, discutimos os algoritmos de maneira independente da linguagem de programação. O nosso objetivo é que cada aluno seja capaz de desenvolver os próprios programas na linguagem mais familiar. Disponibilizamos os códigos desenvolvidos em Fortran90 em um repositório [⁵], porém encorajamos que os leitores desenvolvam suas próprias versões.

Esse artigo está organizado da seguinte maneira. Na seção 2, introduzimos alguns aspectos teóricos e mostramos como adimensionalizar e discretizar o sistema físico, afim de abordá-lo numericamente. Na seção 3, apresentamos o método da tentativa, útil para potenciais simétricos, e o aplicamos a dois sistemas com solução analítica, o poço infinito e o oscilador harmônico. Na seção 4, apresentamos o método do encontro e sua aplicação ao potencial de Lennard-Jonnes, o qual não possui solução analítica. Na seção 5, apresentamos o método matricial, o qual considera uma abordagem de problema de autovalor-autovetor, e o aplicamos ao átomo de Hidrogênio. Na seção 6, apresentamos o método variacional de Monte Carlo (VMC), ilustrando a sua aplicação no átomo de Hélio. Por fim, na seção 7, apresentamos algumas conclusões.

2. Preliminares

2.1. A equação de Schrödinger unidimensional

Estamos interessados na equação de Schrödinger independente do tempo para uma partícula em uma dimensão [⁶]

(1) $$H\varphi (x)=-\frac{{\hslash }^2}{2m}\frac{d^2\varphi (x)}{d{x}^2}+V(x)\varphi (x)=E\varphi (x),$$

onde H é o operador Hamiltoniano do sistema, m é a massa da partícula e ℏ a constante de Planck reduzida. A energia E é constante e corresponde ao autovalor da função de onda ϕ(x). Resolver a equação (1) implica em obter simultaneamente a energia E e a função de onda ϕ(x).

Uma propriedade muito importante da equação de Schrödinger é a linearidade: a soma de soluções da equação também é solução. Assim, podemos expressar a solução geral como a combinação linear

(2) $$\varphi (x)={\displaystyle \sum _n{c}_n}{\psi }_n(x),$$

onde c_n são coeficientes complexos. Ao escrever a equação (2) nessa forma, estamos supondo que o potencial V(x) é confinante. Assim, as energias E_n podem assumir apenas alguns valores discretos, ou seja, são quantizadas. Nesse caso, as funções ψ_n(x) são chamadas de autoestados, com as autoenergias associadas E_n. Os autoestados têm a propriedade de ortonormalidade [⁶],

(3) $${\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}d}x{\psi }_n^{*}(x){\psi }_m(x)={\delta }_{nm},$$

onde δ_nm é a delta de Kronecker.

Além disso, o conjunto de autofunções ψ_n(x) deve ser completo, ou seja, deve ser capaz de gerar o espaço de funções que contém ϕ(x). A completeza do conjunto é assumida na teoria [⁷], sendo uma hipótese lógica, uma vez que, se fosse falsa, existiriam estados físicos indescritíveis.

A interpretação probabilística de ϕ(x) nos diz que |ϕ(x)|²dx é a probabilidade de encontrarmos a partícula entre x e x+dx. Como a soma de todas as probabilidades deve ser 1, temos que

(4) $${\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}d}x|\varphi (x){|}^2=1.$$

Essa condição também é chamada de normalização da função de onda, e vamos utilizá-la ao nosso favor quando introduzirmos os métodos numéricos. Utilizando as equações (2) a (4), vemos que a normalização implica na seguinte condição para os coeficientes c_i,

(5) $${\displaystyle \sum _n|}{c}_n{|}^2=1.$$

Cada coeficiente complexo c_n pode ser obtido projetando a componente do n-ésimo autoestado ψ_n(x) presente em ϕ(x),

(6) $$\begin{array}{lll}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}d}x{\psi }_n{(x)}^{*}\varphi (x)\hfill & =\hfill & {\displaystyle \sum _m{c}_m}\left({\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}d}x{\psi }_n{(x)}^{*}{\psi }_m(x)\right)\hfill \\ \hfill & =\hfill & {\displaystyle \sum _m{c}_m}{\delta }_{nm}=c_n.\hfill \end{array}$$

Ao invés de lidar com a equação (1) e encontrar uma solução geral, nosso objetivo nesse trabalho é resolver a equação

(7) $$H{\psi }_n(x)=-\frac{{\hslash }^2}{2m}\frac{d^2{\psi }_n(x)}{d{x}^2}+V(x){\psi }_n(x)=E_n{\psi }_n(x)$$

para os vários autoestados ψ_n(x). O motivo é que, em posse das autofunções ψ_n(x), e seus autovalores associados E_n, podemos usar a equação (2) para expressar a função de onda total.

2.2. Propriedades dos autoestados

Mais do que apenas objetos matemáticos, as soluções da equação de Schrödinger devem representar sistemas físicos. Quantidades provenientes das autofunções que podem ser medidas em experimentos devem ser “bem comportadas”. Assim, para ser uma solução aceitável, ψ_n(x) e dψ_n(x) devem ser: (i) finitas, (ii) unicamente valoradas, e (iii) contínuas.

Sobre a condição (i), como |ψ_n(x)|² está associado com uma densidade de probabilidade, então ψ_n(x) deve ser finita para evitarmos situações não-físicas. Para dψ_n(x)/dx ser finita, então ψ_n(x) deve ser contínua [condição (iii)], como esperaríamos de um sistema físico.

A condição (ii) vai ser satisfeita sem necessidade de imposição para os sistemas nesse trabalho. Porém, problemas onde há uma simetria cilíndrica ou esférica requerem atenção para satisfazer esse requerimento.

A necessidade de dψ_n(x)/dx ser contínua [condição (iii)] pode ser vista reescrevendo a equação (7) como:

(8) $$\frac{d^2{\psi }_n(x)}{d{x}^2}=\frac{2m}{{\hslash }^2}(V(x)-E){\psi }_n(x).$$

Para V(x) e E finitos, então d²ψ_n(x)/dx² deve ser finito também, logo dψ_n(x)/dx deve ser contínua. O leitor familiarizado com o poço infinito sabe que a derivada primeira da função de onda é descontínua nas paredes do poço, consequência de V(x) ser divergente nesses pontos. Isso é uma idealização que facilita a matemática envolvida, mas devemos lembrar que um potencial impenetrável, não importando quão grande seja a energia, não é um sistema físico.

2.3. A simetria de paridade

Dentre a vasta gama de potenciais descrevendo sistemas físicos, alguns potenciais de interesse possuem uma propriedade muito interessante e útil: são funções pares,

(9) $$V(x)=V(-x).$$

Para ver a consequência dessa propriedade nas soluções, vamos supor que conhecemos uma solução ϕ(x) para a equação de Schrödinger independente do tempo, equação (1). Em seguida, fazemos a mudança de variáveis x → −x,

(10) $$-\frac{{\hslash }^2}{2m}\frac{d^2\varphi (-x)}{d{(-x)}^2}+V(-x)\varphi (-x)=E\varphi (-x).$$

Usando a equação (9) (e d²/d(−x)² → d²/dx²),

(11) $$-\frac{{\hslash }^2}{2m}\frac{d^2\varphi (-x)}{d{x}^2}+V(x)\varphi (-x)=E\varphi (-x),$$

ou seja, ϕ(−x) também é solução da equação de Schrödinger independente do tempo.

Lembrando da linearidade da equação (1), podemos construir duas combinações a partir das soluções ϕ(x) e ϕ(−x), as quais também serão soluções:

(12) $$\begin{array}{lll}{\varphi }_{+}(x)\hfill & =\hfill & \varphi (x)+\varphi (-x),\hfill \\ {\varphi }_{-}(x)\hfill & =\hfill & \varphi (x)-\varphi (-x).\hfill \end{array}$$

Essas duas funções têm paridade bem definida,

(13) $$\begin{array}{lll}{\varphi }_{+}(x)\hfill & =\hfill & +{\varphi }_{+}(-x),\hfill \\ {\varphi }_{-}(x)\hfill & =\hfill & -{\varphi }_{-}(-x),\hfill \end{array}$$

ou seja, ϕ+(x) é par e ϕ−(x) é ímpar. A conclusão é que se o potencial é par, então os autoestados terão paridade bem definida, i.e., serão pares ou ímpares.

Na Figura 1 exibimos um exemplo de uma função cosseno, que é par, e uma função seno, ímpar. Mostramos também o comportamento dessas funções próximo à origem, pois isso será essencial para o método da tentativa, seção 3.

Figura 1
Exemplo de uma função (a) ímpar (seno) e (b) par (cosseno) e seus comportamentos próximos à origem.

2.4. Penetração da região classicamente proibida

Sabemos que a Mecânica Quântica traz fenômenos novos que contradizem a intuição da Mecânica Clássica. Um deles é a penetração da região classicamente proibida, ou seja, o fato de partículas poderem ser encontradas em regiões onde sua energia total é menor que a potencial, E ＜ V (x), o que resultaria em um momento imaginário [⁸]. Para exemplificar esse fenômeno, consideremos o potencial degrau [¹],

(14) $$V(x)=\left\{\begin{array}{ll}0\hfill & x\le 0,\hfill \\ V_0\hfill & x>0,\hfill \end{array}\right.$$

com V₀ ＞ 0. A penetração da região proibida se expressa no fato de ψ(x) não se anular na região x ＞ 0, sendo portanto possível observar a partícula nessa parte do espaço. A Figura 2 exemplifica o comportamento da função de onda nessa região. Note que, embora a partícula adentre a região proibida, a amplitude de probabilidade de encontrá-la decresce exponencialmente até se anular [¹]. O quanto a partícula é capaz de entrar nessa região depende da relação entre a energia total E e o potencial V₀: quanto menor for a energia em relação a V₀, menor será a penetração, embora ela sempre ocorra. Esse fenômeno é intrínseco da Mecânica Quântica, e acontece em qualquer situação em que há uma região classicamente proibida. Logo, também ocorrerá no estudo de potenciais confinantes nas próximas seções.

Figura 2
Penetração da região classicamente proibida no potencial degrau para energias (a) E = 0, 67 V₀ e (b) E = 0, 97 V₀. Note que na situação (a) a penetração é menor que no caso (b), onde a energia da partícula se aproxima de V₀.

2.5. A equação de Schrödinger tridimensional para potenciais esfericamente simétricos

Até agora discutimos apenas uma partícula em uma dimensão. Essa abordagem pode parecer um tanto limitada, mas nessa seção mostramos que é possível transformar o problema de duas partículas em 3D interagindo via um potencial esfericamente simétrico em um problema similar ao de uma partícula em 1D. O procedimento exposto a seguir é similar ao apresentado em livros introdutórios de MQ ao resolver o átomo de hidrogênio. Alunos interessados em mais detalhes do que os fornecidos aqui podem consultar as referências [¹, ⁶].

Uma situação muito comum é quando temos duas partículas de massa m₁ e m₂ submetidas a um potencial V(x), onde r é a distância entre as partículas. Para esses casos, podemos descrever o sistema em termos de uma única partícula e r passa a ser a distância da origem [⁶, ⁹]. A massa μ da partícula é dada por

(15) $$\frac{1}{\mu }=\frac{1}{m_1}+\frac{1}{m_2},$$

chamada de massa reduzida do sistema [⁸].

Podemos escrever a equação de Schrödinger para a partícula de massa reduzida,

(16) $$\left[-\frac{{\hslash }^2}{2\mu }{\nabla }^2+V(r)\right]\psi (\mathit{r})=E\psi (\mathit{r}).$$

Note que a equação (16) é uma equação diferencial parcial em três dimensões, sendo mais complexa de ser resolvida tanto analiticamente quanto numericamente. Contudo, podemos utilizar a simetria esférica do problema e buscar soluções do tipo

(17) $${\psi }_{n\ell m}(\mathit{r})=R_{n\ell }(r)Y_{\ell }^m(\theta ,\varphi ),$$

onde (r, θ, ϕ) são as coordenadas esféricas usuais, n é chamado de número quântico principal e ℓ e m são os números quânticos associados ao momento angular total e sua projeção, conservados em potenciais esfericamente simétricos [⁶].

Substituindo a equação (17) na (16) e aplicando o método da separação de variáveis [¹], reconhecemos as funções $Y_{\ell }^m(\theta ,\varphi )$ como os harmônicos esféricos [⁶, ⁷]. Um comentário que merece ênfase é que os harmônicos esféricos são solução da equação angular de todos potenciais esfericamente simétricos, e não apenas do átomo de hidrogênio. O que é particular do átomo de hidrogênio é que as autoenergias não dependem do número quântico ℓ, particularidade do potencial de Coulomb proporcional a 1/r.

O procedimento de separação de variáveis também fornece uma equação para as funções radiais R_nℓ(r) [¹],

(18) $$\frac{1}{r^2}\frac{d}{dr}\left(r^2\frac{d{R}_{n\ell }(r)}{dr}\right)+\frac{2\mu }{{\hslash }^2}\left[E_{nl}-V(r)-\frac{\ell (\ell +1)}{r^2}\right]R_{n\ell }(r)=0.$$

Em seguida, fazemos a transformação de variáveis u_nℓ(r) = rR_nℓ(r), comumente chamada de função radial reduzida, obtendo a equação diferencial

(19) $$-\frac{{\hslash }^2}{2\mu }\frac{d^2{u}_{n\ell }(r)}{d{r}^2}+\left[\text{}V(r)+\frac{{\hslash }^2\ell (\ell +1)}{2\mu r^2}\text{}\right]u_{n\ell }(r)=E_{n\ell }{u}_{n\ell }(r).$$

Para fazer a conexão com o caso unidimensional, identificamos o termo entre colchetes como um potencial efetivo,

(20) $$V_{\text{eff}}(r)=V(r)+\frac{{\hslash }^2\ell (\ell +1)}{2\mu r^2}.$$

O segundo termo do lado direito da equação existe para qualquer ℓ ≠ 0, é repulsivo, depende do momento angular e aumenta rapidamente próximo da origem, logo recebe o nome de barreira centrífuga. Uma situação similar ocorre na Mecânica Clássica em potenciais centrais [⁸].

Portanto, mesmo se tratando de duas partículas, temos uma equação de Schrödinger unidimensional para a função u_nℓ(r), a qual pode ser resolvida com os métodos desenvolvidos neste artigo. Ressaltamos que após obter as soluções numéricas para u_nℓ(r), ainda é necessário dividí-las por r para obter as soluções radiais, dado que u_nℓ(r) = rR_nℓ(r).

2.6. Soluções reais ou complexas?

Ao encontrar soluções analíticas, empregar funções reais ou complexas é meramente uma questão de conveniência. Por exemplo, se vamos escolher utilizar cos(kx) e sin(kx) ou exp(±ikx) vai depender do nosso sistema de coordenadas, potencial e condições de contorno do problema. Para esse exemplo específico, sabemos que as duas bases são intercambiáveis se usarmos as relações:

(21) $$\begin{array}{lll}\cos (kx)\hfill & =\hfill & \frac{e^{ikx}+e^{-ikx}}{2},\hfill \\ \sin (kx)\hfill & =\hfill & \frac{e^{ikx}-e^{-ikx}}{2i}.\hfill \end{array}$$

A situação não é tão simples se considerarmos uma abordagem numérica. Elaborar um código com quantidades complexas, apesar de possível, tem dificuldades associadas. Não são todas linguagens que possuem suporte nativo a variáveis complexas, logo o usuário precisa implementar as funcionalidades não existentes. Mesmo se todas operações necessárias estiverem disponíveis, temos que tomar cuidado com operações que requerem o complexo conjugado da quantidade. Em muitos casos, é mais conveniente separar a parte real e imaginária da solução e escrever duas equações a serem resolvidas, cada uma apenas com quantidades reais.

Por esses motivos, sempre que possível, gostaríamos de trabalhar apenas com soluções reais no nosso código. A pergunta que aparece é se ao considerarmos apenas soluções reais estamos perdendo generalidade.

Suponha que encontramos uma solução ϕ da equação de Schrödinger independente do tempo, equação (1). Sabemos que ϕ* também é solução dessa equação. Pela propriedade de linearidade, combinações lineares de soluções também são soluções. Assim, vamos compor duas combinações:

(22) $$\begin{array}{lll}{\varphi }_1\hfill & =\hfill & \varphi +{\varphi }^{*},\hfill \\ {\varphi }_2\hfill & =\hfill & i(\varphi -{\varphi }^{*}).\hfill \end{array}$$

As soluções ϕ₁ e ϕ₂ são reais e satisfazem a equação (1). Logo, dada uma solução complexa da equação de Schrödinger independente do tempo, sempre podemos construir duas soluções reais (claro que se ϕ já era puramente real, ϕ₂ = 0). A conclusão é que felizmente podemos trabalhar apenas com variáveis reais no nosso programa e não estamos deixando de encontrar soluções da equação.

2.7. Adimensionalização

A equação de Schrödinger depende da constante reduzida de Planck ℏ ∼ 10⁻³⁴ Js e da massa da partícula (para um elétron, por exemplo, é da ordem de 10⁻³¹ kg). Esses valores são muito pequenos quando comparados com os que lidamos no dia-a-dia. Isso não interfere no trabalho analítico, onde os valores numéricos são substituídos apenas na última etapa dos cálculos. Todavia, nas simulações numéricas precisamos empregar esses valores desde o primeiro momento. Contudo, a representação de números reais no computador é uma aproximação, que depende de uma precisão fixada previamente pelo número de bits [⁹, ¹⁰]. Devido a isso, a representação de números reais em uma máquina tem uma precisão finita. Portanto, é uma tarefa difícil representar quantidades tão pequenas no computador. Mais ainda, realizar um número muito grande de operações com números pequenos faz com que os erros de truncamento se acumulem.

Logo, é útil abandonar as unidades, transformando todas as grandezas do problema em adimensionais, e recuperar as dimensões físicas após termos obtido a solução do problema. Para isso, escolhemos uma grandeza de interesse, como por exemplo um comprimento característico do sistema ℓ, e reescalamos todas as distâncias:

(23) $$\tilde{x}=x/\ell .$$

Nessa seção utilizaremos o símbolo ˜ para denotar quantidades adimensionais. Aplicando essa mudança de escala na equação de Schrödinger, equação (1), aparece a seguinte combinação de unidades:

(24) $$\frac{{\hslash }^2}{m}\frac{d^2}{d{x}^2}=\left(\frac{{\hslash }^2}{m{\ell }^2}\right)\frac{d^2}{d{\tilde{x}}^2}.$$

O termo em parênteses tem dimensão de energia, logo é um candidato para reescalar todas grandezas com unidade de energia,

(25) $$\begin{array}{lll}\tilde{E}\hfill & =\hfill & E/ϵ,\hfill \\ \tilde{V}\hfill & =\hfill & V/ϵ,\hfill \end{array}$$

onde ∈ = ℏ²/(mℓ²).

Aplicando as equações (23) e (25) em (1) temos:

(26) $$-\frac{1}{2}\frac{d^2\tilde{\varphi }(\tilde{x})}{d{\tilde{x}}^2}+\tilde{V}(\tilde{x})\tilde{\varphi }(\tilde{x})=\tilde{E}\tilde{\varphi }(\tilde{x}),$$

onde $\tilde{\varphi }={\ell }^{1/2}\varphi$. A equação (26) é a que queremos resolver numericamente, visto que ela não contém, a princípio, números muito grandes ou pequenos.

Comparando a equação (1) e (26) vemos que poderíamos ter adotado o conhecido “ℏ = m = 1” e obteríamos o mesmo resultado. A vantagem dessa abordagem passo-a-passo é que as equações (23) e (25) deixam claro o que deve ser feito para recuperar as dimensões ao final da simulação. Note que dessa maneira, podemos utilizar a mesma simulação para obter resultados em regimes diferentes, de escalas de energia e comprimento, apenas alterando ℓ e ϵ.

Para evitar uma notação muito carregada, abandonamos os símbolos de ˜ daqui para frente, deixando implícito que todas equações são adimensionais. A exceção é a seção 6, onde recuperamos as unidades físicas explicitamente para obter a energia do átomo de Hélio.

2.8. Discretização

Em cálculos analíticos estamos acostumados com um contínuo de posições. Embora existam métodos numéricos que usem essa base, como o método variacional de Monte Carlo que veremos na seção 6, muitas abordagens discretizam o espaço, seções 3 a 5.

Ao invés de encontrar uma solução para um contínuo de posições de x = 0 até x = L, estamos interessados na função de onda apenas em um conjunto discreto de pontos:

(27) $$x_i=i\text{Δ}x,\text{Δ}x=L/N,i=0,1,\dots ,N.$$

Logo, resolver o problema corresponde a encontrar os valores de ψ(x_i) ≡ ψ_i para os N + 1 pontos do intervalo. Computar o potencial nesses pontos é imediato, V (x_i) ≡ V_i, porém derivadas da função de onda requerem mais trabalho.

Para isso, escrevemos a série de Taylor da função de onda ao redor do ponto x_i+1,

(28) $${\psi }_{i+1}\approx {\psi }_i+{\psi }_i^{\prime }\text{Δ}x+\frac{1}{2}{\psi }_i^{\prime \prime }(\text{Δ}{x}^2),$$

e do ponto x_i−1,

(29) $${\psi }_{i-1}\approx {\psi }_i-{\psi }_i^{\prime }\text{Δ}x+\frac{1}{2}{\psi }_i^{\prime \prime }{(\text{Δ}x)}^2.$$

Subtraindo a equação (29) da (28) obtemos a discretização simétrica da primeira derivada [⁴, ¹¹],

(30) $${\psi }^{\prime }(x_i)\equiv {\psi }_i^{\prime }=\frac{{\psi }_{i+1}-{\psi }_{i-1}}{2\text{Δ}x}+𝒪[{(\text{Δ}x)}^2].$$

Essa equação supõe que sabemos o valor da função de onda no ponto à esquerda (ψ_i−1) e à direita (ψ_i+1). No entanto, existem situações em que temos apenas o valor para um dos vizinhos. Logo, podemos utilizar as equações (28) e (29), menos precisas que a (30), para obter:

(31) $${\psi }_i^{\prime }=\frac{{\psi }_{i+1}-{\psi }_i}{\text{Δ}x}+𝒪[(\text{Δ}x)],$$

(32) $${\psi }_i^{\prime }=\frac{{\psi }_i-{\psi }_{i-1}}{\text{Δ}x}+𝒪[(\text{Δ}x)].$$

Somando as equações (29) e (28) obtemos a discretização simétrica da segunda derivada [⁴],

(33) $${\psi }^{\prime \prime }(x_i)\equiv {\psi }_i^{\prime \prime }=\frac{{\psi }_{i+1}-2{\psi }_i+{\psi }_{i-1}}{{(\text{Δ}x)}^2}+𝒪[{(\text{Δ}x)}^2].$$

Substituindo a equação (33) em (26), obtemos uma equação de Schrödinger discretizada para o ponto x_i,

(34) $$-\frac{1}{2}\left[\frac{{\psi }_{i+1}-2{\psi }_i+{\psi }_{i-1}}{{(\text{Δ}x)}^2}\right]\simeq (E-V_i){\psi }_i.$$

Isolando ψ_i+1, temos

(35) $${\psi }_{i+1}\simeq 2{\psi }_i-{\psi }_{i-1}-2{(\text{Δ}x)}^2(E-V_i){\psi }_i,$$

onde o símbolo ≃ nos lembra que esse não é um resultado exato, mas sim uma aproximação sujeita à precisão da expressão que utilizamos para a derivada segunda. No restante do texto, apesar de não usarmos o símbolo por simplicidade, sabemos que temos uma situação similar se alguma discretização de derivada está envolvida. Se soubermos a função de onda em dois pontos consecutivos (ψ_i−1 e ψ_i), podemos utilizar a equação (35) para obter o ponto seguinte (ψ_i+1). Repetindo o procedimento, podemos obter todos os ψ_i do intervalo de interesse. A equação (35) não é a única opção para obtermos os valores de ψ_i dados dois valores consecutivos da função de onda; um procedimento análogo pode ser feito com o método de Numerov [¹²]. Equações do tipo da equação (35) são conhecidas como algoritmos integradores, pois se “integrarmos” uma equação diferencial obtemos a sua solução.

Temos duas dificuldades para a aplicação imediata dessa equação. Primeiramente, devemos saber o valor da função de onda em dois pontos consecutivos. Segundo, a equação de Schrödinger é uma equação de autovalor-autovetor, portanto a energia E também precisa ser determinada. Os métodos das seções 3 a 5 se baseiam na equação (35) e resolvem essas dificuldades apresentadas.

3. O Método da Tentativa

Iniciamos agora o estudo dos métodos numéricos para a obtenção dos autoestados de sistemas quânticos ligados. Nosso interesse será voltado exclusivamente para potenciais que confinam a partícula em uma região finita do espaço. Em problemas desse tipo precisamos obter ao mesmo tempo a energia E_n e sua autofunção ψ_n(x) associada. Como precisamos buscar essas duas incógnitas simultaneamente, muitos métodos propõem uma forma de teste, onde escolhemos uma energia e verificamos se a autofunção correspondente satisfaz os requisitos da seção 2.2.

A primeira abordagem que vamos apresentar neste trabalho recebe o nome de método da tentativa, e se aplica para potenciais pares. Como visto na seção 2.3, nessa situação as autofunções serão obrigatoriamente pares ou ímpares. Assim, vamos utilizar a equação (35) e buscar por funções que possuem paridade bem definida, pois sabemos que essa é uma consequência da teoria.

Primeiramente, precisamos do valor da função de onda em dois pontos consecutivos para utilizar a equação (35), e para isso vamos utilizar o fato das autofunções serem pares ou ímpares. Como vimos na seção 2.3, essas funções possuem um comportamento bem definido no ponto de simetria (no caso, x = 0) dado por

(36) $$\psi (0)=\left\{\begin{array}{ll}\text{constante}\hfill & \text{(funções pares)},\hfill \\ 0\hfill & \text{(funções ímpares)},\hfill \end{array}\right.$$

conforme ilustrado na Figura 1. Já a derivada primeira nas proximidades do ponto de simetria implica em

(37) $${\psi }^{\prime }(0)=\left\{\begin{array}{ll}0\hfill & \text{(funções pares)},\hfill \\ \text{constante}\hfill & \text{(funções ímpares)}.\hfill \end{array}\right.$$

O fato de termos duas constantes indeterminadas nas equações (36) e (37) pode parecer um problema a princípio. No entanto, podemos usar a propriedade de linearidade da equação de Schrödinger: se encontramos uma solução, então uma constante vezes essa solução também resolve a equação. Podemos tomar valores arbitrários para essas duas contantes indeterminadas e, após obter a solução, normalizamos a função de onda de acordo com a equação (4).

Assim, as equações (36) e (37) podem ser utilizadas para iniciar a autofunção em dois pontos da nossa simulação,

(38) $$\begin{array}{lll}{\psi }_0\hfill & =\hfill & 1\text{e}{\psi }_1=1\text{(funções pares)},\hfill \\ {\psi }_0\hfill & =\hfill & 0\text{e}{\psi }_1=1\text{(funções ímpares)}.\hfill \end{array}$$

A substituição desses valores na equação (31) mostra que a derivada primeira obedece a equação (37).

Podemos novamente fazer uso da simetria do sistema, encontrando a solução em apenas em metade do espaço. Se estivermos interessados no intervalo −L ⩽ x ⩽ L, podemos encontrar a solução apenas na região 0 ⩽ x ⩽ L (ou −L ⩽ x ⩽ 0). Uma vez obtida a autofunção em uma metade do intervalo, podemos encontrar a autofunção completa utilizando a equação (13). Isso é útil pois iniciamos a autofunção no ponto de simetria (x = 0), além de ser um algoritmo mais eficiente.

Ao mesmo tempo que buscamos a autofunção, também temos que encontrar a autoenergia. Até agora discutimos apenas como inicializar a autofunção no espaço. Todavia, para utilizar a equação (35), precisamos também de um valor para E. Como não conhecemos as autoenergias de antemão, vamos encontrá-las por tentativa e erro. Começamos com um valor E de nossa escolha e utilizamos a equação (35) para encontrar todos {ψ_i}. Se E não corresponder a nenhuma autoenergia do sistema, então as condições da seção 2.2 não serão satisfeitas.

Em particular, como estamos interessados em potenciais confinantes, teremos regiões em que E ＜ V (x), ou seja, regiões classicamente proibidas. Como vimos na seção 2.4, e ilustrado na Figura 2, os sistemas quânticos são capazes de adentrar essas regiões, porém as autofunções decrescem exponencialmente nelas. Essa característica, contudo, somente vai ocorrer para as autoenergias E_n do sistema. Para energias diferentes, a função de onda diverge na região classicamente proibida, como mostra a Figura 3.

Figura 3
Representação esquemática da função de onda obtida com a equação (35) quando a energia E difere da energia E₀ do estado fundamental do sistema. Na figura (a) temos uma energia maior que a verdadeira (E ＞ E₀), e vemos que a função de onda diverge para o infinito negativo, enquanto em (b) temos uma energia menor que a verdadeira (E ＜ E₀) e vemos a divergência para o infinito positivo. Logo, o sinal da divergência muda quando passamos pela energia verdadeira.

Este comportamento da divergência para as regiões |x| ＞ L é o ponto-chave para a determinação da energia no método da tentativa. Como o sinal da divergência muda quando passamos de E ＜ E₀ para E ＞ E₀, temos como nos aproximar iterativamente do valor de E₀ utilizando um análogo do método da bisseção [¹⁰].

Começamos com um valor tentativa da energia E_T, obtemos a função de onda utilizando a equação (35) e salvamos o sinal da divergência na região proibida (+∞ ou −∞). Realizamos um incremento na energia, E_T ← E_T + dE e repetimos o processo até o sinal da divergência mudar entre duas iterações consecutivas, indicando que a energia E₀ está entre as duas energias empregadas. Nesse momento, poderíamos parar o nosso algoritmo e teríamos uma estimativa para E₀ com erro de dE. No entanto, podemos obter um valor mais preciso com uma simples modificação. Quando o sinal da divergência mudar, fazemos dE ← −dE/2, ou seja, ao invés de aumentar a energia tentativa devemos diminuí-la, e o fator 1/2 é para centrarmos a nova tentativa entre os dois últimos valores da energia empregados. Repetimos esse processo até que dE seja tão pequeno quanto desejarmos, momento em que teremos encontrado a solução numérica do problema. Uma vez que obtivermos a energia correta, basta normalizar a autofunção utilizando a equação (4) para obtermos o estado físico.

Podemos resumir o algoritmo do método da tentativa da seguinte maneira:

1. Escolha L que determina o intervalo de integração 0 ＜ x ＜ L, o número de pontos N para a sua discretização (ou equivalentemente o espaçamento Δx entre os pontos) e uma tolerância ε para o valor da energia.

2. Escolha os valores iniciais da energia tentativa E_T e incremento dE.

3. Escolha a paridade da autofunção, e inicie a função de onda (ψ₀ e ψ₁) usando a equação (38).

4. Vamos utilizar a variável ultima_div para salvar o sinal da divergência no passo anterior. Para o primeiro passo, podemos tomar ultima_div como sendo positiva se dE ＞ 0.

5. Utilize a equação (35) e o valor atual de E_T para obter todos os {ψ_i} no intervalo escolhido. Na região proibida a função de onda vai divergir. Salve o sinal dessa divergência em uma variável nova_div.

6. Analise se a divergência mudou de direção:

   • (a)

     Se (ultima_div)×(nova_div)＞, a divergência não mudou de sinal. Não altere dE. Faça nova_div←ultima_div.

   • (b)

     Se (ultima_div)×(nova_div)＜, a divergência mudou de sinal, ou seja, passamos pela energia verdadeira. Tome dE ← −dE/2 e nova_div← ultima_div

7. Compare dE e a tolerância ε:

   • (a)

     Se |dE| ＜ ε, obtivemos a energia verdadeira dentro da precisão escolhida. Utilize a equação (13) para obter a autofunção no intervalo completo −L ⩽ x ⩽ L. Normalize a solução encontrada de acordo com a equação (4).

   • (b)

     Caso contrário, tome E_T ← E_T+dE e retorne ao passo 5.

Vamos exemplificar o método aplicando-o a potenciais que possuem solução analítica. Isso nos permitirá testar a implementação do código e a sua precisão.

3.1. Poço infinito

Consideremos o poço como primeiro exemplo. Nas unidades da seção 2.7 [ℓ = L e ∈ = ℏ²/(mL²)], podemos escrever:

(39) $$V(x)=\left\{\begin{array}{ll}0\hfill & \text{se}-1\le x\le 1,\hfill \\ V_0\hfill & \text{caso contrário}\text{.}\hfill \end{array}\right.$$

Esse potencial é par com relação a x = 0 [V (x) = V (−x)], logo vamos resolver o problema com relação a esse ponto. Chamamos de poço infinito o caso em que V₀ → ∞, porém o método do chute também funciona para valores finitos de V₀. Como precisamos analisar a divergência da função de onda na região classicamente proibida (fora do poço), escolhemos o intervalo de integração ligeiramente maior que o poço, 0 ≤ x ≤ 1, 01. A solução analítica deste problema é conhecida [⁴], sendo as autofunções dadas por:

(40) $${\psi }_n(x)=\left\{\begin{array}{ll}\cos (n\pi x/2)\hfill & \text{com}n=1,3,5,\dots ,\hfill \\ \sin (n\pi x/2)\hfill & \text{com}n=2,4,6,\dots ,\hfill \end{array}\right.$$

e as autoenergias,

(41) $$E_n={\pi }^2{n}^2/8.$$

Ao desenvolver programas que buscam soluções numéricas é desejável, sempre que possível, comparar os resultados com problemas de solução analítica conhecida. Isso permite verificar se o programa está correto e também fica evidente a precisão da técnica. Assim, vamos aplicar o método da tentativa ao poço infinito. Os resultados obtidos para as energias dos cinco primeiros estados ligados do sistema estão expostos na Tabela 1.

Tabela 1
Energia do estado fundamental e dos quatro primeiros estados excitados do poço infinito obtidos através do método da tentativa. Foram utilizados 10000 pontos na discretização do intervalo 0 ≤ x ≤ 1. A integração foi feita até o limite |dE| ≤ 0, 001. O poço infinito foi simulado utilizando-se V₀ = 1 × 10¹⁰. Comparamos os valores obtidos com os valores teóricos da equação (40).

n	E Num	E Teo	|ETeo − ENum
1	1,23379	1,23370	0,00009
2	4,9350	4,9348	0,0002
3	11,105	11,103	0,002
4	19,740	19,739	0,001
5	30,845	30,843	0,002

Vemos que os valores obtidos numericamente para a energia E^Num concordam bem com os valores teóricos E^Teo da equação (41). No entanto, as divergências tendem a aumentar para os estados excitados. Isso se deve principalmente ao espaçamento finito Δx, o qual vai fornecer soluções menos precisas quão maior for o número de nós nas autofunções, e ao fato do poço não ser de fato finito (nas soluções da tabela 1 utilizamos V₀ = 1 × 10¹⁰), sendo uma representação do modelo teórico no computador.

As autofunções dos três primeiros estados são apresentadas na Figura 4. O número de nós aumenta conforme incrementamos o número quântico n, como esperado. Além disso, vemos que as autofunções numéricas concordam com as teóricas da equação (40), com o estado fundamental sendo uma função cosseno, e o primeiro estado excitado uma função seno.

Figura 4
Autofunções do estado fundamental e dos dois primeiros estados excitados do poço infinito obtidos através do método da tentativa. A solução foi obtida na região de integração, 0 ＜ x ＜ 1, e é representada por linhas cheias. Na região −1 ≤ x ＜ 0, ou seja, x → −x, obtemos a solução através da simetria de paridade, equação (13), representadas pelas linhas pontilhadas.

3.2. Oscilador Harmônico

Em seguida, aplicamos o método da tentativa para o oscilador harmônico,

(42) $$V(x)=\frac{1}{2}{x}^2.$$

As soluções analíticas para os autoestados e autoenergias são dadas por [⁷]

(43) $$\begin{array}{lll}{\psi }_n(x)\hfill & =\hfill & \sqrt{\frac{1}{2^{n}n!\sqrt{\pi }}}{H}_n(x)e^{-\frac{x^2}{2}},n=0,1,\dots ,\hfill \\ E_n\hfill & =\hfill & n+\frac{1}{2},\hfill \end{array}$$

onde adimensionalizamos o problema utilizando ℓ = [ℏ/(mω)]^1/2 e ϵ = ℏω, como descrito na seção 2.7. As soluções analíticas dependem de funções especiais (polinômios de Hermite, H_n(x), representados na Figura 5 [⁷]), e a interpretação das expressões analíticas requer familiaridade com os H_n(x). Todavia, podemos aplicar a mesma metodologia da seção anterior para resolver o problema numericamente.

Figura 5
Representação dos cinco primeiros polinômios de Hermite.

Para utilizar a equação (35), vamos novamente nos valer da simetria do potencial com relação a x = 0. Queremos as soluções do sistema no intervalo −L ⩽ x ⩽ L, com L definido pelo usuário do programa. Contudo, temos uma situação diferente do poço infinito, equação (39), em que o potencial é infinito para x ＞ 1 ou x ＜ −1. No caso do oscilador harmônico, equação (42), vemos que o potencial é finito em todo intervalo de interesse (com o valor máximo de L²/2 para x = ±L). Um potencial infinito não permite a penetração da região classicamente proibida, com o autoestado tendo que se anular imediatamente nos pontos x = 1 ou x = −1, como vemos na Figura (4). O fato do potencial harmônico ser finito permite que haja penetração da região classicamente proibida. Assim, a escolha do intervalo relacionado a L é importante se queremos observar o autoestado se anular para grandes |x|, pois valores pequenos de L podem não ser suficientes (ver Figura 2). Temos que escolher um L suficientemente grande para que a penetração da região proibida ocorra e a autofunção decresça exponencialmente até valores negligenciáveis, supondo que empregamos a autoenergia correta. Para os primeiros autoestados verificamos que L = 5 é suficiente para que isso ocorra.

Apresentamos as energias obtidas na Tabela 2, e as autofunções correspondentes na Figura 6. Novamente, temos uma boa concordância com os resultados analíticos, equação (43). As autofunções da Figura 6 correspondem ao resultado esperado: polinômios de Hermite (Figura 5) modulados por uma gaussiana.

Tabela 2
Energia do estado fundamental e dos quatro primeiros estados excitados do oscilador harmônico obtidos através do método da tentativa. Foram utilizados 10 000 pontos na discretização e L = 5. Foi considerado o critério de convergência |dE| ⩽ 0, 00001. Comparamos os resultados numéricos com os valores analíticos, equação (43).

n	E Num	E Teo	$|E_n^{\text{num}}-E_n^{\text{teo}}|$
0	0,5001	0,5000	0,0001
1	1,500006	1,500000	0,000005
2	2,5003	2,5000	0,003
3	3,500006	3,500000	0,000006
4	4,5005	4,5000	0,0005

Figura 6
Autofunções do oscilador harmônico obtidas numericamente com o método da tentativa: o estado fundamental e os dois primeiros estados excitados. As soluções foram obtidas na região 0 ⩽ x ⩽ 5, e são representadas por linhas cheias. Na região −5 ⩽ x ＜ 0 construímos as soluções através da equação (13), representadas pelas linhas pontilhadas.

4. O Método do Encontro

O método da seção anterior se aplica a potenciais com paridade bem definida. Todavia, existem potenciais de interesse sem essa propriedade, com um método mais geral se fazendo necessário. Nessa seção descrevemos uma abordagem para potenciais unidimensionais quaisquer, a qual é chamada de método do encontro [⁴, ¹³].

Para tratar o problema devemos considerar propriedades gerais que as autofunções respeitam (vistas na seção 2.2): continuidade da autofunção e de sua derivada primeira e ser de quadrado integrável. A última implica que a autofunção deve se anular nos extremos do intervalo. Mais especificamente, como focamos em potenciais confinantes, sabemos da seção 2.4 que as autofunções decrescem exponencialmente na região classicamente proibida, sendo essencialmente nulas após determinado ponto.

Diferentemente da seção 3, em geral, não teremos potenciais com paridade bem definida e a equação (13) não será válida. Isso não significa que o método do encontro falha para potenciais com paridade bem definida, apenas que estamos desenvolvendo um método que se aplica a uma classe mais geral de potenciais. Logo, não há um ponto privilegiado do espaço, que na seção anterior era o ponto x = 0, e devemos buscar um novo modo de iniciar as autofunções.

O módulo quadrado da função de onda está relacionado com a densidade de probabilidade de encontrar a partícula ao redor de uma posição. Logo, esperamos que |ψ(x)|² seja apreciável próximo ao mínimo do potencial, e que conforme nos afastamos desse mínimo a função de onda decresça até se anular para distâncias muito grandes.

Um modo de inicializar as autofunções abrindo mão da paridade e garantindo que essas sejam de quadrado integrável é forçando-as a se anularem na região classicamente proibida. Para isso, podemos inciar ψ_i = 0 em um ponto da região classicamente proibida (de preferência longe da região permitida), o que impede que ela seja divergente nessa região, como era o caso da Figura 3. Alocamos um valor pequeno ${\psi }_i^{\prime }=\pm h$ para a sua derivada, garantindo que a autofunção cresça quando nos aproximamos da região de interesse, a região classicamente permitida. Dessa forma, a autofunção possui um valor não-nulo na região de relevante e decresce fora dessa, se anulando nos extremos do intervalo.

Essa prescrição vale para pontos à direita e à esquerda da região de interesse. Assim, no método do encontro, consideramos duas funções de onda ψ_e(x) e ψ_d(x). Iniciamos ψ_e(x) à esquerda da região classicamente permitida, com uma derivada pequena e positiva, para que ψ_e(x) seja crescente conforme adentramos a região permitida. Do mesmo modo, ψ_d(x) é iniciada à direita com derivada pequena e negativa, para que seja decrescente conforme nos afastamos da região permitida.

Estamos interessados no caso em que ambas funções coincidem e correspondem à autofunção verdadeira. Assim como na seção anterior, não sabemos a priori a autoenergia que devemos empregar no algoritmo para o cálculo da autofunção. Assim, quando integramos ψ_e(x) para valores crescentes de x (“da esquerda para direita”) passamos pela região permitida classicamente e adentramos novamente na região proibida, onde a função de onda diverge para energias diferentes das autoenergias, conforme discutido na seção 3. O mesmo acontece para ψ_d(x) quando integramos para valores decrescentes de x (“da direita para esquerda”). Podemos contornar o problema da divergência escolhendo um ponto de encontro x_m em que fazemos as duas funções se encontrarem. Assim, a solução ψ(x) do sistema é construída a partir das duas funções:

(44) $$\psi (x)=\left\{\begin{array}{ll}{\psi }_e(x)\hfill & x_e⩽x⩽x_m,\hfill \\ {\psi }_d(x)\hfill & x_m<x⩽x_d.\hfill \end{array}\right.$$

Escolhendo uma energia arbitrária para começar o algoritmo, e utilizando as considerações acima para a inicialização das funções de onda, podemos obter ψ_e(x) através da equação (35). Para ψ_d, precisamos “integrar da direita para a esquerda”, ou seja,

(45) $${\psi }_{i-1}=2{\psi }_i-{\psi }_{i+1}-2{(\text{Δ}x)}^2(E-V_i){\psi }_i.$$

Muito provavelmente as duas funções de onda terão valores diferentes em x_m, como ilustrado na Figura 7. Estamos interessados em normalizar o autoestado somente ao final do procedimento, logo temos liberdade para multiplicar soluções da equação de Schrödinger por constantes, afinal elas também são soluções. Assim, para garantir a continuidade da autofunção, podemos reescalar uma das funções para que os valores se igualem em x_m. Na Figura 7, mostramos o resultado de multiplicar ψ_e(x) por ψ_d(x_m)/ψ_e(x_m).

Figura 7
Ilustração da aplicação do método do encontro ao oscilador harmônico, com x_m = −1. As funções ψ_e(x) e ψ_d(x) (curvas contínuas), obtidas com as equações (35) e (45), diferem no ponto de encontro x_m. Aproveitamos a linearidade da equação de Schrödinger para construir a solução ψ_d(x_m)ψ_e(x)/ψ_e(x_m) (curva pontilhada) a qual, por construção, coincide com ψ_d(x_m) no ponto de encontro. Para energias menores (a) ou maiores (b) que E₀, a derivada não é contínua em x_m, enquanto que uma derivada contínua (c) indica que encontramos uma autoenergia do sistema.

Temos também que garantir a continuidade da derivada primeira de ψ(x). Por construção, o único ponto em que essa pode ser descontínua é no ponto de encontro. Esse é o cerne do método do encontro: ψ′(x) só será contínua para as autoenergias do sistema. Assim, podemos variar a energia empregada nas equações (35) e (45) até encontrar um valor que faça ψ′(x_m) ser contínua. Comparamos as inclinações das funções de onda no ponto de encontro x_m e quando a diferença entre as duas for menor que uma tolerância de nossa escolha, teremos obtido uma solução. Para comparar as inclinações, consideramos a diferença entre as derivadas no ponto de encontro. Assim como na seção anterior, quando essa mudar de sinal teremos passado pelo valor da energia verdadeira.

Com isso, podemos resumir o algoritmo do método do encontro como:

1. Escolha os dois limites do intervalo, pontos x_e (limite esquerdo) e x_d (limite direito), o número de pontos N para a discretização (ou equivalentemente o espaçamento Δx entre os pontos), o ponto de encontro x_m e uma tolerância ε para a descontinuidade da derivada.

2. Escolha os valores iniciais da energia tentativa E_T e incremento dE.

3. Inicie as duas funções de onda ψ_e(x) e ψ_d(x):

   (46)$$\left\{\begin{array}{ll}{\psi }_e(x):{\psi }_0=0\hfill & \text{e}{\psi }_1=\text{Δ}x,\hfill \\ {\psi }_d(x):{\psi }_{N-1}=\text{Δ}x\hfill & \text{e}{\psi }_N=0.\hfill \end{array}\right.$$

   Ao invés de Δx, pode-se usar outra quantidade pequena.

4. Vamos utilizar a variável ultima_dif para salvar o sinal da diferença entre as derivadas. Para o primeiro passo, podemos tomar ultima_dif como sendo positiva se dE ＞ 0.

5. Utilize a equação (35) para obter ψ_e(x) e a equação (45) para ψ_d(x).

6. Multiplique toda a função ψ_e(x) pelo número ψ_d(x_m)/ψ_e(x_m).

7. Utilize a equação (30) para calcular as duas derivadas ${\psi }_e^{\prime }(x_m)$ e ${\psi }_d^{\prime }(x_m)$, e calcule a diferença entre elas $\text{Δ}={\psi }_e^{\prime }(x_m)-{\psi }_d^{\prime }(x_m)$. Salve o sinal de Δ em nova_dif.

8. Caso o módulo da diferença seja menor que a tolerância, |Δ| ＜ ε, encontramos a energia correta dentro da precisão especificada. Escreva a autofunção ψ(x) do sistema seguindo a equação (44) e normalize-a usando a equação (4).

9. Analise se a diferença mudou de sinal:

   • (a)

     Se (ultima_dif)×(nova_dif)＞, a diferença não mudou de sinal. Não altere dE. Faça nova_dif← ultima_dif.

   • (b)

     Se (ultima_dif)×(nova_dif)＜ 0, a diferença mudou de sinal, ou seja, passamos pela energia verdadeira. Tome dE ← −dE/2 e nova_dif← ultima_dif.

10. Tome E_T ← E_T + dE e volte ao passo 5.

Para exemplificar, vamos considerar o potencial de Lennard-Jones [⁴],

(47) $$V(x)=4\alpha \left[{\left(\frac{\sigma }{x}\right)}^{12}-{\left(\frac{\sigma }{x}\right)}^6\right],$$

ilustrado na Figura 8. Apesar de ser empregado em várias áreas da Física e Química, esse potencial não possui solução analítica. Seguindo a prescrição da seção 2.7 para tornar o problema adimensional, consideramos $\overline{x}=x/\sigma$ e $\overline{E}=Em{\sigma }^2/{\hslash }^2$. A profundidade α do potencial determina o número de estados ligados existentes, ou seja, potenciais “rasos” possuem menos estados ligados. Isso ocorre pois os níveis de energia possuem um espaçamento entre si. Logo, caso o potencial não seja profundo o suficiente, esse espaçamento pode não ser possível, não havendo estados excitados ou mesmo estados ligados. Assim, escolhemos $\overline{\alpha }=\alpha m{\sigma }^2/{\hslash }^2=200$ para garantir que tenhamos um potencial “profundo” o suficiente para a existência de estados excitados. É importante ressaltar que estados de espalhamento existem para qualquer profundidade do potencial, mas o nosso foco nesse artigo é exclusivamente lidar com estados ligados.

Figura 8
O potencial de Lennard-Jones para σ = 1 e α = 200. Ele possui uma forte repulsão a curtas distâncias e uma fraca atração para grandes valores de x. Sua profundidade é dada por −α e seu mínimo é em 2^1/6σ. A linha azul em x_m indica o ponto de encontro adotado no algoritmo.

Em seguida temos que escolher o ponto de encontro x_m, o qual possui interferência direta no resultado obtido. É conveniente escolher um ponto dentro da região classicamente permitida, onde as duas funções de onda possuem um valor relativamente grande, o que vai garantir maior precisão no cálculo da diferença das derivadas. É por esse motivo que não tomamos x_m próximo a x_e ou x_d, situação na qual seria necessário um maior número de pontos na integração da equação (35).

Esperamos que a função de onda seja apreciável próximo ao mínimo do potencial. Inspeção da Figura 8 mostra que a escolha x_m = 1,40 é razoável. Assim, somos capazes de obter cinco estados ligados, cujas energias são listadas na Tabela 3. Nossos resultados mostram o comportamento esperado: quanto mais excitado o estado, maior sua energia. As três primeiras autofunções encontram-se na Figura 9 e também apresentam o comportamento esperado, com um aumento do número de nós para estados excitados.

Tabela 3
Energias obtidas para os cinco estados ligados obtidos do potencial de Lennard-Jones com o método do encontro. Integramos o sistema no intervalo x_e = 0,5 ⩽ x ⩽ x_d = 5,5. Foi escolhido α = 200, x_m = 1,4 e foram utilizados 10 000 pontos na integração numérica. A diferença mínima entre as derivadas (critério de parada do método) foi escolhida como Δ = 1 × 10⁻⁶.

n	En
1	−150,642
2	−77,210
3	−32,004
4	−8,865
5	−0,836

Figura 9
Autofunções dos três primeiros estados ligados do potencial de Lennard-Jones obtidos com o método do encontro. As linhas cheias denotam ψ_e(x) e as linhas pontilhadas ψ_d(x). Elas se encontram em x_m para formar a solução no intervalo inteiro. Note que para a energia correta temos uma autofunção ψ_n(x) contínua, e o mesmo para sua derivada.

5. Uma Abordagem Matricial

Os métodos das seções anteriores abordam a equação de Schrödinger sob a perspectiva da solução numérica de equações diferenciais. Contudo, a equação (7) é uma equação de autovalor-autovetor e podemos utilizar esse fato para obter a solução do sistema. Para isso, consideremos a equação (34) e a discretização do espaço em N + 2 pontos {x₀, x₁, x₂, . . . , x_N, x_N+1} igualmente espaçados de Δx. Podemos escrever um sistema de N equações:

(48) $$\begin{array}{l}-\frac{1}{2}\left[\frac{{\psi }_2-2{\psi }_1+{\psi }_0}{{(\text{Δ}x)}^2}\right]+V_1{\psi }_1=E{\psi }_1,\hfill \\ -\frac{1}{2}\left[\frac{{\psi }_3-2{\psi }_2+{\psi }_1}{{(\text{Δ}x)}^2}\right]+V_2{\psi }_2=E{\psi }_2,\hfill \\ ⋮\hfill \\ -\frac{1}{2}\left[\frac{{\psi }_{N+1}-2{\psi }_N+{\psi }_{N-1}}{{(\text{Δ}x)}^2}\right]+V_N{\psi }_N=E{\psi }_N.\hfill \end{array}$$

A razão pela qual não podemos escrever as equações para os pontos i = 0 e N + 1 é porque elas iriam depender de i = −1 e i = N + 2, respectivamente, os quais estão fora do intervalo considerado. Resolvemos esse problema tomando as condições de contorno ψ₀ = ψ_N+1 = 0, ou seja, a função de onda se anula nos extremos do intervalo.

Aplicando as condições de contorno no sistema de equações (48), podemos colocar o problema na forma matricial

(49) $$\mathscr{H}\psi =E\psi ,$$

onde ψ = (ψ₁, ψ₂, . . . , ψ_N)^T é o vetor coluna que representa a autofunção e

(50) $$\mathscr{H}=\frac{1}{2{(\text{Δ}x)}^2}\left(\begin{array}{ccccccc}{\tilde{H}}_1& -1& 0& 0& 0& \dots & 0\\ -1& {\tilde{H}}_2& -1& 0& 0& \dots & 0\\ 0& -1& {\tilde{H}}_3& -1& 0& \dots & 0\\ ⋮& ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& 0& 0& \dots & {\tilde{H}}_N\end{array}\right),$$

é a matriz que representa o operador Hamiltoniano do sistema, com ${\tilde{H}}_i=2+2{(\text{Δ}x)}^2{V}_i$.

Diagonalizando a matriz da equação (50) resolvemos o problema de autovalor-autovetor diretamente. Uma das grandes vantagens desse método é que ele se aplica para um potencial geral. As únicas escolhas necessárias são o intervalo de interesse da autofunção, compatível com as condições de contorno ψ₀ = ψ_N+1 = 0, e o número de pontos em que o intervalo será particionado. Fixados esses parâmetros, aplicamos um método de diagonalização [¹⁰], ou utilizamos uma biblioteca para tal.

Note que o Hamiltoniano da equação (50) é uma matriz tridiagonal simétrica, um tipo especial de matriz para a qual rotinas especiais de diagonalização existem. Em particular, uma rotina de interesse para o problema em questão é a dstemr da biblioteca LAPACK [¹⁴]. Essa rotina calcula os autovalores e (caso seja desejado) os autovetores de uma matriz tridiagonal simétrica e real, que é o nosso caso. Podemos controlar quantos autoestados do espectro queremos calcular através de variáveis da rotina. Uma lista completa dos parâmetros da rotina bem como a sua descrição podem ser obtidas na documentação da biblioteca [¹⁵]. No repositório [⁵] apresentamos nosso código em Fortran90, mas existem implementações que utilizam essa mesma rotina em outras linguagens, como C e Python [¹⁶].

Para exemplificar o método, vamos considerar o átomo de hidrogênio. A princípio trata-se de um problema de duas partículas interagentes, próton e elétron, em três dimensões. No entanto, como o potencial de Coulomb é esfericamente simétrico, podemos utilizar o que foi discutido na seção 2.5 para transformar esse problema no equivalente de uma partícula em 1D.

Seguindo a prescrição da seção 2.7 para tornar o problema adimensional, identificamos a escala de comprimento típica do sistema como sendo a_μ = m_ea₀/μ, onde μ = m_pm_e/(m_p + m_e) é a massa reduzida, m_p é a massa do próton, m_e é a massa do elétron e a₀ é o raio de Bohr. Já a escala típica de energia é dada por ${\hslash }^2/(\mu a_{\mu }^2)$. Com essas considerações, estamos interessados em resolver a equação (19) com o potencial eletrostático,

(51) $$-\frac{1}{2}\frac{d^2{u}_{n\ell }(r)}{d{r}^2}+\left[-\frac{1}{r}+\frac{\ell (\ell +1)}{2r^2}\right]u_{n\ell }(r)=E_{n\ell }{u}_{n\ell }(r).$$

As autoenergias, conhecidas analiticamente [⁷], são dadas por

(52) $$E_n=-\frac{1}{2n^2}.$$

O potencial de Coulomb é particular no sentido que as energias dependem somente do número quântico principal. No entanto, os autoestados dependem do valor de ℓ = 0, 1, . . . , n − 1, de maneira que precisamos dar como entrada do programa o número quântico ℓ associado com o autoestado que estamos buscando.

Devemos salientar que a equação (51) é para a função radial reduzida u_nℓ(r), a qual está relacionada com a função radial do átomo de hidrogênio por R_nℓ(r) = u_nℓ(r)/r. Como teríamos uma divergência se u_nℓ(r → 0) fosse finita, tomamos como uma das condições de contorno u_nℓ(0) = 0. A outra condição é que para grandes valores de r, a função radial reduzida deve se anular, pois estamos buscando um estado ligado. Embora as condições de contorno adotadas estejam corretas, uma justificação rigorosa requer uma análise mais detalhada. Para o leitor interessado, discutimos mais o assunto no Apêndice A. Assim vemos que esse problema obedece as condições de contorno utilizadas acima quando introduzimos a abordagem matricial, u₀ = 0 e u_N+1 = 0, usando a notação compacta u_nℓ(r_i) = u_i.

A matriz relevante para esse sistema é a mesma que a equação (50), substituindo-se Δx por Δr e com

(53) $${\tilde{H}}_n=2+2{(\text{Δ}r)}^2\left(-\frac{1}{r_i}+\frac{\ell (\ell +1)}{2r_i^2}\right).$$

As expressões analíticas das autofunções são conhecidas, embora sejam dadas em termos de funções especiais [⁷]. Resolvendo o problema com o método matricial podemos obter tais autofunções, bem como as energias do sistema. Para isso, tomamos o intervalo 0 ⩽ r ⩽ 30 com N =10 000 pontos e diagonalizamos a matriz. As energias obtidas encontram-se na Tabela 4.

Tabela 4
Comparação entre as energias obtidas para o átomo de Hidrogênio através do método matricial $E_n^{\text{num}}$ com os valores analíticos conhecidos $E_n^{\text{teo}}$, equação (52). Consideramos os três estados de mais baixa energia com ℓ = 0 e ℓ = 1. Discretizamos o sistema na região 0 ⩽ r ⩽ 30 utilizando 10 000 pontos.

n	ℓ	$E_n^{\text{num}}$	$E_n^{\text{teo}}$	$\mathrm{|}{E}_n^{\text{num}}-E_n^{\text{teo}}\mathrm{|}$
1	0	−0,49999887	−0,50000000	1,1 × 10−6
2	0	−0,12499992	−0,12500000	8,0 × 1010−8
3	0	−0,05542365	−0,05555555	1,3 × 10−4
2	1	−0,12500002	−0,12500000	2,0 × 10−8
3	1	−0,05547129	−0,05555555	8,4 × 10−4

Os valores de energia obtidos numericamente apresentam um bom acordo com os valores conhecidos, com divergências aparecendo apenas a partir da ordem de 10⁻⁴ para o segundo estado excitado. Note que os níveis de energia são degenerados, com as energias dos estados n = 2 e 3 independentes do valor de ℓ, como esperado. Podemos também obter as autofunções do sistema, expostas na Figura 10. Vemos que o número de nós das autofunções é dado por n − 1 − ℓ, utilizando a convenção em que os zeros em r = 0 e r → ∞ não são contabilizados. Esse resultado, bem como as autofunções, são conhecidos da literatura, mas sua obtenção através dos métodos dessa seção é relativamente simples e permite a compreensão das implicações de forma que as particularidades das funções especiais presentes nas soluções analíticas não se fazem necessárias.

Figura 10
Densidades radiais reduzidas de probabilidade |u_nℓ(r)|² = |rR_nℓ(r)|² obtidas para os primeiros estados ligados do átomo de Hidrogênio através do método matricial. Utilizamos N =10 000 pontos no intervalo 0 ⩽ r ⩽ 30.

Embora o átomo de hidrogênio seja a princípio um problema de duas partículas em três dimensões, utilizar coordenadas do centro de massa e um potencial esfericamente simétrico fez com que a equação a ser resolvida fosse similar ao problema de uma partícula em uma dimensão. Podemos generalizar esse método para sistemas verdadeiramente bidimensionais e tridimensionais. Temos que discretizar o espaço, agora em uma rede, e construir a função de onda e o operador Hamiltoniano de modo análogo. No entanto, o custo computacional torna essa abordagem inviável, e outros tipos de métodos são necessários.

6. O Método Variacional de Monte Carlo

Os métodos abordados até agora são voltados para solucionar a equação de Schrödinger para uma partícula em uma dimensão. Para duas partículas em potenciais esfericamente simétricos, mostramos que a solução numérica é similar se considerarmos coordenadas do centro de massa e a função de onda radial reduzida, ver seção 2.5 e aplicação para o átomo de hidrogênio, seção 5.

Para sistemas com mais de duas partículas tal redução não é possível, e os métodos apresentados até aqui não podem ser aplicados. Além disso, em três dimensões, caso o potencial não seja central [ou seja, tenhamos um potencial V (r, θ, ϕ)] temos que resolver a equação (16), uma equação diferencial parcial, com uma abordagem distinta da adotada. Por fim, mesmo em problemas unidimensionais, a obtenção dos autoestados pode ser dificultada pela complexidade do potencial. Por exemplo, quando estudamos o dêuteron (um isotopo do hidrogênio), podemos utilizar os métodos da seção 2.5 para obter equações unidimensionais acopladas. Contudo, como o potencial inclui graus de liberdade de spin, isospin, e momento angular, não temos a situação mais simples da equação (19), e os métodos desenvolvidos até agora não se aplicam. Logo, queremos apresentar um método capaz de estimar a energia de sistemas tridimensionais com muitas partículas cujo custo computacional escala mais favoravelmente do que abordagens similares às apresentadas até agora.

O método que vamos apresentar é chamado de método variacional, e nos fornece um valor máximo para a energia do sistema. Para obtê-lo, lembramos que a equação (2) implica que qualquer função de onda pode ser escrita como uma combinação linear das autofunções do sistema, desde que satisfaça as condições de contorno do problema. Assim, se calcularmos o valor esperado da energia E_V do sistema para uma função de onda qualquer ϕ_T(r) a qual chamamos de função de onda teste, temos [⁶]

(54) $$\mathrm{\langle }H\mathrm{\rangle }\equiv E_V=\frac{\mathrm{\langle }{\varphi }_T\left|H\right|{\varphi }_T\mathrm{\rangle }}{\mathrm{\langle }{\varphi }_T|{\varphi }_T\mathrm{\rangle }},$$

onde o denominador garante que o resultado mesmo para estados não-normalizados. Utilizando a expansão da equação (2) temos

(55) $$E_V={\displaystyle \sum }_n{\displaystyle \sum }_m{c}_n^{*}{c}_m\left(\int dx{\psi }_n^{*}(x)H{\psi }_m(x)\right),$$

sendo H o operador Hamiltoniano. Sabemos da equação (7) que Hψ_n = E_nψ_n, logo

(56) $$E_V={\displaystyle \sum _n{\displaystyle \sum _m{c}_n^{*}}}{c}_m{E}_m\left({\displaystyle \int d}x{\psi }_n^{*}(x){\psi }_m(x)\right),$$

utilizando a ortonormalidade das autofunções, equação (3),

(57) $$E_V={\displaystyle \sum _n{\displaystyle \sum _m{c}_n^{*}}}{c}_m{E}_m{\delta }_{mn}={\displaystyle \sum _n{c}_n^{*}}{c}_n{E}_n.$$

Separamos a soma em dois termos,

(58) $$E_V=|c_0{|}^2{E}_0+\left({\displaystyle \sum _{n=1}^{\infty }|}{c}_n{|}^2{E}_n\right)⩾E_0,$$

a qual é o cerne do princípio variacional [⁶]. A equação (58) mostra que temos um limite superior para a energia obtida, ou seja, E_V não pode ser menor que o valor da energia do estado fundamental. Além disso, quanto mais a função ϕ_T se aproximar do verdadeiro estado fundamental, menores se tornaram os coeficientes c_n com n ≠ 0 e mais nos aproximaremos da verdadeira energia E₀.

Como queremos estimar a energia de um sistema o qual não necessariamente conhecemos a solução analítica, temos que propor uma função de onda teste ϕ_T. Sabemos que ϕ_T não é a solução exata do problema, mas queremos que se aproxime dessa. A escolha da função de onda teste é feita levando-se em consideração as características do potencial, suas condições de contorno e o máximo de informação física que conhecemos sobre o problema.

Além disso, podemos adicionar parâmetros a essa função de onda teste e variá-los, buscando o conjunto de parâmetros que fornece a menor energia. Assim, podemos estimar a energia do estado fundamental propondo uma função de onda teste ϕ_T (r, {α_i}), onde {α_i} é um conjunto de parâmetros. Ao calcularmos a energia associada com essa função de onda teste encontraremos uma energia variacional que depende dos parâmetros α_i, E_V(α_i). Podemos minimizá-la com relação aos parâmetros α_i afim de nos aproximarmos da energia verdadeira, assim, garantindo que a função de onda teste utilizada se aproxima o máximo possível do autoestado verdadeiro.

Vamos exemplificar a escolha de uma função de onda teste para a aplicação do método buscando a energia do estado fundamental do átomo de Hélio, a qual, experimentalmente sabemos ser E_He = −79, 02 eV. O Hamiltoniano do sistema é

(59) $$H=-\frac{{\hslash }^2}{2m_e^2}({\nabla }_1^2+{\nabla }_2^2)-\frac{Z{e}^2}{r_1}-\frac{Z{e}^2}{r_2}+\frac{e^2}{|{\mathit{r}}_1-{\mathit{r}}_2|},$$

onde m_e e e são a massa e a carga do elétron, Z = 2 é o número atômico do Hélio, e o núcleo está fixo na origem. Os operadores ${\nabla }_i^2$ atuam apenas nas coordenadas para partícula i. O último termo desse Hamiltoniano representa a interação elétron-elétron, e dificulta a solução analítica do problema. Contudo, os demais termos formam a soma de dois Hamiltonianos de átomos de um elétron [⁷],

(60) $$\begin{array}{lll}H\hfill & =\hfill & \left(-\frac{{\hslash }^2}{2m_e^2}{\nabla }_1^2-\frac{Z{e}^2}{r_1}\right)+\left(-\frac{{\hslash }^2}{2m_e^2}{\nabla }_2^2-\frac{Z{e}^2}{r_2}\right)+\frac{e^2}{|{\mathit{r}}_1-{\mathit{r}}_2|}\hfill \\ \hfill & =\hfill & H_1^{\text{átomo}1e^{-}}+H_2^{\text{átomo}1e^{-}}+\frac{e^2}{|{\mathit{r}}_1-{\mathit{r}}_2|}.\hfill \end{array}$$

Logo, sabemos que a solução do sistema possui semelhanças com a de átomos de um elétron, sendo o caso particular do átomo de Hidrogênio (Z = 1) discutido na seção 5.

Teremos desvios da solução válida para átomos de um elétron devido a interação elétron-elétron, que depende da separação entre eles. Em primeira aproximação desconsiderarmos essa interação e propomos uma função de onda teste que seja o produto de duas autofunções do estado fundamental de átomos de um elétron com Z = 2 [⁶]. Podemos melhorar essa aproximação permitindo que a constante do decaimento exponencial seja um parâmetro variacional α tal que

(61) $${\varphi }_T({\mathit{r}}_1,{\mathit{r}}_2,\alpha )=A\exp \left[-\alpha (r_1+r_2)\right],$$

onde A é uma constante de normalização. Estamos interessados em encontrar o valor de α que fornece a menor energia variacional.

Os resultados numéricos apresentados nessa seção foram obtidos com um programa que emprega grandezas adimensionais, como descrito na seção 2.7. Para isso, escrevemos a energia como sendo $\tilde{E}=E/ϵ$ e utilizamos $ϵ={\hslash }^2/a_0^2{m}_e$, onde o raio de Bohr a₀ é uma escala de comprimento muito apropriada para sistemas atômicos. Como conhecemos a energia do sistema nas unidades de elétron-volt (eV) queremos escrever ϵ nas mesmas unidades. Para isso, multiplicamos por c² o numerador e denominador, $ϵ={(\hslash c)}^2/a_0^2{m}_e{c}^2$, onde c é a velocidade da luz. Sabendo que ℏc = 197, 327 MeV fm, m_ec² = 0,511 MeV e a₀ = 0, 529Å ¹, podemos escrever ϵ = 139, 14 eV. Para recuperar as unidades físicas basta multiplicar a energia adimensional obtida em nosso programa por ϵ ($(E[\text{eV}]=ϵ[\text{eV}]\tilde{E})$), o que é ilustrado na Tabela 5.

Tabela 5
Energias adimensionais (${\tilde{E}}_V$) e com dimensão de eV (E_V) obtidas com o método VMC para o estado fundamental do átomo de Hélio. Os valores do parâmetro variacional α utilizados são próximos do mínimo obtido nas Figuras 11 e 12.

a 0 α	${\tilde{E}}_V$	EV [eV]
1,55	-2,802(3)	-76,72(9)
1,60	-2,829(3)	-77,01(9)
1,70	-2,846(3)	-77,47(8)
1,80	-2,836(3)	-77,21(7)
1,90	-2,793(3)	-76,02(8)

Podemos calcular o valor da energia através das integrais que aparecem na equação (54). Todavia, temos integrais multidimensionais. Mais especificamente, se tivermos N partículas em um sistema d-dimensional, para calcular a energia variacional temos dN integrais no numerador e denominador da equação (54). Para sistemas físicos de interesse, esse problema pode se tornar extremamente complexo e inviável. Contudo, podemos novamente recorrer à física computacional, e calcular essas integrais através de métodos numéricos. Particularmente, uma maneira eficiente de amostrar integrais multidimensionais é chamada de método de Monte Carlo, e sua aplicação no princípio variacional é conhecida como método variacional de Monte Carlo (VMC) [¹⁷, ¹⁸].

Para aplicar o método VMC, vamos escrever a energia variacional da equação (54) no espaço das posições, de maneira que as integrais são feitas em d^3N R, onde temos N partículas em 3 dimensões, e R = {r₁, . . . , r_N representa o vetor posição de todas as partículas do sistema (no exemplo do átomo de Hélio, R = {r₁, r₂}). Assim, temos

(62) $$E_V=\frac{{\displaystyle \int d^{3N}}R{\varphi }_T^{*}(R)H{\varphi }_T(R)}{{\displaystyle \int d^{3N}}R|{\varphi }_T(R){|}^2}.$$

Multiplicamos e dividimos o integrando do numerador por ϕ_T (R),

(63) $$E_V=\frac{{\displaystyle \int d^{3N}}R{\varphi }_T^{*}(R){\varphi }_T(R)\frac{H{\varphi }_T(R)}{{\varphi }_T(R)}}{{\displaystyle \int d^{3N}}R|{\varphi }_T(R){|}^2}.$$

Com isso, a energia variacional pode ser escrita como [¹⁷]

(64) $$E_V={\displaystyle \int d^{3N}}RP(R)E_L(R),$$

com uma energia local

(65) $$E_L(R)=\frac{H{\varphi }_T(R)}{{\varphi }_T(R)},$$

e uma densidade de probabilidade

(66) $$P(R)=\frac{|{\varphi }_T(R){|}^2}{{\displaystyle \int d^{3N}}R|{\varphi }_T(R){|}^2}.$$

Note que essa função é normalizada e não-negativa, logo pode ser interpretada como uma densidade de probabilidade.

Uma das grandes vantagens do método variacional é calcular a integral da equação (64) por amostragem. Para isso, pensamos em P(R) como um indicador de quais configurações dão maior contribuição para o valor da integral. Fisicamente sabemos que essa interpretação está correta, pois P(R) é proporcional a |ϕ_T (R)|², logo teremos uma probabilidade maior de encontrar o sistema em regiões onde o quadrado da função de onda é maior. Assim, usamos P(R) para amostrar um conjunto de configurações {R₁,R₂, . . . ,R_n}. Já a energia local E_L(R_i) recebe esse nome por ter dimensões de energia e depender da configuração R_i onde é calculada. Nesse esquema, uma vez amostrado um conjunto de configurações, podemos estimar o valor da energia variacional como a média do valor de E_L,

(67) $${\overline{E}}_V=\frac{1}{n}{\displaystyle \sum _{i=1}^n}{E}_L(R_i).$$

Quanto mais configurações n amostrarmos, mais o valor da equação (67) se aproxima do valor real da equação (54), com os dois sendo iguais no limite n → ∞. Note que não fizemos nenhuma aproximação na integral original para obter o resultado da equação (67), sendo portando um resultado exato. A única fonte de incertezas existente é o tamanho finito da amostra n, ou seja, temos apenas incerteza de natureza estatística. Tomamos como incerteza dos nossos resultados o erro padrão [¹⁹],

(68) $${\sigma }_V^2=\frac{1}{n(n-1)}{\displaystyle \sum _{i=1}^n(E_L(}{R}_i)-{\overline{E}}_V{)}^2.$$

Com isso, o valor da energia variacional estimado com o método VMC é

(69) $$E_V={\overline{E}}_V\pm {\sigma }_V.$$

Após a apresentação do método, focamos em como amostrar um conjunto de configurações da densidade de probabilidade dada pela equação (60). O método escolhido é o algoritmo de Metropolis [²⁰] e para aplicá-lo precisamos sortear 3N + 1 números aleatórios distribuídos uniformemente no intervalo de 0 a 1 em cada passo. Tipicamente todas linguagens de programação possuem intrinsecamente esse tipo de gerador de números pseudo-aleatórios, ou bibliotecas amplamente divulgadas com essa funcionalidade. Além disso, precisamos de um parâmetro Δ, o qual vai controlar o deslocamento das posições durante a simulação. Podemos resumir o algoritmo de Metropolis como:

1. Comece de uma configuração qualquer R_i (com i = 1).

2. Sorteie 3N números aleatórios de uma distribuição uniforme entre 0 e 1. Subtraia de cada um deles 1/2, de tal forma que eles fiquem distribuídos entre −1/2 e 1/2. Guarde o resultado em um vetor ξ com 3N entradas.

3. Vamos propor uma atualização de cada uma das componentes do vetor R_i em um novo vetor R_novo, através de

   $$R_{\text{novo}}=R_i+\text{Δ}\xi .$$

4. Agora decidimos se aceitamos ou não essa mudança no vetor posição. Para isso, calculamos a função P(R) nas configurações R_i e R_novo e sorteamos um número aleatório η (de uma distribuição uniforme entre 0 e 1).

   • (a)

     Caso P(R_novo)/P(R_i) ⩾ η aceite a mudança, R_i+1 = R_novo.

   • (b)

     Caso P(R_novo)/P(R_i) ＜ η rejeite a mudança, R_i+1 = R_i.

5. Tome i ← i+1. Se o número total de amostras for atingido, pare. Caso contrário, vá para o passo 2.

Uma das grandes vantagens do algoritmo de Metropolis é que não precisamos saber a normalização da função de onda teste, pois ela aparece tanto no numerador quanto no denominador de P(R_novo)/P(R_i). Isso pode parecer trivial nos exemplos mais simples, porém calcular a normalização em um sistema de muitos-corpos interagentes não é uma tarefa fácil.

Se P(R_novo) ⩾ P(R_i), sempre iremos aceitar a mudança. Isso garante que estamos favorecendo as configurações onde |ϕ_T (R_novo)|² é expressivo, e sabemos que o sistema tem uma probabilidade maior de ser encontrado nessa configuração. No entanto, se P(R_novo) ＜ P(R_i) ainda há chance de aceitarmos essa mudança. Isso é necessário para que possamos escapar de máximos locais e assim possamos amostrar o espaço inteiro.

O parâmetro Δ é de nossa escolha. Ele influência o tamanho da mudança proposta para o novo vetor R_novo. Valores de Δ muito pequenos fazem com que R_novo seja muito semelhante a R_i, o que faz com que P(R_novo)/P(R_i) seja próximo de um, ou seja, aceitamos praticamente todas as mudanças. Valores muito grandes de Δ causam um deslocamento abrupto, o que faz com que poucas mudanças sejam aceitas, resultando na repetição da configuração R_i. Em ambos os casos, a correlação entre as configurações geradas é grande. Escolhemos Δ afim de aceitar metade das propostas de mudança, ou seja, mantemos a taxa de aceitação entre 40% e 60%, via tentativa e erro. Como o algoritmo de Metropolis gera a configuração i + 1 a partir da configuração i do sistema, essas duas são correlacionadas. Métodos para descorrelacionar a amostragem devem ser aplicados, e estão amplamente disponíveis na literatura [¹⁷, ¹⁹].

Agora vamos aplicar o método VMC ao átomo de Hélio, com a função de onda teste da equação (61). Vamos fazer simulações independentes para diferentes valores de α, afim de encontrar o valor de α que minimize a energia [¹⁸]. Testando valores no intervalo 0 ＜ α a₀ ＜ 5 espaçados de 0,05 obtemos o resultado da Figura 11. Note a presença de um mínimo de energia, na região 1, 3 ＜ α a₀ ＜ 2,0. Devemos buscar o valor de α que resulta na mínima energia para termos a melhor aproximação possível do estado fundamental com essa função de onda teste. Apresentamos os valores na vizinhança do mínimo na Tabela 5 e na Figura 12.

Figura 11
Energia do estado fundamental do átomo de Hélio como função do parâmetro variacional α obtida com o método VMC para valores no intervalo 0 ＜ α a₀ ＜ 5. Cada energia variacional foi computada a partir de 100 000 amostras obtidas com o método de Metropolis.

Figura 12
Mesmos resultados da Figura 11, mas focando no intervalo 1, 3 ＜ α a₀ ＜ 2 onde podemos ver o mínimo da energia com relação ao parâmetro variacional.

Assim, a nossa estimativa para a energia do estado fundamental do átomo de Hélio é de:

(70) $$E_0^{\text{He}}=(-77,47\pm 0,08)\text{eV}.$$

O valor obtido é maior que a energia experimental do Hélio, −79, 02 eV, como esperávamos do princípio variacional. Embora a função de onda teste utilizada nesse caso seja relativamente simples, equação (61), conseguimos um resultado que difere da ordem de 2% do valor experimental. Caso não conhecêssemos o valor real da energia do Hélio, teríamos um limite superior de −77, 47 eV garantido pelo princípio variacional, ou seja, a energia do estado fundamental só pode ser menor, ou igual, ao valor obtido com o método variacional. Podemos nos aproximar mais da energia experimental melhorando a função de onda teste. Uma possibilidade seria a inclusão de termos que levem em conta a interação entre os elétrons. Note que embora o método variacional nos permita estimar a energia do estado fundamental do sistema, ele não fornece uma medida ou critério de quão próximo dessa estamos.

Nos últimos anos, arquiteturas de processamento paralelo tornaram-se mais comuns e acessíveis, impulsionadas por avanços em processadores com múltiplos núcleos, placas de vídeo e serviços de computação em nuvem. Hardware mais acessível e a crescente popularidade de linguagens adequadas para programação em paralelo contribuíram para que pesquisadores aproveitem o paralelismo para realizar cálculos científicos complexos. O método VMC é considerado trivialmente paralelo [²¹] porque seus cálculos podem ser divididos em tarefas independentes. Cada amostra de Monte Carlo pode ser processada de forma isolada, sem a necessidade de comunicação entre os processos. Isso torna o método ideal para computação paralela, onde várias amostras podem ser calculadas simultaneamente em diversos núcleos, acelerando significativamente as simulações.

7. Conclusão

Nesse trabalho introduzimos quatro abordagens numéricas para encontrar soluções da equação de Schrödinger independente do tempo. Apresentamos cada uma delas a partir dos princípios em que se baseiam, com a finalidade de estudar sistemas físicos. Mostramos que o método da tentativa é uma alternativa de relativamente simples implementação para potenciais simétricos. Para potenciais mais gerais, sem essa simetria, podemos aplicar o método do encontro ou o matricial. O primeiro tem parâmetros de nossa escolha (ponto de encontro e diferença mínima entre as derivadas) que podem influenciar diretamente o resultado. Além disso, em situações onde a região classicamente proibida é expressiva, dificuldades numéricas aparecem devido a integração nessa região levando a divergências, como visto na seção 3.

O método matricial é uma alternativa interessante para resolver o problema na forma de autovalor-autovetor, e pode ser estendido para mais dimensões. Todavia, a busca numérica de autovalores e autovetores é um problema computacionalmente custoso, e o método deve ser aplicado com uma escolha razoável para o tamanho da matriz. Em mais de uma dimensão as soluções da equação de Schrödinger podem ser de difícil obtenção mesmo com métodos numéricos. Logo, o método VMC é uma boa alternativa para estimar a energia do estado fundamental. Duas grandes diferenças entre esse método e os apresentados anteriormente são que trabalhamos no espaço contínuo, sem a necessidade de discretizar o sistema em uma rede, e que estamos restritos ao estado fundamental (embora seja possível obter estados excitados com modificações do método).

Em suma, vimos que não há uma abordagem única para obter soluções numéricas da equação de Schrödinger independente do tempo. A escolha do método mais adequado deve levar em conta características do potencial em questão, número de partículas e dimensão, além do poder computacional disponível. Conhecer os métodos apresentados na literatura, bem como suas vantagens e desvantagens, é imprescindível para resolver problemas de Física Computacional eficientemente.

Agradecimentos

Esse trabalho foi financiado pelos processos n° 2013/07276-1, 2014/50857-8 e 2023/04451-9, Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP), e pelo processo n° 8887.821526/2023-00 da Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES).

A. A equação radial e suas condições de contorno

Na seção 5, abordamos o potencial de Coulomb, o qual é esfericamente simétrico, e utilizamos as condições de contorno para a função de onda radial reduzida: u_nℓ(r → 0) = 0 e u_nℓ(r → ∞) = 0. Nesse apêndice, fornecemos uma explicação mais detalhada dessas escolhas baseando-se nas Refs. [²², ²³].

Sejam duas soluções da equação (19) dadas por $u_{n_1,{\ell }_1}(r)$ e $u_{n_2,{\ell }_2}(r)$ que, por simplicidade, vamos denotar u₁ e u₂. Queremos que o operador Hamiltoniano do sistema seja hermitiano com respeito a essas funções, ou seja:

(A.1) $${\displaystyle \underset{0}{\overset{\infty }{\int }}d}r{u}_1^{*}(H{u}_2)={\left[{\displaystyle \underset{0}{\overset{\infty }{\int }}d}r{u}_2^{*}(H{u}_1)\right]}^{*}={\displaystyle \underset{0}{\overset{\infty }{\int }}d}r{(H{u}_1)}^{*}{u}_2.$$

Escrevendo explicitamente a aplicação do Hamiltoniano nas funções de onda radiais reduzidas, temos:

(A.2) $${\displaystyle {\int }_0^{\infty }d}r\left\{-\frac{{\hslash }^2}{2\mu }{u}_1^{*}{u^{″}}_2+\left[V(r)+\frac{{\hslash }^2\ell (\ell +1)}{2\mu r^2}\right]u_1^{*}{u}_2\right\}={\displaystyle {\int }_0^{\infty }d}r\left\{-\frac{{\hslash }^2}{2\mu }{u}_1^{*\prime \prime }{u}_2+\left[V(r)+\frac{{\hslash }^2\ell (\ell +1)}{2\mu r^2}\right]u_1^{*}{u}_2\right\}.$$

Simplificando,

(A.3) $${\displaystyle \underset{0}{\overset{\infty }{\int }}d}r\left(u_1^{*}{u^{″}}_2-u_1^{*\prime \prime }{u}_2\right)={\displaystyle \underset{0}{\overset{\infty }{\int }}d}r\frac{d}{dr}\left(u_1^{*}{u^{\prime }}_2-u_1^{*\prime }{u}_2\right)=0,$$

de onde temos a condição para que o operador Hamiltoniano do sistema seja hermitiano:

(A.4) $${[u_1^{*}(r){u^{\prime }}_2(r)-u_1^{*\prime }(r)u_2(r)]}_0^{\infty }=0.$$

Se estamos considerando estados ligados, como é o caso desse artigo, então a função de onda não pode se estender até o infinito, de modo que u(r → ∞) = 0, o que anula o limite superior da equação (A.4). Caso estivéssemos lidando com estados de espalhamento por potenciais de alcance finito, u(r → ∞) ∝ e^ikr(com $k=\sqrt{2\mu E}/\hslash$, também zerando o limite superior da equação (A.4).

O limite inferior da equação (A.4) é satisfeito se

(A.5) $$u(r\to 0)=C,$$

com C constante. Como a função de onda radial é dada por R(r) = u(r)/r, ela diverge na origem para C ̸= 0. Isso por si só não é uma restrição, dado que a probabilidade de encontrarmos uma partícula entre r e r+dr é dada por |R(r)|²r²dr = |u(r)|²dr, a qual é finita para r → 0.

O problema de escolhermos C não-nulo é que isso leva a uma divergência da energia cinética na origem. O problema fica evidente se aplicarmos o operador na equação (17) com ℓ = m = 0,

(A.6) $$-\frac{{\hslash }^2}{2\mu }{\nabla }^2\left[\frac{C}{r}{Y}_0^0(\theta ,\varphi )\right]=-\frac{{\hslash }^2}{2\mu }\frac{C}{\sqrt{4\pi }}{\nabla }^2\left(\frac{1}{r}\right).$$

Como [²]

(A.7) $${\nabla }^2\left(\frac{1}{r}\right)=-4\text{π}{\delta }^3(\mathit{r}),$$

onde δ³(r) é a função delta de Dirac, temos uma divergência na origem se C ̸= 0. O único jeito da equação de Schrödinger ser satisfeita na origem nessa caso é se o termo de potencial cancelar exatamente essa divergência [²⁴]. Como para todos outros potenciais esse não é o caso, devemos tomar C = 0 para garantir que o operador H seja hermitiano, ou seja, u(r → 0) = 0.

References

• ^[1] R.M. Eisberg e R. Resnick, Quantum Physics of Atoms, Molecules, Solids, Nuclei, and Particles (Wiley, Nova Iorque, 1985), v. 1.
• ^[2] E. Butkov, Mathematical Physics (Addison-Wesley Publishing Company, Boston, 1968).
• ^[3] P.B. Visscher, Computers in Physics 5, 596 (1991).
• ^[4] N.J. Giordano e H. Nakanishi, Computational Physics (Pearson/Prentice Hall, Upper Saddle River, 2006).
• ^[5] R. Tonhon, Repositório do Github, disponível em:https://github.com/rafaeltonhon/eqschrtimeindep
  » https://github.com/rafaeltonhon/eqschrtimeindep
• ^[6] D.J. Griffiths, Introduction to Quantum Mechanics (Cambridge University Press, Cambridge, 2017).
• ^[7] J.J. Sakurai e J. Napolitano, Modern Quantum Mechanics (Addison-Wesley, Boston, 2011).
• ^[8] S.T. Thornton e J.B. Marion, Classical Dynamics of Particles and Systems (Brooks/Cole, Belmont, 2007).
• ^[9] M. Macêdo-Lima e L. Madeira, Revista Brasileira de Ensino de Física 45, e20230079 (2023).
• ^[10] W.H. Press, B.P. Flannery, S.A. Teukolsky e W.T. Vetterling, Numerical Recipes in FORTRAN 77: The Art of Scientific Computing (Cambridge University Press, Cambridge, 1992), v. 1.
• ^[11] M.E.J. Newman, Computational Physics (CreateSpace Independent Publishing Platform, Washington DC, 2013).
• ^[12] L.B. Carvalho, W.C. Santos e E.A. Correa, Revista Brasileira de Ensino de Física 41, e20180359 (2019).
• ^[13] F. Caruso e V. Oguri, Revista Brasileira de Ensino de Física 36, 2310 (2014).
• ^[14] E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Dongarra, J.D. Croz, A. Greenbaum, S. Hammarling, A. McKenney e D. Sorensen, LAPACK Users’ Guide (Society for Industrial and Applied Mathematics, Philadelphia, 1999).
• ^[15] LAPACK: dstemr, disponível em:https://www.netlib.org/lapack/explore-html/d4/dec/group__stemr_ga71cbf49a0387762afa39582e6abbe466.html#ga71cbf49a0387762afa39582e6abbe466, acessado em: 01/08/2024.
  » https://www.netlib.org/lapack/explore-html/d4/dec/group__stemr_ga71cbf49a0387762afa39582e6abbe466.html#ga71cbf49a0387762afa39582e6abbe466
• ^[16] scipy.linalg.eigh_tridiagonal - SciPy v1.13.0 Manual, disponível em:https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.eigh_tridiagonal.html, acessado em: 01/08/2024.
  » https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.eigh_tridiagonal.html
• ^[17] R. Canabrava e S.A. Vitiello, Revista Brasileira de Ensino de Física 25, 35 (2003).
• ^[18] M.A. Reis e S.A. Vitiello, Revista Brasileira de Ensino de Física 28, 45 (2006).
• ^[19] I.G. Hughes e T.P.A. Hase, Measurements and their Uncertainties: A practical guide to modern error analysis (Oxford University Press, Oxford, 2010).
• ^[20] N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller e TellerE., The Journal of Chemical Physics 21, 1087 (1953).
• ^[21] W.M.C. Foulkes, L. Mitas, R.J. Needs e G. Rajagopal, Reviews of Modern Physics 73, 33 (2001).
• ^[22] R. Shankar, Principles of Quantum Mechanics (Springer, Nova, Iorque, 2012).
• ^[23] E. Merzbacher, Quantum Mechanics (Wiley, Nova, Iorque, 1998).
• ^[24] F.A.B. Coutinho e M. Amaku, European Journal of Physics 30, 1015 (2009).

• 1
  Valores consultados emhttps://physics.nist.gov/cuu/Constants/ em 31/07/2024.

Datas de Publicação

• Publicação nesta coleção
  22 Nov 2024
• Data do Fascículo
  2024

Histórico

• Recebido
  13 Ago 2024
• Revisado
  02 Out 2024
• Aceito
  04 Out 2024