Ising

1. Introduction

Le modèle d'Ising ([1]) a été introduit pour expliquer la physique des matériaux ferromagnétiques, en particulier la transition de phase ferromagnétique/paramagnétique se produisant à la température de Curie.

Ce document présente une simulation numérique du modèle d'ising à une dimension, par la méthode de Monte-Carlo Metropolis. Comme la solution exacte est connue, une comparaison avec celle-ci pourra être faite.

2. Modèle d'Ising

2.a. Définition

Le ferromagnétisme est causé par les moments magnétiques de spin des électrons. Dans le modèle d'Ising, ces moments magnétiques sont supposés tous alignés dans une direction. Soit S_i le spin d'un atome dans cette direction, qui peut prendre les valeurs -1/2 ou 1/2. Le spin d'un atome peut interagir avec les spins des atomes voisins du réseau cristallin. Il est par ailleurs sensible au champ magnétique extérieur dont la composante dans la direction des spins est noté B.

L'énergie du système de N spins s'écrit :

E = - g μ_{B} B \sum_{i = 0}^{N - 1} S_{i} - J \sum_{i, j} S_{i} S_{j} (1)

Le premier terme représente l'interaction entre les spins et le champ extérieur. μ_B est le magnéton de Bohr et g le facteur de Landé. Cette énergie est minimisée lorsque tous les spins sont alignés avec le champ, c'est-à-dire lorsqu'ils sont tous positifs.

Le second terme représente l'interaction entre les spins des atomes voisins. La somme porte sur les paires d'atomes voisins. J est une constante positive appelée constante de couplage. Cette énergie est minimisée lorsque les spins voisins sont dans le même sens, donc lorsque tous les spins sont dans le même sens. Ce terme de couplage entre spins explique la possibilité d'une aimantation dans le matériau ferromagnétique, en l'absence de champ extérieur.

Le moment magnétique total du système de spins est :

M = g μ_{B} \sum_{i = 0}^{N - 1} S_{i} (2)

Le système est en contact avec un thermostat de température T. La probabilité d'une configuration particulière des spins parmi toutes les configurations possibles est proportionnelle au facteur de Boltzmann :

exp (- \frac{E_{μ}}{k_{B} T}) = e^{- β E_{μ}} (3)

La valeur moyenne (l'espérance) du moment magnétique est donc :

< M > = \frac{\sum_{μ} M_{μ} e^{- β E_{μ}}}{\sum_{μ} e^{- β E_{μ}}} (4)

où μ désigne une configuration des spins.

2.b. Modèle d'Ising à une dimension

Les spins sont disposés sur une droite, régulièrement espacés. Le spin i a deux voisins avec lequel il interagit, le spin i-1 et le spin i+1. On utilise une condition limite périodique, qui consiste à considérer que le spin 0 interagit avec le spin N-1 et le spin 1, et que le spin N-1 interagit avec le spin N-2 et le spin 0. Cela revient à placer les spins sur un cercle.

L'énergie s'écrit :

E = - \sum_{i = 0}^{N - 1} g μ_{B} B S_{i} + J S_{i} S_{i + 1} (5)

L'expression exacte du moment magnétique moyen est obtenue par le calcul ([1]) :

< M > = \frac{N g μ_{B}}{2} \frac{s h (β \frac{g μ_{B} B}{2})}{\sqrt{{s h}^{2} (β \frac{g μ_{B} B}{2}) + e^{- β J}}} (6)

3. Méthode de Monte-Carlo Metropolis

3.a. Algorithme

On applique l'algorithme de Metropolis au modèle d'Ising à une dimension. On commence par introduire des grandeurs réduites (sans dimensions). La constante de couplage est prise égale à 1 et l'énergie s'écrit :

E = - \sum_{i = 0}^{N - 1} B S_{i} + S_{i} S_{i + 1} (7)

Pour simplifier, les valeurs possibles des spins sont -1 et 1 (au lieu de 1/2). Tout cela revient à faire les substitutions suivantes :

\begin{matrix} \frac{1}{2} g μ_{B} B \to B & (8) \\ \frac{1}{4} J \to 1 & (9) \end{matrix}

L'algorithme consiste à construire une chaîne d'états de spins. La construction d'un état à partir du précédent se fait de la manière suivante :

Sélection d'un spin au hasard et calcul de la variation d'énergie ΔE qui résulterait du changement de signe de ce spin.
Si ΔE<=0, le spin est changé de signe.
Si ΔE>0, un nombre aléatoire x est tiré avec une probabilité uniforme dans l'intervalle [0,1[. Si $x < e^{- β Δ E}$ , le spin est changé. Dans le cas contraire, le spin n'est pas changé.

Pour chaque nouvelle configuration, on calcule le moment magnétique totale M_k. Lorsque P configurations de spin ont été obtenues de cette manière, on calcule une estimation de la valeur moyenne par :

M_{P} = \frac{1}{P} \sum_{k = 1}^{P} M_{k} (10)

Voyons comment est modifiée l'énergie du système lorsqu'un spin change de signe. Le terme d'interaction avec le champ extérieur est changé de +2BS où S est l'état du spin avant le changement. Pour le terme d'interaction entre spins voisins, il y a trois cas suivant l'état des spins avant le changement :

Les deux spins voisins S_i-1 et S_i+1 sont de même signe que S_i. La variation d'énergie est +4.
Les deux spins voisins sont de même signe mais de signe opposé à S_i. La variation d'énergie est -4.
Les deux spins voisins sont de signes opposés. La variation d'énergie est nulle.

Finalement, il y a 6 variations d'énergie totale possibles :

2 B + 4, 2 B - 4, 2 B, - 2 B + 4, - 2 B - 4, - 2 B (11)

Le facteur e^-βΔE de ces 6 termes devra être calculé au préalable et stocké dans une liste.

3.b. Implémentation

On définit une classe dont le constructeur prend en argument le nombre de spins. Les spins sont initialisés à 1, ce qui correspond à une aimantation complète.

ising1D.py

import numpy
import numpy.random
import random
import math

class Ising1D:
    def __init__(self,N):
        self.N = N
        self.spin = numpy.ones(N)
        self.E = 0
        self.M = N

La première fonction permet d'affecter des valeurs au champ magnétique extérieur et à la température. Les facteurs de Boltzmann sont calculés et stockés.

    def constantes(self,B,T):
        beta = 1.0/T
        self.dE = [2*B+4,2*B-4,2*B,-2*B+4,-2*B-4,-2*B]
        self.A = []
        for k in range(6):
            self.A.append(math.exp(-beta*self.dE[k]))

La fonction suivante renvoit l'état des deux spins voisins, en tenant compte des conditions limites périodiques :

    def voisins(self,i):
        if i==0:
            s1 = self.spin[self.N-1]
        else:
            s1 = self.spin[i-1]
        if i==self.N-1:
            s3 = self.spin[0]
        else:
            s3 = self.spin[i+1]
        return (s1,s3)

La fonction suivante détermine la variation d'énergie et le facteur de Botzmann pour le basculement d'un spin :

    def deltaE(self,i):
        s2 = self.spin[i]
        (s1,s3) = self.voisins(i)
        if s2==1:
            p = 0
        else:
            p = 3
        if (s1,s2,s3)==(1,1,1) or (s1,s2,s3)==(-1,-1,-1):
            pass
        elif (s1,s2,s3)==(1,-1,1) or (s1,s2,s3)==(-1,1,-1):
            p += 1
        else:
            p += 2
        return (self.dE[p],self.A[p])

La fonction suivante effectue un pas de l'algorithme de Metropolis.

    def _metropolis(self):
        i = random.randint(0,self.N)
        (dE,A)= self.deltaE(i)
        if dE<=0:
            self.spin[i] = -self.spin[i]
            self.E += dE
            self.M += 2*self.spin[i]
        else:
            x = random.random()
            if x<A:
                self.spin[i] = -self.spin[i]
                self.E += dE
                self.M += 2*self.spin[i]

La fonction suivante effectue une boucle de calcul et calcule la moyenne du moment magnétique et l'écart-type. Le nombre d'itérations est donné en argument. Une itération est constituée de N changement (ou tentative) de spin, ce qui fait en moyenne un changement pour chaque spin du système.

    def _boucle(self,n):
        m = numpy.zeros(n)
        for k in range(n):
            m[k] = self.M
            for i in range(self.N):
                self.metropolis()
        return (m,numpy.mean(m),numpy.std(m))

L'écart-type permet de mesurer les fluctuations autour de la valeur moyenne. En le divisant par la racine carrée du nombre de tirages, on obtient l'écart-type de l'estimation de la moyenne.

L'efficacité des deux fonctions précédentes peut être améliorée en générant les nombres aléatoires avant la boucle. On utilise pour cela les fonction numpy.random.randint et numpy.random.random_sample.

Voici les deux nouvelles fonctions :

    def metropolis(self,i,x):
        (dE,A)= self.deltaE(i)
        if dE<=0:
            self.spin[i] = -self.spin[i]
            self.E += dE
            self.M += 2*self.spin[i]
        else:
            if x<A:
                self.spin[i] = -self.spin[i]
                self.E += dE
                self.M += 2*self.spin[i]
                
    def boucle(self,n):
        m = numpy.zeros(n)
        tab_i = numpy.random.randint(0,self.N,size=n*self.N)
        tab_x = numpy.random.random_sample(size=n*self.N)
        j = 0
        for k in range(n):
            m[k] = self.M
            for i in range(self.N):
                self.metropolis(tab_i[j],tab_x[j])
                j += 1
        return (m,numpy.mean(m),numpy.std(m))

3.c. Résultats

import numpy
import random
import math
from matplotlib.pyplot import *
import ising1D

Voyons l'évolution du moment magnétique pour une température donnée, à champ magnétique nul :

N=1000
ising = ising1D.Ising1D(N)
ising.constantes(0.0,1.0)
(m,M,dM)=ising.boucle(100)
figure(figsize=(10,6))
plot(m)
xlabel('iteration')
ylabel('M')
axis([0,100,-N,N])
grid()

figA.pdf

Le système évolue vers un état d'équilibre où le moment magnétique fluctue autour de zéro. On augmente la température d'un facteur 10 :

ising.constantes(0.0,10.0)
(m,M,dM)=ising.boucle(100)
figure(figsize=(10,6))
plot(m)
xlabel('iteration')
ylabel('M')
axis([0,100,-N,N])
grid()

figB.pdf

L'amplitude des fluctuations est plus faible. Voyons une température très faible :

ising.constantes(0.0,0.1)
(m,M,dM)=ising.boucle(100)
figure(figsize=(10,6))
plot(m)
xlabel('iteration')
ylabel('M')
axis([0,100,-N,N])
grid()

figC.pdf

Les fluctuations autour de zéro sont plus amples, mais il n'apparait pas de momemt magnétique. Effectivement, le modèle d'Ising à une dimension ne présente pas de transition vers le ferromagnétisme et ne peut donc expliquer l'aimantation permanente.

On revient à la température la plus élevée et on ajoute un champ magnétique :

ising.constantes(5.0,10.0)
(m,M,dM)=ising.boucle(100)
figure(figsize=(10,6))
plot(m)
xlabel('iteration')
ylabel('M')
axis([0,100,-N,N])
grid()

figD.pdf

Une simulation complète consiste à faire varier un paramètre, ici le champ magnétique. Pour chaque valeur du paramètre, on procède à une mise en équilibre du système puis on calcule les moyennes et écart-types sur un nombre suffisant d'itérations. On trace aussi la courbe donnée par la relation (6) pour comparaison.

bmax = 20.0
B = numpy.arange(0.0,bmax,bmax/10)
n = B.size
listM = numpy.zeros(n)
listDM = numpy.zeros(n)
T = 10.0
N=1000
ising = ising1D.Ising1D(N)
ising.constantes(0.0,T)
ising.boucle(50)
for k in range(n):
    ising.constantes(B[k],T)
    ising.boucle(10)
    (m,M,dM)=ising.boucle(100)
    listM[k] = M
    listDM[k] = dM
figure(figsize=(8,6))
errorbar(B,listM,yerr=listDM,fmt=None)
axis([0,bmax,0,N])
ylabel("M")
xlabel("B")
title("T=10")
grid()
def moment(B,T):
    sh = math.sinh(B/(1*T))
    return N*sh/math.sqrt(sh*sh+math.exp(-4.0/(T)))
B = numpy.arange(0.0,bmax,bmax/100)
M = numpy.zeros(B.size)
for k in range(B.size):
    M[k] = moment(B[k],T)
plot(B,M)

figE.pdf

Le moment magnétique total augmente avec le champ jusqu'à la saturation, qui est obtenue lorsque tous les spins sont dans le sens du champ.

On refait les mêmes calculs avec une température plus faible :

bmax = 20.0
B = numpy.arange(0.0,bmax,bmax/10)
n = B.size
listM = numpy.zeros(n)
listDM = numpy.zeros(n)
T = 5.0
N=1000
ising = ising1D.Ising1D(N)
ising.constantes(0.0,T)
ising.boucle(50)
for k in range(n):
    ising.constantes(B[k],T)
    ising.boucle(10)
    (m,M,dM)=ising.boucle(100)
    listM[k] = M
    listDM[k] = dM
figure(figsize=(8,6))
errorbar(B,listM,yerr=listDM,fmt=None)
axis([0,bmax,0,N])
ylabel("M")
xlabel("B")
title("T=5")
grid()
B = numpy.arange(0.0,bmax,bmax/100)
M = numpy.zeros(B.size)
for k in range(B.size):
    M[k] = moment(B[k],T)
plot(B,M)

figF.pdf

L'augmentation du moment magnétique est plus rapide, car l'effet de l'agitation thermique est moindre.

Modèle d'Ising à une dimension