Problema de Störmer: um estudo computacional sobre a formação do cinturão de radiação de Van Allen

Störmer problem: a computational study on the formation of the Van Allen radiation belt

Resumo

Neste trabalho, apresentamos o problema de Störmer no contexto da mecânica clássica: o movimento de partículas carregadas não relativísticas sob a ação do campo magnético dipolar da Terra. A partir da formulação lagrangiana em coordenadas cilíndricas, obtemos as equações de movimento e identificamos as constantes de movimento do sistema. O sistema resultante é não integrável, e recorremos ao método simplético de Störmer-Verlet de segunda ordem para integrar numericamente as equações de movimento. Três casos são analisados: o caso equatorial, em que a partícula permanece restrita ao plano $Z=0$; o caso tridimensional completo; e o caso vinculado a uma superfície esférica. No caso equatorial, mostramos que o potencial efetivo possui um mínimo e um máximo locais, e que partículas com energia abaixo do máximo ficam confinadas em órbitas oscilatórias. No caso tridimensional, a dinâmica apresenta sensibilidade às condições iniciais e trajetórias que reproduzem a estrutura toroidal do cinturão de radiação de Van Allen. No caso esférico, o sistema torna-se integrável e admite solução analítica em funções elípticas de Jacobi, permitindo validar a precisão do integrador numérico.

Palavras-chave:
Problema de Störmer; cinturão de Van Allen; integração simplética; mecânica clássica

Abstract

We present the Störmer problem in the framework of classical mechanics: the motion of non-relativistic charged particles under the Earth’s dipolar magnetic field. Starting from the Lagrangian formulation in cylindrical coordinates, we derive the equations of motion and identify the conserved quantities. The resulting system is non-integrable, and we employ the second-order symplectic Störmer-Verlet method to numerically integrate the equations of motion. Three cases are analyzed: the equatorial case, where the particle is restricted to the $Z=0$ plane; the full three-dimensional case; and the case where the particle is constrained to a spherical surface. In the equatorial case, we show that the effective potential has a local minimum and maximum, and that particles with energy below the maximum are confined to oscillatory orbits. In the three-dimensional case, the dynamics exhibits sensitivity to initial conditions and trajectories that reproduce the toroidal structure of the Van Allen radiation belt. In the spherical case, the system becomes integrable and admits an analytical solution in Jacobi elliptic functions, which we use to validate the numerical integrator.

Keywords:
Störmer problem; Van Allen radiation belt; symplectic integration; classical mechanics

1. Introdução

1.1. As auroras e o problema de Störmer

As auroras polares são registradas em relatos desde a Antiguidade, porém uma explicação física para o fenômeno só surgiu no final do século XIX. Em 1896, o físico norueguês Kristian Birkeland propôs que as auroras seriam causadas por feixes de partículas carregadas emitidas pelo Sol e deflectidas pelo campo magnético da Terra em direção aos polos [¹]. Para testar essa hipótese, Birkeland construiu um dispositivo experimental conhecido como terrella: uma esfera magnetizada, com um eletroímã interno, colocada dentro de uma câmara de vácuo. Ao disparar raios catódicos contra a terrella, Birkeland observou que os elétrons se concentravam em regiões anulares ao redor dos polos magnéticos da esfera, reproduzindo em laboratório o padrão geográfico das auroras [¹].

Quem se dedicou a tratar o problema de forma matemática rigorosa foi Carl Störmer (1874–1957), matemático e físico norueguês contemporâneo de Birkeland. Störmer formulou as equações de movimento de partículas carregadas sob a ação de um campo magnético dipolar e dedicou-se, ao longo de décadas, a calcular numericamente milhares de trajetórias de partículas, identificando padrões de confinamento e espalhamento que dependiam das condições iniciais [², ³]. Os cálculos eram feitos à mão, com o auxílio de assistentes, e os resultados foram compilados em seu livro The Polar Aurora, publicado em 1955 [⁴]. Os métodos numéricos que Störmer desenvolveu para integrar as equações de movimento, baseados em fórmulas de diferenças finitas com múltiplos passos, são, de certa forma, os precursores dos integradores simpléticos modernos, como veremos na Seção 3.

Em paralelo ao trabalho de Störmer, o estudo dos raios cósmicos trouxe evidências independentes da influência do campo geomagnético sobre partículas carregadas. No início da década de 1930, Clay e Berlage [⁵] e, logo em seguida, Compton [⁶] observaram que a intensidade dos raios cósmicos que atingem a superfície terrestre depende da latitude geográfica: a intensidade é menor próxima ao equador e maior em latitudes elevadas. Esse efeito de latitude tem uma explicação direta no contexto do problema de Störmer. O campo magnético dipolar é mais intenso na região equatorial e deflete as partículas de menor energia, que não conseguem penetrar até a superfície. Em latitudes altas, onde as linhas do campo magnético do dipolo convergem e se aproximam da vertical, as partículas encontram menos resistência magnética e atingem a atmosfera com maior facilidade. A análise quantitativa de Störmer já previa esse comportamento, e a confirmação experimental do efeito de latitude fortaleceu a validade do modelo dipolar.

1.2. A descoberta dos cinturões de Van Allen

A confirmação mais direta das previsões de Störmer veio em 1958, no contexto da corrida espacial. O satélite Explorer 1, lançado em janeiro daquele ano, carregava um detector Geiger projetado pelo físico James Van Allen e sua equipe na Universidade de Iowa. Os dados transmitidos pelo satélite revelaram algo inesperado: em determinadas altitudes, o contador Geiger saturava, indicando níveis de radiação muito acima do previsto. A princípio, houve dúvida sobre se o instrumento havia falhado. Medições complementares do satélite Explorer 3, lançado em março de 1958, confirmaram que a saturação era real [⁷]. O que os detectores registravam era a presença de grandes quantidades de prótons e elétrons aprisionados pelo campo magnético terrestre, concentrados em duas regiões toroidais ao redor da Terra, exatamente o tipo de confinamento que Störmer havia previsto teoricamente décadas antes.

Essas regiões receberam o nome de cinturões de radiação de Van Allen. O cinturão interno, situado entre aproximadamente 1.000 e 6.000 km de altitude, é composto predominantemente por prótons de alta energia. O cinturão externo, entre cerca de 13.000 e 60.000 km, contém majoritariamente elétrons. A existência desses cinturões tem consequências práticas: satélites e sondas que atravessam ou orbitam nessas regiões ficam expostos a níveis elevados de radiação, que podem danificar componentes eletrônicos e representar riscos para astronautas [⁸, ⁹, ¹⁰]. A Anomalia Magnética do Atlântico Sul, uma região onde o cinturão interno se aproxima da superfície terrestre devido a irregularidades no campo geomagnético, é um exemplo concreto desse problema: satélites em órbita baixa sofrem maior incidência de falhas eletrônicas ao passarem sobre essaregião [⁸].

1.3. Contexto e objetivos deste trabalho

Do ponto de vista da mecânica analítica, o problema de Störmer é um sistema hamiltoniano com três graus de liberdade, nas coordenadas cilíndricas $(\rho ,\varphi ,Z)$. A simetria axial torna a coordenada $\varphi$ cíclica e conserva o momento canonicamente conjugado a $\varphi$, que denotamos $c_2$; com isso, a dinâmica se reduz ao plano meridional $(\rho ,Z)$. O sistema dispõe, contudo, de apenas duas constantes de movimento independentes, a energia e $c_2$, uma a menos do que as três que a integrabilidade no sentido de Liouville exigiria [¹¹, ¹²]. Esse fato, demonstrado rigorosamente por Dilão e Alves-Pires, significa que não existe uma terceira constante de movimento independente que permita reduzir o sistema a quadraturas. A consequência prática é que as equações de movimento só podem ser resolvidas numericamente, e que certas trajetórias apresentam sensibilidade exponencial às condições iniciais, ou seja, comportamento caótico.

O objetivo deste artigo é apresentar o problema de Störmer com a formulação matemática completa e a implementação computacional, em um nível adequado para estudantes de graduação em física. O problema reúne, num mesmo contexto físico, as formulações lagrangiana e hamiltoniana, a identificação de constantes de movimento por simetria, a análise de potencial efetivo e a integração numérica de sistemas dinâmicos. O método numérico escolhido, o integrador de Störmer-Verlet, tem, além disso, uma conexão histórica com o próprio problema: como discutiremos na Seção 3, Störmer foi um dos primeiros a usar esquemas de diferenças finitas desse tipo, precisamente para calcular as trajetórias de partículas no campo dipolar.

Na Seção 2, derivamos as equações de movimento a partir da força de Lorentz e da formulação lagrangiana, e introduzimos as grandezas adimensionais utilizadas na simulação. Na Seção 3, apresentamos o método de Störmer-Verlet, sua origem histórica e suas propriedades como integrador simplético. Na Seção 4, discutimos os resultados numéricos para o caso equatorial, o caso tridimensional e o caso vinculado à esfera. Na Seção 5, apresentamos as conclusões.

2. Formulação do Problema

2.1. Equações de movimento

Consideremos uma partícula não relativística de massa $m$ e carga $q$ movendo-se sob a ação de um campo magnético $\mathit{B}$. A equação de movimento, dada pela força de Lorentz, é

(1) $$m\ddot{\mathit{r}}=q\left(\dot{\mathit{r}}\times \mathit{B}\right).$$

Modelamos o campo magnético terrestre como o campo gerado por um dipolo magnético ${\mathit{\mu }}_z=\mu {\mathit{e}}_z$ posicionado no centro da Terra, com o eixo do dipolo alinhado ao eixo de rotação. O potencial vetor associado a esse dipolo é

(2) $$\mathit{A}=\frac{{\mu }_0}{4\pi }\frac{{\mathit{\mu }}_z\times {\mathit{e}}_r}{r^2}=M_z\frac{1}{r^3}\left(-y{\mathit{e}}_x+x{\mathit{e}}_y\right),$$

onde ${\mu }_0$ é a permeabilidade magnética do vácuo, $M_z={\mu }_0\mu /\left(4\pi \right)$ é o momento de dipolo magnético efetivo e $r=\sqrt{x^2+y^2+z^2}$. O campo magnético é então obtido por $\mathit{B}=\nabla \times \mathit{A}$.

Substituindo o campo resultante na (equação 1), obtemos as equações de movimento em coordenadas cartesianas:

(3) $$\{\begin{array}{cc}\ddot{x}=3\alpha \frac{z}{r^5}\left(\dot{y}z-\dot{z}y\right)-\alpha \dot{y}\frac{1}{r^3}\hfill & \\ \ddot{y}=-3\alpha \frac{z}{r^5}\left(\dot{x}z-\dot{z}x\right)+\alpha \dot{x}\frac{1}{r^3}\hfill & \\ \ddot{z}=3\alpha \frac{z}{r^5}\left(\dot{x}y-\dot{y}x\right)\hfill & \end{array}$$

onde a constante $\alpha =q{M}_z/m$ condensa as propriedades da partícula e do dipolo. Para o campo magnético terrestre, o momento de dipolo efetivo vale $M_z\approx \mathrm{8,2}\times {10}^{\text{15}}\text{T}\cdot {\text{m}}^3$, e os valores de $\alpha$ para elétrons e prótons são

(4) $$\alpha =\{\begin{array}{cc}-\mathrm{1,44}\times {10}^{\text{27}}{\text{m}}^3/\text{s}\hfill & \text{(elétrons)}\hfill \\ \mathrm{7,87}\times {10}^{\text{23}}{\text{m}}^3/\text{s}\hfill & \text{(prótons)}\hfill \end{array}$$

2.2. Adimensionalização

Para trabalhar com grandezas de ordem unitária, realizamos uma mudança de escala usando o raio terrestre $r_0=\mathrm{6\; 378\; 136}\text{m}$ como unidade de comprimento. Definimos as coordenadas adimensionais $X=x/r_0$, $Y=y/r_0$ e $Z=z/r_0$, de modo que $R=\sqrt{X^2+Y^2+Z^2}$ é a distância do centro da Terra à partícula, em unidades de raios terrestres. As equações de movimento tornam-se

(5) $$\{\begin{array}{cc}\ddot{X}=3{\alpha }_1\frac{Z}{R^5}\left(\dot{Y}Z-\dot{Z}Y\right)-{\alpha }_1\dot{Y}\frac{1}{R^3}\hfill & \\ \ddot{Y}=-3{\alpha }_1\frac{Z}{R^5}\left(\dot{X}Z-\dot{Z}X\right)+{\alpha }_1\dot{X}\frac{1}{R^3}\hfill & \\ \ddot{Z}=3{\alpha }_1\frac{Z}{R^5}\left(\dot{X}Y-\dot{Y}X\right)\hfill & \end{array}$$

com a constante adimensionalizada

(6) $${\alpha }_1=\frac{\alpha }{r_0^{\text{3}}}=\{\begin{array}{cc}-\mathrm{5,57}\times {10}^{\text{6}}{\text{s}}^{-1}\hfill & \text{(elétrons)}\hfill \\ \mathrm{3,04}\times {10}^{\text{3}}{\text{s}}^{-1}\hfill & \text{(prótons)}\hfill \end{array}$$

Observe que a diferença de três ordens de grandeza entre ${\alpha }_1$ para elétrons e prótons reflete a diferença de massa entre as duas partículas, e implica que os elétrons realizam oscilações muito mais rápidas.

2.3. Formulações lagrangiana e hamiltoniana

As (equações 5) descrevem completamente o movimento em coordenadas cartesianas, mas não tornam explícitas nem a simetria do problema nem as constantes de movimento associadas a ela. Por isso, passamos a uma formulação lagrangiana em coordenadas adaptadas a essa simetria. O problema possui simetria cilíndrica em torno do eixo do dipolo. É conveniente, portanto, utilizar coordenadas cilíndricas $(\rho ,\varphi ,Z)$, onde $\rho =\sqrt{X^2+Y^2}$. Nessas coordenadas, a lagrangiana do sistema é

(7) $$L=\frac{m}{2}\left({\dot{\rho }}^2+{\rho }^2{\dot{\varphi }}^2+{\dot{Z}}^2\right)+m{\alpha }_1\frac{{\rho }^2}{R^3}\dot{\varphi },$$

onde $R=\sqrt{{\rho }^2+Z^2}$. O último termo é a contribuição do potencial vetor, e acopla o grau de liberdade azimutal ao campo magnético. A derivação completa da lagrangiana a partir do formalismo de partícula em campo eletromagnético pode ser encontrada em [¹³, ¹⁴].

O sistema possui três momentos canônicos: $p_{\rho }=m\dot{\rho }$, $p_Z=m\dot{Z}$ e $p_{\varphi }$. Como a lagrangiana não depende explicitamente de $\varphi$, o momento canônico conjugado a $\varphi$ é conservado:

(8) $$p_{\varphi }=\frac{\partial L}{\partial \dot{\varphi }}=m{\rho }^2\dot{\varphi }+m{\alpha }_1\frac{{\rho }^2}{R^3}=c_2\text{(constante)}.$$

Essa quantidade não é o momento angular mecânico usual, pois inclui uma contribuição do potencial vetor. A constante $c_2$ depende das condições iniciais (por meio de ${\rho }_0$, $Z_0$ e ${\dot{\varphi }}_0$) e permanece fixa ao longo de toda a trajetória.

A hamiltoniana correspondente é

(9) $$H=\frac{1}{2m}\left[p_{\rho }^2+p_Z^2+{\left(\frac{p_{\varphi }}{\rho }-m{\alpha }_1\frac{\rho }{R^3}\right)}^2\right],$$

que coincide com a energia cinética da partícula (a força magnética não realiza trabalho). Como $H$ não depende explicitamente do tempo, a energia total é a segunda constante de movimento do sistema.

As equações de Hamilton fornecem as equações de movimento em coordenadas cilíndricas:

(10) $$\{\begin{array}{cc}\ddot{\rho }=\frac{c_2^{\text{2}}}{m^2{\rho }^3}+3{\alpha }_1^{\text{2}}\frac{{\rho }^3}{R^8}-{\alpha }_1^{\text{2}}\frac{\rho }{R^6}-3\frac{{\alpha }_1{c}_2}{m}\frac{\rho }{R^5}\hfill & \\ \ddot{Z}=3{\alpha }_1^{\text{2}}\frac{{\rho }^{2}Z}{R^8}-3\frac{{\alpha }_1{c}_2}{m}\frac{Z}{R^5}\hfill & \end{array}$$

enquanto a coordenada $\varphi$ evolui de acordo com

(11) $$\dot{\varphi }=\frac{c_2}{m{\rho }^2}-\frac{{\alpha }_1}{R^3}.$$

Essas equações constituem um sistema autônomo, ou seja, não dependem explicitamente do tempo nas variáveis $(\rho ,\dot{\rho },Z,\dot{Z})$; a coordenada $\varphi$ pode ser recuperada a posteriori pela integração da (equação 11).

2.4. Potencial efetivo e regiões proibidas

Para evidenciar o mecanismo de confinamento, separamos as contribuições da hamiltoniana (9) e reescrevemos a energia por unidade de massa na forma $H=T+V_{\text{eff}}$. Os termos que contêm os momentos $p_{\rho }$ e $p_Z$ dão a energia cinética do movimento no plano meridional, $T=\frac{1}{\text{2}}\left({\dot{\rho }}^2+{\dot{Z}}^2\right)$, enquanto o termo restante, que envolve $p_{\varphi }=c_2$ e depende apenas das coordenadas,

(12) $$V_{\text{eff}}(\rho ,Z)=\frac{1}{2m^2}{\left(\frac{c_2}{\rho }-m{\alpha }_1\frac{\rho }{R^3}\right)}^2$$

é o potencial efetivo, manifestamente não negativo por ser um quadrado perfeito, qualquer que seja o sinal de ${\alpha }_1$. Como $V_{\text{eff}}\ge 0$, a condição $T\ge 0$ impõe que

(13) $$V_{\text{eff}}(\rho ,Z)\le E,$$

onde $E$ é a energia total (também por unidade de massa). Os pontos do plano $(\rho ,Z)$ para os quais $V_{\text{eff}}=E$ definem as curvas de energia zero (ou curvas de Hill), que delimitam as regiões acessíveis ao movimento da partícula.

O sistema tem três graus de liberdade e apenas duas constantes de movimento independentes ($E$ e $c_2$); a integrabilidade no sentido de Liouville exigiria uma terceira, que não existe. De fato, Dilão e Alves-Pires [¹¹] demonstraram que o problema de Störmer é não integrável e exibe comportamento caótico para determinadas condições iniciais.

3. Método de Störmer-Verlet

3.1. Origem histórica

O método numérico utilizado neste trabalho foi redescoberto de forma independente ao menos três vezes ao longo de três séculos. A ideia de usar diferenças finitas de segunda ordem para integrar equações diferenciais do tipo $\ddot{q}=f\left(q\right)$ já aparece, de forma implícita, nos Principia de Newton (1687), na sua análise do movimento sob forças centrais [¹⁵]. Newton discretizava o tempo em intervalos iguais e calculava a posição seguinte a partir das duas posições anteriores e da aceleração no ponto atual, obtendo a relação

(14) $$q_{n+1}=2q_n-q_{n-1}+\text{Δ}{t}^{2}f\left(q_n\right),$$

que é a forma mais elementar do que hoje chamamos de método de Störmer-Verlet [¹⁶, ¹⁷].

Na sua forma explícita, esse esquema foi redescoberto e sistematizado por Carl Störmer no início do século XX, precisamente no contexto do problema que leva seu nome. Störmer necessitava de um método numérico para integrar as equações de movimento de partículas carregadas no campo dipolar, e os cálculos eram realizados à mão, linha por linha, em tabelas de diferenças finitas. Störmer utilizou versões com múltiplos passos desse esquema e publicou seus resultados entre 1907 e os anos 1930 [², ⁴]. A eficiência e a estabilidade do método eram questões práticas: erros acumulados significavam semanas de cálculos perdidos.

A segunda redescoberta ocorreu em 1967, quando o físico francês Loup Verlet publicou um artigo sobre simulações de dinâmica molecular com o potencial de Lennard-Jones [¹⁸]. Verlet utilizou exatamente o esquema da (equação 14) para propagar as posições dos átomos ao longo do tempo, e o método passou a ser amplamente adotado na comunidade de dinâmica molecular, onde ficou conhecido simplesmente como “método de Verlet”. Verlet não fez referência ao trabalho de Störmer, e a conexão entre os dois só foi reconhecida posteriormente, quando Hairer, Lubich e Wanner consolidaram a história do método em seu tratado sobre integração geométrica [¹⁶]. Desde então, a denominação Störmer-Verlet passou a ser utilizada para dar crédito a ambas as tradições.

3.2. Formulação do método

A (equação 14) é conveniente quando se dispõe apenas das posições, mas no contexto hamiltoniano é preferível trabalhar com posições e momentos. A forma equivalente mais utilizada em mecânica clássica é a chamada leapfrog (ou velocity Verlet), que intercala atualizações de posição e momento. Dado um sistema hamiltoniano com coordenadas generalizadas $q$ e momentos conjugados $p$, o algoritmo avança a solução de um passo temporal $\text{Δ}t$ segundo o esquema

(15) $$\{\begin{array}{cc}{p}_{n+1/2}=p_n-\frac{\text{Δ}t}{2}\frac{\partial H}{\partial q}\left(q_n\right)\hfill & \\ q_{n+1}=q_n+\text{Δ}t\frac{\partial H}{\partial p}\left(p_{n+1/2}\right)\hfill & \\ p_{n+1}=p_{n+1/2}-\frac{\text{Δ}t}{2}\frac{\partial H}{\partial q}\left(q_{n+1}\right)\hfill & \end{array}$$

O esquema é simétrico: a atualização do momento é dividida em dois meios-passos, um antes e um depois da atualização da posição. O método é explícito, de segunda ordem em $\text{Δ}t$ e reversível no tempo [¹⁶].

3.3. Propriedades simpléticas e conservação da energia

Por que preferir o método de Störmer-Verlet a outros integradores de mesma ordem, como o Runge-Kutta de segunda ordem? A razão é que o Störmer-Verlet preserva a estrutura simplética do espaço de fase. Em termos concretos, o mapa $(q_n,p_n)\mapsto (q_{n+1},p_{n+1})$ definido pelo esquema (15) é uma transformação canônica, ou seja, preserva o volume no espaço de fase e as relações entre coordenadas e momentos que caracterizam a mecânica hamiltoniana.

Uma consequência direta dessa propriedade é o que se chama de análise de erro retroativo (backward error analysis) [¹⁶, ¹⁷]: pode-se demonstrar que a sequência de pontos $(q_n,p_n)$ gerada pelo método de Störmer-Verlet é a solução exata de uma hamiltoniana modificada

(16) $$\tilde{H}=H+\text{Δ}{t}^2{H}_2+\text{Δ}{t}^4{H}_4+\cdots ,$$

onde $H_2,H_4,\dots$ são funções que dependem de $H$ e de suas derivadas. Em outras palavras, o método não resolve exatamente a hamiltoniana original, mas resolve exatamente uma hamiltoniana que difere da original por termos de ordem $\text{Δ}{t}^2$. Como $\tilde{H}$ é uma constante de movimento da dinâmica numérica, a energia $H$ oscila em torno do seu valor verdadeiro com amplitude proporcional a $\text{Δ}{t}^2$, mas não apresenta o desvio secular (drift), o crescimento monotônico e cumulativo do erro de energia ao longo do tempo, que é característico de métodos não simpléticos.

A diferença fica evidente em simulações longas. Com o método de Runge-Kutta de quarta ordem, por exemplo, a energia total pode crescer ou decrescer monotonicamente ao longo do tempo, mesmo para passos temporais pequenos. A taxa desse desvio é lenta, mas em integrações de ${10}^{\text{6}}$ ou mais passos, como as que realizamos neste trabalho, o efeito acumulado pode ser significativo. O método de Störmer-Verlet, por outro lado, mantém a energia confinada em uma faixa estreita, o que o torna adequado para simulações de longa duração em sistemas conservativos [¹⁶].

Vale observar que ser simplético não implica conservar a energia exatamente: a energia oscila, e a amplitude da oscilação diminui com o passo temporal. Para sistemas em que a conservação exata da energia é mais importante do que a preservação da geometria do espaço de fase, existem integradores que conservam a energia por construção, mas esses métodos sacrificam o caráter simplético e podem distorcer outras propriedades qualitativas da dinâmica [¹⁶].

3.4. Implementação

Em nossa implementação, escrita em linguagem C, as variáveis dinâmicas são $(\rho ,p_{\rho },Z,p_Z)$, com $p_{\rho }$ e $p_Z$ definidos na Seção 2.3. As derivadas parciais de $H$ em relação às coordenadas fornecem as forças generalizadas que aparecem nas (equações 10):

(17) $$-\frac{\partial H}{\partial \rho }=\frac{c_2^{\text{2}}}{m^2{\rho }^3}+3{\alpha }_1^{\text{2}}\frac{{\rho }^3}{R^8}-{\alpha }_1^{\text{2}}\frac{\rho }{R^6}-3\frac{{\alpha }_1{c}_2}{m}\frac{\rho }{R^5},$$

(18) $$-\frac{\partial H}{\partial Z}=3{\alpha }_1^{\text{2}}\frac{{\rho }^{2}Z}{R^8}-3\frac{{\alpha }_1{c}_2}{m}\frac{Z}{R^5},$$

e as derivadas em relação aos momentos dão as velocidades, $\dot{\rho }=p_{\rho }/m$ e $\dot{Z}=p_Z/m$. A coordenada $\varphi$ não participa do esquema simplético e é atualizada separadamente a cada passo temporal via (equação 11). A transformação final de coordenadas cilíndricas para cartesianas $(x,y,z)$ é feita a posteriori para a visualização das trajetórias. Todas as simulações apresentadas neste trabalho foram realizadas com o integrador de Störmer-Verlet; não empregamos outros integradores. A comparação qualitativa com métodos não simpléticos, no que se refere à conservação de energia, está discutida na Seção 3.3.

4. Resultados

4.1. Caso equatorial

Consideramos inicialmente o movimento restrito ao plano equatorial, ou seja, $Z\left(t\right)=0$ para todo $t\ge 0$. Essa restrição elimina a equação em $Z$ e simplifica o potencial efetivo: com $R=\rho$, temos

(19) $$V_{\text{eff}}\left(\rho \right)=\frac{c_{20}^{\text{2}}}{2m^2{\rho }^2}-\frac{{\alpha }_1{c}_{20}}{m{\rho }^3}+\frac{{\alpha }_1^{\text{2}}}{2{\rho }^4},$$

onde $c_{20}$ denota o valor da constante $c_2$ no plano equatorial. A equação de movimento para a coordenada radial reduz-se a

(20) $$\ddot{\rho }=\frac{c_{20}^{\text{2}}}{m^2{\rho }^3}-3\frac{{\alpha }_1{c}_{20}}{m{\rho }^4}+2\frac{{\alpha }_1^{\text{2}}}{{\rho }^5},$$

que corresponde a $\ddot{\rho }=-\partial V_{\text{eff}}/\partial \rho$. A coordenada angular é determinada pela relação

(21) $$\dot{\varphi }=\frac{c_{20}}{m{\rho }^2}-\frac{{\alpha }_1}{{\rho }^3}.$$

Quando $c_{20}>0$, a derivada $\partial V_{\text{eff}}/\partial \rho =0$ possui duas raízes positivas:

(22) $${\rho }_1=\frac{m{\alpha }_1}{c_{20}},{\rho }_2=2\frac{m{\alpha }_1}{c_{20}}.$$

A verificação da derivada segunda mostra que ${\rho }_1$ é um ponto de mínimo local e ${\rho }_2$ é um ponto de máximo local do potencial efetivo. O valor do potencial no máximo,

(23) $$V_{\text{max}}=V_{\text{eff}}\left({\rho }_2\right)=\frac{c_{20}^{\text{4}}}{32m^4{\alpha }_1^{\text{2}}},$$

determina a energia limiar para o confinamento: partículas com energia total $E<V_{\text{max}}$ ficam aprisionadas na vizinhança de ${\rho }_1$, oscilando entre dois pontos de retorno. Partículas com energia superior a $V_{\text{max}}$ podem escapar.

A Figura 1 ilustra o potencial efetivo no caso equatorial para $c_{20}>0$. A forma do poço de potencial mostra que o confinamento radial é assimétrico: a barreira é mais íngreme para $\rho <{\rho }_1$ (em direção à Terra) do que para $\rho >{\rho }_1$ (em direção ao espaço exterior).

Figura 1
Potencial efetivo para o caso equatorial com $c_{20}>0$. Os pontos ${\rho }_1$ e ${\rho }_2$ (indicados pelas linhas tracejadas) correspondem ao mínimo e ao máximo locais, respectivamente. Partículas com energia total inferior a $V_{\text{eff}}\left({\rho }_2\right)$ permanecem confinadas entre os pontos de retorno.

Para verificar numericamente essas previsões, simulamos o movimento de prótons no plano equatorial usando o método de Störmer-Verlet com passo temporal $\text{Δ}t=0.0001$. A constante de movimento $c_{20}$ é determinada pelas condições iniciais através da (equação 8), avaliada em $Z=0$ (onde $R=\rho$):

(24) $$c_{20}=m{\rho }_0^{\text{2}}{\dot{\varphi }}_0+m{\alpha }_1\frac{1}{{\rho }_0}.$$

Para o primeiro conjunto de condições iniciais, por exemplo, obtém-se $c_{20}=\mathrm{1,844}\times {10}^{-24}$ (em unidades de kg$\cdot$raios terrestres${}_{}{}^2\cdot$s^-1). Os três conjuntos foram escolhidos para ilustrar diferentes regimes dinâmicos:

• Próton 1: $\rho \left(0\right)=3.0$, $\dot{\rho }\left(0\right)=10.0$, $\varphi \left(0\right)=0.0$, $\dot{\varphi }\left(0\right)=10.0$.

• Próton 2: $\rho \left(0\right)=3.0$, $\dot{\rho }\left(0\right)=80.0$, $\varphi \left(0\right)=3\pi /2$, $\dot{\varphi }\left(0\right)=10.0$.

• Próton 3: $\rho \left(0\right)=3.0$, $\dot{\rho }\left(0\right)=100.0$, $\varphi \left(0\right)=\pi$, $\dot{\varphi }\left(0\right)=10.0$.

Os resultados estão apresentados na Figura 2. Os prótons 1 e 2, com energia total menor do que $V_{\text{max}}$, permanecem confinados em órbitas oscilatórias ao redor da Terra. A amplitude da oscilação radial é maior para o próton 2, cuja velocidade radial inicial é mais alta. Já o próton 3, com energia cinética inicial suficiente para ultrapassar a barreira de potencial, é espalhado e escapa.

Figura 2
Trajetórias de três prótons no plano equatorial, vistas de cima (plano $xy$), com $\text{Δ}t=0.0001$ e ${\alpha }_1=3.04\times {10}^{\text{3}}$s^-1. Condições iniciais: Próton 1 (${\rho }_0=3.0$, ${\dot{\rho }}_0=10.0$, ${\varphi }_0=0$, ${\dot{\varphi }}_0=10.0$); Próton 2 (${\rho }_0=3.0$, ${\dot{\rho }}_0=80.0$, ${\varphi }_0=3\pi /2$, ${\dot{\varphi }}_0=10.0$); Próton 3 (${\rho }_0=3.0$, ${\dot{\rho }}_0=100.0$, ${\varphi }_0=\pi$, ${\dot{\varphi }}_0=10.0$). Os prótons 1 e 2, com energia abaixo de $V_{\text{max}}$, ficam confinados em órbitas ao redor da Terra (disco cinza de raio $R=1$). O próton 3, com energia acima do limiar, escapa. Cada próton é identificado por uma cor: Próton 1 (azul), Próton 2 (verde) e Próton 3 (laranja).

4.2. Caso tridimensional

No caso geral, as equações de movimento (10) devem ser resolvidas simultaneamente para $\rho \left(t\right)$ e $Z\left(t\right)$. A presença da coordenada $Z$ modifica o potencial efetivo (12), que agora depende de duas variáveis. A estrutura de regiões permitidas e proibidas torna-se mais complexa: as curvas de Hill no plano $(\rho ,Z)$ definem faixas ao redor do equador magnético onde a partícula pode se mover.

Para a simulação tridimensional, utilizamos o método de Störmer-Verlet com passo temporal $\text{Δ}t=0.0002$. A constante $c_2$ é calculada a partir das condições iniciais pela (equação 8); para o primeiro conjunto de condições iniciais (Próton 1), obtém-se $c_2=\mathrm{1,776}\times {10}^{-24}$. Três conjuntos de condições iniciais foram utilizados:

• Próton 1: $\rho \left(0\right)=3.0$, $\dot{\rho }\left(0\right)=10.0$, $\varphi \left(0\right)=0.0$, $\dot{\varphi }\left(0\right)=10.0$, $Z\left(0\right)=0.5$, $\dot{Z}\left(0\right)=0.0$.

• Próton 2: $\rho \left(0\right)=3.0$, $\dot{\rho }\left(0\right)=80.0$, $\varphi \left(0\right)=\pi$, $\dot{\varphi }\left(0\right)=10.0$, $Z\left(0\right)=0.5$, $\dot{Z}\left(0\right)=0.0$.

• Próton 3: $\rho \left(0\right)=3.0$, $\dot{\rho }\left(0\right)=100.0$, $\varphi \left(0\right)=5.26$, $\dot{\varphi }\left(0\right)=10.0$, $Z\left(0\right)=0.5$, $\dot{Z}\left(0\right)=0.0$.

A Figura 3 mostra as trajetórias resultantes projetadas no plano meridional $(\rho ,Z)$ e vistas de cima no plano $xy$. Ao contrário do caso equatorial, onde a trajetória confinada é oscilatória, no caso tridimensional a partícula oscila simultaneamente na direção radial e na direção $Z$, percorrendo a região toroidal em torno da Terra sem repetir a mesma trajetória. A partícula permanece confinada pela barreira de potencial, mas sua trajetória no plano meridional preenche densamente a região acessível. A depender das condições iniciais, a órbita pode ser quase-periódica ou caótica, uma consequência direta da não integrabilidade do sistema.

Figura 3
Trajetórias de três prótons no caso tridimensional, com $\text{Δ}t=0.0002$ e ${\alpha }_1=3.04\times {10}^{\text{3}}$s^-1. Painel (a): projeção no plano meridional $(\rho ,Z)$; painel (b): vista de cima no plano $xy$. As cores distinguem os prótons: Próton 1 (azul), Próton 2 (verde) e Próton 3 (laranja); a Terra corresponde à região $R\le 1$, em cinza. Condições iniciais: Próton 1 (${\rho }_0=3.0$, ${\dot{\rho }}_0=10.0$, ${\varphi }_0=0$, ${\dot{\varphi }}_0=10.0$, $Z_0=0.5$, ${\dot{Z}}_0=0$); Próton 2 (${\rho }_0=3.0$, ${\dot{\rho }}_0=80.0$, ${\varphi }_0=\pi$, ${\dot{\varphi }}_0=10.0$, $Z_0=0.5$, ${\dot{Z}}_0=0$); Próton 3 (${\rho }_0=3.0$, ${\dot{\rho }}_0=100.0$, ${\varphi }_0=5.26$, ${\dot{\varphi }}_0=10.0$, $Z_0=0.5$, ${\dot{Z}}_0=0$). O painel meridional mostra que as partículas confinadas oscilam tanto radialmente quanto na direção $Z$, preenchendo a região permitida. A estrutura toroidal resultante é uma representação aproximada do cinturão de radiação de Van Allen.

A região toroidal onde as partículas permanecem confinadas corresponde, em nosso modelo simplificado, ao cinturão de radiação de Van Allen. De fato, na natureza, o cinturão de Van Allen é formado por prótons e elétrons de diferentes energias, e sua estrutura detalhada depende de efeitos que não foram incluídos neste modelo, como variações temporais do campo geomagnético, interações com o vento solar e efeitos relativísticos. Ainda assim, o modelo dipolar reproduz corretamente o mecanismo básico de confinamento magnético.

Vale notar que, durante a oscilação latitudinal, a trajetória da partícula confinada pode atingir a superfície terrestre (representada pela esfera de raio $R=1$ em nossas unidades). Partículas que chegam a altitudes suficientemente baixas interagem com a atmosfera, depositam energia e produzem emissões luminosas: as auroras polares. Esse fenômeno ocorre preferencialmente em latitudes altas, onde as linhas do campo magnético do dipolo convergem. Uma modelagem mais detalhada exigiria efeitos atmosféricos e, possivelmente, uma formulação vinculada à superfície terrestre, como a que apresentamos a seguir.

Do ponto de vista físico, o caráter caótico de parte das trajetórias tem um papel concreto na estrutura do cinturão. Órbitas regulares (quase-periódicas) permanecem restritas a superfícies bem definidas no espaço de fase, enquanto órbitas caóticas exploram de maneira mais uniforme toda a região permitida pelas curvas de Hill, contribuindo para o preenchimento difuso da faixa toroidal observado na Figura 3. Além disso, a sensibilidade às condições iniciais estabelece uma fronteira efetiva entre partículas duradouramente aprisionadas e partículas que, após muitas oscilações, alcançam baixas altitudes e precipitam na atmosfera, ou seja, o mecanismo por trás das auroras. Uma caracterização quantitativa dessa fronteira, por meio de seções de Poincaré ou expoentes de Lyapunov, foge ao escopo deste trabalho e fica como proposta para investigação futura.

4.3. Caso vinculado: movimento sobre umasuperfície esférica

O problema de Störmer na forma livre, discutido nas seções anteriores, é não integrável. Existe, porém, uma variante do problema que admite tratamento analítico: o caso em que a partícula carregada é restrita a se mover sobre uma superfície esférica de raio $R_s$ centrada no dipolo magnético [¹⁹, ²⁰]. Essa formulação, proposta por Cortés e Cortés Poza [¹⁹], modela de forma simplificada a interação de partículas carregadas com a atmosfera terrestre em altitudes fixas. O interesse teórico dessa formulação é que o vínculo holonômico ($R=R_s$) torna integrável um sistema que, sem ele, não o é: ao eliminar um grau de liberdade, o vínculo faz com que as duas constantes de movimento passem a ser suficientes para a integrabilidade.

4.3.1. Formulação sobre a esfera

Com a partícula restrita à esfera $R=R_s$, restam dois graus de liberdade: o ângulo polar $\theta$ (medido a partir do polo norte) e o ângulo azimutal $\phi$. O potencial vetor do dipolo magnético, restrito à superfície, tem apenas componente azimutal:

(25) $$\mathit{A}=-\frac{{\mu }_0\mu \sin \theta }{4\pi R_s^2}{\mathit{e}}_{\phi },$$

onde $\mu$ é o módulo do momento de dipolo magnético, o mesmo que aparece na definição de $M_z$ na (equação 2). A lagrangiana do sistema é

(26) $$L=\frac{1}{\text{2}}M{R}_s^2\left({\dot{\theta }}^2+{\dot{\phi }}^2{\sin }^2\theta \right)-k_s\dot{\phi }{\sin }^2\theta ,$$

onde $M$ é a massa da partícula e $k_s=q{\mu }_0\mu /\left(4\pi R_s\right)$ é uma constante com dimensão de momento angular que condensa as propriedades do dipolo e da partícula. O último termo acopla o campo magnético ao movimento azimutal.

Como no caso livre, a coordenada $\phi$ é cíclica, e o momento canônico conjugado

(27) $$p_{\phi }=M{R}_s^2\dot{\phi }{\sin }^2\theta -k_s{\sin }^2\theta$$

é uma constante de movimento. A hamiltoniana coincide com a energia cinética (a força magnética não realiza trabalho) e fornece a segunda constante de movimento. Com dois graus de liberdade e duas constantes de movimento independentes, o sistema é integrável no sentido de Liouville [¹²].

4.3.2. Potencial efetivo e regimes de movimento

A hamiltoniana pode ser escrita como

(28) $$H=\frac{1}{2M{R}_s^2}\left[p_{\theta }^2+\frac{{\left(p_{\phi }+k_s{\sin }^2\theta \right)}^2}{{\sin }^2\theta }\right]=\frac{p_{\theta }^2}{2M{R}_s^2}+V_{\text{eff}}\left(\theta \right),$$

onde o potencial efetivo é

(29) $$V_{\text{eff}}\left(\theta \right)=\frac{1}{2M{R}_s^2}\frac{{\left(p_{\phi }+k_s{\sin }^2\theta \right)}^2}{{\sin }^2\theta }.$$

A estrutura desse potencial depende da relação entre $\left|p_{\phi }\right|$ e $k_s$. Quando $\left|p_{\phi }\right|<k_s$, o potencial exibe dois poços separados por uma barreira no equador ($\theta =\pi /2$), formando um potencial biestável. Nesse regime, a partícula pode ficar confinada em um dos hemisférios, oscilando entre dois ângulos de retorno em $\theta$ sem cruzar o equador. Quando $\left|p_{\phi }\right|\ge k_s$, a barreira desaparece e a partícula oscila entre os dois hemisférios.

Os polos ($\theta =0$ e $\theta =\pi$) são sempre repulsivos, pois $V_{\text{eff}}\to \infty$ quando $\theta \to 0$ ou $\pi$ (devido ao fator ${\csc }^2\theta$). A partícula jamais atinge os polos com energia finita.

4.3.3. Solução analítica e validação numérica

Com a mudança de variável $z=\cos \theta$, a equação de movimento para $\theta$ se reduz à equação de um oscilador em um potencial quártico [²⁰]:

(30) $${\dot{z}}^2=-b^2\left(z^2-A\right)\left(z^2-B\right),$$

onde $A$ e $B$ são constantes que dependem dos parâmetros adimensionais $a=p_{\phi }/\ell$ e $b=k_s/\ell$, com $\ell =\sqrt{2M{R}_s^{2}K}$ sendo o momento angular característico e $K$ a energia cinética (constante). Piña e Cortés [²⁰] demonstraram que a solução dessa equação se expressa em termos das funções elípticas de Jacobi:

(31) $$z\left(\tau \right)=\{\begin{array}{cc}\sqrt{B}\text{dn}(b\sqrt{B}\tau ,\kappa )\hfill & \text{se}A>0\text{(1 hemisfério)}\hfill \\ \sqrt{B}\text{cn}(b\sqrt{B-A}\tau ,\kappa )\hfill & \text{se}A<0\text{(2 hemisférios)}\hfill \end{array}$$

onde $\tau$ é o tempo adimensional, dn e cn são funções elípticas de Jacobi, e $\kappa$ é o módulo elíptico, determinado por $A$ e $B$.

Essa solução analítica permite validar diretamente o integrador de Störmer-Verlet. Utilizando os mesmos parâmetros do artigo original [²⁰] ($M=2$, $R_s=10$, $k_s=0.5$, em unidades arbitrárias), simulamos numericamente o movimento da partícula sobre a esfera e comparamos as trajetórias com a solução elíptica exata. As Figuras 4–6 mostram essa comparação para três casos representativos que cobrem os diferentes regimes do potencial efetivo (29):

Figura 4
Validação: caso banda no hemisfério norte (${\theta }_0=\pi /3$, $p_{\theta 0}=0$, $p_{\phi }=0.394$, $t=3000$, $\mathrm{1,5}\times {10}^{\text{7}}$ passos). Acima: trajetórias sobre a esfera (azul: Störmer-Verlet, vermelho: analítico). Abaixo: erro de posição normalizado $\Vert {\mathbf{r}}_{\text{SV}}-{\mathbf{r}}_{\text{an}}\Vert /R_s$.

Figura 5
Validação: caso travessia do equador (θ₀ = 0.6 rad, pθ₀ = 0.2525, p_φ = 0.25, t = 3000, 1,5 × 107 passos). A partícula cruza a barreira equatorial do potencial biestável e percorre ambos os hemisférios.

Figura 6
Validação: caso laços (${\theta }_0=\pi /4$, $p_{\theta 0}=0.3$, $p_{\phi }=-0.15$, $t=10000$, $5\times {10}^{\text{7}}$ passos). Quando $-k_s<p_{\phi }<0$, a trajetória exibe laços que não envolvem nenhum dos polos [¹⁹].

• Banda no hemisfério norte (${\theta }_0=\pi /3$, $p_{\theta 0}=0$, $p_{\phi }=0.394$): a partícula oscila em $\theta$ dentro de uma faixa estreita no hemisfério norte, traçando uma banda sobre a esfera. Regime $A>0$, solução via dn.

• Travessia do equador (${\theta }_0=0.6$ rad, $p_{\theta 0}=0.2525$, $p_{\phi }=0.25$): a partícula possui energia suficiente para cruzar a barreira equatorial do potencial biestável e percorre ambos os hemisférios. Regime $A<0$, solução via cn.

• Laços (${\theta }_0=\pi /4$, $p_{\theta 0}=0.3$, $p_{\phi }=-0.15$): quando $-k_s<p_{\phi }<0$, a trajetória exibe laços que não envolvem nenhum dos polos [¹⁹]. Regime $A>0$, solução via dn.

Em cada figura, o painel superior sobrepõe as trajetórias do integrador numérico (azul) e da solução analítica (vermelho) sobre a esfera. O painel inferior mostra o erro de posição normalizado $\Vert {\mathbf{r}}_{\text{SV}}\left(t\right)-{\mathbf{r}}_{\text{an}}\left(t\right)\Vert /R_s$ ao longo da simulação. Nos três casos, esse erro oscila periodicamente na faixa de ${10}^{-8}$ a ${10}^{-7}$, sem crescimento secular. O erro relativo de energia, por sua vez, permanece no nível da precisão aritmética de ponto flutuante ($\sim {10}^{-15}$). Esse comportamento é esperado de um integrador simplético (Seção 3.3): conforme a análise de erro retroativo [¹⁶, ¹⁷], o algoritmo de Störmer-Verlet resolve exatamente uma hamiltoniana modificada $\tilde{H}=H+𝒪\left(\text{Δ}{t}^2\right)$. A órbita numérica segue, portanto, uma trajetória fisicamente consistente, ligeiramente deslocada da original, mas sem acúmulo de erro.

A existência dessa solução analítica revela ainda uma conexão entre o problema de Störmer e a dinâmica de corpos rígidos. Piña e Cortés [²⁰] mostraram que as soluções em termos de dn e cn têm a mesma estrutura matemática das soluções do pião de Euler (rotação livre de um corpo rígido assimétrico) e do pião de Lagrange (corpo rígido simétrico com gravidade). Essa correspondência não é acidental: em ambos os problemas, as equações de movimento se reduzem a uma equação diferencial com potencial quártico em $z=\cos \theta$, cuja solução é dada pelas mesmas classes de funções elípticas. Partindo de um contexto eletromagnético, o problema de Störmer sobre a esfera revela, assim, uma ponte inesperada com a dinâmica de corpos rígidos e pode servir como porta de entrada para as funções elípticas de Jacobi.

5. Conclusões

O problema de Störmer reúne, num mesmo sistema, tópicos que os cursos de graduação costumam tratar em separado: as formulações lagrangiana e hamiltoniana, constantes de movimento, potencial efetivo, integração numérica de sistemas dinâmicos e, no caso vinculado à esfera, funções especiais. Neste artigo, percorremos essas etapas da derivação analítica à simulação com o integrador simplético de Störmer-Verlet, incluindo a validação do método numérico contra uma solução exata.

A análise do caso equatorial mostrou como o potencial efetivo controla o confinamento radial: partículas com energia abaixo do máximo do potencial ficam presas em órbitas oscilatórias, enquanto partículas com energia acima desse limiar escapam. As simulações confirmaram esse critério e permitiram visualizar os dois regimes.

O caso tridimensional acrescentou a oscilação latitudinal, e as simulações mostraram a formação de trajetórias confinadas em uma região toroidal ao redor da Terra, o que corresponde, em nosso modelo idealizado, ao cinturão de radiação de Van Allen. A sensibilidade às condições iniciais, característica da não integrabilidade do sistema, abre espaço para investigações sobre o caráter caótico das trajetórias, cuja caracterização mais precisa, via seções de Poincaré ou expoentes de Lyapunov, fica como proposta para um trabalho futuro.

O caso vinculado à esfera mostrou como a imposição de um vínculo holonômico transforma o sistema não integrável em um sistema integrável, com solução analítica em funções elípticas de Jacobi. A comparação direta entre a solução analítica e a integração numérica pelo método de Störmer-Verlet validou o integrador: ao longo de simulações com até $5\times {10}^{\text{7}}$ passos temporais, o erro de posição permaneceu na faixa de ${10}^{-8}$ a ${10}^{-7}$ do raio da esfera, sem crescimento secular.

O código-fonte das simulações, tanto para o caso livre quanto para o caso esférico, foi implementado em linguagem C e está disponível em repositório público.¹

Agradecimentos

Este trabalho foi desenvolvido no âmbito do Programa Institucional de Bolsas de Iniciação Científica (PIBIC) da Universidade Federal Fluminense, com bolsa financiada pelo CNPq.

Disponibilidade de Dados

Todo o conjunto de dados que dá suporte aos resultados deste estudo foi publicado no próprio artigo.

Referências

• ^[1] K. Birkeland, The Norwegian Aurora Polaris Expedition 1902–1903 (H. Aschehoug & Co., Christiania, 1908).
• ^[2] C. Störmer, Zeitschrift für Astrophysik 1, 237 (1930).
• ^[3] C. Störmer, Terrestrial Magnetism and Atmospheric Electricity 35, 193 (1930).
• ^[4] C. Störmer, The Polar Aurora (Clarendon Press, Oxford, 1955).
• ^[5] J. Clay e H.P. Berlage, Naturwissenschaften 20, 687 (1932).
• ^[6] A.H. Compton, Physical Review 41, 111 (1932).
• ^[7] J.A. Van Allen e L.A. Frank, Nature 183, 430 (1959).
• ^[8] C.I. Underwood, D.J. Brock, P.S. Williams, S. Kim, R. Dilão, P. Ribeiro Santos, M.C. Brito, C.S. Dyer e A.J. Sims, IEEE Transactions on Nuclear Science 41, 2353 (1994).
• ^[9] E.J. Daly, ESA Bulletin 12, 229 (1988).
• ^[10] E.G. Stassinopoulos e J.P. Raymond, Proceedings of the IEEE 76, 1423 (1988).
• ^[11] R. Dilão e R. Alves-Pires, em: Differential Equations, Chaos and Variational Problems, editado por V. Staicu (Springer, Basel, 2007).
• ^[12] J.V. José e E.J. Saletan, Classical Dynamics: A Contemporary Approach (Cambridge University Press, Cambridge, 1998).
• ^[13] H. Goldstein, C. Poole e J. Safko, Classical Mechanics (Addison-Wesley, San Francisco, 2002), 3 ed.
• ^[14] J.R. Taylor, Classical Mechanics (University Science Books, Sausalito, 2005).
• ^[15] I. Newton, Philosophiæ Naturalis Principia Mathematica (Royal Society, Londres, 1687).
• ^[16] E. Hairer, C. Lubich e G. Wanner, Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations (Springer, Berlin, 2006).
• ^[17] E. Hairer, C. Lubich e G. Wanner, Acta Numerica 12, 399 (2003).
• ^[18] L. Verlet, Physical Review 159, 98 (1967).
• ^[19] E. Cortés e D. Cortés Poza, European Journal of Physics 36, 055009 (2015).
• ^[20] E. Piña e E. Cortés, European Journal of Physics 37, 065009 (2016).

• 1
  https://github.com/Shnrqpdr/Stormer-Problem

Editado por

• Editor-Chefe:
  Marcello Ferreira https://orcid.org/0000-0003-4945-3169

Datas de Publicação

• Publicação nesta coleção
  24 Ago 2026
• Data do Fascículo
  2026

Histórico

• Recebido
  11 Abr 2026
• Revisado
  27 Jul 2026
• Aceito
  28 Jul 2026