Modelo Blume-Capel para partículas de spin 1

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 \sum_{< i, j >} σ_{i} σ_{j}$

Com $σ_{i}$ representando o valor do spin da partícula $i$ , podendo ser $\pm 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á

$ℋ_{M a g} = - H \sum_{n = 1}^{N} σ_{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 $- H N$ , onde todas as $N$ partículas estão com o spin alinhado.

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

$ℋ = - J \sum_{< i, j >} σ_{i} σ_{j} - H \sum_{n = 1}^{N} σ_{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 (λ) \propto e^{- β E (λ)}$

com $β = \frac{1}{k_{b} T}$ . Considerando que a soma de todas as probabilidades é 1, deves-se dividir $e^{- β E}$ pela soma de $e^{- β E_{n} (λ)}$ 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 (λ) = \frac{e^{- β E (λ)}}{{\sum_{λ}}^{'} 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 = \frac{1}{N} \sum_{n = 1}^{N} σ_{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} \in {- \frac{1}{2}, \frac{1}{2}}$ , trabalha-se com partículas de spin 1 com $σ_{i} \in {- 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 = \sum_{n = 1}^{N} σ_{n}^{2}$

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 \sum_{< i, j >} σ_{i} σ_{j} - H \sum_{n = 1}^{N} σ_{n} + D \sum_{m = 1}^{N} σ_{m}^{2}$

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 \sum_{< i, j >} σ_{i} σ_{j} + D \sum_{m = 1}^{N} σ_{m}^{2}$

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 / J > > 0$ : 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 / J > > 0$ e $D / J$ coerente: Aqui, analisando-se $P (λ)$ , como $β \propto T^{- 1}$ , 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 $M \approx 0$ 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 $T \approx 0.610 \pm 0.005$ e $D \approx 1.965 \pm 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 $k_{B} = 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^{*} = \frac{k_{B} T}{J};$ $D^{*} = \frac{D}{J};$ $H^{*} = \frac{H}{J} .$

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:

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 $σ \in {- 1, 0, 1} .$
Calcula-se a variação de energia associada à mudança proposta: $Δ E = E_{novo} - E_{atual} .$

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 (i \to j) = P (j) W (j \to i),$

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

\frac{P (λ_{1} \to λ_{2})}{P (λ_{2} \to λ_{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} \to λ_{2}) = m i n (1, e^{\frac{E (λ_{2}) - E (λ_{1})}{T}})

Equivalentemente, gera-se um número aleatório uniforme $r \in [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

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.

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 $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.

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

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

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,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,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

[mecEst-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

[point-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

[1]

[2]

Modelo Blume-Capel para partículas de spin 1

Índice

Modelo de Ising