MC-06 量子 Wang-Landau

このチュートリアルの目的は、量子 Wang-Landau（QWL）アルゴリズムを用いて量子スピン系の熱力学的性質を計算することです。固定温度でシミュレートする経路積分 QMC 法とは異なり、QWL は分配関数の展開次数の空間でランダムウォークを行い、状態密度を直接推定します。そのため、一度のシミュレーションから、再実行することなく全温度領域にわたって熱力学量——エネルギー、比熱、エントロピー、帯磁率——を評価できます。

量子ハイゼンベルクスピン鎖の熱力学

強磁性ハイゼンベルク鎖

パラメータファイル parm6a は、40 サイトの鎖上の $S=\frac{1}{2}$ ハイゼンベルク強磁性体の QWL シミュレーションを設定します。

LATTICE="chain lattice"
MODEL="spin"
local_S = 1/2
L       = 40
CUTOFF  = 500
{J = -1}

CUTOFF はアルゴリズムがサンプリングする分配関数の展開次数の上限を設定します。ここで扱う温度範囲の 40 サイト鎖には 500 で十分です。

コマンドライン

入力ファイルを準備し、qwl コードを実行したあと：

parameter2xml parm6a
qwl parm6a.in.xml

qwl_evaluate を使って、すべての熱力学量について XML プロットファイルを生成します。

qwl_evaluate --T_MIN 0.1 --T_MAX 10 --DELTA_T 0.1 parm6a.task1.out.xml

これにより、物理量ごとに 1 つの XML ファイルが生成されます。

parm6a.task1.plot.energy.xml
parm6a.task1.plot.free_energy.xml
parm6a.task1.plot.entropy.xml
parm6a.task1.plot.specific_heat.xml
parm6a.task1.plot.uniform_structure_factor.xml
parm6a.task1.plot.staggered_structure_factor.xml
parm6a.task1.plot.uniform_susceptibility.xml

plot2text でデータをプレーンテキストとして取り出します。

plot2text parm6a.task1.plot.energy.xml

ツール plot2xmgr（Grace）と plot2gp（Gnuplot）は、それぞれの形式で同等の出力を生成します。解析には次に述べる Python が推奨されます。

Python

スクリプト tutorial6a.py はシミュレーションを設定・実行し、その後 1 回の呼び出しですべての物理量を評価します。

import pyalps
import matplotlib.pyplot as plt

data = pyalps.evaluateQWL(pyalps.getResultFiles(prefix='parm6a'),
                          DELTA_T=0.1, T_MIN=0.1, T_MAX=10.0)

for s in pyalps.flatten(data):
    plt.figure()
    plt.title("Ferromagnetic Heisenberg chain")
    pyalps.plot.plot(s)
    plt.show()

強磁性体（ $J=-1$ ）では、低温に幅の広い比熱のピーク（ショットキー異常に似たもの）が現れ、一様帯磁率は $T\to 0$ で発散するはずです。これは 1 次元において強磁性秩序が絶対零度で発達することと整合します。

反強磁性ハイゼンベルク鎖

反強磁性鎖をシミュレートするには、 $J=-1$ の代わりに $J=1$ とします。パラメータは parm6b に、Python スクリプトは tutorial6b.py にあります。それ以外は強磁性の場合とまったく同じです。

反強磁性体では、一様帯磁率は $T\to 0$ でも有限にとどまります（これは 1 次元反強磁性体のスピン液体基底状態の特徴です）。一方、比熱のピークは位置がずれ、広がり方も異なります。

3 次元ハイゼンベルク反強磁性体

3 次元量子ハイゼンベルク反強磁性体のシミュレーション

パラメータファイル parm6c は、 $4^3=64$ サイトの単純立方格子上の $S=\frac{1}{2}$ ハイゼンベルク反強磁性体の QWL シミュレーションを設定します。 Python スクリプトは tutorial6c.py です。実行と評価は上の鎖の場合とまったく同じです。

スタッガード構造因子 $S(\pi,\pi,\pi)$ は $T\approx 1$ の近くで急激に立ち上がり、反強磁性相関の発現を示すはずです。比熱にはそれに対応するピークが現れ、同じ温度領域でエントロピーは急速に減少します。

臨界点を決めるための有限サイズスケーリング解析

有限サイズスケーリングによれば、臨界点におけるスタッガード構造因子は $S(L) \propto L^{2-\eta}$ でスケールし、 $\eta \approx 0.034$ （3 次元古典ハイゼンベルク普遍性クラス）です。 $S(L)/L^{2-\eta}$ を温度に対してプロットすると、異なる $L$ の曲線が臨界温度 $T_c$ で交わるはずです。

パラメータファイル parm6d（または tutorial6d.py）は、低温での精度を保つために大きめの CUTOFF=1000 を用いて、2 つのシステムサイズ（ $L=4$ と $L=6$ ）で立方格子反強磁性体を実行します。実行後、結果を読み込みます。

import pyalps
import matplotlib.pyplot as plt
import copy

results = pyalps.evaluateQWL(pyalps.getResultFiles(prefix='parm6d'),
                             DELTA_T=0.05, T_MIN=0.5, T_MAX=1.5)

各システムサイズについてスタッガード構造因子を取り出し、 $L^{-(2-\eta)}$ で再スケールして $L$ でラベルを付けます。

eta = 0.034
data = []
for s in pyalps.flatten(results):
    if s.props['ylabel'] == 'Staggered Structure Factor per Site':
        d = copy.deepcopy(s)
        l = s.props['L']
        d.props['label'] = 'L=' + str(l)
        d.y = d.y * pow(float(l), -(2 - eta))  # L^{-(2-eta)} で再スケール
        data.append(d)

続いてスケーリング曲線をプロットします。

plt.figure()
plt.title("Scaling plot for cubic lattice Heisenberg antiferromagnet")
pyalps.plot.plot(data)
plt.legend()
plt.xlabel('Temperature $T/J$')
plt.ylabel('$S(\pi,\pi,\pi)\, L^{-(2-\eta)}$')
plt.show()

異なる $L$ の曲線は $T_c \approx 0.946$ の近くで交わるはずです。

問題

強磁性（ $J=-1$ ）と反強磁性（ $J=1$ ）の鎖について、比熱と一様帯磁率を比較してください。違いが最も顕著なのはどこでしょうか。また、なぜ高温では両者がほとんど同一になるのでしょうか。
鎖のエントロピーは $T=0$ および $T\to\infty$ でいくらになりますか。（ $L=40$ のデータが低温で乱れている場合は $L=8$ でやり直してください。）これは熱力学第三法則と整合しますか。
なぜ一様帯磁率は強磁性体では $T\to 0$ で発散するのに、反強磁性体では有限にとどまるのでしょうか。
なぜ 3 次元反強磁性体のスタッガード構造因子は $T\approx 1$ の近くで増大し始めるのでしょうか。これに伴う他の熱力学的な兆候は何でしょうか。
再スケールした構造因子の曲線 $S(L)\,L^{-(2-\eta)}$ は 1 つの温度で交わりますか。 $T_c$ の見積もりはいくらでしょうか。文献値 $T_c = 0.946$ と比較してください。
$T_c$ をより精密に見積もるにはどうすればよいでしょうか。シミュレーションで何を変える必要がありますか。
立方格子上の量子ハイゼンベルク強磁性体の臨界温度は、反強磁性体と同じになると思いますか。どのように確かめますか。

MC-05 ボゾン MC-07 相転移