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.
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 Si 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 :
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 :
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 :
La valeur moyenne (l'espérance) du moment magnétique est donc :
où μ désigne une configuration des spins.
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 :
L'expression exacte du moment magnétique moyen est obtenue par le calcul ([1]) :
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 :
Pour simplifier, les valeurs possibles des spins sont -1 et 1 (au lieu de 1/2). Tout cela revient à faire les substitutions suivantes :
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 :
Pour chaque nouvelle configuration, on calcule le moment magnétique totale Mk. Lorsque P configurations de spin ont été obtenues de cette manière, on calcule une estimation de la valeur moyenne par :
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 :
Finalement, il y a 6 variations d'énergie totale possibles :
Le facteur e-βΔE de ces 6 termes devra être calculé au préalable et stocké dans une liste.
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.
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))
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.