O problema Fermi-Pasta-Ulam-Tsingou: Equiparticão de energia vista através de simulações computacionais

The Fermi-Pasta-Ulam-Tsingou problem: Equipartition of energy through computational simulations

Resumos

Equipartição e conservação de energia são conceitos difíceis de serem ilustrados em sala de aula. Com essa finalidade, implementamos dois experimentos computacionais do problema Fermi-Pasta-Ulam-Tsigou. Nos experimentos, há uma situação na qual acontece a equipartição de energia e outro em que ela não se realiza. A segunda situação levou ao conhecido paradoxo de Fermi-Pasta-Ulam-Tsigou com grandes implicações para o desenvolvimento da Física. Os experimentos são realizados através de uma implementação numérica “feita em casa” e outra através de uma ferramenta computacional disponível na internet e bastante utilizada pela comunidade científica. As implementações são acompanhadas de uma análise matemática, utilizando-se modos normais.

Palavras-chave:
Equipartição de energia; paradoxo de Fermi-Pasta-Ulam-Tsingou; simulação numérica; dinâmica molecular

Equipartion and conservation of energy are hard concepts to teach or exemplify. With this goal in mind, we implemented two computational experiments of Fermi-Pasta-Ulam-Tsigou. There is an experiment in which the equipartion happens and the other one not. The latter experiment leads to known Fermi-Pasta-Ulam-Tsingou Paradox with great implications to development of Physics. The experiments are carried out through a home-made numerical implementation and a computational research tool available on internet and usually adopted by scientific community. There is a mathematical analysis, using normal modes, to support the computational experiments.

Keywords
Equipartition of energy; Fermi-Pasta-Ulam-Tsigou paradox; numerical simulation; molecular dynamics

1. Introdução

O paradoxo de FPU (Enrico Fermi, John R. Pasta, Stanislaw M. Ulam) resulta de uma análise computacional de um problema hipotético, em que a suposição não foi confirmada pelos resultados [¹]. A análise do problema gerou um paradoxo que começaria a ser respondido somente 10 anos depois, o que ajudou no desenvolvimento das teorias de sólitons e do caos [²]. O intuito original era estudar como um sistema evolui para seu equilíbrio térmico. Para isso, decidiram simular uma corrente de partículas que interagiam entre seus vizinhos apenas, sendo os extremos ligados à paredes fixas. As interações seriam forças elásticas, ou seja, dependentes do deslocamento relativo a seus centros de equilíbrio. As forças teriam um termo linear como a Força de Hooke e mais um termo não-linear, podendo ser quadrático ou cúbico.

Quando o caso foi estudado pela primeira vez em Los Alamos, Estados Unidos, além dos três participantes que relataram o caso, dando origem ao nome do paradoxo, em 1955, tivemos a participação de mais uma pessoa, Mary Tsigou, que ajudou no estudo numérico [³]. Ela implementou o código que seria utilizado. Seu nome não é mencionado no relatório da época, mas, nos padrões atuais para coautoria, ela seria mencionada no artigo que poderia ser escrito. Sendo assim, atualmente, o paradoxo é denominado pela sigla FPUT (Fermi-Pasta-Ulam-Tsigou). Os resultados do relatório nunca vieram a ser publicados em um formato de artigo devido a morte prematura de Fermi, que gerou dúvidas de como ele, sendo um dos autores, “assinaria” o artigo.

A premissa inicial do paradoxo de FPUT consiste no Teorema da Equipartição de Energia. Assim, era esperado que a energia total do sistema fosse distribuída igualmente entre suas várias partículas. No caso em questão, a distribuição de energia entre as partículas pode ser descrita através dos modos normais de vibração das mesmas. Pretendia-se observar a distribuição uniforme de energia entre os diversos modos normais de vibração com o passar do tempo. Isso significaria que o sistema alcançou um equilíbrio térmico e seria uma exemplificação computacional do Teorema de Equipartição de Energia. Caso as forças entre as partículas fossem estritamente lineares, isso não ocorreria, pois a energia alocada em cada modo não conseguiria acessar outros modos. Imaginava-se que uma componente não-linear na força tornaria acessível qualquer modo de vibração, porém não foi o observado.

A princípio, foi observada a tendência do sistema de distribuir a energia. O primeiro modo de vibração antes estimulado, perdeu energia ao longo do tempo, a qual começou a se alocar nos modos de energia mais baixos. Entretanto, por um descuido, deixaram a simulação decorrer por um tempo maior do que era planejado. Ao retornar ao laboratório para corrigir tal erro, se depararam com um resultado inesperado. A energia, que supostamente deveria estar igualmente partilhada entre os modos de vibração, estava quase completamente alocada no primeiro modo de vibração. De fato, somente 3% da energia não estava presente no primeiro modo. Devido a esta observação, deixaram a simulação correr por ainda mais tempo. Notaram, então, que existia um ciclo, no qual a energia saía do primeiro modo de vibração, começava a se distribuir nos modos mais baixos, para, por fim, voltar quase que inteiramente para o primeiro modo de vibração. Contudo, em 2015, relatou-se que o sistema FPUT poderia termalizar – atingir equipartição de energia – pelo menos entre modos normais livres (interação entre três fônons) [⁴].

O presente artigo ilustra, portanto, a equipartição e a conservação de energia, tomando como ponto de partida este paradoxo histórico para a área de física computacional. Acreditamos que os conceitos envolvidos podem ser ricamente explorados com este exemplo através de simulações computacionais e análises de dados extraídos das mesmas, tanto em salas de aula do Ensino Médio como Superior. Primeiramente, desenvolveremos detalhadamente uma explicação sobre o problema físico e sobre modos normais, de maneira a possibilitar a compreensão dos aspectos gerais do problema. Em seguida, exploraremos o Teorema da Equipartição de Energia, o qual nos fornece a expectativa de como o sistema deve se comportar. Por fim, resolveremos o problema através da resolução das equações de Newton, por meio de dois métodos computacionais distintos.

No primeiro método, implementaremos um código em Python, como seria abordado em uma disciplina padrão de “Física Computacional”. Esse método será similar ao empregado por FPUT. O outro procedimento envolverá a resolução do problema com a ajuda de um pacote de simulação de Dinâmica Molecular, que está disponível abertamente na internet. Simulação por Dinâmica Molecular, por sua vez, é um método muito utilizado atualmente para se estudar desde sistemas de física da matéria condensada, até sistemas complexos, como proteínas em biofísica, por exemplo. A explicação teórica e as descrições das metodologias computacionais permitem que professores tenham bom embasamento para abordarem os conceitos e, além disso, demostrem como recursos computacionais podem ser agregados à análise científica. Toda a abordagem utilizada neste trabalho pode ser realizada em computadores pessoais, como os que são disponibilizados em laboratórios de informática escolares, e está detalhadamente descrita para que professores possam usá-la em sala de aula.

2. Descrição Formal do Problema

Temos, então, um sistema unidimensional com N partículas com forças entre seus vizinhos com termos não-lineares. Esses termos podem ser considerados ora quadráticos, ora cúbicos. Sendo x_i o deslocamento da i-ésima partícula em relação a seu centro de equilíbrio, podemos representar o problema pela Figura 1 e, por consequência, escrever a força de interação linear combinada com a interação quadrática como

Figura 1
Esquema do problema FPUT para forças de interação linear. As partículas nos extremos estão fixas às paredes.

(1) $$F_i=m\ddot{x_i}=k_1\{(x_{i+1}-x_i)+(x_{i-1}-x_i)\}+k_2\{{(x_{i+1}-x_i)}^2+{(x_{i-1}-x_i)}^2\},$$

onde x_i±1−x_i = Δx_i representa a distância entre duas partículas consecutivas e $\ddot{x_i}$ é a aceleração da i-ésima partícula numa representação de Newton para a derivada segunda da posição em relação ao tempo. Neste problema, todas as partículas possuem a mesma massa m. Analogamente, podemos escrever o caso de interações lineares e cúbicas na forma

(2) $$F_i=m\ddot{x_i}=k_1\{(x_{i+1}-x_i)+(x_{i-1}-x_i)\}+k_2\{{(x_{i+1}-x_i)}^3+{(x_{i-1}-x_i)}^3\}.$$

Para obtermos as acelerações de cada partícula e conseguirmos a descrição cinética do movimento, podemos considerar a razão da constante elástica do termo linear pela massa como a unidade e definir coeficientes α e β para os termos quadráticos e cúbicos. Esses novos coeficientes definem o quanto esses termos afetam a dinâmica do sistema. Sendo assim, temos as seguintes equações:

(3) $$\ddot{x_i}=(x_{i+1}-x_i)+(x_{i-1}-x_i)+\alpha \{{(x_{i+1}-x_i)}^2+{(x_{i-1}-x_i)}^2\}$$

(4) $$\ddot{x_i}=(x_{i+1}-x_i)+(x_{i-1}-x_i)+\beta \{{(x_{i+1}-x_i)}^3+{(x_{i-1}-x_i)}^3\}$$

Então, considerando-se uma força conservativa, podemos escrevê-la como $F_i=\frac{-\partial U}{\partial x_i}$, em que U é a energia potencial. Para o caso quadrático, a energia potencial será

(5) $$U={\displaystyle \sum _{i=0}^N}{\displaystyle \frac{1}{2}}\{{(x_{i+1}-x_i)}^2+{(x_{i-1}-x_i)}^2\}+{\displaystyle \frac{\alpha }{3}}\{{(x_{i+1}-x_i)}^3+{(x_{i-1}-x_i)}^3\}$$

e, para o caso cúbico,

(6) $$U={\displaystyle \sum _{i=0}^N}{\displaystyle \frac{1}{2}}\{{(x_{i+1}-x_i)}^2+{(x_{i-1}-x_i)}^2\}+{\displaystyle \frac{\beta }{4}}\{{(x_{i+1}-x_i)}^4+{(x_{i-1}-x_i)}^4\}.$$

Supondo que x₀ e x_{N + 1} são as paredes do sistema, a Hamiltoniana, dependendo de cada caso, pode ser escrita como

(7) $$H={\displaystyle \sum _{i=0}^N}{\displaystyle \frac{p_i^2}{2m}}+{\displaystyle \frac{1}{2}}{(x_{i+1}-x_i)}^2+{\displaystyle \frac{\alpha }{3}}{(x_{i+1}-x_i)}^3$$

e

(8) $$H={\displaystyle \sum _{i=0}^N}{\displaystyle \frac{p_i^2}{2m}}+{\displaystyle \frac{1}{2}}{(x_{i+1}-x_i)}^2+{\displaystyle \frac{\beta }{4}}{(x_{i+1}-x_i)}^4,$$

para os casos quadrático e cúbico, respectivamente.

Se possuímos um sistema que realiza pequenas oscilações acopladas, podemos aproximar uma solução particular em que todas as partículas realizam movimentos harmônicos simples com a mesma frequência. A solução mais geral possível seria a soma dessas soluções particulares, chamadas de modos normais de vibração. Portanto, em princípio, sabendo as frequências de vibração e o deslocamento dos modos normais obteríamos a solução do problema.

3. Modos Normais

Desejando obter a solução do problema, é preciso revisar os conceitos de modos normais [⁵, ⁶]. Suponhamos um sistema idêntico ao do paradoxo FPUT, exceto que, agora, iremos considerar apenas termos lineares na força. A equação de movimento de cada partícula então se resume a

(9) $$m\ddot{x_i}=k(x_{i+1}+x_{i-1}-2x_i).$$

Vamos propor uma solução do tipo

(10) $$x_i(t)=A_i{e}^{jwt},$$

onde j é o número imaginário tal que j² = −1. Inserindo-o na equação anterior, temos que

(11) $$-w^2{A}_i=\frac{k}{m}(A_{i+1}+A_{i-1}-2A_i).$$

Estabeleceremos que $w_0^2=\frac{k}{m}$. Vale notar que, como a constante elástica de uma mola possui dimensões, no Sistema Internacional de Unidades, de $\frac{[Kg]}{{[s]}^2}$ (massa por tempo ao quadrado), a constante w₀ tem dimensão de $\frac{1}{[s]}$, a mesma dimensão de uma frequência. Desta forma, podemos reescrever a equação anterior como

(12) $$\frac{2w_0^2-w^2}{w_0^2}=\frac{A_{i+1}+A_{i-1}}{A_i}.$$

Para uma dado valor w, temos que o lado esquerdo da equação deve ser constante, pois w₀ é um valor arbitrário. Portanto, para estar condizente com o lado esquerdo da equação, o lado direito também deve ser constante. Vamos supor, agora, que 0≤w≤2w₀. Isto nos permite escrever a constante do lado esquerdo da equação como 2cos⁡(θ), pois a fração irá variar, no máximo, entre −2 e 2. Assim,

(13) $$\cos (\theta )={\displaystyle \frac{A_{i+1}+A_{i-1}}{2A_i}},$$

(14) $$2\cos (\theta )={\displaystyle \frac{2w_0^2-w^2}{w_0^2}}.$$

Suponhamos que qualquer A_i possa ser escrito como A_i = Bcos⁡(iθ) + Csin⁡(iθ). Desejamos considerar desta forma para provar, por indução, que tal suposição é válida. Podemos determinar o valor de B através de A₀≡B. C, por sua vez, pode ser determinado a partir de A₁ = Bcos⁡(θ) + Csin⁡(θ). Precisamos, por fim, determinar A_{i + 1}. Da equação 13, temos que

(15) $$A_{i+1}=2A_i\cos (\theta )-A_{i-1},$$

(16) $$A_{i+1}=2\cos (\theta )\{B\cos (i\theta )+C\sin (i\theta )\}-\{B\cos ((i-1)\theta )+C\sin ((i-1)\theta )\}$$

e

(17) $$A_{i+1}=B\{2\cos (\theta )\cos (i\theta )-\cos ((i-1)\theta )\}+C\{2\cos (\theta )\sin (i\theta )-\sin ((i-1)\theta )\},$$

mas

(18) $$\cos ((i-1)\theta )=\cos (i\theta -\theta )=\cos (n\theta )\cos (\theta )+\sin (i\theta )\sin (\theta ),$$

e

(19) $$\sin ((i-1)\theta )=\sin (i\theta -\theta )=\sin (i\theta )\cos (\theta )-\sin (\theta )\cos (i\theta ),$$

então,

(20) $$A_{i+1}=B\{\cos (\theta )\cos (i\theta )-\sin (i\theta )\sin (\theta )\}+C\{\cos (\theta )\sin (i\theta )+\sin (\theta )\cos (i\theta )\},$$

o que nos leva a

(21) $$A_{i+1}=B\cos ((i+1)\theta )+C\sin ((i+1)\theta ).$$

Portanto, como A_i é válido para A₀,A₁,A_{i + 1}, provamos, por indução, que ele é válido para qualquer A_i e obtemos a solução

(22) $$x_i=\{B\cos (i\theta )+C\sin (i\theta )\}{e}^{jwt},$$

desde que a suposição sobre a frequência seja respeitada. Para o estudo em questão, possuímos duas condições de contorno, relativas às paredes fixas do sistema. De x₀(t) = 0, temos que B = 0. De x_{N + 1}(t) = 0, temos que $\theta =\frac{m\pi }{N+1}$, para uma solução geral não trivial, sendo m um número inteiro. Desta maneira, resta apenas o termo senoidal para descrever o movimento da partícula: para cada modo normal, podemos associar uma curva senoidal, que limitará os deslocamentos máximos de cada partícula em relação ao seu centro de equilíbrio, assim como fixará a direção do deslocamento relativo das mesmas. Isto é evidente na Figura 2, ao analisarmos o movimento das partículas no seu terceiro modo de vibração. Quando as partículas 5 e 6 estiverem na sua distância máxima para à direita do seu centro de equilíbrio, as partículas 16 e 17 estarão deslocadas esta mesma distância, para à esquerda, de seus respectivos centros.

Figura 2
Representação de partículas e o seu deslocamento relativo ao ponto de equilíbrio para os três primeiros modos normais. As linhas tracejadas, pontilhadas e circulares espaçadas representam o primeiro, segundo e terceiro modo, respectivamente.

Encontrados os θ possíveis, podemos descobrir as frequências associadas a cada modo. Da equação 14, podemos escrever

(23) $$w^2=2w_0^2(1-\cos (\theta )),$$

(24) $$w^2=2w_0^2\left\{1-\cos \left({\displaystyle \frac{m\pi }{N+1}}\right)\right\},$$

(25) $$w^2=2w_0^2\left\{2{\sin }^2\left({\displaystyle \frac{m\pi }{2(N+1)}}\right)\right\},$$

e

(26) $$w=2w_0\left\{\sin \left(\frac{m\pi }{2(N+1)}\right)\right\}.$$

O movimento da partícula é dado como a soma de todas as soluções particulares (modos normais). Como estamos levando em consideração um número discreto de partículas, consequentemente teremos uma quantidade discreta de modos. Logo, a solução mais geral, em que tomamos a parte real da exponencial, será

(27) $$x_i(t)=\sum _{m=1}^N{C}_m\sin \left(\frac{im\pi }{N+1}\right)\cos (w_{m}t),$$

e

(28) $$w_m=2w_0\sin \left(\frac{m\pi }{2(N+1)}\right).$$

A energia de cada modo se divide entre potencial e cinética. Portanto, a energia total, que é a soma da energia alocada em cada modo (descartando termos de ordem maior), será dada por

(29) $$E_{total}=\sum _{i=1}^N\frac{w_i^2{A}_i^2+{\stackrel{.}{A}}_i^2}{2}.$$

Notemos que A_i se refere à amplitude da oscilação de um dado modo normal e ${\stackrel{.}{A}}_i$ se refere à amplitude da velocidade deste mesmo modo. A questão principal, portanto, será como obter estas amplitudes para todo o tempo estabelecido. Se, para cada partícula, sabemos a sua posição, podemos aplicar uma transformada de Fourier no conjunto de posições das partículas e obter, assim, a amplitude de oscilação de todos os modos normais. De maneira análoga, se aplicarmos uma transformada de Fourier no conjunto de velocidades de todas as partículas, obteremos a amplitude de velocidade de todos os modos normais. Então, se criarmos um programa que solucione as equações de movimento e retorne as posições e velocidades para qualquer tempo, podemos obter as amplitudes de oscilação da posição e da velocidade para cada frequência normal e, assim, calcular tanto a energia cinética alocada em cada modo normal, quanto a energia total.

4. Teorema da Equipartição de Energia

Suponhamos que um sistema esteja isolado, ou seja, que a energia interna dele seja constante ao longo do tempo. Para um determinado tempo, é esperado que este sistema atinja equilíbrio térmico, isto é, sozinho, o sistema seja permitido assumir determinadas características intrínsecas. Estas características podem ser, no caso de um gás, a pressão ou a temperatura próprias. Já para o sistema FPUT, uma característica relevante é a contribuição média de cada termo do Hamiltoniano para a energia do sistema.

A mecânica estatística clássica fornece um resultado muito útil quando avaliamos justamente essa contribuição média para a energia interna do sistema, baseada em duas premissas. A energia mencionada na equação 29 pode ser função das coordenadas generalizadas e dos momentos correspondentes

(30) $$E=E(q_i,\mathrm{\dots },q_N,p_i,\mathrm{\dots },p_N).$$

Dependendo de cada situação, a energia pode ser aditiva, a qual pode se separar em funções de uma variável apenas, como

(31) $$E={\epsilon }_i(p_i)+E^{\prime }(q_1,\mathrm{\dots },p_N),$$

onde ε_i é parte da energia dependente apenas da variável p_i, podendo ser uma coordenada ou momento, e o restante da energia total E′ não depende desta variável. A energia também pode se separar em funções quadráticas, seguindo a forma

(32) $${\epsilon }_i(p_i)=a_i{p}_i^2,$$

onde a_i é uma constante. Essas duas características são utilizadas pela mecânica estatística para derivar o Teorema da Equipartição de Energia. Vale notar que o sistema FPUT segue essas duas premissas: ele se divide em 2N funções quadráticas de uma variável apenas, sendo passível da aplicação do teorema.

Com o nosso sistema em equilíbrio térmico, ele assume uma temperatura bem definida, a qual se relaciona com a energia pela constante de Boltzman, K_b. Esse sistema também permite que seja utilizado o ensemble canônico para calcular o valor médio esperado conforme

(33) $$⟨{\epsilon }_i(p_i)⟩=\frac{{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\epsilon }_i(p_i)e^{-\beta E}\mathit{d}{q}_1\mathrm{\dots }\mathit{d}{p}_f}{{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{e}^{-\beta E}\mathit{d}{q}_1\mathrm{\dots }\mathit{d}{p}_f},$$

em que β é relacionado à temperatura por β = (K_bT)⁻¹. Ao utilizar as características da energia, podemos continuar desenvolvendo a equação 33 como

(34) $$⟨{\epsilon }_i(p_i)⟩={\displaystyle \frac{{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\epsilon }_i(p_i)e^{-\beta {\epsilon }_i(p_i)}\mathit{d}{p}_i{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{e}^{-\beta E^{\prime }}\mathit{d}{q}_1\mathrm{\dots }\mathit{d}{p}_f}{{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{e}^{-\beta {\epsilon }_i(p_i)}\mathit{d}{p}_i{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{e}^{-\beta E^{\prime }}\mathit{d}{q}_1\mathrm{\dots }\mathit{d}{p}_f}},$$

(35) $$⟨{\epsilon }_i(p_i)⟩=a_i{\displaystyle \frac{{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{p}_i^2{e}^{-\beta a_i{p}_i^2}\mathit{d}{p}_i}{{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{e}^{-\beta a_i{p}_i^2}\mathit{d}{p}_i}},$$

e

(36) $$⟨{\epsilon }_i(p_i)⟩=\frac{1}{2}{K}_{b}T.$$

Portanto, cada termo quadrático do Hamiltoniano contribui com um valor médio de $\frac{1}{2}{K}_{b}T$ para a energia média do sistema, quando o mesmo se encontrar em equilíbrio térmico [⁷]. Vale notar que termos não quadráticos para osciladores não-harmônicos, desde que aditivos, estão sujeitos à valores médios similares aos correspondentes aos termos quadráticos. A energia interna média depende da potência que a posição e o momento conjugados são elevados na função da energia. Alguns casos não quadráticos são abordados pedagogicamente por Prasanth e co-autores [⁸].

Como demonstrado, se o sistema está em equilíbrio térmico, as coordenadas e momentos generalizados q_i e p_i deveriam contribuir na mesma proporção para a energia total do sistema (se ambos forem termos quadráticos). Visando a solução com apenas termos lineares na força, da Seção ³, em equilíbrio térmico, as amplitudes e suas derivadas temporais, A_i e ${\stackrel{.}{A}}_i$, deveriam contribuir na mesma quantidade para a energia total do sistema, por ambos serem termos quadráticos. Em outras palavras: a energia associada a cada modo normal de vibração deveria ser a mesma, em média.

Buscando esta expectativa da equipartição de energia de um sistema em equilíbrio térmico, analisaremos o problema FPUT, em sua totalidade, através de simulações numéricas. Exemplificaremos situações em que acontecem a equipartição de energia, como previsto pelo teorema apresentado, e uma situação específica que levou ao conhecido paradoxo FPUT, em que a equipartição de energia não acontece.

5. Implementação em Python

Para o caso em estudo, consideremos a equação 1, em que as forças têm uma componente quadrada e que tenhamos N partículas. Podemos criar um vetor de tamanho 2N, no qual os primeiros N termos representam o deslocamento relativo e os últimos N termos representam as velocidades. Assim, o sistema de N equações diferenciais de segunda ordem se torna um sistema de 2N equações diferencias de primeira ordem. Este sistema pode ser descrito da seguinte maneira, com i variando de 2 a (N-1)

(37) $${\stackrel{.}{x}}_1=x_{N+1},$$

(38) $${\stackrel{.}{x}}_i=x_{N+i},$$

(39) $${\stackrel{.}{x}}_N=x_{2N},$$

(40) $${\stackrel{.}{x}}_{N+1}=x_2-2x_1+\alpha \{{(x_2-x_1)}^2+x_1^2\},$$

(41) $${\stackrel{.}{x}}_{N+i}=x_{i+1}+x_{i-1}-2x_i+\alpha \{{(x_{i+1}-x_i)}^2+{(x_{i-1}-x_i)}^2\}$$

e

(42) $${\stackrel{.}{x}}_{2N}=x_{N-1}-2x_N+\alpha \{x_N^2+{(x_{N-1}-x_N)}^2\}.$$

A solução de tal sistema pode ser realizada em qualquer linguagem computacional, mas neste trabalho adotaremos o Python. Os códigos serão fornecidos como material suplementar acessíveis via repositório GitHub [⁹] e os mesmos estão detalhadamente comentados. As integrações numéricas são feitas pela função odeint da biblioteca scipy [¹⁰], adequados para resolver sistemas na forma $\frac{dy}{dt}=f(t,y)$ através do método de Adams (predictor-corretor) [¹¹].

É esperado que, em condições de equilíbrio térmico, o Teorema de Equipartição de Energia ocorra. No programa implementado, se tivermos uma interação não-linear forte ou uma quantidade inicial de energia consideravelmente alta, observaremos que a equipartição de energia acontecerá. Para ilustrar essa situação, executemos o código ajustando a sua condição inicial: todas as partículas possuirão velocidade inicial nula. Além disso, seus deslocamentos estarão apenas no primeiro modo normal, cuja amplitude máxima será de 10. Estas condições têm o efeito prático de tornar a energia do sistema relativamente alta. Com estes resultados, dois gráficos foram gerados. Na Figura 3, temos a representação da quantidade de energia alocada nos primeiros cinco modos de vibração pelo tempo de evolução. Já a Figura 4 mostra a energia média final de cada modo normal. É importante notar que, para se obter a energia média final, foram utilizados os últimos dez valores obtidos de cada modo normal. Desta forma, notamos que a energia é distribuída ao longo do tempo entre os diversos modos normais.

Figura 3
Gráfico das energias armazenadas nos cinco primeiros modos normais: 1° modo – azul, 2° modo – laranja, 3° modo – verde, 4° modo – vermelho e 5° modo – lilás.

Figura 4
Gráfico da energia média final associada a cada modo normal.

Este resultado não foi o encontrado por Fermi e co-autores, pois os mesmos estabeleceram uma energia total inicial menor. Para realizarmos uma simulação similar à realizada, executaremos o código anteriormente mencionado alterando a amplitude máxima do primeiro modo normal para 1, sendo os resultados demonstrados graficamente nas Figuras 5 e 6.

Figura 5
Gráfico das energias armazenadas nos 5 primeiros modos normais: 1° modo – azul; 2° modo – laranja; 3° modo – verde; 4° modo – vermelho; e 5° modo – lilás (no primeiro ciclo).

Figura 6
Gráfico das energias armazenadas nos dois primeiros modos normais: 1° modo – azul e 2° modo – laranja.

Nota-se, na Figura 5, que a energia está, no começo, totalmente no primeiro modo de vibração, devido às condições iniciais das partículas. Ao acompanhar a energia desse modo, percebe-se que a mesma se distribui ao longo do tempo para outros modos normais de valor baixo. Contudo, antes que esta energia possa ser alocada em modos mais altos, a mesma retorna quase inteiramente para o primeiro modo normal. Esta recorrência é exibida na segunda coluna da Tabela 1 e é denominada de Frequência de Recorrência de Primeira Ordem. A primeira coluna da tabela contém a frequência natural do modo normal – no caso, para o primeiro modo – e é calculada diretamente da equação 25. Aqui, ela nos ajuda a fornecer uma dimensão das frequências envolvidas no problema.

Tabela 1
Tabela das frequências para a solução numérica, em Hz.

1º modo normal	1ª ordem	2ª ordem
9,52×10−2	9,62×10−5	4,18×10−6

A Figura 6 representa apenas os modos 1 e 2, de modo a facilitar a compreensão dos dados. É notável, portanto, o fluxo do ciclo de energia partindo do primeiro modo normal e depois se distribuindo para os demais, representados pelo segundo modo normal. Desejamos enfatizar que a quantidade de energia retornada para o modo 1 varia com o tempo. Ou seja, além da recorrência de primeira ordem, temos uma recorrência de segunda ordem com uma frequência menor, como mostrado na Tabela 1.

Na recorrência de segunda ordem, podemos notar que os picos do primeiro modo normal diminuem, perdendo até 40% da sua energia inicial. Depois, retornam basicamente às condições iniciais, com apenas 3% de perda. Para melhor visualização dessa oscilação, temos o gráfico (Figura 7), com os valores dos picos do primeiro modo normal em função do tempo. Neste, esta oscilação é evidenciada e podemos, inclusive, determinar sua frequência. Na Figura 7, a linha vermelha com pontos mostra o valor relativo do pico em relação ao valor da energia total, de forma a facilitar a percepção da quantidade percentual da perda de energia para os demais modos normais.

Figura 7
Gráfico dos picos do primeiro modo normal em função do tempo, gerando uma função “envelope”, representada pela linha azul. A quantidade de energia alocada no primeiro modo normal em relação à energia total é sinalizada pela linha vermelha com pontos.

Ao longo dos anos, com o estudo do paradoxo de FPUT, foi questionada a existência de recorrências de níveis mais altos. Com este questionamento surge, no entanto, a dificuldade de se notar numericamente as alterações conforme aumenta-se a ordem de recorrência. Além disso, para as alterações serem notadas, a simulação deve ocorrer por mais tempo, o que está além da capacidade de um computador pessoal. Apesar de já ter sido comprovada a existência de recorrências de terceira e quarta ordem [¹²], não foi possível demonstrar nesta simulação, demasiadamente curta.

Dois exemplos foram analisados: o primeiro, no qual foi fornecido uma grande quantidade de energia para o primeiro modo normal e outro, em que a energia fornecida não foi tão alta quanto a utilizada neste trabalho. No primeiro caso, vemos uma dispersão da energia por todos os modos normais, enquanto que, no segundo caso, a energia não se equiparte. Podemos explicar este comportamento através de um teorema enunciado por Kolmogorov [¹³], o qual deu surgimento à teoria KAM (Kolmogorov-Arnol’d-Moser). Esta teoria diz respeito ao fato que, se possuímos um sistema oscilatório conservativo não-linear, em que existe uma interação não-linear fraca, essa interação servirá apenas para mudar levemente as frequências de oscilação e introduzir pequenos harmônicos não-lineares. Porém, se possuímos uma interação não-linear forte, ou um sistema com energia alta, é esperado que ele atinja equilíbrio térmico, no caso abordado neste trabalho [¹⁴].

6. Pacote computacional

A solução numérica da seção anterior é relevante e factível apenas para um número pequeno de partículas e passos de integrações. Pesquisadores, atualmente, necessitam lidar com um número da ordem de milhares, ocasionalmente bilhões, e até cerca de 10¹⁰ passos de integrações. Desta forma, a aplicação de métodos numéricos mais simples é o usual em ferramentas computacionais utilizados pelos mesmos [¹¹]. Nesta seção, utilizaremos um pacote computacional que é bastante adotado na comunidade científica, fornecendo uma oportunidade para o estudante utilizar uma ferramenta no seu estado da arte.

Para se obter as posições e velocidades de cada partícula, podemos usar a abordagem de Dinâmica Molecular implementada no LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) [¹⁵]. Esse tipo de abordagem é bastante utilizada atualmente nos campos da Física, Química, Biologia e Engenharias, pois permite uma fácil visualização dos movimentos atômicos e, computacionalmente, há vários pacotes de Dinâmica Molecular disponíveis.

A Dinâmica Molecular, por sua vez, se baseia na possibilidade de descrever a trajetória de um sistema de partículas através da solução das Leis de Newton por meio de computadores [¹⁶] (ou, de maneira equivalente, um sistema Hamiltoniano para mecânica clássica [¹⁷]). No caso do LAMMPS, as integrais são resolvidas por um algoritmo chamado velocity-Verlet [¹⁸]. Apesar da solução por Dinâmica Molecular se basear na solução das leis de Newton, como feito na seção anterior, não será necessário reescrever o conjunto de equação diferenciais, sendo necessário informar apenas o tipo de ligação entre partículas. Para uma melhor compreensão dos conceitos e fundamentos da Dinâmica Molecular, recomenda-se os livros de Mark Tuckerman [¹⁶] e Allen e Tildesley [¹⁹]. Visando a solução do sistema de equações, deve-se conhecer as condições iniciais do mesmo e, para isso, devemos informar a posição inicial de cada partícula. Além disso, é necessário determinar qual o tipo de ligação que essas partículas possuem e com quais elementos as mesmas estão conectadas. No caso abordado neste trabalho, o número de partículas deve ser uma potência de dois (o que permite uma transformação rápida de Fourier no programa de análise), adicionado de duas partículas a mais (que neste caso atuam como paredes), totalizando N partículas. Além disso, estamos considerando molas de constante igual que unem as partículas, portanto temos apenas uma espécie de ligação. Como estas partículas estão ligadas apenas a suas vizinhas, temos (N-1) ligações. Em um outro tipo de situação, com escala mais aprofundada e características mais complexas, essas ligações descreveriam, por exemplo, ligações químicas usando potenciais adequados, os quais podemos encontrar publicados em revistas científicas.

O LAMMPS é processado através de comandos fornecidos por um arquivo de entrada, chamado script. Entretanto, para facilitar a preparação da simulação, foram elaborados dois arquivos editáveis: um arquivo contém as condições iniciais para a simulação e o outro conta com a descrição das forças de interação. O arquivo com as condições iniciais possui duas seções: uma dedicada aos chamados “átomos” e outra para as “ligações”.

Para a sessão “átomos”, o que determina as informações necessárias é o tipo de partícula a ser simulado, sendo nesse caso imprescindível a identificação das partículas e suas coordenadas. Na segunda seção (“ligações”), precisamos identificar os elementos conectados. Além destas informações, na sintaxe do LAMMPS, existe um cabeçalho que contém os limites do nosso espaço de simulação. A necessidade da criação do arquivo relacionado às forças de interação surgiu do fato de nenhum modelo contido no LAMMPS assemelhar-se às forças de interação exigidas pelo nosso estudo de caso. O imprescindível deste arquivo é a descrição da magnitude das forças e a energia de ligação em função da distância da ligação. Todos os arquivos mencionados anteriormente estão escritos a partir de códigos em Python de forma a serem facilmente ajustáveis.

Para podermos iniciar a simulação, foi elaborado um script que condensa os comandos necessários para se obter os resultados. Neste, são especificadas as unidades a serem usadas, o tipo de partícula e a dimensão do espaço. Uma descrição detalhada de cada comando, com seus argumentos possíveis, pode ser encontrada na documentação do LAMMPS, em inglês. Há comandos que determinam o tipo de integração a ser adotado – aqui velocity-Verlet. Para podermos obter resultados comparáveis com a seção anterior, precisamos estabelecer como condição de contorno um sistema físico com paredes fixas. Assim, apenas um sub-grupo de partículas que não contém os elementos das bordas terá seu movimento integrado, excluindo as partículas das extremidades do sistema. No final do script é solicitado para o LAMMPS gerar um arquivo com a posição e velocidade de todos as partículas para cada iteração do simulador.

Na posse deste arquivo, foi escrito um código em Python para podermos analisar os resultados e compará-los à solução numérica. Basicamente, o que realmente é necessário na análise é agrupar as informações de cada partícula e ajustar a posição de cada um para o deslocamento relativo. Depois disso, como realizado na solução numérica, executamos uma transformada de Fourier e, por fim, conseguimos comparar os resultados obtidos anteriormente na seção, com a solução numérica.

Podemos notar que os gráficos das Figuras 8 e 9 são muito similares aos apresentados na solução numérica. Porém, há uma distinção entre eles, a saber, o número de iterações que é realizado. No caso de Dinâmica Molecular, este número é função do timestep (passo de integração). Na natureza, o tempo é uma variável contínua, já em computadores, a variável tempo é discreta – passa em passos. Para haver conservação da energia é necessário que haja uma continuidade temporal. Em um computador, aproximamos da realidade através da diminuição do timestep [¹⁶]. Por um lado, timesteps menores nos levam à conservação de energia melhores, mas também um número de iterações maiores, com processamentos mais longos, o que implica no aumento dos custos computacionais.

Figura 8
Representação das simulações de Dinâmica Molecular com dois timesteps diferentes, mostrando as energias armazenadas nos cinco primeiros modos normais: 1° modo – azul; 2° modo – laranja; 3° modo – verde; 4° modo – vermelho; e 5° modo – lilás (no primeiro ciclo).

Figura 9
Representação das simulações de Dinâmica Molecular com dois timesteps diferentes, mostrando as energias armazenadas nos dois primeiros modos normais: 1° modo – azul e 2° modo – laranja.

Quanto maior o passo de integração numa simulação de Dinâmica Molecular, mais impreciso se torna a conservação de energia. Isso pode ser percebido pela maior participação dos modos normais, além do primeiro, na alocação de energia, assim como nas frequências de recorrência que foram calculadas. Tomando as frequências obtidas pela solução numérica como referência, os dados para cada frequência, assim como suas concordâncias associadas, são apresentados na Tabela 2. Em Dinâmica Molecular, a conservação de energia depende do tamanho do timestep. Longos passos de integração levam a uma pior conservação da energia [¹⁶].

Tabela 2
Tabela de frequências, em Hz, para dois timesteps diferentes e valores relativos percentuais Δf_t comparados aos da solução numérica, sendo o subscrito _t indicativo de timestep 0,1 ou 0,5.

Timestep	1ª ordem	Δf0,1	2ª ordem	Δf0,5
0.1	9.55×10−5	99	4.34×10−6	96
0.5	8.25×10−5	86	3.75×10−6	89

A partir das Figuras 8 e 9, podemos mensurar as frequências f_DM de 1^a e 2^a ordem para os timesteps de 0,1 e 0,5. Os resultados dessas determinações são mostrados na Tabela 2, junto com os valores de concordância, Δf, entre a solução por Dinâmica Molecular f_DM e a solução numérica f_SN

(43) $$\mathrm{\Delta }f=100\times \left(1-\frac{|f_{DM}-f_{SN}|}{f_{SN}}\right).$$

Notamos que há uma redução relativa das frequências para a solução numérica e, quanto maior o timestep, maior ainda é a redução da frequência. Os valores de 2^a ordem são mais afetados pois necessitam de processamentos mais longos para se obter precisão similar que os de 1^a ordem. Assim, a conservação de energia para processamentos mais longos é diretamente afetada pela grandeza do timestep, como discutido anteriormente.

7. Exercícios pedagógicos

Nossos programas e scripts estão disponíveis através de um repositório Github [⁹]. Para processar programas em Python, será necessário a instalação de um pacote Python. Indicamos o pacote gratuito Anaconda disponível em https://www.anaconda.com/. O pacote LAMMPS está disponível gratuitamente em https://lammps.sandia.gov/download.html. Os exercícios propostos abaixo podem necessitar de adaptações dos códigos de programação em Python, as quais os alunos são aptos a realizar. Nossa intenção é que os exercícios permitam desde uma simples repetição dos resultados apresentados, até uma investigação inicial mais profunda, e que possam gerar novas abordagens para os professores de Ensino Médio e Superior.

1. Processe os programas e obtenha as frequências observadas no artigo.

2. Faça as alterações necessárias (como comentado no script do LAMMPS) para realizar as simulações para um caso cúbico. Explique o porquê de certos modos de energia não serem excitados.

3. Estabeleça a relação entre a energia depositada inicialmente no primeiro modo normal e a distribuição nos modos normais com o passar do tempo (observe que o valor praticamente se conserva).

4. Adotamos no artigo o Sistema Internacional de Unidades. Para isso, tomamos como exemplo a massa da partícula como 1 kg e a unidade de tempo como o segundo. Porém, utilizar unidades reduzidas não acarretará em alterações. Sugere-se explorar a transformação de unidades para sistemas não convencionais.

5. Se utilizar o LAMMPS, pode-se imprimir as posições de forma a criar um filme das partículas em movimento. Para isso, pode-se usar um visualizador de partículas e moléculas como o Ovito, disponível em https://www.ovito.org. Tente relacionar a visualização com o que é lido, por exemplo, dos gráficos 1.

6. Para qual valor de α máximo, no código em Python para solução numérica, ainda conseguimos observar as recorrências de FPUT?

7. Examine a diferença entre posições de partículas vizinhas x_{i + 1}−x_i (é possível detectar a presença de estruturas propagantes que são, na verdade, sólitons de Zabusky-Kruskal [²⁰]).

References

• ^1.
  E. Fermi, J. Pasta e S. Ulam, Los Alamos Internal Report, University of California, California (1955).
• ^2.
  G. Gallavotti, The Fermi-Pasta-Ulam problem: a status report (Springer, Berlin, 2007), v. 728.
• ^3.
  T. Dauxois, Phys. Today 61, 55 (2008).
• ^4.
  M. Onorato, L. Vozella, D. Proment e Y. Lvov, Proceedings of the National Academy of Sciences 14, 4208 (2015).
• ^5.
  G. Arfken e H. Weber, Mathematical methods for physicists (Elsevier, Amsterdam, 1999).
• ^5.
  W. Boyce e R. DiPrima, Elementary differential equations and boundary value problems (John Wiley & Sons, New York,2012).
• ^7.
  F. Reif, Fundamentals of Statistical and Thermal Physics (Waveland Press, New York, 2009).
• ^8.
  P. Prasanth, P. Nishanth e K. Udayanandan, Physics Education 35, 1 (2019).
• ^9.
  D. Lucero, scrpits, disponível em: https://github.com/davi-lucero/codes-art, acessado em 03/12/2020.
  » https://github.com/davi-lucero/codes-art
• ^10.
  P. Virtanen, R. Gommers, T. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Brightet al., Nature methods 17, 261 (2020).
• ^11.
  R. Landau, M. Paez e C. Bordeianu, Computational physics: Problem solving with Python (John Wiley & Sons,Weinheim, 2015).
• ^12.
  S. Pace e D. Campbell, Chaos: An Interdisciplinary Journal of Nonlinear Science 29, 023132 (2019).
• ^13.
  A. Kolmogorov, Doklad. Acad. Nauk USSR, 98, 527 (1954).
• ^14.
  G. Walker e J. Ford, Phys. Rev. 188, 416 (1969).
• ^15.
  S. Plimpton, Journal of Computational Physics 117, 1 (1995).
• ^16.
  M. Tuckerman, Statistical mechanics: theory and molecular simulation (Oxford University Press, Oxford, 2010).
• ^17.
  H. Goldstein, C. Poole e J. Safko, Classical mechanics (Addison Weslwy, San Francisco, 2002).
• ^18.
  M. Newman, Computational physics (University of Michigan, Michigan, 2013).
• ^19.
  M. Allen e D. Tildesley, Computer simulation of liquids (Oxford University Press, Oxford, 2017).
• ^20.
  N. Zabusky e M. Kruskal, Physical Review Letters 15, 240 (1965).

Datas de Publicação

• Publicação nesta coleção
  03 Mar 2021
• Data do Fascículo
  2021

Histórico

• Recebido
  03 Dez 2020
• Revisado
  01 Fev 2021
• Aceito
  01 Fev 2021