MODÈLE D’ISING EN DIMENSION 1 ET 2 - CORRIGÉ
On a rédigé les programmes demandés en pseudo-code, les variables des programmes
apparaissant en gras.
1. D’après la propriété (i), les spins des électrons ne peuvent être observés que dans
les positions +1 et −1 relativement à une direction ~u ; et la probabilité donnée par
l’équation (1) est bien une probabilité sur l’ensemble des fonctions σ : E → {±1}. Si
deux électrons voisins e et e0 ont le même spin dans une configuration σ, ceci contribue
à un facteur eβ dans la probabilité de σ, tandis que s’ils ont un spin différent, ceci
contribue à un facteur e−β dans la probabilité de σ. Comme e−β < eβ , la propriété
(ii) est bien retranscrite. Finalement, si un électron e pris individuellement a un spin
positif, ceci contribue à un facteur eα dans la probabilité de σ, tandis qu’un spin
négatif apporte une contribution e−α . Supposons α > 0. Alors, eα > e−α , donc un
électron a plus tendance à avoir un spin positif, ce qui correspond à la propriété (iii).
~ ui. On a ainsi montré
En particulier, le coefficient α doit avoir le même signe que hB|~
que les trois propriétés étaient retranscrites par l’équation (1).
Si α tend vers l’infini, on s’attend à ce que les configurations avec des spins négatifs
aient une probabilité beaucoup plus petite que les configurations avec des spins positifs
à la place des spins négatifs, et à la limite, la probabilité Pα,β = P+∞,β charge avec
poids 1 la configuration σ : E → {±1} avec seulement des spins positifs : σ(e) = +1
pour tout e ∈ E.
Si β tend vers l’infini, alors les configurations avec des électrons voisins de spins
différents deviennent très rares, et à la limite, la probabilité Pα,β = Pα,+∞ charge
uniquement les configurations σ+ et σ− constantes avec que des spins positifs ou que
des spins négatifs (ou pour être plus précis, les configurations constantes sur chaque
composante connexe de E). Finalement, si α = β = 0, alors Pα,β est la configuration
uniforme sur l’ensemble des configurations : P0,0 [σ] = 2−card E pour tout σ.
2. Les 4 configurations de spins sur [[1, 2]] ⊂ Z sont ++, +−, −+ et −−, et elles ont
probabilités respectives :
e2α+β e−β e−β e−2α+β
; ; ; .
Z(α, β) Z(α, β) Z(α, β) Z(α, β)
La constante de normalisation Z(α, β) est donc égale à e2α+β + e−2α+β + 2e−β .
Les 16 configurations de spins sur [[1, 2]]2 ⊂ Z2 et leurs poids Pα,β [σ] Z(α, β) sont
données par le tableau suivant :
+ + + + + + + −
+ + e4α+4β + − e2α − + e2α + + e2α
− + + + + − − −
+ + e2α − − 1 + − 1 + + 1
− + − + −4β + − −4β − − −2α
− + 1 + − e − + e − + e
− − −2α − − − −
+ − e −
+
− e−2α +
− − e−2α − − e−4α+4β
La fonction de partition vaut donc Z(α, β) = e4α+4β +e−4α+4β +2e−4β +4+4e2α +4e−2α .
1
3. Notons que P (x, y) = eβ si xy = 1, et P (x, y) = e−β si xy = −1 ; ainsi, P (x, y) = eβ xy .
Soit (σ(1), . . . , σ(n)) une configuration de spins de taille n ; calculons la probabilité de
cette configuration en tant que trajectoire de la chaîne de Markov décrite :
P[σ(1), . . . , σ(n)] = P[σ(1)]P (σ(1), σ(2)) P (σ(2), σ(3)) · · · P (σ(n − 1)σ(n))
1 1
= eβ(σ(1)σ(2)+σ(2)σ(3)+···+σ(n−1)σ(n))
2 (2 cosh β)n−1
en utilisant la remarque précédente. La probabilité prend bien la forme de la formule
(2) avec α = 0, et ceci démontre la proposition 1, avec pour valeur de la fonction de
partition Zn (0, β) = 2(2 cosh β)n−1 .
On en déduit un programme qui calcule en temps linéaire une configuration de taille
n sous la loi P0,β :
Configuration(n,beta):
res = liste vide
etat = tirage d’une variable de Bernoulli de parametre 1/2
adjoindre etat a res
pour i entre 2 et n:
mult = tirage d’une variable de Bernoulli de
parametre exp(beta)/(2 cosh(beta))
etat = etat*mult
adjoindre etat a res
retourner la liste res
(Par variable de Bernoulli de paramètre p, on entend ici une variable qui prend les
valeurs ±1 avec probabilité p pour la valeur +1.)
4. La matrice de transfert a la propriété suivante : T (x, y) = eαy+βxy pour tous signes x, y.
Le poids d’une configuration σ s’écrit donc :
Pn Pn−1
Zn (α, β) Pα,β [σ] = eα i=1 σ(i)+β i=1 σ(i)σ(i+1)
n
Y
= eασ(1) eασ(i+1)+βσ(i)σ(i+1)
i=1
= Vσ(1) Tσ(1)σ(2) · · · Tσ(n−1)σ(n) .
La somme sur toutes les configurations de ces quantités est
Zn (α, β) = (V T n−1 )+1 + (V T n−1 )−1 = V T n−1 W.
Les valeurs propres de T sont les racines de (X − eα+β )(X − e−α+β ) − e−2β = 0, il s’agit
donc de
q p
λ± = e cosh α ± e2β cosh2 α − 2 sinh 2β = eβ cosh α ± e2β sinh2 α + e−2β ,
β
en utilisant cosh2 α = sinh2 α + 1 pour la seconde égalité. Comme elles sont distinctes,
T est diagonisable et il existe une matrice inversible P telle que
n−1
λ+ 0
Zn (α, β) = V P 0 λ− P −1 W.
2
Ceci prouve que Zn (α, β) est une combinaison linéaire de (λ+ )n−1 et de (λ− )n−1 . Pour
identifier les coefficients de cette combinaison, on écrit
Z1 (α, β) = eα + e−α = a+ + a− ;
Z2 (α, β) = e2α+β + e−2α+β + 2e−β = a+ λ+ + a− λ− .
La résolution du système donne bien les valeurs de a+ et a− :
eβ sinh2 α + e−β
a± = cosh α ± p .
e2β sinh2 α + e−2β
La connaissance de la fonction de partition et donc des probabilités de chaque configu-
ration ne fournit pas d’algorithme efficace de construction de configuration aléatoire,
car il faudrait pour cela calculer a priori toutes les probabilités des 2n configurations,
ce qui est très long.
5. L’espérance sous Pα,β de etMn s’écrit :
X Pn
Eα,β [etMn ] = et i=1 σ(i) Pα,β [σ]
σ∈{±1}E
1 X Pn Pn Pn−1
= et i=1 σ(i) α
e i=1 σ(i)+β i=1 σ(i)σ(i+1)
Zn (α, β)
σ∈{±1}E
Zn (α + t, β)
= .
Zn (α, β)
Cette fonction de t est une combinaison linéaire finie de fonctions t 7→ ect , elle est donc
clairement infiniment dérivable. Sa dérivée s’écrit :
∂Eα,β [etMn ] X ∂ etM (σ) X
= Pα,β [σ] = Pα,β [σ] etM (σ) Mn (σ)
∂t σ
∂t σ
où la somme porte sur les configurations de taille n. En particulier, pour t = 0,
∂Eα,β [etMn ] X
= Pα,β [σ] M (σ) = Eα,β [Mn ].
∂t t=0 σ
Le terme de gauche de cette expression est aussi
1 ∂Zn (α, β) ∂(log Zn )(α, β)
= .
Zn (α, β) ∂α ∂α
D’après la question précédente,
n−1
a− λ−
log Zn (α, β) = (n − 1) log λ+ + log a+ + log 1 + n−1
a+ λ+
n−1
a− λ−
1 ∂λ+ (α, β) 1 ∂a+ (α, β) ∂ log 1 +
(α, β) n−1
a+ λ+
Eα,β [Mn ] = (n − 1) + + .
λ+ ∂α a+ ∂α ∂α
Le terme log 1 + a− λn−1 n−1
− /a+ λ+ tend exponentiellement vers 0 lorsque n tend vers
l’infini, et pour (α, β) restant dans une partie compacte, on voit facilement que c’est
aussi le cas de ses dérivées. De même, le terme a1+ ∂a+∂α(α,β)
reste un O(1) lorsque (α, β)
3
reste dans une partie compacte. Par conséquent,
2β
eβ sinh α + √e 2βcosh α2 sinh−2β
α
e sinh α+e
Eα,β [Mn ] = (n − 1) p + O(1)
eβ cosh α + e2β sinh2 α + e−2β
eβ sinh α
=np + O(1)
e2β sinh2 α + e−2β
Lorsque β = 0, c’est-à-dire que les interactions entre les spins sont nulles, on obtient
α −α
une magnétisation moyenne par électron de cosh sinh α
α
= eeα −e +e−α
, ce qui est logique pour
eα
des variables de Bernoulli ±1 de paramètre eα +e−α . Lorsque α = 0, on obtient une
magnétisation moyenne nulle, ce qui découle aussi dans cette situation de la symétrie
du système.
6. Pour illustrer le théorème central limite, on peut au minimum dresser un histogramme
des magnétisations sur un grand nombre N de configurations aléatoires :
Histogramme(N,n,alpha,beta):
res = liste vide
m = exp(beta)*sinh(alpha)/
sqrt(exp(2*beta)*(sinh(alpha))^2 + exp(-2*beta))
var = exp(-beta)*cosh(alpha)/
(exp(2*beta)*(sinh(alpha))^2 + exp(-2*beta))^(3/2)
pour i entre 1 et N:
conf = Configuration(n,alpha,beta)
magnet = somme(conf)
adjoindre (magnet - n*m)/sqrt(n*var) a res
retourner un histogramme des valeurs de res
Un autre programme utile est un estimateur de la variance de la magnétisation, ce qui
permet de tester la formule du théorème 4. On s’attend à une variance proche de la
valeur
n e−β cosh α
.
(e2β sinh2 α + e−2β )3/2
On peut utiliser l’estimateur sans biais
N
1 X (i)
(M − M )2 ,
N − 1 i=1
(i)
où les M sont les magnétisations de N tirages indépendants de configurations, et
1 N
M = n i=1 M (i) est leur moyenne empirique. Ainsi :
P
EstimateurVariance(N,n,alpha,beta):
res = liste vide
var = exp(-beta)*cosh(alpha)/
(exp(2*beta)*(sinh(alpha))^2 + exp(-2*beta))^(3/2)
pour i entre 1 et N:
conf = Configuration(n,alpha,beta)
adjoindre somme(conf) a res
moyenne = somme(res)/N
estimateur = 0
pour i entre 1 et N:
4
ajouter (res[i]-moyenne)^2/(N-1) a estimateur
afficher estimateur/n et var
permet de comparer la valeur théorique et une valeur empirique de la variance.
z Mn√−nm
Pour la preuve du théorème 4, détaillons le calcul de log Eα,β [e n ], avec z variable
complexe. Pour éviter les singularités du logarithme complexe, il est préférable d’évaluer
z Mn√−nm
directement Eα,β [e n ], qui est une fonction bien définie et entière en tant que
combinaison linéaire de fonctions z 7→ ecz :
!
√ β
h Mn −nm i
z √n z e sinh α
Eα,β e = Fn √ , α, β exp −z n p
n e2β sinh2 α + e−2β
!
√ β
z e sinh α
= Zn √ + α, β Zn (α, β)−1 exp −z n p .
n e2β sinh2 α + e−2β
Les paramètres α et β étant fixés, la fonction de partition Zn (α, β) a pour développe-
ment asymptotique
Zn (α, β) = a+ (α, β)(λ+ (α, β))n−1 (1 + o(1)).
Si l’on remplace α par α + √zn avec z variable complexe fixée, les expressions de
a+ , λ+ , a− , λ− avec des racines carrées restent bien définies pour n assez grand (on
reste proche de l’axe réel R∗+ , autour duquel la fonction racine se prolonge de façon
holomorphe). L’estimée précédente reste donc vraie avec un o(1) uniforme pour z dans
une partie compacte de C :
n−1
z z z
Zn α + √ , β = a+ α + √ , β λ+ α + √ , β (1 + o(1))
n n n
n
z
= a+ (α, β) λ+ α + √ , β (1 + o(1))
n
Par conséquent,
!n
√z , β)
!
h
z Mn√−nm
i λ+ (α + n √ eβ sinh α
Eα,β e n = exp −z n p (1 + o(1))
λ+ (α, β) e2β sinh2 α + e−2β
!!n
λ+ (α + √z , β)
n z eβ sinh α
= exp − √ p (1 + o(1))
λ+ (α, β) n e2β sinh2 α + e−2β
avec un reste uniforme sur les parties compactes de C. On écrit alors le développement
de Taylor à l’ordre 2 en √zn du terme qui est élevé à la puissance n. On obtient après
calcul
2 n
h Mn −nm i
z √n e−β cosh α z2 z
Eα,β e = 1 + 2β 2 −2β
+o (1 + o(1))
(e sinh α + e ) 2n 3/2 n
e−β cosh α z2
= exp (1 + o(1)),
(e2β sinh2 α + e−2β )3/2 2
cette expression asymptotique étant toujours valable pour z variable complexe, n as-
sez grand et avec un reste uniforme sur les parties compactes du plan complexe. En
particulier, on peut poser z = it, et on reconnaît alors la fonction caractéristique d’une
5
variable gaussienne. Le critère de Lévy permet de conclure quand à la convergence en
loi de Mn√−nm
n
vers un variable gaussienne.
7. La loi de σ (k+1) conditionnellement à (σ (0) , . . . , σ (k) ) est
P[σ (k+1) = τ |σ (0) , . . . , σ (k) ] = P[Fek ,uk (σ (k) ) = τ |σ (0) , . . . , σ (k) ]
1 X
= p (k) 1(τ (e)=+1, τ (e0 )=σ(k) (e0 )) + (1 − pe,σ(k) ) 1(τ (e)=−1, τ (e0 )=σ(k) (e0 ))
card E e∈E e,σ
en utilisant l’indépendance de (ek , uk ) par rapport à la tribu engendrée par les confi-
gurations σ (0) , . . . , σ (k) , et où
(k) (k)
eα+β(P (e,σ )−N (e,σ ))
pe,σ(k) = .
2 cosh(α + β(P (e, σ (k) ) − N (e, σ (k) )))
On obtient bien une fonction P (σ (k) , τ ) qui ne dépend que des états σ (k) et τ , donc une
chaîne de Markov. Lorsque E = [[1, 2]], le graphe des transitions de cette chaîne est :
e−α−β
eα+β 4 cosh(α+β) e−α−β eα−β
2 cosh(α+β) 4 cosh(α+β) + 4 cosh(α−β)
++ −+
α+β
e
4 cosh(α+β)
e−α−β eα+β e−α+β eα−β
4 cosh(α+β) 4 cosh(α+β) 4 cosh(α−β) 4 cosh(α−β)
e−α+β
4 cosh(α−β)
+− −−
e−α−β eα−β e−α+β
4 cosh(α+β) + 4 cosh(α−β) eα−β
2 cosh(α−β)
4 cosh(α−β)
Il est clair sur le graphe que la chaîne est irréductible et apériodique (il y a à chaque
fois une probabilité non nulle de rester sur la même configuration) ; et ceci reste vrai
pour tout ensemble d’états E ⊂ Zd . La réversibilité de la mesure d’Ising Pα,β pour
cette chaîne de Markov découle des arguments du texte : avec σ et σ 0 voisins dans
le graphe de la chaîne de Markov, c’est-à-dire ne différant qu’en un site e, on sépare
dans Pα,β [σ] P (σ, σ 0 ) les termes dépendant du site e et les autres termes, qui seront les
mêmes dans Pα,β [σ 0 ] P (σ 0 , σ) car σ et σ 0 ne différent qu’au niveau du site e. Supposant
par exemple σ(e) = +1 et σ 0 (e) = −1, les termes dépendant de e dans Pα,β sont
σ(e)σ(e0 )
P
eασ(e)+β e0 ∼e = eα+β(P (e,σ)−N (e,σ)) ,
e−α−β(P (e,σ)−N (e,σ))
et P (σ, σ 0 ) = 2 cosh(α+β(P (e,σ)−N (e,σ)))
. Ainsi,
1
Pα,β [σ] P (σ, σ 0 ) = × termes ne dépendant pas de e,
2 cosh(α + β(P (e, σ) − N (e, σ)))
6
et symétriquement,
1
Pα,β [σ 0 ] P (σ 0 , σ) = × termes ne dépendant pas de e.
2 cosh(α + β(P (e, σ 0 ) − N (e, σ 0 )))
Comme P (e, σ) = P (e, σ 0 ) et N (e, σ) = N (e, σ 0 ), la réversibilité est bien démontrée.
Par suite, la mesure d’Ising est invariante par la chaîne de Markov, car
!
X X X
(P × P )(σ) = P[σ 0 ] P (σ 0 , σ) = P[σ] P (σ, σ 0 ) = P[σ] P (σ, σ 0 ) = P[σ].
σ0 σ0 σ0
On a donc une chaîne de Markov irréductible apériodique qui admet une mesure de
probabilité invariante, donc qui est récurrente positive. Les deux théorèmes ergodiques
assurent alors que :
(a) La loi de la configuration σ (n) de la chaîne de Markov tend vers la loi d’Ising :
∀σ ∈ {±1}E , lim P[σ (n) = σ] = Pα,β [σ].
n→∞
(b) Les lois empiriques de la chaîne de Markov convergent également vers la loi d’Ising,
au sens presque sûr :
n
1X
E
∀σ ∈ {±1} , lim 1σ(n) =σ = Pα,β [σ].
n→∞ n
i=1
8. D’après la question précédente, pour simuler une configuration sous la loi d’Ising, on
peut utiliser la chaîne de Markov et regarder un état σ (k) avec k grand. En pratique,
sur un domaine carré de taille n × n, on s’attend à ce que la configuration aient une
probabilité proche de la loi d’Ising après que les sites ont tous été tirés par les variables
aléatoires (e1 , e2 , . . . , ek ) ; c’est le cas si k & n2 log n2 . D’autre part, on peut prendre
pour configuration initiale n’importe quelle configuration ; un choix possible est de
partir d’une configuration aléatoire où chaque site suit une loi de Bernoulli indépen-
dante des autres et de paramètre 1/2. On a séparé le programme en les composantes
suivantes :
% creation de la configuration initiale
ConfigurationInitiale(n):
M = matrice de taille n*n
pour i de 1 a n:
pour j de 1 a n:
M[i,j] = tirage d’une variable de Bernoulli de parametre 1/2
retourner M
% calcul du nombre de voisins positifs et negatifs d’un site
% dans une configuration (attention aux effets de bords !)
VoisinPositif(M,e):
i,j = coordonnees de e
n = taille de M
res = 0
si i > 1 et M[i-1,j] = 1:
augmenter res de 1
si i < n et M[i+1,j] = 1:
augmenter res de 1
7
si j > 1 et M[i,j-1] = 1:
augmenter res de 1
si j < n and M[i,j+1] = 1:
augmenter res de 1
retourner res
VoisinNegatif(M,e):
retourner VoisinPositif(-M,e)
% transition de la chaine de Markov associee a un site e et a une quantite u
% (on effectue les transitions sur la matrice en place)
Transition(M,e,u,alpha,beta):
i,j = coordonnees de e
P = VoisinPositif(M,e)
N = VoisinNegatif(M,e)
t = alpha + beta*(P-N)
si u < exp(t)/(2*cosh(t)):
M[i,j] = 1
sinon:
M[i,j] = -1
%Configuration aleatoire apres k etapes de la chaine de Markov
ConfigurationMarkov(k,n,alpha,beta):
res = ConfigurationInitiale(n)
pour i entre 1 et k:
x = tirage d’un entier uniformement entre 1 et n
y = tirage d’un entier uniformement entre 1 et n
e = (x,y)
u = tirage d’une variable uniforme sur [0,1]
Transition(res,e,u,alpha,beta)
retourner res
%Configuration aleatoire sous une loi tres proche de la loi d’Ising
Ising(n,alpha,beta):
retourner ConfigurationMarkov(3*(n^2)*log(n),n,alpha,beta)