468
Revista Brasileira de Ensino de F
sica, vol. 22, no. 4, Dezembro, 2000
Simulac~ao Graca do Comportamento
Qu^antico da Cadeia Monoat^omica
(Graphical simulation of the quantum behavior of the atomic chain)
Luciano Terra Peixoto e Newton Eloi Oliveira de Azevedo
Departamento de F
sica, Universidade Federal do Esp
rito Santo
29060-900 Vit
oria, ES, Brasil
Recebido em 20 de Abril, 2000. Aceito em 14 de Julho, 2000.
No presente trabalho exploramos o uso da computac~ao graca para simulac~ao do comportamento
qu^antico de uma cadeia at^omica. Preliminarmente, completando trabalho anterior, desenvolvemos
um programa em C++ para visualizac~ao do comportamento classico de uma cadeia diat^omica, tanto
nos modos acusticos quanto nos opticos. Trabalhamos ent~ao com tratamento qu^antico, em relac~ao
inicialmente a um unico oscilador harm^onico, e depois a um conjunto de osciladores acoplados
analogo a cadeia at^omica na aproximac~ao harm^onica. Uma vez superado o problema de estabilidade
que tais algortmos apresentam, o objetivo de simulac~ao em relac~ao ao oscilador unico foi atingido
sem restric~oes. No caso dos osciladores acoplados, a necessidade de trabalhar com func~oes de tantas
variaveis quanto o numero de osciladores nos leva a um algortimo com multiplos loops aninhados e
a exig^encia de memoria e de velocidade de processamento incompatveis, nos computadores atuais,
com a visualizaca~o do comportamento do sistema atraves da animac~ao de imagens. Recorremos
ent~ao a transformac~ao de coordenadas para modos normais do sistema, o que evita a diculdade
apontada, mas com perda parcial da ideia de simulac~ao.
In this work we explore the use of computer graphics to simulate the quantum behavior of an
atomic chain. Preliminarly, concluding a previous work, we developed a C++ program to visualise
the classical behavior of a diatomic chain, in order to observe this behavior for acoustical and optical
modes. Then we turned to the quantum approach, rst for a single harmonic oscilator, and then
for a set of coupled oscilators analog to the atomic chain in the harmonic approximation. Once
surmounted the stability problem posed by such algorithms, the simulation for the single oscilator
have been achieved without any restriction. In the case of the coupled oscilators, the need to work
with functions of so many variables as the number of atoms in the chain, lead us to an algorithm
with multiple nested loops and to memory and speed requirements incompatibles in the existing
computers with the visualization of the behavior of the system trough images animation. We then
resorted to the coordinate transformation to normal coordinates of the system, a procedure which
avoids the problem but with partial loss of the idea of simulation.
I
Introdu
c~
ao
Em um trabalho anterior descrevemos o desenvolvimento de um programa de computador, designado \Simulac~ao de um Laboratorio para estudo do Comportamento da Cadeia At^omica", servindo para visualizac~ao
do comportamento de um conjunto de osciladores acoplados na abordagem da fsica classica (PEIXOTO e
GOMES, 1998). O termo simulac~ao no caso signica
que o programa se encarrega de encontrar a posic~ao e a
velocidade de cada atomo em cada instante por soluca~o
numerica da equac~ao de Newton na forma diferencial, a
partir de posic~oes e velocidades escolhidas para um ins Estudante
de Fsica, bolsista de Iniciac~ao Cientca (CNPq)
tante inicial. Neste sentido, tivemos oportunidade de
observar peculiaridades do comportamento do sistema
ate ent~ao insuspeitas. Se tais peculiaridades correspondem ou n~ao ao comportamento real da cadeia e uma
quest~ao em aberto.
No presente trabalho descrevemos o complemento
daquele programa, primeiramente nos referindo ao
comportamento da cadeia diat^omica na aproximac~ao
harm^onica. A cadeia foi representada por uma
sequ^encia horizontal de crculos desenhados na tela
do computador, capazes de oscilar verticalmente com
a passagem de excitac~oes transversais, ou horizontalmente com a passagem de excitac~oes longitudinais. A
469
Luciano Terra Peixoto e Newton El
oi Oliveira de Azevedo
cadeia at^omica e tomada como modelo unidimensional
de um solido, para observac~ao e estudo de seu comportamento vibracional. Reescrevemos aqui a equac~ao
de Newton e a express~ao de denic~ao da velocidade na
forma discreta, tal como utilizadas naquele trabalho:
vp t +
Æ
Æ
2 = vp t 2 + ! Æ[yp (t)
+yp (t) 2yp (t)]
(1 a)
2
0
+1
1
yp (t) = yp (t Æ) + Ævp t
Æ
(1 b)
2 :
A velocidade e a posic~ao do p-esimo atomo em um instante s~ao aqui obtidas da velocidade e da posic~ao do
mesmo atomo e de atomos vizinhos em instantes anteriores; Æ representa o incremento de tempo e ! e
a frequ^encia natural de oscilac~ao, dada pela raz quadrada da raz~ao k=m, sendo k a constante de mola e
m a massa do atomo. A equac~ao de Newton acima
esta adaptada para a situac~ao em que cada atomo interage apenas com os dois vizinhos mais proximos (ASHCROFT, 1976). Notar que as velocidades e as posic~oes
s~ao consideradas em instantes alternados, o que constitui um recurso para tornar o algoritmo rapido o suciente e livre de instabilidades (HOCKNEY, 1970, PEIXOTO e GOMES, 1998).
No caso da cadeia diat^omica o par de equac~oes
acima se desdobra em dois pares, um valido para a velocidade vp e a coordenada yp de um tipo de atomo,
envolvendo as coordenadas yp0 e yp0 dos atomos vizinhos que s~ao do segundo tipo, e o outro par valido
para a velocidade vp0 e a coordenada yp0 do segundo tipo
de atomo, envolvendo as coordenadas yp e yp dos
atomos vizinhos que s~ao do primeiro tipo. Os dois pares de equac~oes tambem se diferenciam pelo valor de ! ,
isso porque os dois tipos de atomos tem massas diferentes. Tudo que e necessario fazer e duplicar o programa
original, com escolhas apropriadas de ! e com o necessario cuidado no trato das condic~oes iniciais. Elas
devem corresponder a presenca em t=0 de modos normais acusticos ou modos normais oticos.
0
+1
1
+1
1
0
0
II
O tratamento qu^
antico do
oscilador harm^
onico
A equac~ao de Schrodinger dependente do tempo para
um oscilador harm^onico de massa m e constante elastica
k, com afastamento em relac~ao a posic~ao de equilbrio
designado por y, e
@ (y; t)
=
i~
@t
1
2m @y + 2 ky (y; t)
~2 @ 2
2
2
(2)
Aqui as variaveis y e t tornam-se discretas, e a equac~ao
diferencial torna-se uma equac~ao de diferencas nitas,
i~
Æ
y; t + delta
2
y; t 2Æ
~
= 2ma
[(y+a; t)+(y a; t) 2(y; t)]+ 12 ky (y; t);
(3)
onde Æ e o incremento no tempo e a o incremento em
y, tal que ym = ma. Esta equac~ao se desdobra em
duas quando consideramos separadamente a parte real
R(y; t) e a parte imaginaria I (y; t) da func~ao de onda
(y; t): Obtemos ent~ao ,
2
2
2
R+m = Rm + c[Im+1 + Im
I +m = Im c[Rm+1 + Rm
1
2Im ] dm Im ; (4 a)
2
2Rm]+ dm Rm ; (4 b)
onde introduzimos ndices superiores + e - para representar o valor da func~ao em t + Æ=2 e t Æ=2, respectivamente, de forma que sem o ndice o valor da
func~ao se refere ao instante t. Tambem introduzimos os
par^ametros c = ~Æ=2ma e d = kÆa =2~, de modo que
a escolha de valores de k, m, Æ e a se reduzem a escolha
de valores de c e d. As express~oes (4) s~ao aplicadas dentro da condica~o de contorno periodica, tal que, sendo
N o numero de atomos da cadeia, ent~ao N = e
m = N para m = 1. Calculadas as partes real e
imaginaria da func~ao de onda, o que vamos representar
na tela do computador sera, instante por instante, a
curva da distribuic~ao de probabilidade, dada por
2
1
2
2
+1
1
1
2
Pm = Rm
+ Im2 :
(5)
Podemos fazer o programa desenhar pontos na tela do
computador, cujos afastamentos verticais em relac~ao
a uma linha horizontal v~ao representar probabilidades
n~ao nulas do atomo ocupar as posic~oes horizontais correspondentes.
Na implementac~ao do calculo das express~oes (4-a) e
(4-b), e desejavel um algortmo sucientemente rapido
para que, ao se imprimir toda a sequ^encia de crculos
na tela do computador em ciclos sucessivos do processamento, se tenha impress~ao de continuidade de movimento da curva da distribuic~ao de probabilidade. Neste
sentido e apropriada a escolha dos chamados
, para os quais a variavel dependente (Rm ou
Im ) so aparece no lado direito da express~ao quando
valida para um instante anterior. Tais metodos apresentam entretanto o problema da instabilidade, por causa da magnitude do erro numerico - da ordem do
proprio incremento do tempo - a que est~ao sujeitos.
Metodos implcitos, por outro lado, podem apresentar
erro de segunda ordem no incremento do tempo e s~ao
algor
tmos
expl
citos
470
Revista Brasileira de Ensino de F
sica, vol. 22, no. 4, Dezembro, 2000
estaveis, mas s~ao lentos, porque um sistema de equaco~es
algebricas lineares tem que ser resolvida em cada ciclo de processamento (BUTCHER, 1994). E possvel
contornar o problema utilizando um algortmo explcito
em que as partes real e imaginaria da func~ao de onda
s~ao dadas em instantes de tempo alternados (HOCKNEY, 1970, VISSCHER, 1991). Ele esta discutido no
Ap^endice.
Neste algortmo os valores de R na express~ao (4)
s~ao denidos, por exemplo, nos instantes 0, Æ, 2Æ..., e
os de I nos instantes Æ=2, 3Æ=2, 5Æ=2... Para ocorrer a
conservac~ao de probabilidade, e necessario reescrever a
express~ao (5) na forma
2
Pm = Rm
+Im+ Im
(para P e R nos instantes 0, Æ, 2Æ...)
(6 a)
Pm = Rm + Rm + Im (para P e I nos instantes
Æ/2, 3Æ/2, 5Æ/2...):
(6 b)
Em (6-a), para valores de P e R tomados no instante
t = Æ, entram os valores de I e I nos instantes 3Æ=2 e
Æ=2, respectivamente, e em (6-b), para os valores de P
e I tomados no instante t = Æ=2, entram os valores de
R e R nos instantes t = Æ e t = 0; respectivamente.
E facil vericar que (5) deve ser substituda pelas duas
express~oes acima. Basta tomar uma depend^encia temporal na forma exp(iwt) em e substituir a parte real
R por cos(wt) e a parte imaginaria I por sen(wt) e,
tendo em conta os instantes de denic~ao dos varios termos, vericar que o cumprimento da condic~ao
X
Pm = 1
(7)
+
2
+
+
m
ca garantida instante por instante.
Devemos escolher valores iniciais para os varios Rm
e Im tais que representem a parte real no instante t = 0
e a parte imaginaria no instante t = Æ=2. Trabalhamos
com curvas gaussianas, dadas por
p 2
Rm (0) = Ae ym= cos(!0)
(8 a)
Æ
(
2 )
p2
ym = 2 sen ! Æ
(8 b)
2 = Ae
2
onde a constante A e escolhida em func~ao do tamanho
da gura na tela do computador. As duas curvas entram como R e I em (4), para o calculo de R e
I nos instantes t = Æ e t = 3Æ=2, respectivamente, os
quais, por sua vez, retornam como R e I em (4) para
o calculo de R e de I nos instantes t = 2Æ e t = 5Æ=2;
respectivamente, e assim por diante.
Im
(
)
+
+
III
Osciladores acoplados
Consideremos aqui N osciladores de massa m, situados
sobre o eixo x, acoplados por molas id^enticas, de constante elastica k, separados pela dist^ancia b. A posic~ao
de equilbrio do n-esimo oscilador sera dada por x = nb:
Cada oscilador movimenta-se na direc~ao y, na nossa
descrica~o dando saltos entre pontos vizinhos sobre uma
linha de M pontos separadas pela dist^ancia a. Ou seja,
representando o deslocamento do n-esimo oscilador em
relac~ao a respectiva posic~ao de equilbrio por yn, cada
variavel yn tera M possveis valores dados por yn = ma.
A func~ao de onda (y ; y ; :::yN ), fornece a probabilidade de ocorrer um dado conjunto de valores para as
coordenadas y , y ,...yN :
1
1
2
2
P (y1 ; y2 ; :::yN ) = j(y1; y2 ; :::yN )j2 :
(9)
A func~ao de N partculas (y ; y ; :::yN ) n~ao e separavel em produto de func~oes de onda de uma
partcula porque os osciladores s~ao acoplados. Para
obter a curva de distribuic~ao de probabilidade para o
i-esimo oscilador e necessario calcular a express~ao
1
P (yi ) =
XX X
y1 y2
:::
yN
2
P (y1 ; 2y2; :::yN )
(10)
em que ca excludo o somatorio sobre yi . Isso signica que o algortmo de simulac~ao devera ter um total
de N 1 loops aninhados. Sendo N e M numeros
da ordem de dezenas, isso implicara em velocidade de
processamento incompatvel com a intenca~o de visualizac~ao de movimento contnuo da curva de distribuic~ao
de probabilidade para cada oscilador. Alem disso, poderemos enfrentar o problema de limitaca~o de memoria
do computador, pois precisaremos estocar M N valores
para cada variavel P , R e I (cada uma delas com N
ndices, cada ndice podendo tomar M valores).
Sem essas limitac~oes, para observar um modo normal bastaria introduzir no programa a condic~ao inicial
correspondente a sua excitac~ao, e observar a resposta
do sistema. Com elas, a observac~ao do modo normal so
sera possvel se levarmos em conta que basta fazer a simulac~ao para um unico oscilador e aproveitar o fato de
que a curva que representa a distribuic~ao de probabilidade em dado instante para este oscilador vale tambem
para os demais osciladores, apenas introduzindo-se a
relac~ao de fase apropriada ao comprimento de onda do
modo normal que se deseja observar. O preco a pagar e
a perda parcial da ideia de simulac~ao, porquanto estaremos introduzindo no programa informac~oes referentes
ao conhecimento previo do comportamento do sistema.
A transformac~ao das coordenadas cartesianas yn
para as coordenadas normais um e a transformac~ao inversa podem ser escritas na forma (ASHCROFT, 1976)
um =
r X
2 N
N
n=0
yn senqm xn
(11 a)
Luciano Terra Peixoto e Newton El
oi Oliveira de Azevedo
yn =
r X
2 N
N
m=0
um senqm xn
471
(11 b)
onde xn = nb: Os vetores de onda qm podem ser determinados de modo a satisfazer a condic~ao de contorno periodica, tal que yn N = yn , o que implica em
qm = (2=Nb)(m N=2), com m variando entre 0 e
N para que os valores de q quem compreendidos entre =b e =b (primeira zona de Brillouin). A express~ao (11-b) nos diz que o afastamento do n-esimo
oscilador em relac~ao a respectiva posic~ao de equilbrio
e obtido pela superposic~ao dos afastamentos provocados pelos varios modos normais, cada um deles uma
onda senoidal de amplitude um e comprimento de onda
m = Nb=(m N=2).
+
IV
Figura 1. Aspecto da tela do computador em um estagio do
processamento do programa para visualizac~ao do comportamento de uma cadeia diat^omica.
Resultados e Discuss~
oes
Os programas foram codicados com a linguagem C
(Press
, 1992). A Fig. 1 mostra o aspecto da
tela do computador em um estagio da execuca~o do
codigo para simulac~ao do comportamento de modos
normais acustico e optico de uma cadeia diat^omica
com 60 atomos. Os dois tipos de atomos s~ao indicados por crculos em preto e em branco. Ambos os
modos mostrados na gura tem comprimento de onda
igual a metade do comprimento da cadeia, cando evidente a oscilac~ao de dois atomos vizinhos com defasagem de 180Æ no caso do modo optico e sem defasagem no caso do modo acustico. E sugerido ao aluno
medir as duas frequ^encias com o cron^ometro de seu
relogio de pulso (contando 10 ou mais oscilac~oes) e vericar, a partir da raz~ao entre as massas dos atomos
informada na tela do computador, que a raz~ao entre as
frequ^encias satisfaz conhecida express~ao valida para a
cadeia diat^omica no limite de grandes comprimentos de
onda (ver, por exemplo, ASHCROFT, 1976). Podemos
dizer que aqui os objetivos da simulac~ao s~ao atingidos
sem restric~oes, sendo possvel simular modos normais de
diferentes comprimentos de onda, em cada caso para os
dois modos de oscilac~ao da rede.
A Fig. 2 mostra o aspecto da tela em um instante durante a execuc~ao do programa para simulac~ao
do comportamento qu^antico de um unico oscilador
harm^onico. Os crculos na tela do computador indicam
pontos da curva P (x), representando a probabilidade
do oscilador ser encontrado na posic~ao x. Novamente
podemos dizer que os objetivos da simulac~ao s~ao atingidos sem restric~oes, percebendo-se o que o pulso que
representa P (x) se desloca para a esquerda e para a direita, com velocidade maxima no centro do espaco de
oscilac~ao.
++
et al
Figura 2. Aspecto da tela do computador em um estagio do
processamento do programa para visualizac~ao do comportamento qu^antico de um oscilador harm^onico.
Na Fig. 3 mostramos o aspecto da tela do computador em um instante da simulac~ao do programa para o
modo normal com comprimento de onda igual ao comprimento da cadeia. A gura, no caso, tem um carater
tridimensional. A curva de distribuica~o de probabilidade, que chamamos de P (x), e representada por um
pulso que oscila horizontalmente na tela do computador. Ela corresponde a um atomo que oscila na mesma
direc~ao. Diferentes atomos est~ao dispostos segundo a
direc~ao y. O modo normal se propaga na direca~o y,
provocando o deslocamento transversal dos pulsos de
probabilidade de uma forma que lembra o serpentear
de uma cobra. Aqui se deve ter em conta que a simulac~ao, pelo menos no sentido que temos empregado
neste trabalho, so ocorre para um unico oscilador, sendo
a posic~ao do pulso para cada um dos demais osciladores
determinada de modo a se cumprir a relac~ao de fases
correspondentes ao modo normal com = comprimento
da cadeia.
472
Revista Brasileira de Ensino de F
sica, vol. 22, no. 4, Dezembro, 2000
Uma vers~ao experimental do programa para simulac~ao pode ser solicitado pelo endereco eletr^onico
[email protected].
obter a equac~ao de evoluc~ao temporal, consideremos a
representaca~o matricial correspondente para a equac~ao
(A.3)
R
I+
1 0
A 1
=
e substituindo (A.5) em (A.4):
R
I
(A 5)
R
= 1 AA A1
: (A 6)
I
Agora temos uma equac~ao de evoluc~ao temporal, em
que a matriz de evoluc~ao tem autovalores
R+
I+
2
=
(A 7)
1 11 A A 14 A 1
A matriz sera unitaria se os autovalores tiverem modulo
1, o que vai acontecer se jAj 2, e esta sera a condic~ao
de estabilidade para o sistema.
Consideremos a atuac~ao de H sobre a func~ao de
onda do estado fundamental do oscilador harm^onico,
de forma gaussiana, que corresponde a escolha que zemos anteriormente para representar o estado inicial do
sistema, dada por (LEIGHTON, 1959):
2
Figura 3. Figura 3. Aspecto do modo normal correspondente a = comprimento da cadeia.
Ap^endice
Um dos problemas com o presente algortmo e o de
estabilidade, isto e, o de se conseguir impedir que os
crculos na tela do computador passem a oscilar com
amplitudes cada vez maiores ate desaparecerem e isso
provocar a interrupc~ao do programa. Ja discutimos o
criterio de estabilidade para o tratamento classico da
cadeia linear em um trabalho anterior (PEIXOTO E
GOMES, 1998). Vamos aqui discut-lo para o tratamento qu^antico. Basicamente adaptamos os argumentos de Visscher (VISSCHER, 1991), validos para o comportamento qu^antico de uma partcula livre. Discretizando a equac~ao de Schrodinger,
@ (y; t)
i~
= H (ymt)
(A 1)
@t
separando a parte real da parte imaginaria e utilizando
a notac~ao adotada acima, obtemos
+ HÆ
I (t = 0; Æ; 2Æ; :::)
~
R+ = R
I+ = I
HÆ
(A 2)
(t = Æ=2; 3Æ=2; 5Æ=2; :::) (A 3)
Uma representac~ao matricial para (A.2), com A =
HÆ=~, e,
~
R
R
= 10 A1
(A 4)
I
que n~ao e uma equac~ao de evoluc~ao temporal, pois do
lado esquerdo temos a parte real da func~ao de onda
considerada em um instante posterior ao da parte imaginaria, enquanto o reverso ocorre no lado direito. Para
R+
I+
+
1 2
2
(y) = B exp[ y =2 ]:
2
(A:8)
2
onde B e uma constante de normalizac~ao e = mw=h:
Aplicando o operador H obtemos, como na express~ao
(3),
~
H
(y; t) = 2ma ((y + a; t)
~
2
2
a2
+(y a; t) 2(y; t)) + 2~ km (y; t)
2
(A 9)
onde, como antes, y = ma: Substituindo (A.8) em (A.9)
e multiplicando ambos os membros pelo incremento de
tempo Æ, encontramos,
HÆ
~
(y; t) = fc[1 e
a2 =22
cos(a m= )
2
2
+a m = ]g(y; t)
onde c e aqui a mesma constante introduzida na sec~ao
II. O termo entre chaves e aquele acima chamado de A,
e assim, levando em conta que o valor maximo de m e
M e representando a largura do pacote gaussiano por
= m a, a condic~ao de estabilidade ca
4
2
2
0
cj1 e
cos(M=m ) + (M=m ) j 2
p
Por exemplo, para M = 30 e mo = 3, encontramos
c 0:02 como condic~ao para o sistema car estavel.
1=2m2
0
2
0
2 2
0
Luciano Terra Peixoto e Newton El
oi Oliveira de Azevedo
References
[1] L. T. Peixoto e K.Q. Gomes, Rev. Bras. Ens. Fs., 20
111 (1998).
[2] N. W. Ashcroft e N. D. Mermin, Solid
Philadelphia, Holt-Saunders, 1976.
State Physics.
[3] J. C. Butcher, Comp. Phys. 8(4), 411 (1994).
473
[4] B. P. Flannery, W. H. Press, S. A. Teukolsky, W. T.
Vetterling, Numerical Recipes in C. Cambridge, Cambridge University Press, 1992.
[5] P. B. Visscher, Comp. Phys. nov/dec, 596, 1991.
[6] R. W. Hockney; Methods Comput. Phys., 9, 136 (1970).
[7] Robert B. Leighton, Principles of Modern Physics. New
York, McGraw-Hill, 1959.