8.3. 【例題】イジング模型の平均場近似¶

8.3.1. 方程式 ¶

イジング模型は強磁性相転移を記述する最も簡単な模型です。ハミルトニアンは

$\mathcal{H} = -J \sum_{\langle i, j \rangle} \sigma_i \sigma_j$

で与えられます。ここで、 $\sigma_i$ は格子点 $i$ にあるスピンの状態を表し、+1または-1を取ります（+1が上向き、-1が下向き）。 $\langle i, j \rangle$ は隣り合う格子点の組み合わせに関する和を意味します。

この模型は多体模型なので、分配関数を厳密に評価することは、特別な場合を除いてできません。そこで、統計力学の教科書では、平均場近似を適用して一体模型に置き換えて物理量を計算します。

平均場近似では、平均磁化 $m = \langle \sigma_i \rangle$ は以下の自己無撞着方程式を解くことで決定できます：

(1)¶ $m = \tanh (\beta z J m)$

ここで、 $\beta=k_\mathrm{B}T$ は逆温度、 $z$ は最近接格子点の数を表します。この方程式は $T<T_\mathrm{c} \equiv zJ$ を満たす低温領域で非ゼロの解を持ちます（つまり自発磁化が生じる）。具体的には、 $T = 0$ で $m= \pm1$ 、 $T \simeq T_\mathrm{c}$ で $m=\pm\sqrt{3(T_\mathrm{c}-T)}$ になることが解析的に示せます。しかし、中間の温度領域では解析的に解けません。数値計算を使ってこの方程式を解き、自発磁化の温度変化を求めます。

磁化 $m$ が求まったら、エントロピーと比熱はそれぞれ

$\frac{S}{Nk_\mathrm{B}} = -\frac{1+m}{2} \ln \frac{1+m}{2} -\frac{1-m}{2} \ln \frac{1-m}{2}$

$\frac{C}{Nk_\mathrm{B}} = \frac{1-m^2}{4[1-(1-m^2)/t]} \left( \ln \frac{1+m}{1-m} \right)^2$

から計算できます。ここで、 $t \equiv T/T_\mathrm{c}$ です。

8.3.2. 実装 ¶

ソースコード ising_mf.py

方程式を解くにはSciPyの scipy.optimize.root_scalar 関数を使います。 scipy.optimize.root_scalarの使い方で説明したように、まず解きたい関数 f を定義します。 (1) 式を $f(m)=0$ の形に変形して $f(m)$ を返すようにします。

def func(m, zj, beta, h):
    return m - np.tanh(beta * (zj * m + h))

第二引数以下はパラメーターです。パラメーターの数値はroot_find関数の args 引数に与えます。

方程式の解を求める関数は以下のように書けます。

def find_solution(zj, t, h, m_min=1e-6, m_max=1, verbose=False):
    beta = np.inf if t==0 else 1./t
    sol = optimize.root_scalar(func, args=(zj, beta, h), method='brentq', bracket=(m_min, m_max))
    if verbose:
        print(sol)
    return sol.root

温度 t から逆温度 beta=1/t の変換では、t=0 の場合に beta=numpy.inf をセットすることで、t=0 も扱えるようにしています。非線型方程式の解法として、ここでは``method='brentq'`` によりBrent法を指定しています。 Brent法は挟み込むアルゴリズムですので、解の範囲を bracket 引数に与えます。 $m\neq 0$ の解を選ぶような値をデフォルトとしていますが、外部からも変更できるようにしています。磁場 h がある場合には調整する必要があります。 verbose は詳細な結果を出力するオプションです。計算結果がおかしい場合などに verbose=True として詳細な出力を得るようにします。

あとはこの関数を呼び出せば自己無撞着方程式を満たす解が得られます。

def main():
    t = 0.1
    zj = 1.0
    h = 0.0
    m = find_solution(zj, t, h, verbose=True)
    s = calc_entropy(m)
    c = calc_specific_heat(m, zj, t)
    print(f"zJ = {zj}")
    print(f"T  = {t}")
    print(f"h  = {h}")
    print(f"m  = {m}")
    print(f"S  = {s}")
    print(f"C  = {c}")

比熱とエントロピーも計算していますが、特別な計算はしていないので、解説は省きます。ソースコードを見てください。

8.3.3. 出力結果 ¶

$ python ising_mf.py
      converged: True
           flag: 'converged'
 function_calls: 4
     iterations: 3
           root: 0.9999999958776924
zJ =  1.0
T  = 0.1
h  = 0.0
m  = 0.999999995878
S  = 4.32842296841e-08
C  = 8.2446158185e-07

最初の5行は verbose=True としたことによる、詳細の出力です。 m はほぼ1で完全に強磁性状態になっていることが分かります。

パラメーターを変化させて正しい解が得られていることを確認してください。

8.3.4. 温度変化 ¶

さて方程式を正しく解けていることが確認出来たら、温度変化のグラフを作ってみましょう。そのために、先ほど作った関数 find_solution は使いまわして、main関数だけを新たに作ります。

ソースコード ising_mf_t.py

t_dep.pyの冒頭で先ほどのファイルをインポートします。

from ising_mf import find_solution, calc_entropy, calc_specific_heat

これにより、ising_my.pyの中にある find_solution 関数を使えるようになります。ising_mf.pyとt_dep.pyは必ず同じディレクトリにおいてください。そうしないとising_mf.pyを見つけられずにエラーがでます。

あとは温度を変化させながら find_solution を繰り返し呼び出せば完了です。

def main():
    h = 0
    tmin = 0
    tmax = 1
    n = 201
    for t in np.linspace(tmin, tmax, n):
        try:
            m = find_solution(1., t, h)
            s = calc_entropy(m)
            c = calc_specific_heat(m, 1., t)
        except ValueError:
            print(f"# {t:.4f}")
        else:
            print(f"{t:.4f} {m:.5e} {s:.5e} {c:.5e}")

解が見つからなかった場合のために、例外処理をしてあります。 tryブロックでエラー（例外）が発生したら、該当する例外に対応するexceptブロックが実行されます。もし、該当するexceptブロックが存在しない場合には、プログラムが終了します。エラーがなくtryブロックが終了した場合には、elseブロックが実行されます。例外処理については公式ドキュメントのほか、文献 [2] を参照してください。