MC-08 量子相転移

このチュートリアルの目的は、量子モンテカルロデータの有限サイズスケーリングを用いて、二次元スピン模型における量子相転移を特定し、その性質を明らかにすることです。古典相転移とは異なり、量子相転移は $T=0$ で起こり、熱ゆらぎではなく量子ゆらぎによって駆動されます。ここではビンダーキュムラントとスピン剛性を用いて臨界結合定数を決定し、続いて相関長指数 $\nu$ 、異常次元 $\eta$ 、動的臨界指数 $z$ を求めて、この転移がどの普遍性クラスに属するかを明らかにします。この模型には量子臨界点が二つ存在することが分かるので、その両方を調べます。

模型は、正方格子上に配置されたスピン $\frac{1}{2}$ ハイゼンベルグ梯子から成ります。各梯子は水平方向に伸びており、梯子内の格子点はレッグに沿って結合 $J_0$ で、ラングを介して結合 $J_1$ でつながっています。隣り合う梯子は $J_2$ によって垂直方向に結合しています。

o-J0-o-J0-o-J0-o   <- レッグ
|         |
J1        J1        <- ラング（梯子内）
|         |
o-J0-o-J0-o-J0-o   <- レッグ
|         |
J2        J2        <- 梯子間結合
|         |
o-J0-o-J0-o-J0-o   <- レッグ
|         |
J1        J1        <- ラング（梯子内）
|         |
o-J0-o-J0-o-J0-o   <- レッグ

パラメータ L は各レッグに沿った格子点数を指定し、W = L/2 は行数を指定するので、梯子は全部で W/2 本になります。このチュートリアルでは $J_0=J_1=1$ とし、梯子間結合 $J_2$ を変化させます（模型の図については Wenzel and Janke, Phys. Rev. B 79, 014410 (2009) の Fig. 1 も参照してください）。二次元ハイゼンベルグ模型では有限温度で相転移は起こりませんが（マーミン・ワグナーの定理）、 $T=0$ では異なる基底状態の間の転移が起こり得ます。

バックグラウンドでの計算

複数のシステムサイズにわたる有限サイズスケーリングの計算は、以下で行う単一システムのスキャンよりもはるかに時間がかかります。チュートリアルの残りの部分を進めている間にバックグラウンドで実行されるよう、今すぐ開始しておきましょう。

コマンドライン

parm8b をダウンロードして、次を実行します。

parameter2xml parm8b
loop parm8b.in.xml &

Python

先へ進む前に、別のターミナルまたはバックグラウンドプロセスとして tutorial8b.py の前半部分（セットアップと pyalps.runApplication の呼び出し）を実行してください。

異なる相を同定する

まず、二つの単純な極限、すなわち非結合の梯子（ $J_2=0$ ）と等方的な正方格子（ $J_2=1$ ）を考えます。非結合の梯子は短距離相関をもつ基底状態を示し、有限のスピンギャップをもちます。これがスピン液体相です。一方、正方格子は有限のスタッガード磁化を伴う長距離秩序を示します。これが反強磁性ネール相です。

この二つの相を調べる分かりやすい方法は、磁化率 $\chi$ を見ることです。両方の場合について $8\times 8$ のシステムをさまざまな温度でシミュレーションし、結果を比較しましょう。非結合の梯子では、スピンギャップのために低温で帯磁率が活性化型の振る舞いを示します。一方、正方格子では $T\to 0$ で有限の定数に近づきます。なお、有限系では十分低い温度になると有限サイズギャップのために $\chi$ は最終的にゼロに向かいますが、ここではそれは関心の対象ではありません。パラメータファイル parm8a を使ってコマンドラインからシミュレーションを実行します。

parameter2xml parm8a
loop parm8a.in.xml

あるいは Python スクリプト tutorial8a.py を使います。

import pyalps
import matplotlib.pyplot as plt
import pyalps.plot
import numpy as np
import pyalps.fit_wrapper as fw

parms = []
for j2 in [0.,1.]:
    for t in [0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0]:
        parms.append(
            { 
              'LATTICE'        : "coupled ladders", 
              'LATTICE_LIBRARY': 'lattices.xml',  # カスタム格子定義へのパス
              'MODEL_LIBRARY'  : 'models.xml',    # カスタムモデル定義へのパス
              'local_S'        : 0.5,
              'ALGORITHM'      : 'loop',
              'SEED'           : 0,
              'T'              : t,
              'J0'             : 1 ,
              'J1'             : 1,
              'J2'             : j2,
              'THERMALIZATION' : 5000,
              'SWEEPS'         : 50000, 
              'MODEL'          : "spin",
              'L'              : 8,
              'W'              : 4
            }
    )
    
input_file = pyalps.writeInputFiles('parm8a',parms)
pyalps.runApplication('loop',input_file)

$J_2=0$ の場合、スピンギャップの値は磁化率の有限温度での振る舞いから見積もることができます（Phys. Rev. B 50, 13515 (1994) で導出されています）。

\chi = \frac{A}{\sqrt{T}} \exp\!\left(-\frac{\Delta}{T}\right),

ここで $A$ とスピンギャップ $\Delta$ はフィッティングパラメータです。 $T\leq1$ のデータをフィットして、スピンギャップの推定値を求めましょう。Python でこの解析を行う例を以下に示します。

data = pyalps.loadMeasurements(pyalps.getResultFiles(pattern='parm8a.task*.out.h5'),['Staggered Susceptibility','Susceptibility'])
susc1=pyalps.collectXY(data,x='T',y='Susceptibility', foreach=['J2'])

lines = []
gap_j0 = None
for data in susc1:
    pars = [fw.Parameter(1), fw.Parameter(1)]
    data.y= data.y[data.x < 1]
    data.x= data.x[data.x < 1]
    f = lambda self, x, pars: (pars[0]()/np.sqrt(x))*np.exp(-pars[1]()/x)
    fw.fit(None, f, pars, np.array([v.mean for v in data.y]), data.x)
    prefactor = pars[0].get()
    gap = pars[1].get()
    if data.props['J2'] == 0.0:
        gap_j0 = gap
    
    lines += plt.plot(data.x, f(None, data.x, pars))
    lines[-1].set_label('$J_2=%.4s$: $\\chi = \\frac{%.4s}{\\sqrt{T}}\\exp\\left(\\frac{-%.4s}{T}\\right)$' % (data.props['J2'], prefactor,gap))

フィットによってギャップ $\Delta$ と前置因子 $A$ が得られます。得られた曲線を以下に示します。

plt.figure()
pyalps.plot.plot(susc1)
plt.xlabel(r'$T$')
plt.ylabel(r'$\chi$')
plt.title('spin gap $\\Delta \\approx$ %.4s' % gap_j0)
plt.legend()
plt.show()

$J_2=0$ では帯磁率は低温で急激に減少し、フィット曲線はデータによく追随して、スピンギャップ $\Delta \approx 0.5$ が得られるはずです。 $J_2=1$ では、正方格子にギャップが存在しないことを反映して、帯磁率は $T \to 0$ で有限の定数に近づくはずです。

相転移を特定する

$J_2=0$ と $J_2=1$ で異なる二つの相を同定したので、それらを隔てる量子相転移が（少なくとも一つ）存在するはずです。先ほど開始したバックグラウンドの計算では、システムサイズ $L=8,10,12,16$ について、逆温度 $\beta=2L$ で結合定数の範囲 $J_2 \in [0.2,0.4]$ をスキャンしています（パラメータの全体については parm8b / tutorial8b.py を参照してください）。計算が終了したら、以下の手順で結果を読み込んで解析しましょう。

逆温度の選び方

シミュレーションは厳密な $T=0$ ではなく、有限の逆温度 $\beta = 2L$ で実行されます。 $\beta \propto L$ という選択は恣意的なものではありません。この量子相転移では動的臨界指数が $z=1$ であり、これは時間と空間が同じようにスケールすることを意味します。 $\beta = 2L$ を保つことで、経路積分の虚時間方向の広がりがシステムサイズと同じ比率で大きくなり、シミュレーションが熱励起ではなく基底状態の物理をサンプリングすることが保証されます。 $\rho_s L$ の曲線が一点で交わるのを観察することで、この選択を暗黙のうちに検証することになります。この交差がきれいに現れるのは、 $z=1$ かつ $\beta \propto L^z$ の場合だけです。

また、ループアルゴリズムの計算時間は $\beta$ に線形にスケールすることにも注意してください。有限温度の量子モンテカルロアルゴリズムは時空体積 $\beta L^d$ に対して線形より良くスケールすることはできないので、このスケーリングは量子臨界点においてさえ最適です。

結果が有限温度の影響を受けていないことを確認するには、 $\beta=2L$ を $\beta=4L$ に変えて同じシミュレーションを繰り返します。スピン剛性とビンダーキュムラントは変化しないはずです。 $\beta=L/4$ も試して、系が基底状態から遠く離れている場合に結果がどのように悪化するかを観察してみましょう。

スタッガード磁化、ビンダーキュムラント、スピン剛性

古典モンテカルロのチュートリアルと同様に、二つの物理量を用いて相転移点を特定します。

一つ目は、反強磁性相の秩序変数であるスタッガード磁化 $m_s$ のビンダーキュムラントです。

U_4 = \frac{\langle m_s^4\rangle}{\langle m_s^2\rangle^2}.

二つ目はスピン剛性（Wenzel and Janke, Phys. Rev. B 79, 014410 (2009)）で、この模型ではビンダーキュムラントの交差よりも有限サイズ補正が小さくなります。

\rho_s = \frac{3}{4\beta} \langle w_x^2 + w_y^2\rangle,

ここで $w_x$ 、 $w_y$ はワールドラインの $x$ 方向および $y$ 方向の巻き付き数です。量子臨界点では $\rho_s \propto L^{d-2-z}$ となります。ここで $d$ は空間次元、 $z$ は動的臨界指数です。 $z=1$ の場合、組み合わせ $\rho_s L$ は臨界点で無次元となるため、異なる $L$ に対する曲線はすべて $J_2^c$ で交わります。どちらか一方が発散するのではなく、両方の物理量が一点で交わるという事実は、この転移が一次相転移ではなく連続相転移であることを示しています。

物理量の読み込みとプロットは次のように行います。

data = pyalps.loadMeasurements(pyalps.getResultFiles(pattern='parm8b.task*.out.h5'),['Binder Ratio of Staggered Magnetization','Stiffness'])

binder=pyalps.collectXY(data,x='J2',y='Binder Ratio of Staggered Magnetization', foreach=['L'])
stiffness =pyalps.collectXY(data,x='J2',y='Stiffness', foreach=['L'])

for q in stiffness:
    q.y = q.y*q.props['L']

plt.figure()
pyalps.plot.plot(stiffness)
plt.xlabel(r'$J_2$')
plt.ylabel(r'Stiffness $\rho_s L$')
plt.title('coupled ladders')

plt.figure()
pyalps.plot.plot(binder)
plt.xlabel(r'$J_2$')
plt.ylabel(r'$g(m_s)$')
plt.title('coupled ladders')
plt.show()

異なる $L$ に対する $\rho_s L$ の曲線は、左側でゼロから立ち上がり、ほぼ一点で交差して、 $J_2^c \approx 0.30$ – $0.31$ という最初の推定値を与えるはずです。ビンダーキュムラントの曲線も同じ領域で交差しますが、有限サイズ補正がより大きいため、これらのシステムサイズでは交点があまり鋭くありません。どちらか一方が発散するのではなく両方の物理量が交差することは、この転移が連続的であることを裏付けています。

臨界指数の推定

ここまでで量子臨界点 $J_2^c$ のおおまかな推定値が得られました。古典の場合と同様に、臨界指数を求めるには $J_2^c$ をより精密に決定する必要があります。

そのためには、parm8d および tutorial8d.py で設定されているように、より大きなシステムサイズについて $J_2$ のより細かいグリッドで計算を実行します。これらのシミュレーションは CPU 時間を多く消費するため、演習として残しておきます。異なるシステムサイズについてビンダーキュムラント $U_4$ と規格化したスピン剛性 $\rho_s L$ をプロットしましょう。その交点から $J_2^c$ のより精密な推定値が得られます。 $\nu$ を求めるには、これらの量の $J_2$ に関する微分を $J_2^c$ で評価したものが、システムサイズに対してどのようにスケールするかを考えます。これらの微分は原理的にはモンテカルロで直接測定できますが、このチュートリアルでは、細かい $J_2$ グリッドを利用した数値微分で計算すれば十分です。

両方の量について異なるシステムサイズで数値微分を行い、 $J_2^c$ での値をシステムサイズの関数としてプロットしましょう。データはべき則に従ってスケールするはずです。

\left.\frac{dU_4}{dJ_2}\right|_{J_2^c} \propto L\left.\frac{d\rho_s}{dJ_2}\right|_{J_2^c} \propto L^{1/\nu}.

べき則フィットから $\nu$ を求めましょう。どちらの量からも、三次元古典ハイゼンベルグ模型の値である $\nu \approx 0.71$ に近い、互いに整合する推定値が得られるはずです。

演習： 古典の場合と同様に、データコラプスを行うことで推定値の精度を視覚的に判断できます。 $U_4$ と $\rho_s L$ のスケーリング形は、古典相転移の例におけるビンダーキュムラントのものと同一です。

$z$ と $\nu$ のほかに、最後の独立な指数 $\eta$ は、臨界点におけるスタッガード帯磁率の有限サイズスケーリングから得られます。

\chi_s(J_2^c) \sim L^{2-\eta}.

ここでは一様帯磁率 $\chi$ ではなく、秩序変数のゆらぎに対応する帯磁率であるスタッガード帯磁率 $\chi_s$ を用いなければならないことに注意してください。 $J_2^c$ における $\chi_s$ をシステムサイズの関数としてプロットし、その傾きから $\eta$ を求めましょう。期待される値は $\eta \approx 0.034$ とゼロに近く、これは臨界点でスタッガード帯磁率がほぼ $L^2$ に比例して増大することを意味します。

このチュートリアルで扱う二つの量子相転移は、いずれも三次元古典ハイゼンベルグ模型の有限温度転移の普遍性クラスに属します。

二つ目の量子臨界点を特定する

$J_2$ が非常に大きい場合、梯子間結合が梯子内の結合 $J_0 = J_1 = 1$ を上回り、系は $J_2$ ボンドに沿った実効的な一重項二量体へと再編成されます。これも再びギャップをもつスピン液体相なので、相図は $J_2$ を 0 から大きな値へ増やすにつれて、スピン液体 → ネール反強磁性 → スピン液体という構造をもちます。したがって、ネール秩序が壊れる二つ目の量子臨界点 $J_2^{c_2}$ が存在するはずです。

パラメータファイル parm8c またはスクリプト tutorial8c.py を用いて、パラメータ領域 $J_2 \in [1.8, 2.1]$ で有限サイズスケーリング解析を繰り返します。これらは parm8b と同じシステムサイズと $\beta=2L$ を使いますが、より高い $J_2$ の範囲をスキャンします。

コマンドライン

parameter2xml parm8c
loop parm8c.in.xml

Python

import pyalps
import matplotlib.pyplot as plt
import pyalps.plot
import numpy as np

parms = []
for l in [8,10,12,16]:
    for j2 in [1.8,1.85,1.9,1.95,2.,2.05,2.1]:
        parms.append(
            {
              'LATTICE'        : "coupled ladders",
              'LATTICE_LIBRARY': 'lattices.xml',
              'MODEL_LIBRARY'  : 'models.xml',
              'local_S'        : 0.5,
              'ALGORITHM'      : 'loop',
              'SEED'           : 0,
              'BETA'           : 2*l,
              'J0'             : 1,
              'J1'             : 1,
              'J2'             : j2,
              'THERMALIZATION' : 5000,
              'SWEEPS'         : 50000,
              'MODEL'          : "spin",
              'L'              : l,
              'W'              : l//2
            }
        )

input_file = pyalps.writeInputFiles('parm8c', parms)
pyalps.runApplication('loop', input_file)

先ほどと同じ解析コマンドを使い、ファイルパターンの parm8b を parm8c に置き換えて、スピン剛性 $\rho_s L$ とビンダーキュムラントを読み込んでプロットしましょう。

$\rho_s L$ の曲線は再び一点で交差し、今度はその位置は $J_2^{c_2} \approx 1.91$ 付近になるはずです。ビンダーキュムラントの交差も同じ領域に現れるはずです。この二つ目の転移は最初の転移と同じ普遍性クラス、すなわち三次元古典ハイゼンベルグ模型に属するので、同じ臨界指数 $\nu$ と $\eta$ が当てはまります。

課題

$J_2=0$ と $J_2=1$ について磁化率をプロットしてください。低温での振る舞いはどのように異なり、そこから基底状態について何が分かりますか。
$J_2=0$ について、 $\chi$ を $A/\sqrt{T}\exp(-\Delta/T)$ にフィットしてスピンギャップ $\Delta$ を求めてください。得られた推定値は文献の値（Phys. Rev. Lett. 73, 886 (1994); Phys. Rev. Lett. 77, 1865 (1996)）と比べてどうですか。
$\rho_s L$ とビンダーキュムラントの曲線の交差から、最初の量子臨界点 $J_2^c$ の推定値はいくらになりますか。
$\beta=2L$ を $\beta=4L$ に変えて計算を繰り返してください。スピン剛性とビンダーキュムラントは影響を受けますか。次に $\beta=L/4$ を試してください。結果は目に見えて変わりますか。このことから、これらのシミュレーションにおける温度の役割について何が言えますか。
べき則スケーリング $dU_4/dJ_2|_{J_2^c} \propto L^{1/\nu}$ から、 $\nu$ の値はいくらになりますか。
スタッガード帯磁率の有限サイズスケーリング $\chi_s(J_2^c) \sim L^{2-\eta}$ から、 $\eta$ の値はいくらになりますか。
得られた $\nu$ と $\eta$ の推定値を、Phys. Rev. B 65, 144520 (2002) で報告されている三次元古典ハイゼンベルグ普遍性クラスの値と比較してください。
（発展）求めた $J_2^c$ と $\nu$ を用いて、 $U_4$ と $\rho_s L$ のデータコラプスを行ってください。コラプスの良さは $J_2^c$ の正確な値にどの程度敏感ですか。
$J_2 \in [1.8, 2.1]$ をスキャンして二つ目の量子臨界点を特定してください。得られた推定値は Wenzel and Janke, Phys. Rev. B 79, 014410 (2009) の高精度な結果と比べてどうですか。

MC-07 相転移