Modelo de Potts -- 2D

Introdução

O Modelo de Potts

O modelo de Potts (Potts, 1951) pode ser descrito como uma generalização do modelo de Ising (Ising, 1925) para um sistema de N spins e Q-estados (Q > 2).

Originalmente, o problema proposto por Domb era o de compreender o modelo de Ising como um sistema de spins, em que os spins podem ser paralelos ou antiparalelos. Desse modo, a generalização seria considerar que os spins estivessem sobre um plano, em que cada spin estivesse apontando para direções diferentes, igualmente espaçadas por um ângulo $\mathrm{\Theta }$ definido por Q:

${\mathrm{\Theta }}_n=\frac{2\pi n}{Q}$, em que n = (0, 1, 2, ..., Q-1).

Representação dos spins para $Q=2$.

Representação dos spins para $Q=3$.

Representação dos spins para $Q=4$.

A energia de cada interação entre dois spins é dada pelo potencial $V(s_i,s_j)=-J\delta (s_i,s_j)$, em que $J$ é a constante de interação entre os dois spins . Esse tipo de sistema tem o Hamiltoniano na forma:

${\mathscr{H}}=-J\sum _{⟨i,j⟩}\delta (s_i,s_j)$

Para o caso em que o modelo de Potts é equivalente ao modelo de Ising, temos que Q = 2, então:

${\mathrm{\Theta }}_n=\pi n$, n = (0, 1, 2, ..., Q-1)

$$

Objetivos

Neste trabalho, primeiro buscamos analisar o comportamento do sistema utilizando dois algoritmos diferentes: o Algoritmo de Metropolis-Hastings e o Algoritmo de Banho Térmico, em Q = 2. Para isso, realizamos medidas de energia e magnetização média, além da susceptibilidade magnética e calor específico. Depois, com FSS, aplicamos o método de U4 para analisar o comportamento do sistema em diferentes temperaturas, para valores diferentes de L. Por fim, procuramos analisar a susceptibilidade magnética e o calor específico bruto, também com L de diferentes valores. Além disso, expandimos a análise computacional para $Q=6$ e $Q=7$ a fim de investigar a mudança na natureza da transição de fase, observando o surgimento de transições de primeira ordem em contrapartida às transições contínuas observadas para $Q\le 4$.

Metodologia

Para realizar este trabalho, utilizamos uma rede quadrada 2D de comprimento L com N=L*L spins. A rede possui condições de contorno periódicas e vizinhança de von Neumann para primeiros vizinhos. A constante de interação é $J=1$ e $\beta =1$.

Método de Monte Carlo

No método de Monte Carlo, começamos sorteando aleatoriamente os spins da rede. Em cada passo de Monte Carlo (MCS), sorteamos spins aleatoriamente N vezes, de modo que cada spin tenha uma chance 1/N de ser sorteado. Quando sorteado, um novo valor (diferente do anterior) para esse spin é escolhido e aceito com base em critérios estabelecidos por cada algoritmo. Esse processo acontece em 1 passo de Monte Carlo e número total de passos de Monte Carlo para cada algoritmo foi escolhido separadamente, levando em conta número de passos do estado transiente em cada algoritmo.

Algoritmo de Metropolis-Hasting

(texto retirado de : )

O primeiro algoritmo utilizado para gerar as configurações do sistema foi o algoritmo de Metropolis-Hasting. O algoritmo escolhe repetidamente um novo estado para o sistema $\nu$ e aceitando ou rejeitando ele de acordo com uma probabilidade de aceitação $A(\mu \to \nu )$ de transitar de um estado antigo $\mu$ para o novo estado $\nu$. O algoritmo que iremos descrever utiliza a dinâmica de inversão única de spins.

Temos que a condição de balanceamento detalhado é dada por :

$\frac{A(\mu \to \nu )}{A(\nu \to \mu )}=e^{-\beta \mathrm{\Delta }E},(3)$

onde $\mathrm{\Delta }E=E_{\nu }-E_{\mu }$ é a diferença de energia entre o novo e o antigo estado.

Vamos supor que tenhamos os estados $\mu$ e $\nu$ e que temos a relação de energias: $E_{\mu }<E_{\nu }$. Então, a maior das duas chances de aceitação é $A(\nu \to \mu )$, portanto iremos igualar essa probabilidade a 1. Para que $(3)$ seja respeitada, iremos definir o valor de $A(\mu \to \nu )$ como $e^{-\beta \mathrm{\Delta }E}$. Temos, assim, o algoritmo de Metropolis-Hasting:

$A(\mu \to \nu )=\{\begin{array}{l}{e}^{-\beta \mathrm{\Delta }E},\text{se }\mathrm{\Delta }E>0\\ \\ 1,\text{caso contrario}.\end{array}$

Dessa forma, sempre que tivermos um estado cuja energia seja menor do que a do estado atual, iremos aceitar a transição, mas se a energia for maior, teremos uma pequena probabilidade de trocarmos de estado.

Algoritmo de Banho Térmico

(texto retirado de : )

O algoritmo de Metropolis-Hasting para inversão única de spins é eficaz para o modelo de Potts em baixos valores de $Q$ ou temperaturas acima da temperatura crítica, entretanto para valores altos de $Q$ ou baixas temperaturas o algoritmo falha em convergir o sistema rapidamente para o estado estacionário.

Considerando um caso onde $Q=100$ e um spin que possui 4 vizinhos, se todos os vizinhos do spin possuem valores diferentes uns do outro e do próprio spin, poderá levar em média $100/4=25$ passos de Monte Carlo para sortear um novo valor de $q$ que tem a transição aceita, e dessa forma o algoritmo irá demorar mais tempo para alcançar a configuração de equilíbrio do sistema. A dificuldade de aceitar transições é maior ainda para baixas temperaturas, onde a probabilidade de transicionar para um novo estado tem um peso maior para qualquer $q$ diferente dos spins da vizinhança, dessa maneira poderá demorar 96 passos para gerar um spin que seja igual a algum spin da vizinhança e realizar a transição. Para contornar este problema podemos utilizar o algoritmo de banho térmico.

O algoritmo de banho térmico, assim como o metropolis, é um algoritmo em que mudamos um spin por vez. O algoritmo segue as seguintes etapas: primeiro escolhemos um spin na rede ($i$), e independente do seu valor atual, escolhemos um novo valor para $s_i$. Esse novo valor será aceito, ou não, de acordo com os pesos de Boltzmann. Temos que no algoritmo de Banho Térmico nós atribuímos um valor $n$, entre 1 e $q$, ao spin com uma probabilidade

$A(\mu \to \nu )=\frac{e^{-\beta E_{\nu }}}{\sum _{m=1}^q{e}^{-\beta E_m}},$

onde $E_n$ é a energia do sistema quando $s_i=n$ e o somatório é dado em todas energias possíveis. Pelo fato desse algoritmo permitir que o spin assuma qualquer valor ele satisfaz a condição de ergodicidade. Temos que $A(\mu \to \nu )=a_{\nu }$ e $A(\nu \to \mu )=a_{\mu }$, assim, pela descrição do algoritmo temos

$\frac{A(\mu \to \nu )}{A(\nu \to \mu )}=\frac{a_{\nu }}{a_{\mu }}=\frac{e^{-\beta E_{\nu }}}{\sum _{m=1}^q{e}^{-\beta E_m}}\times \frac{\sum _{m=1}^q{e}^{-\beta E_m}}{e^{-\beta E_{\mu }}}=e^{-\beta \mathrm{\Delta }E}.$

Ou seja, o algoritmo de banho térmico respeita a condição de balanço detalhado.

Veremos que para valores pequenos de $q$ (como $q=2$, que é o modelo de Ising), o algoritmo de Metropolis é mais eficiente. Porém, para valores altos de $q$ ou altas temperaturas o algoritmo de banho térmico é mais eficiente.

O Parâmetro de Ordem e Transições de Fase

No modelo de Ising (Q = 2), a magnetização média do sistema atua como um parâmetro de ordem natural. Contudo, para o modelo de Potts com Q > 2, calcular a média escalar dos spins perde o sentido físico, uma vez que os estados não representam intensidades, mas sim simetrias discretas. Portanto, para Q >= 3, substituímos a magnetização clássica por um Parâmetro de Ordem (m) baseado no estado majoritário da rede:

$m=\frac{Q{N}_{max}-N}{N(Q-1)}$

Onde N_max é o número de sítios no estado mais populoso e N é o número total de sítios. Este parâmetro vai a zero na fase desordenada e a 1 na fase perfeitamente ordenada. Adicionalmente, implementamos análises de histograma de energia no estado estacionário para investigar a ordem da transição de fase. Enquanto transições de segunda ordem (Q <= 4) apresentam flutuações em torno de um único pico de energia, transições de primeira ordem (Q > 4) devem apresentar uma distribuição bimodal na temperatura crítica, indicando a coexistência de fases (ordenada e desordenada).

Resultados

Comparando os Algoritmos em Q=2

Aqui se encontram os resultado para Q=2. As simulações foram realizadas para MCS = 50.000 passos. O transiente para o banho térmico foi de 22.000 passos, enquanto para o algoritmo de Metropolis foram de 10.000 passos, garantindo que estávamos no estado transiente em ambos os algoritmos ao tirarmos as medidas. Os valores de energia e magnetização usados foram obtidos a partir da média dos valores de energia e magnetização no estado estacionário, respectivamente.

Medidas para Q = 2 e L = 64 utilizando os algoritmos de Metropolis-Hastings e Banho Térmico.

Aqui observamos que o modelo realmente se aproxima do modelo de ISING quando Q=2. Primeiramente, vemos que a temperatura crítica do Potts Q=2 se aproxima à temperatura crítica do modelo de ISING, que sabemos que é $T_C\approx 1.1346$. Além disso, vemos que na energia média os dois algoritmos apresentam resultados praticamente idênticos, com exceção de T=0.7. No gráfico de magnetização média, ambos os algoritmos apresentam grandes variações num período inicial de T, mas logo após atingir a temperatura crítica, os algoritmos se tornam quase indistinguíveis e mantém valores muito próximos até T = 2. Por fim, para o calor específico e a susceptibilidade magnética também vemos que os dois algoritmos se assemelham na maior parte do intervalo da temperatura, com exceção dos picos apresentados em cada gráfico que vemos que os algoritmos se diferenciam próximos à região da temperatura crítica.

Discrepâncias com o Banho Térmico

Além disso, produzimos alguns gif's que mostram o estado da rede em snapshots tiradas a cada 100 passos, a fim de compreender os pontos discrepantes nos nossos gráficos. Observamos que esses pontos apareceram aleatoriamente em diferentes T's (menores que $T_C$) nas simulações que realizamos (que não estão expostas neste texto). Abaixo, podemos comparar o estado da rede ao longo do tempo em T = 0.7 com estado da rede em T = 0.6 e T=0.8, ambos com o algoritmo de Banho Térmico:

Snapshots da rede para T= 0.6, 0.7, 0.8 com Banho Térmico.

Observamos que, diferentemente do que acontece em T=0.6 e T=0.8, em T=0.7 o sistema atinge um estado em que os dois grandes domínios de spins coexistem. Isso esclarece o porquê da energia média em 0.7 ser um pouco maior com o banho térmico do que com o metropolis nesse ponto. Além disso, também entendemos a magnetização média ser $\approx 0.5$ em T=0.7, ao invés de 1 ou 0 como o que é observado para as simulações com $T<T_C$. Da mesma forma, do gráfico da susceptibilidade magnética podemos entender o que acontece em T=0.7. Como a susceptibilidade magnética é alta, a nossa hipótese é de que isso acontece aleatoriamente, e já que o algoritmo de Banho Térmico atua localmente, ele não consegue desmanchar completamente esses domínios se o sistema atinge esse ponto em $T<T_C$. acreditamos que isso pode ser contornado aumentando o tempo total de simulação, esperando que eventualmente um desses domínios tome conta do sistema.

Não foi possível colocar aqui os gif's para T's maiores, mas podemos dizer que para T próximo da temperatura crítica, os spins que são a maioria na rede trocam: em dado momento são os spins 0, enquanto em outros momentos são os spins 1 que estão em maioria. Para temperaturas muito altas, o sistema tem aproximadamente o mesmo número de spins de cada tipo, de modo que a magnetização média é próxima de 0.

Ademais, inicialmente havíamos usado o intervalo de transiente para o algoritmo de Metropolis-Hastings maior que o do Banho Térmico, pois percebemos que o estado estacionário demorava mais para ser atingido no primeiro algoritmo. Contudo, depois de perceber estados como o encontrado em T=0.7 no Banho Térmico, optamos por aumentar o transiente para esse algoritmo, mas ressaltamos que (com exceção do T=0.7) o Banho Térmico atingiu o estado estacionário mais rápido.

Comparação entre os casos Q <= 4

Aqui, comparamos os resultados obtidos para Q = 2, 3 e 4. Analisamos as grandezas termodinâmicas do sistema (Energia Média, Parâmetro de Ordem de Potts, Susceptibilidade Magnética e Calor Específico) utilizando o algoritmo de Metropolis-Hastings para uma rede de tamanho L = 64, no intervalo de temperatura de 0.5 <= T <= 2.0.

Nos gráficos, observamos que o aumento do número de estados Q desloca o ponto de transição de fase (temperatura crítica) para valores menores: temos $T_c\approx 1.13$ para Q = 2, $T_c\approx 1.0$ para Q = 3 e $T_c\approx 0.9$ para Q = 4.

Na análise do Parâmetro de Ordem, vemos claramente a quebra de simetria do sistema. Em baixas temperaturas, o parâmetro tende a 1, indicando que a grande maioria dos spins da rede se alinhou em um único estado majoritário (fase ordenada). À medida que a temperatura cruza a temperatura crítica, a agitação térmica destrói essa ordem e o parâmetro cai rapidamente para zero, caracterizando a fase desordenada, onde todos os Q estados têm populações equivalentes. É importante notar que, para todos esses três casos, a transição ocorre de maneira contínua (transição de segunda ordem).

O comportamento das flutuações corrobora essa continuidade. No gráfico de calor específico, os picos nos pontos de transição de fase se tornam cada vez mais estreitos e acentuados à medida que Q aumenta. Isso nos indica que a energia necessária para o sistema realinhar seus spins e destruir a ordem global cresce significativamente com o número de estados possíveis. A susceptibilidade magnética acompanha esse padrão, divergindo próximo à temperatura crítica devido ao surgimento de domínios de spins de todos os tamanhos, uma marca registrada de sistemas críticos. A energia do sistema reflete essa física com uma curvatura que se torna cada vez mais íngreme perto da temperatura crítica para valores maiores de Q, embora ainda sem apresentar a descontinuidade abrupta característica de transições de primeira ordem.

FSS

O Finite Size Scaling (FSS) é uma ferramenta computacional que nos permite pegar uma rede de tamanho finito (como L = 64) e descobrir como ela se comportaria no limite termodinâmico quando L tende a infinito. Em Q = 2, sabemos que a susceptibilidade e o calor específico deveriam ir para o infinito quando T = Tc. Contudo, como a nossa rede tem tamanho finito, observamos que em Q = 2 não é isso que acontece. Portanto, usamos essa ferramenta nas simulações a fim de verificar estes resultados já conhecidos na literatura.

U4

O cumulante de Binder, ou U4, é uma ferramenta estatística que mede o 'quarto momento' das flutuações do parâmetro de ordem da rede. Matematicamente, ele é definido pela relação entre a média do momento de quarta ordem e o quadrado do momento de segunda ordem do parâmetro de ordem $m$:

$U_4=1-\frac{⟨m^4⟩}{3⟨m^2{⟩}^2}$

Fisicamente, essa grandeza quantifica o formato da distribuição de probabilidade das flutuações, indicando se o sistema está bem ordenado ou completamente caótico. Quando a temperatura é baixa (sistema perfeitamente ordenado), a medida para todos os tamanhos de rede tende a 2/3. Em temperaturas altas (o sistema é um caos desordenado), a medida vai a zero.

Exatamente na transição de fase, o valor do cumulante não depende do tamanho L da rede, o que significa que as curvas plotadas para redes de tamanhos diferentes irão convergir e se cruzar em um único ponto: a temperatura crítica.

Para expandir a validação da ferramenta, aplicamos o mesmo método para estados de Potts maiores. A figura abaixo demonstra a eficácia do U4 em localizar o ponto crítico para Q=3, onde o cruzamento das curvas ocorre exatamente na temperatura teórica esperada ($T_c\approx 0.9950$), validando o método para transições contínuas antes de aplicá-lo em regimes mais complexos.

Aprofundando a análise para Q=4, que representa o caso marginal exato entre transições contínuas e de primeira ordem no modelo de Potts 2D, a ferramenta U4 novamente demonstra grande precisão. O cruzamento ocorre na temperatura esperada (Tc ≈ 0.9102). Notavelmente, a queda do parâmetro de ordem logo após o ponto crítico é muito mais abrupta, exibindo um mergulho profundo para a rede maior (L=64), o que reflete as fortes flutuações e correções logarítmicas características deste regime limite.

Susceptibilidade Magnética

A susceptibilidade magnética mede o quão 'sensível' o sistema é. Se fosse aplicado um campo magnético minúsculo, o quanto o sistema responderia mudando sua magnetização? Estatisticamente, isso é calculado através das flutuações (variância) da magnetização. O que acontece é que longe da temperatura crítica, os spins estão ou muito ordenados ou muito caóticos. Em ambos os casos, as flutuações globais são pequenas, então a susceptibilidade é baixa. Perto da temperatura crítica, o sistema não está nem muito ordenado, nem muito caótico, então flutuam violentamente. Por isso, a susceptibilidade é altíssima e diverge para o infinito na temperatura crítica. Como na simulação a rede é finita, o comprimento de correlação (tamanho de domínios de spins que flutuam muito) tenta crescer, mas 'bate nas paredes' que tem tamanho finito L da rede. É por isso que o pico não vai propriamente a infinito.

A teoria diz que, muito perto do ponto crítico, a física do sistema 'esquece' os detalhes microscópicos da rede e passa a depender somente da proporção entre o tamanho do sistema e o comprimento de correlação. Nós sabemos que no ponto crítico o pico máximo da susceptibilidade deve crescer como $L^{\gamma /\nu }$. Então se nós pegarmos o eixo Y e dividirmos por $L^{\gamma /\nu }$, estamos achatando os picos na mesma proporção. Da mesma forma, estivamos o eixo X multiplicando por $L^{1/\nu }$.

Calor Específico

O calor específico mede a a variância de energia. Perto da temperatura crítica, como os spins estão virando em blocos massivos, a energia total do sistema oscila drasticamente. A variância da energia sobe, criando o pico do calor específico. A teoria nos diz que o expoente crítico que governa o calor específico é o $\alpha$. Para a imensa maioria dos modelos 3D, $\alpha$ é um número positivo, fazendo o pico crescer muito rápido com L. Mas, incrivelmente, para o modelo de Potts Q=2 (Ising) em duas dimensões, $\alpha =0$. Quando isso acontece, a divergência não é uma lei de potência normal, mas sim logarítmica. É por isso que, quando olhamos o gráfico de $C_v$ o pico da rede de 64 é maior que o da rede 16, mas ele cresce de forma mais suave, proporcional ao $\ln (L)$.

A Transição de Primeira Ordem em Q > 4

Conforme previsto analiticamente para o modelo de Potts em duas dimensões, o comportamento crítico do sistema muda drasticamente para Q > 4, passando de uma transição contínua para uma transição de primeira ordem. Evitamos simular o caso Q = 5, pois ele apresenta uma transição fracamente de primeira ordem cujo comprimento de correlação excede significativamente o tamanho finito da nossa rede, comportando-se pseudo-continuamente nas simulações. Focamos a análise em Q = 6 e Q = 7.

Como o sistema passa por uma transição abrupta de primeira ordem, o comportamento do cumulante de Binder ($U_4$) muda radicalmente em relação ao regime contínuo. Em vez de convergirem para um ponto de cruzamento fixo e positivo perto de 2/3, as curvas sofrem uma deformação severa na região crítica, desenvolvendo um poço pronunciado que pode atingir valores negativos.

Esse comportamento do cumulante é o reflexo direto da coexistência de fases na temperatura crítica. A flutuação do sistema entre estados macroscópicos muito distintos (ordenado e desordenado) altera drasticamente a distribuição estatística do parâmetro de ordem, fazendo com que o momento de quarta ordem ($⟨m^4⟩$) cresça desproporcionalmente em relação ao quadrado do segundo momento ($⟨m^2{⟩}^2$). Conforme demonstrado nos gráficos, esse efeito é discreto em Q = 6 devido ao seu caráter fracamente de primeira ordem, mas torna-se violento e nitidamente negativo em Q = 7, onde a profundidade do poço diverge com o aumento do tamanho da rede L, servindo como uma evidência numérica robusta da mudança de regime termodinâmico.

A prova definitiva da mudança da natureza da transição se encontra na análise de energia. Em uma transição de primeira ordem, há liberação de calor latente, o que implica que o sistema pode coexistir em duas fases distintas exatamente na temperatura crítica.

Como pode ser observado na figura acima, a rede de Q = 7 exibe uma clara distribuição bimodal, em que a energia do sistema flutua entre dois estados macroscópicos bem definidos (um de menor energia, ordenado, e outro de maior energia, desordenado). Isso confirma o forte caráter de primeira ordem da transição. Curiosamente, para Q = 6, observamos apenas um pico alargado e assimétrico. Isso ocorre devido a efeitos de tamanho finito: como a transição em Q = 6 é fracamente de primeira ordem, o comprimento de correlação do sistema ainda excede o tamanho da nossa rede (L = 64), fazendo com que ele exiba um comportamento pseudo-contínuo. Para observar a separação bimodal nítida em Q = 6, seria necessária a simulação de redes significativamente maiores.

Conclusões

A partir da simulação do modelo de Potts em duas dimensões, foi possível validar a eficácia e as diferenças entre os algoritmos de Metropolis e Banho Térmico. A comparação evidenciou que, enquanto o algoritmo de Metropolis é adequado para valores baixos de Q (como Q = 2) e altas temperaturas, o Banho Térmico se mostrou mais vantajoso para contornar a baixa taxa de aceitação em regimes de alta frustração (Q elevados ou baixas temperaturas), convergindo mais rapidamente para o equilíbrio. Contudo, observamos alguns problemas nas simulações com o Banho Térmico em Q's e T's baixos, mas que acreditamos que podem ser contornados se o tempo de simulação fosse maior.

A análise feita para Q <= 4 confirmou as expectativas teóricas, mostrando que o aumento do número de estados desloca a temperatura crítica para valores menores e torna as curvas de transição de fase mais abruptas. O modelo Q = 2 reproduziu o que esperávamos do modelo de Ising, com a temperatura crítica convergindo para o valor analítico já conhecido.

Além disso, a aplicação do Finite Size Scaling contornou as limitações da rede finita. O cruzamento do U4 confirmou o ponto crítico independentemente de L. O colapso dos dados de susceptibilidade magnética validou a relação de escala com os expoentes críticos, enquanto a análise do calor específico para Q = 2 mostrou a divergência logarítmica característica em vez de uma lei de potência padrão.

Por fim, a expansão da análise para Q = 6 e Q = 7 demonstrou a limitação da magnetização escalar em descrever quebras de simetria múltiplas, exigindo o uso de um parâmetro de ordem baseado na densidade populacional de estados. Mais importante, as simulações capturaram com sucesso a mudança na natureza termodinâmica do sistema: enquanto o histograma bimodal nítido em Q = 7 evidenciou a coexistência de fases característica de transições de primeira ordem, o pico único e alargado em Q = 6 ilustrou de forma empírica as limitações de tamanho finito em transições fracamente de primeira ordem. Em conjunto, esses resultados confirmam que o modelo de Potts bidimensional exibe transições de primeira ordem para Q > 4, rompendo com o regime de transições contínuas observado para Q <= 4. De maneira geral, os métodos numéricos empregados se mostraram de acordo com o esperado pela literatura.

Código

Para gerar estas simulações, foi produzido um código em Python3. Utilizamos as bibliotecas seguintes bibliotecas: Numpy, em que usamos os gerador de números aleatórios e as arrays; Numba, para otimizar as funções e acelerar o código; e Matplotlib, que utilizamos para gerar os códigos.

Link: https://colab.research.google.com/drive/1V7GsP-fEIkhU29QCPvh_UcZKIXtupR7p?usp=sharing