Dinâmica Molecular de Partículas-Planares Interagentes

Grupo: Antônio Carlan, Gabriela Pereira, Renan Soares e Victor Gandara

Objetivo: O trabalho tem como objetivo descrever a dinâmica molecular de um conjunto de partículas-planares infinitas paralelas. o modelo simplificado para estudar o comportamento de muitos corpos interagentes ao longo de um único eixo espacial é usado para feixes de partículas carregadas e também para sistemas autogravitantes. O estudo desse tipo de interação é importante para ter insights de sistemas de maior dimensionalidades. Entre os fatores mais interessantes de serem analisados no espaço de fase da distribuição de partículas, caracterizado com posição no eixo x e momento ou velocidade no eixo y, temos a formação de aréolas no espaço de fase. Essas aréolas ocorrem quando distribuições iniciais acabam fornecendo energia superior a presente na média inicial de partículas devido a oscilações na posição do núcleo (parte da distribuição do espaço de fase que contém a grande maioria das partículas) e, por conta dessa formação, o feixe tende a perder qualidade e até número de partículas quando há contato com as paredes do acelerador.

O modelo

É possível criar um modelo simplificado de interação de partículas imaginando que todas elas são compostas por planos de mesma massa e sem espessura, paralelos entre si e com liberdade de movimento relevante no eixo perpendicular às superfícies das mesmas.

Tal modelo pode ser utilizado para estudar a dinâmica de um número elevado de partículas e conta com algumas particularidades simplificadoras de complexidade numérica. Alguns estudos utilizam esse modelo para estudar o comportamento de feixes carregados com a menor dimensionalidade possível e também para sistemas auto-gravitantes. As duas aplicações citadas serão as demonstradas a seguir mas o modelo é versátil para diferentes condições iniciais e de contorno.

Descrição do Modelo

Disposição espacial

Um número "N" de partículas-plano de mesma massa é disposto ao longo do eixo x com uma função de distribuição n(x) e com todas as partículas paralelas ao plano yz. No caso do modelo auto-gravitante, movimentos ao longo do eixo z não têm relevância prática mas para o modelo de um feixe de planos é definido que todos os planos se movem na direção z com mesma velocidade.

A consideração de movimento ao longo do eixo z para o modelo de feixe é importante para a análise do comportamento do espaço de fases ao longo do comprimento do acelerador de partículas. Aqui, podemos imaginar estudar casos em que o feixe gerado inicialmente no espaço de fase terá um comportamento constante no tempo para cada secção transversal em xy para uma coordenada em z pré-estabelecida.

Ausência de Efeitos de Colisão

O modelo leva em consideração que todas os planos possuem a mesma massa e que as colisões não levam a perda de energia no sistema.

Dessa forma podemos prever o que ocorre quando a posição de dois planos coincide pois há conservação no momento linear em x e o coeficiente de restituição (razão entre as velocidades de afastamento e aproximação dos corpos) é igual a 1.

${\begin{cases}mv_{1}-mv_{2}=-mv_{1}'+mv_{2}'\\\\{\frac {v_{1}'+v_{2}'}{v_{1}+v_{2}}}=1\end{cases}}$

O sistema de equações resulta em: $v_{1}'=v_{2}\ \ \ ;\ \ \ v_{2}'=v_{1}$

Dessa forma, é irrelevante ponderar se houve de fato uma colisão com troca de velocidades ou se as partículas simplesmente cruzam sem interferir de forma relevante. A última consideração é a levada em conta no ponto de vista computacional pois exclui a necessidade de atualizar a velocidade das partículas que passaram pela mesma posição no eixo x. Devido a igualdade de massas e também ao fato de que os encontros entre duas partículas se dão de forma sempre frontal, podemos de fato desconsiderar efeitos de colisão durante a modelagem da integração temporal do feixe.

Interação entre as partículas

Para analisar a interação entre as partículas é preciso entender como uma única partícula afeta seu entorno.

No caso das partículas-plano carregadas eletricamente podemos considerar os efeitos de campo elétrico de um isolante plano infinito:

Pela lei de Gauss, para uma dada distribuição de carga com densidade superficial $\sigma$:

$\oint {\overrightarrow {E}}.d{\overrightarrow {S}}={\frac {Q_{int}}{\epsilon _{0}}}$

$E.A-(-E.A)={\frac {\sigma .A}{\epsilon _{0}}}$

$E={\frac {\sigma }{2\epsilon _{0}}}$

O mesmo pode ser feito para um plano com densidade de massa distribuída uniformemente, para o caso do sistema auto-gravitante, com mudança na constante física multiplicativa e no sinal.

O fato importante a ser destacado aqui é a total independência do campo gerado pela partícula e a distância. Dessa forma, a interação entre as partículas em qualquer sistema que seguir esse modelo não terá dependência da distância entre elas e sim da posição relativa entre as partículas consideradas.

Essa consideração é razoável visto que as interações eletrostáticas de longo alcance possuem pouca influência relativa entre as cargas e reduz bastante o esforço computacional quando simulações são realizadas. Diferentes de sistemas de partículas que possuem força resultante dependente do quadrado da distância e que devem ser calculadas para todos os pares de partículas, aqui as forças sentidas pela interação eletrostática entre as partículas-plano depende exclusivamente da diferença do número de partículas de um lado e do outro, sendo necessário para realização dos cálculos numéricos apenas a implementação de uma função de ordenamento.

Interferência externa do sistema

Para o caso dos planos no sistema auto-gravitante não é necessária ação externa alguma para manter o equilíbrio em um sistema dinâmico estável. Sempre é levado em conta que todas as partículas têm velocidade em x inferior a velocidade de escape do restante do sistema e, por conta disso e das forças sempre atrativas de interação, o sistema não expande indefinidamente.

Já no caso de planos carregados eletricamente as forças de interação entre os planos serão sempre repulsivas e a energia total do sistema é positiva. Sem uma interação externa, mesmo um sistema composto por um número limitado de partículas de fraca interação expandiria indefinidamente. Por causa disso, nesse caso é considerado que há uma força colimadora no feixe de planos proporcional a distância da origem. A origem dessa força externa vem de modelos de maior dimensionalidade de aceleradores de partícula que contam com um campo magnético externo capaz de colimar um feixe de cargas radialmente, por não existir um paralelo direto para o modelo unidimensional consideramos apenas uma força restauradora que impede a expansão indefinida do sistema como um todo. Também é possível ajustar condições de contorno simulando as paredes externas do acelerador como placas condutoras infinitas e com potencial definido. Esse tipo de consideração pode ser feita para sistemas que não são simétricos e impede que as cargas se afastem demasiadamente da origem.

Comportamento de uma Partícula no Sistema

Graças a caracterização das partículas temos que o campo gerado por uma partícula-plano pode ser caracterizado (para ambos os sistemas citados) da seguinte forma:

$E_{i}(x)=-\lambda \Theta [-(x-x_{i})]+\lambda \Theta (x-x_{i})$

Em que a constante $\lambda$ pode ser tanto positiva (cargas positivas) quanto negativa (cargas negativas ou interação gravitacional) e a função $\Theta$ equivale a função degrau de Heaviside:

$\Theta (x)={\begin{cases}0,\ \ x<0\\1,\ \ x>0\end{cases}}$

Dessa forma, para exemplificação, um sistema de 5 partículas-plano teria uma interação da seguinte forma:

Aqui, F seria a força de interação entre duas partículas no caso de planos carregados. É possível notar que a força de interação varia de duas unidades dessa constante para cada partícula contada em uma direção e que a força máxima sentida (pelas partículas-plano localizadas nas extremidades) é proporcional a $(N-1)$.

É possível notar que para uma distribuição de partículas-plano do tipo:

$n(x)=\sum _{i=1}^{N}\delta (x_{i})$

com o índice $\ \ ''i''\$ sendo a ordem da partícula contada do sentido negativo para o positivo do eixo x, teremos uma força resultante de interação para uma partícula que pode ser definida por

$F_{i}=[-(N-1)+2(i-1)]F$

$F_{i}=(2i-N-1)F$

Essa força é suficiente para manter o sistema auto-gravitante coeso mas para as partículas carregadas devemos adicionar a força restauradora descrita anteriormente:

$F_{R,i}^{grav.}=(2i-N-1)F$

$F_{R,i}^{el.}=(2i-N-1)F-\kappa x_{i}$

Conservação da Energia do Sistema

Ao realizar uma simulação numérica com um grande número de entes que representam um sistema físico é importante que seja feita uma verificação sobre o quão apropriados os resultados da simulação de fato são. Uma das formas mais comuns de verificar se uma simulação numérica está correta na dinâmica molecular de vários corpos e efetuar periodicamente a soma da energia do sistema.

No caso de sistemas como os descritos, não há nenhum incremento de energia proveniente de uma fonte externa. Mesmo para o sistema eletricamente carregado, a força restauradora é potencial e não muda a energia total do sistema dada uma distribuição inicial de posição e momento.

Para todos os efeitos, a energia dos sistemas estudados pode ser dividida em três partes:

Energia Cinética

Todas as partículas da distribuição escolhida terão massa e velocidade, de forma que a energia cinética do sistema será dada por:

$En_{cin}={\frac {m}{2}}\sum _{i=1}^{N}{\dot {x_{i}}}^{2}$

Energia Potencial Restauradora

Observada no sistema eletricamente carregado citado, pode ser descrita como:

$En_{pot}={\frac {\kappa }{2}}\sum _{i=1}^{N}x_{i}^{2}$

Energia Potencial Intrínseca

Energia contida na interação entre as partículas. A energia de uma partícula depende do potencial na sua posição da sua carga (ou massa).

Como potenciais não são absolutos, é estabelecido um Potencial nulo de referência nas condições de contorno do sistema. Como pode ser visto no gif com as partículas-plano carregadas em azul, são colocadas placas condutoras planas nas posições $x=L$ e $x=-L$ (caso haja simetria no espaço de fase do sistema) e é estabelecido que os potenciais nessas posições são $V(-L)=V(L)=0$.

Uma maneira mais simples de determinar o potencial das N partículas é considerar que há nas posições de contorno, $x_{0}$ e $x_{N+1}$, partículas não interagentes.

A diferença de potencial entre duas posições vizinhas pode ser determinada por:

$\Delta V_{i}=(x_{i}-x_{i-1})E_{i}$

com $E_{i}$ sendo o módulo do campo entre a i-ésima e a (i - 1)-ésima partícula, dado por:

$E_{i}=(2i-N-1)\lambda$

O que resulta na energia de uma partícula específica (k-ésima) descrita por:

$En_{k}=\alpha \sum _{i=1}^{k}(x_{i}-x_{i-1})(2i-N-1)$

Aqui foram aglutinadas as possíveis constantes multiplicativas em uma constante $\alpha$ que pode ser escolhida para se adequar ao sistema estudado e a energia intrínseca do sistema é dada pela soma da energia de interação de cada uma das partículas. Como o sistema é dividido em intervalos com variação de potencial discretos, podemos escrever a energia intrínseca da seguinte forma:

$En_{int}=\alpha \sum _{k=1}^{N}\sum _{i=1}^{k}(x_{i}-x_{i-1})(2i-N-1)$

A soma dessas três energias deve se manter constante (ou com pouca variação relativa significativa, pelo menos) durante a simulação:

$E_{total}=En_{cin}+En_{pot}+En_{int}=cte.$

Modelagem Computacional

Criando a Distribuição do Sistema

O primeiro passo para criar a simulação das N partículas é definir as coordenadas fundamentais das partículas em si (posição e velocidade).

A simulação que será utilizada de exemplo é simétrica no espaço de fase, ou seja, para cada partícula com coordenadas $(x_{i},y_{i})$ existirá uma partícula espelhada com coordenadas $(-x_{i},-y_{i})$.

Esse artifício é utilizado para que o centro de massa mantenha-se na origem e que o momento linear do sistema em x seja sempre nulo. Com isso é garantido que haverá sempre, também, a simetria par entre as energias das partículas espelhadas e a simetria ímpar entre as forças resultantes.

Com isso garantido, basta fazer o cálculo numérico de uma metade da distribuição levando em conta a presença da outra metade, porém não tornando necessário que os cálculos sejam repetidos pois o comportamento de cada partícula espelhada torna-se perfeitamente previsível.

# Importando as bibliotecas necessárias 
from math import *
import matplotlib.pyplot as plt
import numpy as np

A=[]     #criamos uma lista vazia

N =  1000       # Determinamos o número de partículas utilizado

for i in range(int(N)):

    A.append([np.random.rand() ,2*(np.random.rand() - 1/2)]) #valores aleatórios de x entre 0 e 1 e de v entre -1 e 1 são adicionados


A = np.array(A)
X= np.array([])     #Arrays usados para exemplificar
V= np.array([])     #distribuição inicial

for i in range (len(A)):    #simetrização da distribuição
    X=np.append(X,A[i][0])
    V=np.append(V,A[i][1])
    V=np.append(V,-A[i][1])
    X=np.append(X,-A[i][0])
    
plt.xlabel('posição')             #plot da distribuição inicial
plt.ylabel('momento linear')
plt.title('Espaço de fase inicial')
plt.grid(True)

plt.scatter(X,V)

Função de Indexação

Em outras modelagens como a interação de cargas em sistemas bidimensionais ou de massas em tridimensionais é necessário calcular a força de interação entre partículas em função das distâncias. Na modelagem proposta não é necessário calcular distâncias mas sim saber o ordenamento das partículas ao longo do eixo $x$.

No entanto, é muito comum em simulações que utilizam métodos de integração utilizar vetores que possuam apenas uma coluna e $2N$ linhas, sendo a primeira metade de linhas correspondentes às posições e a outra metade correspondente às velocidades.

Para não precisar atualizar as linhas de posição e velocidade do vetor 2N-dimensional, utiliza-se uma função de indexação que será responsável por indicar qual é a ordem da partícula correspondente a uma determinada linha em um dado estado.

Assim, a posição da partícula é atualizada no tempo e seu índice também, porém a posição dela no vetor de fase nunca se modifica. Esse fato é importante para estudos que pretendem analisar a trajetória de um único corpo no tempo, como é feito nos diagramas de Poincaré, por exemplo.

teste = np.array([2,6,3,9,0,1,4])   #criado um vetor de teste  para que seus valores sejam indexados

def indexador(arr):                 #define-se a função
    
    ordenador = []                  #é criada uma lista vazia
    
    for i in range(len(arr)):
        
        ordenador.append([abs(arr[i]),i]) # para cada valor criamos uma lista com o valor absoluto e a ordem atual
        
    ordenador.sort()                      # é feito o ordenamento da lista 'ordenador' de acordo com o primeiro elemento
    
    for i in range(len(arr)):             

        ordenador[i].append(i)            # é adicionado um terceiro valor para cada elemento de ordenador com o valor da sua nova ordem na sequência
        
    for i in range(len(arr)):
        
        a = ordenador[i][0]
        b = ordenador[i][1]             # É feita a troca do segundo com o primeiro valor de cada elemento de ordenador para...
        ordenador[i][0] = b
        ordenador[i][1] = a
        
    ordenador.sort()                    # reordenar os elementos de acordo com a ordem antiga
        
        
    lista = []
    
    for i in range(len(arr)):
        
        lista.append(ordenador[i][2])   #por fim guarda-se a ordem correspondente a cada valor da lista colocada como variável 
        
    indexador = np.array(lista)  
        
    return indexador                     # e essa ordem é retornada como saída da função

print(indexador(teste))

Função de movimento e aceleração das partículas

Para simular como um estado evolui no tempo, é necessário que tenhamos as derivadas temporais de todas as coordenadas do vetor de estado ($y$).

Para as coordenadas espaciais x ($y[i],\ 0\leq i<N$) temos que:

${\dot {x_{i}}}=y_{[i+N]}$

Para as coordenadas de velocidade v ($y[i],\ N\leq i<2N-1$) teremos que usar a relação de forças que foi descrita no modelo. Algumas simplificações comuns serão feitas como considerar constantes físicas unitárias ($m=q=\kappa =1$) e relações de proporcionalidade com o número de partículas do sistema para mantê-lo com um comportamento mais suave ($F=1/N$). O programa será escrito levando em conta o modelo de cargas elétricas:

${\dot {v_{i}}}={\ddot {x_{i}}}={\frac {1}{m}}{\Big [}-\kappa y_{[i-N]}+(2i-N-1)F{\Big ]}$

A relação acima pode ser descrita de forma diferente se levarmos em conta que a simulação é simétrica em x e que todas as posições iniciais são positivas.

A partícula que recebe índice $1$ possui $N$ partículas à esquerda (x negativo) e $(N-1)$ partículas a direita (x positivo), o que resulta uma força de interação de $F$ para a direita. A de índice 2 terá força de interação $3F$, a de índice 3 terá $5F$ e assim por diante.

Antes de reescrever a equação da aceleração deve-se levar em consideração de uma das partículas distribuídas no espaço de fase assumirem coordenadas negativas (caso tenham proximidade da origem e velocidade negativa). Novamente usando o argumento de simetria, sabe-se que caso uma partícula assuma uma posição negativa, a partícula espelhada assumirá a posição simétrica em relação a origem. Dessa forma basta calcular a força de interação sobre a partícula espelhada e inverter o sinal para a partícula original. Levando tudo isso em conta teremos a seguinte equação da aceleração:

${\ddot {x_{i}}}={\frac {1}{m}}{\Big [}-\kappa y_{[i-N]}+{\frac {y_{[i-N]}}{|y_{[i-N]}|}}(2i-1)F{\Big ]}$

def deriv(t, y):

    p = len(y)
    k = indexador(y[0:int(p/2)])
    sinal=abs(y[0:int(p/2)])/y[0:int(p/2)]
    
    dydt = np.zeros(int(p))

    dydt[0:int(p/2)] = y[int(p/2):int(p)] #definindo as derivadas temporais de posição
    
                                                         #definindo as derivadas temporais de velocidade
    dydt[int(p/2):int(p)] = -y[0:int(p/2)] + ((np.ones(int(p/2)) + 2*k )/(p))*sinal
    
    deriv = dydt[0:int(p)]
    
    return deriv

Função de Energia do estado

def Energia_total(y):     # Para calcular a energia precisamos apenas do vetor de estado y do sistema de partículas utilizado

    l = int(len(y)/2)     # a variável l representa a quantidade de valores que teremos 
    pot = np.zeros(l+1)   # cria-se um array com l + 1 dimensões para representar as l partículas e os respectivos potenciais 

    r = abs(np.array(y[0:l]))       # cria-se um array r com as coordenadas de posição do estado y
    r = np.append(r,4)   # adiciona-se uma coordenada para representar o contorno externo do experimento
    r.sort()              # Faz-se o ordenamento desse array por distância da origem
    
    pot_0 = 0             # O potencial do contorno deve ser nulo 
    
    for i in range(l,0,-1):     # Loop decrescente que define o potencial de posição para cada partícula 
        
        pot[i] = pot_0
        
        pot_0 = pot_0 + (1/2)*(i/l)*abs((r[i] - r[i-1]))
        
    pot[0] = pot_0


    v = np.array(y[l:-1])       # Cria-se um array com os valores de velocidade
    
    E_total = 0


    for i in range(l-1):        #Soma-se as energias de cada partícula

      E_total = E_total + pot[i] + 0.5*v[i]**2 + 0.5*r[i]**2


    return 2*E_total, r[0:l+1], pot[0:l+1]

Evolução da Distribuição no Tempo

Para simular a evolução temporal da distribuição inicial escolhida, precisamos de um método de integração numérica no tempo. O método escolhido deverá partir do estado da distribuição e da função de derivada temporal atribuída para estimar como a evolução ocorre.

from scipy.integrate import solve_ivp  #integrador numérico escolhido

def nova_fase(y,i):      #função que define a função do espaço de fase em cada instante de tempo, adaptada para o integrador numérico escolhido
    
    fase = np.array([])
    
    for j in range(len(y)):
            
        fase = np.append(fase,y[j][i])
        
    return fase


Fase = []   # Cria uma lista para reservar a evolução do espaço de fase no tempo
T = np.array([])     # Cria uma lista que registra os tempos de marcação do integrador numérico
E = np.array([])     # Cria uma lista que registra a energia da distribuição no tempo
E_pot = []           # Cria uma lista que registra o comportamento da energia potencial ao longo de x no tempo
X_pot = []           # O análogo ao anterior para a posição
tempo = 0.0  # o tempo começa em 0
delta_t = 1   # o intervalo de tempo é 0.1

X = np.random.rand(4000)
V = np.random.rand(4000) - np.ones(4000)*0.5
y = np.append(X,2*V)    #define-se a fase inicial

for j in range(50):

    sol = solve_ivp(deriv, (tempo,tempo + delta_t), y , rtol = 0.00001)     # integrador numérico realiza a integração

    T= np.append(T, sol.t)

    tempo = tempo + delta_t
    for i in range(len(sol.t)):       #o integrador possui um número de vetores intermediários ao vetor da solução final sol.y[-1]
        
        
        y = nova_fase(sol.y,i)                              # aqui são registradas as funções definidas para cada estado de integração
        Fase.append(y)

        Ei =  Energia_total(y)

        E = np.append(E, Ei[0])

        E_pot.append(Ei[2])
        X_pot.append(Ei[1])
    
    y = nova_fase(sol.y, -1)                     # A distribuição atual é atualizada para nova integração

Animações

Foi utilizado o código abaixo para gerar os gifs da evolução temporal

from matplotlib.animation import FuncAnimation
from IPython import display

# código base para montar gifs da evolução temporal 
fig, (ax1, ax2, ax3) = plt.subplots(3,1 ,figsize = (6,14))    #cria uma figura de três linhas

lines = []                                                    #objeto padrão da biblioteca FuncAnimation


def animate(frame):


    # As entradas são limpas para evitar superposição de imagens nos frames
    ax1.cla()              
    ax2.cla()
    ax3.cla()

    # Plot do espaço de fases
    ax1.scatter(Fase[frame][0:4000], Fase[frame][4000:8000], color = 'k', s  =5)
    ax1.scatter(-Fase[frame][0:4000], -Fase[frame][4000:8000], color = 'k', s = 5)
    ax1.set_xlabel('$posição$')
    ax1.set_ylabel('$velocidade$')
    ax1.set_xlim(-4,4)
    ax1.set_ylim(-2.5,2.5)   

    # Plot da distribuição de Potencial
    ax2.plot(X_pot[frame],E_pot[frame], color = 'grey')
    ax2.plot(-np.flip(X_pot[frame]),np.flip(E_pot[frame]), color = 'grey')
    ax2.scatter(X_pot[frame],E_pot[frame], color = 'b', s  =5)
    ax2.scatter(-np.flip(X_pot[frame]),np.flip(E_pot[frame]), color = 'b', s  =5)
    ax2.set_ylabel("potencial")
    ax2.set_xlabel("posição")

    # Plot da energia total no tempo
    ax3.plot(T[0:frame],E[0:frame], color = 'k')
    ax3.set_ylabel("Energia")
    ax3.set_xlabel("tempo")
    ax3.set_xlim(0, T[-1])
    ax3.set_ylim(E[0]- 10, E[0] + 10)
    
    return lines

# Execução do programa e geração do gif inline   
anima = FuncAnimation(fig, animate,frames=len(T), interval=20, blit=True)
video = anima.to_html5_video()
html = display.HTML(video)
display.display(html)

plt.close()

Referências

CARLAN JR, Antônio H.; PAKTER, Renato. Core-halo boundary in a sheet beam model. Physics of Plasmas, v. 28, n. 11, p. 113103, 2021.

R. L. Gluckstern, Phys. Rev. Lett. 73, 1247 (1994).

C. Chen and R. C. Davidson, Phys. Rev. Lett. 72, 2195 (1994).