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

${\mathscr{H}}=-J\sum _{<i,j>}{\sigma }_i{\sigma }_j$

Com ${\sigma }_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á

${{\mathscr{H}}}_{Mag}=-H\sum _{n=1}^N{\sigma }_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á:

${\mathscr{H}}=-J\sum _{<i,j>}{\sigma }_i{\sigma }_j-H\sum _{n=1}^N{\sigma }_n$

Agora, consideremos o ensemble canonico, isso é, o conjunto de todos os estados $\lambda$ 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$ é:

$P(\lambda )\propto e^{-\beta E(\lambda )}$

com $\beta =\frac{1}{k_{b}T}$. Considerando que a soma de todas as probabilidades é 1, deves-se dividir $e^{-\beta E}$ pela soma de $e^{-\beta E_n(\lambda )}$ para todos os estados $\lambda$ ao qual o sistema pode ser encontrado. Importante ressaltar que soma é sobre $\lambda$ 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(\lambda )=\frac{e^{-\beta E(\lambda )}}{{\sum _{\lambda }}^{\prime }{e}^{-\beta E({\lambda }^{\prime })}}$

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{\sigma }_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 ${\sigma }_i\in \{-\frac{1}{2},\frac{1}{2}\}$, trabalha-se com partículas de spin 1 com ${\sigma }_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 ${\mathscr{H}}$. 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, definido como:

$Q=\sum _{n=1}^N{\sigma }_n^2$

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

${\mathscr{H}}=-J\sum _{<i,j>}{\sigma }_i{\sigma }_j-H\sum _{n=1}^N{\sigma }_n+D\sum _{m=1}^N{\sigma }_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:

${\mathscr{H}}=-J\sum _{<i,j>}{\sigma }_i{\sigma }_j+D\sum _{m=1}^N{\sigma }_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(\lambda )$, como $\beta \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(\lambda )$, estados onde a energia do sistema seja menor que zero serão favorecidos, pois estados positivos terão uma probabilidade menor que 1 ($e^{-{\mathscr{H}}/T<1}$) e estados de energia menor que zero terão uma probabilidade maior ($e^{-{-}{|}{\mathscr{H}}{|}/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$

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:

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 $\sigma \in \{-1,0,1\}.$
3. Calcula-se a variação de energia associada à mudança proposta: $\mathrm{\Delta }E=E_{\text{novo}}-E_{\text{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(\lambda )$ definida anteriormente, chega-se que:

$\frac{P({\lambda }_1\to {\lambda }_2)}{P({\lambda }_2\to {\lambda }_1)}=e^{-\mathrm{\Delta }E/T}.$

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

$P({\lambda }_1\to {\lambda }_2)=min(1,e^{\frac{E({\lambda }_2)-E({\lambda }_1)}{T}})$

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

$r<e^{-\mathrm{\Delta }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 $\lambda$. 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 ${\lambda }_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

[IMAGEM]

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. Para tornar mais claro o terceiro caso, onde o que importa é que $<M>$ e $<Q>$ sejam diferentes, também será desenhado o gráfico do módulo da diferença entre eles.

[IMAGEM]

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.

[IMAGEM]

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 diferentes estados. Porém, não foi possível determinar com certeza o TCP, pois, enquanto a borda entre o estado estático com spins fica claro quando comparado com os outros estados, o estado sem spin e com as partículas aleatórias acabam tendo uma divisão mais gradual.

Mesmo assim, o código aqui descrito mostra claramente uma implementação do modelo de Blume-Capel capaz de simular os três estados de maneira funcional.