Modelo Blume-Capel para partículas de spin 1: mudanças entre as edições

De Física Computacional
Ir para navegação Ir para pesquisar
Misalocin (discussão | contribs)
m Ajuste de sinal
Misalocin (discussão | contribs)
Formatação
 
(14 revisões intermediárias pelo mesmo usuário não estão sendo mostradas)
Linha 1: Linha 1:
== <span style="color: red;"> Essa página está em construção </span> ==  
== Modelo de Ising ==


Considere um grupo de partículas com spin <math>1/2</math> que interagem com seus vizinhos mais próximos. Pode-se calcular a energia do sistema através do Hamiltoniano <ref name="mecEst" />


== Teoria ==
<center>
 
=== Modelo de Ising ===
 
Considere um grupo de partículas com spin <math>1/2</math> que interagem com seus vizinhos mais próximos. Pode-se calcular a energia do sistema através do Hamiltoniano [CITATION]
 
<math>
<math>
H=-J\sum_{<i,j>} \sigma_i \sigma_j
\mathcal{H}=-J\sum_{<i,j>} \sigma_i \sigma_j
</math>
</math>
</center>


Com <math>\sigma_i</math> representando o valor do spin da partícula <math>i</math>, podendo ser <math>\pm 1/2</math>. O somatório representa a energia potencial das partículas vizinhas, onde J representa o comportamento e a propensão ao alinhamento do sistema, podendo ser qualquer número real, no qual valores positivos de <math>J</math> diminuindo a energia quando os spins estão alinhados e, valores negativos a aumentando. Assim, quanto mais longe de 0 $J</math> estiver, mais o sistema tende a um extremo.
Com <math>\sigma_i</math> representando o valor do spin da partícula <math>i</math>, podendo ser <math>\pm 1/2</math>. O somatório representa a energia potencial das partículas vizinhas, onde J representa o comportamento e a propensão ao alinhamento do sistema, podendo ser qualquer número real, no qual valores positivos de <math>J</math> diminuindo a energia quando os spins estão alinhados e, valores negativos a aumentando. Assim, quanto mais longe de 0 $J</math> estiver, mais o sistema tende a um extremo.
Linha 16: Linha 13:
Consideremos agora a aplicação de um campo magnético <math>H</math> externo, paralelo a direção do spin. Esse campo irá adicionar uma energia ao sistema que será
Consideremos agora a aplicação de um campo magnético <math>H</math> externo, paralelo a direção do spin. Esse campo irá adicionar uma energia ao sistema que será


<center>
<math>
<math>
H_{Mag}=-H\sum_{n=1}^N\sigma_n
\mathcal{H}_{Mag}=-H\sum_{n=1}^N\sigma_n
</math>
</math>
</center>


O termo negativo da equação surge pois a energia total do sistema será menor quando os spins estiverem alinhados com campo, com o mínimo de energia sendo <math>-HN</math>, onde todas as <math>N</math> partículas estão com o spin alinhado.
O termo negativo da equação surge pois a energia total do sistema será menor quando os spins estiverem alinhados com campo, com o mínimo de energia sendo <math>-HN</math>, onde todas as <math>N</math> partículas estão com o spin alinhado.


Assim, a energia total do sistema será[CITATION]:
Assim, a energia total do sistema será<ref name="mecEst" />:


<center>
<math>
<math>
H=-J\sum_{<i,j>} \sigma_i \sigma_j -H\sum_{n=1}^N\sigma_n
\mathcal{H}=-J\sum_{<i,j>} \sigma_i \sigma_j -H\sum_{n=1}^N\sigma_n
</math>
</math>
</center>


Agora, consideremos o ensamble canônico, isso é, o conjunto de todos os estados <math>\lambda</math> que o sistema pode estar. Se a energia do sistema for medida com todas as partículas dentro dele estiverem dispostas aleatoriamente, a probabilidade de encontrar o sistema na energia <math>E</math> é: [CITATION]
Agora, consideremos o ensemble canonico, isso é, o conjunto de todos os estados <math>\lambda</math> que o sistema pode estar. Se a energia do sistema for medida com todas as partículas dentro dele estiverem dispostas aleatoriamente, a probabilidade de encontrar o sistema na energia <math>E</math> é: <ref name="mecEst" />


<center>
<math>
<math>
P(\lambda)\propto e^{-\beta E(\lambda)}
P(\lambda)\propto e^{-\beta E(\lambda)}
</math>
</math>
</center>


com <math>\beta=\frac{1}{k_b T}</math>. Considerando que a soma de todas as probabilidades é 1, deves-se dividir <math>e^{-\beta E}</math> pela soma de <math>e^{-\beta E_n(\lambda)}</math> para todos os estados <math>\lambda</math> ao qual o sistema pode ser encontrado.
com <math>\beta=\frac{1}{k_b T}</math>. Considerando que a soma de todas as probabilidades é 1, deves-se dividir <math>e^{-\beta E}</math> pela soma de <math>e^{-\beta E_n(\lambda)}</math> para todos os estados <math>\lambda</math> ao qual o sistema pode ser encontrado.
Importante ressaltar que soma é sobre <math>\lambda</math> e não sobre todas as energias, pois muitos estados podem gerar a mesma energia, então não é feita uma soma simples por todo o sistema.


Assim, o resultado é:
Assim, o resultado é:


<center>
<math>
<math>
P(\lambda)=\frac{e^{-\beta E(\lambda)}}{\sum_{\lambda_n} e^{-\beta E(\lambda_n)}}
P(\lambda)=\frac{e^{-\beta E(\lambda)}}{\sum_\lambda' e^{-\beta E(\lambda')}}
</math>
</math>
</center>


Outra coisa importante ainda do modelo de Ising é a magnetização, que representa o quâo ordenado o sistema é à partir do valor médio de todos os spins:
Outra coisa importante ainda do modelo de Ising é a magnetização, que representa o quão ordenado o sistema é à partir do valor médio de todos os spins:


<center>
<math>
<math>
M=\frac{1}{N}\sum_{n=1}^{N}\sigma_n
M=\frac{1}{N}\sum_{n=1}^{N}\sigma_n
</math>
</math>
</center>


== Modelo de Blume-Capel ==
== Modelo de Blume-Capel ==


=== Ennergia do sistema ===
=== Energia do sistema ===


Esse modelo é uma extensão do modelo de Ising onde, ao invés de se trabalhar com partículas de spin <math>1/2</math> onde <math>\sigma_i\in \{-\frac{1}{2},\frac{1}{2}\}</math>, trabalha-se com partículas de spin 1 com <math>\sigma_i\in \{-1,0,1\}</math>.
Esse modelo é uma extensão do modelo de Ising onde, ao invés de se trabalhar com partículas de spin <math>1/2</math> onde <math>\sigma_i\in \{-\frac{1}{2},\frac{1}{2}\}</math>, trabalha-se com partículas de spin 1 com <math>\sigma_i\in \{-1,0,1\}</math>.


Essa mudança de um para três estados não muda de forma significativa o que foi visto anteriormente, mas adiciona um terceiro termo no hamiltoniano <math>H</math>. Como aqui partículas podem ter spin nulo, define-se uma nova medida que está relacionada a quantidade de spins diferentes de zero, o mmomento de quadripolo magnético[CITATION], definido como:
Essa mudança de um para três estados não muda de forma significativa o que foi visto anteriormente, mas adiciona um terceiro termo no hamiltoniano <math>\mathcal{H}</math>. Como aqui partículas podem ter spin nulo, define-se uma nova medida que está relacionada a quantidade de spins diferentes de zero, o momento de quadrupolo magnético<ref name="point" />, definido como:


<center>
<math>
<math>
Q= \sum_{n=1}^N\sigma_n^2
Q= \sum_{n=1}^N\sigma_n^2
</math>
</math>
</center>


Esse momento também é adicionado a energia, pois sair do estado <math>\sigma = 0</math> adiciona energia aos sistema, isso faz com que o hamiltoniano completo seja:
Esse momento também é adicionado a energia, pois sair do estado <math>\sigma = 0</math> adiciona energia aos sistema, isso faz com que o hamiltoniano completo seja:


<center>
<math>
<math>
H= - J\sum_{<i,j>} \sigma_i \sigma_j - H\sum_{n=1}^N\sigma_n + D\sum_{m=1}^N\sigma_m^2
\mathcal{H}= - J\sum_{<i,j>} \sigma_i \sigma_j - H\sum_{n=1}^N\sigma_n + D\sum_{m=1}^N\sigma_m^2
</math>
</math>
</center>


No modelo de Blume-Capel, o quadrado do momento de quadripolo pode ser substituido por um símples módulo, já que os spins serão -1,0 ou 1.
Onde o parâmetro <math>D</math> é chamado de campo cristalino. Nessa implementação do modelo de Blume-Capel, o quadrado do momento de quadrupolo pode ser substituido por um módulo, já que os spins serão -1,0 ou 1.




=== Ponto tri-crítico ===
=== Ponto tri-crítico ===


O ponto mais interessante do modelo Blume-Capel é o ponto tri-crítico.
Analisa-se o sistema para um dado valor de J e considerando o sistema sem a presença de um campo magnético, de forma que a energia total seja:
 
Analisa-se o sitema para um dado valor de J e considerando o sistema sem a presença de um campo magnético, de forma que a energia total seja:


<center>
<math>
<math>
H= - J\sum_{<i,j>} \sigma_i \sigma_j + D\sum_{m=1}^N\sigma_m^2
\mathcal{H}= - J\sum_{<i,j>} \sigma_i \sigma_j + D\sum_{m=1}^N\sigma_m^2
</math>
</math>
</center>


Analisando esse sistema simplificado, considerando tempos muito grandes e dividindo tudo por J, o sistema tenderá a três possíveis estados de energia mínima:
Analisando esse sistema simplificado, considerando tempos muito grandes e dividindo tudo por J, o sistema tenderá a três possíveis estados de energia mínima:


* <math>D/J >> 0</math>: Aqui, o spin de cada partícula tem um peso grande sobre a energia total, o sistema tenderá a um estado onde grande parte dos spins são 0, logo, <math>M</math> e <math>Q</math> são próximos de 0.
* <math>D/J \gg 0</math>: Aqui, o spin de cada partícula tem um peso grande sobre a energia total, o sistema tenderá a um estado onde grande parte dos spins são 0, logo, <math>M</math> e <math>Q</math> são próximos de 0.
 
* <math>T/J \gg 0</math> e <math>D/J</math> coerente: Aqui, analisando-se <math>P(\lambda)</math>, como <math>\beta\propto T^{-1}</math>, beta é aproximadamente 0 e a probabilidade de se encontrar o sistema em dada configuração é a mesma para todas as configurações. Nesse caso, é mais provável que o sistema esteja com <math>M\approx 0</math> e <math>Q > 0 </math> com as partículas com spins aleatórios.
 
* <math>T/J </math> e <math>D/J</math> coerentes: Aqui, analisando novamente <math>P(\lambda)</math>, estados onde a energia do sistema seja menor que zero serão favorecidos, pois estados positivos terão uma probabilidade menor que 1 (<math>e^{-\mathcal{H}/T}<1</math>) e estados de energia menor que zero terão uma probabilidade maior (<math>e^{-(|\mathcal{H}|)/T}>1</math>), isso fará com que os spins tendam a se alinhar e <math>M</math> e <math>Q</math> serão próximos de 1
 


* <math>T/J >> 0</math> e <math>D/J</math> coerente: Aqui, analisando-se <math>P(\lambda)</math>, como <math>\beta\propto T^{-1}</math>, beta é approximadamente 0 e a probabilidade de se encontrar o sistema em dada configuração é a mesma para todas as configurações. Nesse caso, é mais provável que o sistema esteja com <math>M\approx 0</math> e <math>Q > 0 </math> com as partículas com spins aleatórios.
Considerando que existem 3 regiões e se é possível mudar de uma para a outra, existe um ponto de interseção dos três estados. Esse ponto é chamado de ponto tri-crítico(TCP). Os valores de T e D desse ponto ainda não possuem um valor exato, mas aproximações numéricas estipulam que eles devem ser <math>T\approx 0.610\pm0.005</math> e <math>D\approx 1.965\pm0.001</math> <ref name="point" />


* <math>T/J </math> e <math>D/J</math> coerentes: Aqui, analisando novamente <math>P(\lambda)</math>, estados onde a energia do sistema seja menor que zero serão favorecidos, pois estados positivos terão uma probabilidade menor que 1 (<math>e^{-\tfrac{H}{T}}<1</math>) e estados de energia menor que zero terão uma probabilidade maior (<math>e^{\frac{|H|}{T}}>1</math>), isso fará com que os spins tendam a se alinhar e <math>M</math> e <math>Q</math> serão próximos de 1
=== Unidades Naturais ===
Antes de começar, é conveniente introduzir o conceito de unidades reduzidas. Como apenas razões entre os parâmetros físicos aparecerão nas equações finais, pode-se escolher uma constante como unidade de energia e utilizar unidades reduzidas, usando <math>k_B = 1</math> e tomando <math>J</math> como unidade de energia. Dessa forma, todas as grandezas passam a ser medidas em unidades de <math>J</math>, sendo comum trabalhar com as variáveis adimensionais


<center>
<math>T^* = \frac{k_B T}{J};</math>
<math>D^* = \frac{D}{J};</math>
<math>H^* = \frac{H}{J}.</math>
</center>


Considerando que exitem 3 regiões e se é possível mudar de uma para a outra, existe um ponto de interseção dos três estados. Esse ponto é chamado de ponto tri-crítico. Os valores de T e D desse ponto ainda não possuem um valor exato, mas aproximações numéricas estipuilam que eles devem ser <math>T\approx 0.610\pm0.005</math> e <math>D\approx 1.965\pm0.001</math> [CITATION]
Ao longo deste trabalho, os asteriscos serão omitidos, e todas as grandezas são entendidas como estando em unidades reduzidas.


==  Método Monte-Carlo ==
==  Método Monte-Carlo ==


Tendo definido todas as informações do modelo, pode-se agora
=== Algoritmo ===
 
Como o número de estados cresce exponencialmente com o número de partículas, não é possível calcular essa distribuição diretamente. Em vez disso, constrói-se uma sequência de estados através de atualizações locais dos spins.
 
O procedimento consiste em repetir os seguintes passos:
 
# Escolhe-se aleatoriamente uma partícula.
# Propõe-se um novo valor de spin para essa partícula, escolhido uniformemente entre os três valores possíveis <math>\sigma \in \{-1, 0, 1\}.</math>
# Calcula-se a variação de energia associada à mudança proposta: <math>\Delta E = E_{\text{novo}} - E_{\text{atual}}.</math>
 
Sabendo a variação de energia, pode-se então decidir se o sistema será ou não atualizado.
 
=== Probabilidade de transição ===
 
Partindo do balanço detalhado
 
<center>
<math>
P(i) W(i \rightarrow j) = P(j) W(j \rightarrow i),
</math>
</center>
 
E a probabilidade do sistema estar no estado <math>P(\lambda)</math> definida anteriormente, chega-se que:
 
<center><math>
\frac{P(\lambda_1 \to \lambda_2)}{P(\lambda_2 \to \lambda_1)} = e^{-\Delta E/T}.
</math></center>
 
uma das maneiras de se definir <math>P</math> que resolvem esse sistema é o algoritmo de Metropolis <ref name="mecEst" />, aqui, a probabilidade de transição é:
 
<center><math>
P(\lambda_1 \to \lambda_2)=min\left(1,e^{\frac{E(\lambda_2)-E(\lambda_1)}{T}}\right)
</math></center>
 
Equivalentemente, gera-se um número aleatório uniforme <math>r \in [0,1)</math>. A atualização é aceita se
 
<center><math>
r < e^{-\Delta E/T}.
</math></center>
 
essa solução garante,também, que a solução da cadeia de Markov seja exatamente a distribuição de Boltzmann.
 
Quando o spin aleatório sorteado for igual ao atual, o sistema deve também guardar essa mudança, isso se dá para que o sistema respeite o balanço detalhado e não assuma que a probabilidade de permanência em um estado seja nula.
 
==== Passo de Monte Carlo ====
 
Considere um sistema contendo <math>N</math> partículas em uma configuração <math>\lambda</math>. Em algoritmos de Monte Carlo, um passo elementar consiste na seleção de uma partícula e na tentativa de atualização de seu estado segundo a dinâmica escolhida (por exemplo, o algoritmo de Metropolis).
 
Como o número de partículas depende do tamanho do sistema, utilizar apenas o número de passos elementares dificulta a comparação entre simulações realizadas em redes de diferentes dimensões. Por esse motivo, define-se uma unidade de tempo computacional normalizada pelo tamanho do sistema.
 
O Passo de Monte Carlo (MCS) corresponde a NNN tentativas de atualização elementares. Dessa forma, após um MCS, cada partícula terá sido selecionada, em média, uma vez para possível atualização. Essa definição permite comparar de maneira consistente a evolução temporal e as grandezas medidas em sistemas de diferentes tamanhos.
 
== Implementação Computacional ==
 
Começa-se criando uma matriz <code>state</code>. Como esse código será executado diversas vezes com loops longos, diversos pontos devem ser otimizado para velocidade ao invés de legibilidade. O primeiro deles é o uso de uma tabela de referência <code>lookUp</code>. Ela é gerada uma vez sempre que o tamanho da grade é alterado e possui 2 índices, o primeiro (<code>i</code>) sendo a partícula que se quer saber a posição das células adjacentes, e o segundo é um valor de 0 a 3 que guarda, em ordem, o índice da partícula imediatamente acima, abaixo, a esquerda e a direita da partícula <code>i</code>.
 
Isso também permite que seja montado o sistema como periódico, onde a primeira e a última linha seriam vizinhas, assim como a primeira e última coluna. Isso é útil para diminuir os erros que surgem nas bordas, onde as partículas não teriam vizinhas nas quatro direções.
 
<syntaxhighlight lang="python">
def genLookup(side):
    index = np.arange(side)
   
    up    = (index - side) % (side**2)
    down  = (index + side) % (side**2)
    left  = np.where(index      % side == 0, index - 1 + side, index - 1)
    right = np.where((index + 1) % side == 0, index + 1 - side, index + 1)
 
    return np.stack((up, down, left, right))
</syntaxhighlight>
 
Após isso, é pode-se iniciar o sistema com uma matriz 1D de N partículas, onde N é o quadrado do tamanho desejado de um lado da grade. Usa-se matrizes 1D ao invés de 2D para maximizar a velocidade do código.
 
Com o sistema iniciado, se implementa o algoritmo de Metropolis em uma função que, quando aplicada sobre o estado <math>\lambda_i</math>, evolui-rá ele para outro estado com uma quantidade definida de passos, podendo ser um passo fundamental, um MCS, ou qualquer outro valor conveniente.
 
<syntaxhighlight lang="python">
def timeSteps(state, lookup, steps = 1, J = 1, D = 1, T = 1):
    N = state.size
    varE = 0.0
 
    # Repete a função para a quantidade de passos dada
    for i in range(steps):
   
        # A cada passo, escolhe uma partícula e um spin aleatórios
        index = np.random.randint(0, N)
        oldSpin = state[index]
        newSpin = np.random.randint(-1, 2)
   
        # avança a iteração se os spins forem iguais
        if newSpin == oldSpin:
            continue
       
        # Soma o estado dos vizinhos através da tabela de referência
        neiSum = np.sum(state[lookup[:, index]])
 
        # Calcula a mudança de energia do passo
        dE = -J * (newSpin - oldSpin) * neiSum + D * (newSpin**2 - oldSpin**2)
       
        # Verifica se o passo foi ou não aceito
        if np.random.random() < np.exp(-dE / T):
            state[index] = newSpin
            varE += dE       
           
    # Retorna o estado depois de todos os passos e a variação de energia total   
    return state, varE
</syntaxhighlight>
 
Com todas essas etapas montadas, pode-se gerar a primeira simulação. Abaixo, pode-se analisar uma simulação com valores de D e T fixos e um número arbitrário de MCS
 
[[file:Evolução blume capel.gif|thumb|400px|center|Simulação da evolução de um sistema com parâmetros aleatórios]]
 
Aqui, pode-se observar a evolução de um único sistema, mas o grande trunfo do modelo de Blume-Capel é analisar os estados do sistema, especificamente o TCP. Para isso, uma única simulação não é suficiente, mas será necessário ver o estado de equilíbrio do sistema em um gráfico de Temperatura (T)/Campo cristalino (D) e analisar a magnetização média e o momento de quadrupolo médio.
 
Como foi dito anteriormente, essas duas medidas serão iguais em dois casos e diferentes em um. Então o gráfico Temperatura (T)/Campo cristalino (D) irá mostrar <math><M></math> e <math><Q></math>, além de um terceiro gráfico mostrando o módulo da diferença entre eles.
 
[[file:Blume Capel espaco de fase.gif|thumb|800px|center|Evolução da magnetização, momento de quadrupolo e a diferença de ambos em uma grande área, simulações utilizada para identificar a posição dos estados]]
 
Desenhando esses três gráficos, pode-se perceber que, onde o TCP teórico está próximo dos três estados de equilíbrio teóricos. Mas para avaliar se ele realmente é o TCP, é necessário verificar duas propriedades.
 
# Nele, o sistema não deve chegar a um estado de equilíbrio único, mas sim alternar entre os conhecidos;
# A vizinhança dele deve ser composta por simulações que levam a um estado bem definido.
 
Com isso, para avaliar o sistema, foram geradas várias simulações com diferentes valores de <math>T</math> e <math>D</math> dem volta do valor teórico para encontrar qual é o valor que se chega nessa simulação e, se for possível encontrá-lo, qual é o tamanho do ponto.
 
[[File:Blume Capel espaco de fase zoom.gif|thumb|800px|center|Espaços de fases Temperatura em y e D em x aproximado onde na região do PTC. A esquerda a evolução da magnetização do sistema, no meio a evolução do momento de quadrupolo e, a direita, o absoluto da diferença entre os dois. Um ponto vermelho localiza-se na previsão teórica do TCP, com três pontos azuis distribuídos ao seu redor mostrando pontos de simulação.]]
 
Definindo os pontos a serem simulados separadamente, pode-se animar a evolução do sistema com os devidos parâmetros.
 
{| style="width:100%"
|-
| [[File:Blume Capel Estado1.gif|thumb|300px|center|Estados onde Q é maior que E, com as partículas oscilando entre os três estados aleatoriamente.]]
| [[File:Blume Capel Estado2.gif|thumb|300px|center|Estados com Q e M próximos de 1, com a simulação estando de maneira quase homogênea em um único spin diferente de 0.]]
| [[File:Blume Capel Estado3.gif|thumb|300px|center|Estados com Q e M estão próximos de 0, com a maior parte das partículas estando com spin igual a 0.]]
|}
Além disso, pode-se também ver a evolução do estado no TCP:
 
{| style="width:100%"
|-
| [[File:Blume Capel PTC1.png|thumb|250px|center|Momento em que é observado o estado com Q e M baixos.]]
| [[File:Blume Capel PTC2.png|thumb|250px|center|Momento em que é observado o estado com Q e M altos.]]
| [[File:Blume Capel PTC3.png|thumb|250px|center|Momento em que é observado o estado com Q muito maior que M.]]
| [[File:Blume Capel PTC.gif|thumb|200px|center|Animação total onde os gráficos foram extraídos, clique para ver animado (aviso de luzes piscantes).]]
|}
 
== Conclusão ==
 
Esse trabalho deixou claro a existência de três estados possíveis, além de mostrar que em uma interface de estados o sistema muda voluntariamente entre os três diferentes estados. pode-se ver, parando o gráfico em alguns pontos específicos, que ele passou pelos três estados.
 
Assim, o código aqui descrito mostra claramente uma implementação do modelo de Blume-Capel capaz de simular os três estados e um ponto onde os três oscilam entre si aleatoriamente. É importante notar que a variação entre os estados não é tão clara no TCP, pois os estados com M próximo de 0 acabam sendo semelhantes devido ao número de alterações que ocorrem a todo MCS.


https://journals.aps.org/prb/abstract/10.1103/PhysRevB.33.1717
Mas o algoritmo aqui descrito está disponível na página [[Código:Modelo de Blume-Capel]] feito em Jupyter Notebook, é capaz de prever de maneira satisfatória a evolução de um sistema à partir das variáveis supridas à ele.  


https://fiscomp.if.ufrgs.br/index.php?title=Modelo_de_Blume-Capel_bidimensional
== Referências ==
<references>
<ref name="mecEst">Pathria, R. K.; Beale, Paul D.; ''Statistical Mechanics'', third edition; Elsevier Ltd, 2011; secções 3.4, 12.3, 16.2A</ref>
<ref name="point">Beale, Paul D.; ''Finite-size scaling study of the two-dimensional Blume-Capel model''; APS Journals, 1986; DOI https://doi.org/10.1103/PhysRevB.33.1717</ref>
</references>

Edição atual tal como às 15h23min de 1 de junho de 2026

Modelo de Ising

Considere um grupo de partículas com spin 1/2 que interagem com seus vizinhos mais próximos. Pode-se calcular a energia do sistema através do Hamiltoniano [1]

=J<i,j>σiσj

Com σi representando o valor do spin da partícula i, podendo ser ±1/2. O somatório representa a energia potencial das partículas vizinhas, onde J representa o comportamento e a propensão ao alinhamento do sistema, podendo ser qualquer número real, no qual valores positivos de J diminuindo a energia quando os spins estão alinhados e, valores negativos a aumentando. Assim, quanto mais longe de 0 $J</math> estiver, mais o sistema tende a um extremo.

Consideremos agora a aplicação de um campo magnético H externo, paralelo a direção do spin. Esse campo irá adicionar uma energia ao sistema que será

Mag=Hn=1Nσn

O termo negativo da equação surge pois a energia total do sistema será menor quando os spins estiverem alinhados com campo, com o mínimo de energia sendo HN, onde todas as N partículas estão com o spin alinhado.

Assim, a energia total do sistema será[1]:

=J<i,j>σiσjHn=1Nσn

Agora, consideremos o ensemble canonico, isso é, o conjunto de todos os estados λ que o sistema pode estar. Se a energia do sistema for medida com todas as partículas dentro dele estiverem dispostas aleatoriamente, a probabilidade de encontrar o sistema na energia E é: [1]

P(λ)eβE(λ)

com β=1kbT. Considerando que a soma de todas as probabilidades é 1, deves-se dividir eβE pela soma de eβEn(λ) para todos os estados λ ao qual o sistema pode ser encontrado. Importante ressaltar que soma é sobre λ e não sobre todas as energias, pois muitos estados podem gerar a mesma energia, então não é feita uma soma simples por todo o sistema.

Assim, o resultado é:

P(λ)=eβE(λ)λeβE(λ)

Outra coisa importante ainda do modelo de Ising é a magnetização, que representa o quão ordenado o sistema é à partir do valor médio de todos os spins:

M=1Nn=1Nσn

Modelo de Blume-Capel

Energia do sistema

Esse modelo é uma extensão do modelo de Ising onde, ao invés de se trabalhar com partículas de spin 1/2 onde σi{12,12}, trabalha-se com partículas de spin 1 com σi{1,0,1}.

Essa mudança de um para três estados não muda de forma significativa o que foi visto anteriormente, mas adiciona um terceiro termo no hamiltoniano . Como aqui partículas podem ter spin nulo, define-se uma nova medida que está relacionada a quantidade de spins diferentes de zero, o momento de quadrupolo magnético[2], definido como:

Q=n=1Nσn2

Esse momento também é adicionado a energia, pois sair do estado σ=0 adiciona energia aos sistema, isso faz com que o hamiltoniano completo seja:

=J<i,j>σiσjHn=1Nσn+Dm=1Nσm2

Onde o parâmetro D é chamado de campo cristalino. Nessa implementação do modelo de Blume-Capel, o quadrado do momento de quadrupolo pode ser substituido por um módulo, já que os spins serão -1,0 ou 1.


Ponto tri-crítico

Analisa-se o sistema para um dado valor de J e considerando o sistema sem a presença de um campo magnético, de forma que a energia total seja:

=J<i,j>σiσj+Dm=1Nσm2

Analisando esse sistema simplificado, considerando tempos muito grandes e dividindo tudo por J, o sistema tenderá a três possíveis estados de energia mínima:

  • D/J0: Aqui, o spin de cada partícula tem um peso grande sobre a energia total, o sistema tenderá a um estado onde grande parte dos spins são 0, logo, M e Q são próximos de 0.
  • T/J0 e D/J coerente: Aqui, analisando-se P(λ), como βT1, beta é aproximadamente 0 e a probabilidade de se encontrar o sistema em dada configuração é a mesma para todas as configurações. Nesse caso, é mais provável que o sistema esteja com M0 e Q>0 com as partículas com spins aleatórios.
  • T/J e D/J coerentes: Aqui, analisando novamente P(λ), estados onde a energia do sistema seja menor que zero serão favorecidos, pois estados positivos terão uma probabilidade menor que 1 (e/T<1) e estados de energia menor que zero terão uma probabilidade maior (e(||)/T>1), isso fará com que os spins tendam a se alinhar e M e Q serão próximos de 1


Considerando que existem 3 regiões e se é possível mudar de uma para a outra, existe um ponto de interseção dos três estados. Esse ponto é chamado de ponto tri-crítico(TCP). Os valores de T e D desse ponto ainda não possuem um valor exato, mas aproximações numéricas estipulam que eles devem ser T0.610±0.005 e D1.965±0.001 [2]

Unidades Naturais

Antes de começar, é conveniente introduzir o conceito de unidades reduzidas. Como apenas razões entre os parâmetros físicos aparecerão nas equações finais, pode-se escolher uma constante como unidade de energia e utilizar unidades reduzidas, usando kB=1 e tomando J como unidade de energia. Dessa forma, todas as grandezas passam a ser medidas em unidades de J, sendo comum trabalhar com as variáveis adimensionais

T*=kBTJ; D*=DJ; H*=HJ.

Ao longo deste trabalho, os asteriscos serão omitidos, e todas as grandezas são entendidas como estando em unidades reduzidas.

Método Monte-Carlo

Algoritmo

Como o número de estados cresce exponencialmente com o número de partículas, não é possível calcular essa distribuição diretamente. Em vez disso, constrói-se uma sequência de estados através de atualizações locais dos spins.

O procedimento consiste em repetir os seguintes passos:

  1. Escolhe-se aleatoriamente uma partícula.
  2. Propõe-se um novo valor de spin para essa partícula, escolhido uniformemente entre os três valores possíveis σ{1,0,1}.
  3. Calcula-se a variação de energia associada à mudança proposta: ΔE=EnovoEatual.

Sabendo a variação de energia, pode-se então decidir se o sistema será ou não atualizado.

Probabilidade de transição

Partindo do balanço detalhado

P(i)W(ij)=P(j)W(ji),

E a probabilidade do sistema estar no estado P(λ) definida anteriormente, chega-se que:

P(λ1λ2)P(λ2λ1)=eΔE/T.

uma das maneiras de se definir P que resolvem esse sistema é o algoritmo de Metropolis [1], aqui, a probabilidade de transição é:

P(λ1λ2)=min(1,eE(λ2)E(λ1)T)

Equivalentemente, gera-se um número aleatório uniforme r[0,1). A atualização é aceita se

r<eΔE/T.

essa solução garante,também, que a solução da cadeia de Markov seja exatamente a distribuição de Boltzmann.

Quando o spin aleatório sorteado for igual ao atual, o sistema deve também guardar essa mudança, isso se dá para que o sistema respeite o balanço detalhado e não assuma que a probabilidade de permanência em um estado seja nula.

Passo de Monte Carlo

Considere um sistema contendo N partículas em uma configuração λ. Em algoritmos de Monte Carlo, um passo elementar consiste na seleção de uma partícula e na tentativa de atualização de seu estado segundo a dinâmica escolhida (por exemplo, o algoritmo de Metropolis).

Como o número de partículas depende do tamanho do sistema, utilizar apenas o número de passos elementares dificulta a comparação entre simulações realizadas em redes de diferentes dimensões. Por esse motivo, define-se uma unidade de tempo computacional normalizada pelo tamanho do sistema.

O Passo de Monte Carlo (MCS) corresponde a NNN tentativas de atualização elementares. Dessa forma, após um MCS, cada partícula terá sido selecionada, em média, uma vez para possível atualização. Essa definição permite comparar de maneira consistente a evolução temporal e as grandezas medidas em sistemas de diferentes tamanhos.

Implementação Computacional

Começa-se criando uma matriz state. Como esse código será executado diversas vezes com loops longos, diversos pontos devem ser otimizado para velocidade ao invés de legibilidade. O primeiro deles é o uso de uma tabela de referência lookUp. Ela é gerada uma vez sempre que o tamanho da grade é alterado e possui 2 índices, o primeiro (i) sendo a partícula que se quer saber a posição das células adjacentes, e o segundo é um valor de 0 a 3 que guarda, em ordem, o índice da partícula imediatamente acima, abaixo, a esquerda e a direita da partícula i.

Isso também permite que seja montado o sistema como periódico, onde a primeira e a última linha seriam vizinhas, assim como a primeira e última coluna. Isso é útil para diminuir os erros que surgem nas bordas, onde as partículas não teriam vizinhas nas quatro direções.

def genLookup(side):
    index = np.arange(side)
    
    up    = (index - side) % (side**2)
    down  = (index + side) % (side**2)
    left  = np.where(index       % side == 0, index - 1 + side, index - 1)
    right = np.where((index + 1) % side == 0, index + 1 - side, index + 1)

    return np.stack((up, down, left, right))

Após isso, é pode-se iniciar o sistema com uma matriz 1D de N partículas, onde N é o quadrado do tamanho desejado de um lado da grade. Usa-se matrizes 1D ao invés de 2D para maximizar a velocidade do código.

Com o sistema iniciado, se implementa o algoritmo de Metropolis em uma função que, quando aplicada sobre o estado λi, evolui-rá ele para outro estado com uma quantidade definida de passos, podendo ser um passo fundamental, um MCS, ou qualquer outro valor conveniente.

def timeSteps(state, lookup, steps = 1, J = 1, D = 1, T = 1):
    N = state.size
    varE = 0.0

    # Repete a função para a quantidade de passos dada
    for i in range(steps):
    
        # A cada passo, escolhe uma partícula e um spin aleatórios
        index = np.random.randint(0, N)
        oldSpin = state[index]
        newSpin = np.random.randint(-1, 2)
    
        # avança a iteração se os spins forem iguais
        if newSpin == oldSpin:
            continue
        
        # Soma o estado dos vizinhos através da tabela de referência
        neiSum = np.sum(state[lookup[:, index]])

        # Calcula a mudança de energia do passo
        dE = -J * (newSpin - oldSpin) * neiSum + D * (newSpin**2 - oldSpin**2)
        
        # Verifica se o passo foi ou não aceito
        if np.random.random() < np.exp(-dE / T):
            state[index] = newSpin
            varE += dE         
            
    # Retorna o estado depois de todos os passos e a variação de energia total    
    return state, varE

Com todas essas etapas montadas, pode-se gerar a primeira simulação. Abaixo, pode-se analisar uma simulação com valores de D e T fixos e um número arbitrário de MCS

Simulação da evolução de um sistema com parâmetros aleatórios

Aqui, pode-se observar a evolução de um único sistema, mas o grande trunfo do modelo de Blume-Capel é analisar os estados do sistema, especificamente o TCP. Para isso, uma única simulação não é suficiente, mas será necessário ver o estado de equilíbrio do sistema em um gráfico de Temperatura (T)/Campo cristalino (D) e analisar a magnetização média e o momento de quadrupolo médio.

Como foi dito anteriormente, essas duas medidas serão iguais em dois casos e diferentes em um. Então o gráfico Temperatura (T)/Campo cristalino (D) irá mostrar <M> e <Q>, além de um terceiro gráfico mostrando o módulo da diferença entre eles.

Evolução da magnetização, momento de quadrupolo e a diferença de ambos em uma grande área, simulações utilizada para identificar a posição dos estados

Desenhando esses três gráficos, pode-se perceber que, onde o TCP teórico está próximo dos três estados de equilíbrio teóricos. Mas para avaliar se ele realmente é o TCP, é necessário verificar duas propriedades.

  1. Nele, o sistema não deve chegar a um estado de equilíbrio único, mas sim alternar entre os conhecidos;
  2. A vizinhança dele deve ser composta por simulações que levam a um estado bem definido.

Com isso, para avaliar o sistema, foram geradas várias simulações com diferentes valores de T e D dem volta do valor teórico para encontrar qual é o valor que se chega nessa simulação e, se for possível encontrá-lo, qual é o tamanho do ponto.

Espaços de fases Temperatura em y e D em x aproximado onde na região do PTC. A esquerda a evolução da magnetização do sistema, no meio a evolução do momento de quadrupolo e, a direita, o absoluto da diferença entre os dois. Um ponto vermelho localiza-se na previsão teórica do TCP, com três pontos azuis distribuídos ao seu redor mostrando pontos de simulação.

Definindo os pontos a serem simulados separadamente, pode-se animar a evolução do sistema com os devidos parâmetros.

Estados onde Q é maior que E, com as partículas oscilando entre os três estados aleatoriamente.
Estados com Q e M próximos de 1, com a simulação estando de maneira quase homogênea em um único spin diferente de 0.
Estados com Q e M estão próximos de 0, com a maior parte das partículas estando com spin igual a 0.

Além disso, pode-se também ver a evolução do estado no TCP:

Momento em que é observado o estado com Q e M baixos.
Momento em que é observado o estado com Q e M altos.
Momento em que é observado o estado com Q muito maior que M.
Animação total onde os gráficos foram extraídos, clique para ver animado (aviso de luzes piscantes).

Conclusão

Esse trabalho deixou claro a existência de três estados possíveis, além de mostrar que em uma interface de estados o sistema muda voluntariamente entre os três diferentes estados. pode-se ver, parando o gráfico em alguns pontos específicos, que ele passou pelos três estados.

Assim, o código aqui descrito mostra claramente uma implementação do modelo de Blume-Capel capaz de simular os três estados e um ponto onde os três oscilam entre si aleatoriamente. É importante notar que a variação entre os estados não é tão clara no TCP, pois os estados com M próximo de 0 acabam sendo semelhantes devido ao número de alterações que ocorrem a todo MCS.

Mas o algoritmo aqui descrito está disponível na página Código:Modelo de Blume-Capel feito em Jupyter Notebook, é capaz de prever de maneira satisfatória a evolução de um sistema à partir das variáveis supridas à ele.

Referências

  1. 1,0 1,1 1,2 1,3 Pathria, R. K.; Beale, Paul D.; Statistical Mechanics, third edition; Elsevier Ltd, 2011; secções 3.4, 12.3, 16.2A
  2. 2,0 2,1 Beale, Paul D.; Finite-size scaling study of the two-dimensional Blume-Capel model; APS Journals, 1986; DOI https://doi.org/10.1103/PhysRevB.33.1717