イジング模型

と書くことができます。<i,j>は最近接のスピンの和だけを取るということですが、少し分かりづらいと思うのでプログラムの一部を実際に書き下してみましょう。

ここでは見やすさのため Pyhton を用いることにします。スピンが存在するサイトの数は N×M で、相互作用の結合定数を J>0 としましょう。

さてエネルギーは、

E = 0

for i in range(N):

for j in range(M):

E = E -J (s[i, j] s[(i+1)%N, j] + s[i, j] s[i, (j+1)%M])   ← エネルギー計算

と書けます。 つまり、あるサイトにあるスピン s[i,j]  は（x,y方向それぞれの方向の）そのすぐ隣のスピンとのみ相互作用するということです。

%は余りを計算する演算で、今の場合は x(y) 方向に N(M) 個のサイトがあるので、端のスピンを考える場合は逆側の端のスピンと相互作用するよう導入しています。

s[N,j] s[N+1,j] を例にすると、s[N+1,j] は存在しないサイトなので%の演算を使って s[0,j] として、s[N,j] s[0,j] にするということです。

(Python の配列は通常 0 から番号が振られます)

 これを周期境界条件といい、感覚的には x 軸や y 軸として直線でなく円になっている形を思い浮かべればO.K.です。

さて、上のエネルギーの表式を見てみるとスピンが揃っているとエネルギーが下がります。実際、-J×1×1 < 0 ですし、-J×(-1)×(-1) <0 ですね。

普通、世界はエネルギーが低い状態を好みますので（人類の経験則です）、この系はスピンが揃う状況を好みます。

スピンが揃った状況は明らかに系として向きがありますので強磁性の状態（つまり磁石の状態）と解釈できます。

ここから、この模型は磁性の議論ができる模型になっていることが分かると思います。

実は上の議論は、ゼロ温度の場合に正しい記述です。

温度がある場合は揺らぎがあり、配位として（統計力学を学んでいないと理解しづらいですが）エネルギーが一番低いもの以外も一定の確率で存在します。

つまり、様々なエネルギー状態の配位が確率分布に従って生成され、そのミクロな情報の平均としてマクロな状態が記述されるわけです。

結果として実現する状態はゼロ温度とは異なってくることがあります。これをプログラムしてみましょう。温度を T で導入します。

通常、数値計算ではメトロポリス法が便利で、N_step はマルコフ過程を作る際の step 数とします。具体的には、

s_new = s   ← スピンの更新のため新しい配列を用意

for i in range(N_step)

j = int((N-1)*random.random())    ← 乱数を使ってサイトを一つ選ぶ

k = int((M-1)*random.random())

s_new = s.copy()

s_new[j,k] = - s[j,k]   ← 選んだサイトのスピンを反転 

E_new = 0   ← 一つスピンを反転させた時のエネルギー計算

for i in range(N):

for j in range(M):

E_new = E_new -J (s_new[i, j] s_new[(i+1)%N, j] + s_new[i, j] s_new[i, (j+1)%M])

if E_new < E:   ← スピン反転後の方がエネルギーが下がった場合、そのスピンを受け入れる（sを更新する）

s = s_new.copy()

elif E_new > E:   ← スピン反転後の方がエネルギーが上がった場合

prob1 = math.e**(-(E_new-E)/T)

prob2 = random.random()

if prob1 < prob2:   ← 確率でスピンを受け入れる（sを更新する）

s = s_new

となります。最後の if 文のところが温度揺らぎを考慮する部分で、exp(-(E_new-E)/T) は Gibbs 分布と呼ばれます。

この過程を繰り返していくことでもっともらしいスピン配位が得られ、統計和を取ることで磁化などが計算できるようになります。

（この過程は本物の時間変化では無いことに注意）

少し式を見て構造を考えて見ましょう。

まず、温度が高くなると確率 exp(-(E_new-E)/T) が1に近づくので、スピンの反転が受け入れられやすくなります。

つまり、高温ではスピンはランダムな方向を向き、系としての方向はなくなります（強磁性でなくなる）。

一方、温度が低いと確率 exp(-(E_new-E)/T) が0に近づくのでエネルギーが一番低い配位が支配的に成るためスピンは揃います。

つまり、ある温度において強磁性から常時性への相転移が存在しうるわけです。

（実際に相転移するかや、相転移の次数などは系の空間の次元などに依存しますが、ここでは省略します）

本当はマルコフ鎖等の議論や外場の導入なども確認する必要がありますが、とりあえずは上記の記載を理解できるよう勉強してみてください。

2次元イジング模型のプログラムを組んだのち、実際の卒研テーマに入っていきましょう。

実際に上記の内容で一次元イジング模型を計算すると以下のような図が得られます。

横軸は温度、縦軸はスピンの期待値です（N=20の計算）。