ED-05 ED Phase Transition

次近接相互作用を持つハイゼンベルク鎖の臨界点

このチュートリアルでは、ハイゼンベルク鎖を扱った ED-04 チュートリアルの最後の部分を引き継ぎます。ハミルトニアンに次近接結合項を加えます。

H = J_1 \sum_{\langle i,j \rangle} \mathbf{S}_i \cdot \mathbf{S}_j \; + \; J_2 \sum_{\langle\langle i,j \rangle\rangle} \mathbf{S}_i \cdot \mathbf{S}_j ,

ここで最初の和は同じ周期鎖上の最近接ボンドについて、2番目の和は次近接ボンドについて取られます。これは有名な $J_1$ – $J_2$ （あるいは Majumdar-Ghosh）鎖です。

$J_2 = 0$ の極限では、この模型は Bethe 仮説で解ける臨界ハイゼンベルク鎖に帰着します。 $J_2/J_1=0.5$ でもこの模型は厳密に解けます：これは Majumdar-Ghosh 点であり、そこでは基底状態が最近接一重項の厳密な直積になります（C.K. Majumdar and D.K. Ghosh, J. Math. Phys. 10, 1388 (1969)）。

|\Psi\rangle = \left(|\uparrow\rangle_1 |\downarrow\rangle_2 - |\downarrow\rangle_1 |\uparrow\rangle_2\right) (|\uparrow\rangle_3 |\downarrow\rangle_4 - |\downarrow\rangle_3 |\uparrow\rangle_4) (|\uparrow\rangle_5 |\downarrow\rangle_6 - |\downarrow\rangle_5 |\uparrow\rangle_6)

この二量体化した二重縮退の基底状態は、 $J_2=0$ におけるギャップレスで並進不変な基底状態とは定性的に異なります。これはもちろん、ある中間の $J_2/J_1 \in (0,1/2)$ に相転移があることを示しています——数値的にはこれは $J_2/J_1 \approx 0.2411$ 付近で起こることが知られています。

このチュートリアルの前半では、結合を調節するにつれてスペクトル、特に異なる対称性セクターにおけるギャップがどのように変化するかを見ることで、臨界点の位置を特定します。後半では、臨界鎖の CFT 内容を改めて検討します。解析的には、臨界点における模型はハイゼンベルク鎖と同じ CFT で記述されますが、対数的に消える有限サイズ補正を引き起こす辺境演算子の重みはゼロになるため、スケーリング次元をはるかに精度良く求めることができることが示せます。

パラメータ

パラメータ	意味	値
`LATTICE`	最近接・次近接ボンドを持つ周期鎖	`nnn chain lattice`
`MODEL`	量子スピン模型	`spin`
`local_S`	各サイトのスピン量子数	1/2
`J`	最近接結合 $J_1$	1
`J1`	次近接結合 $J_2$ （注：ALPS の `nnn chain lattice` はこのボンドパラメータを `J1` と名付けており、上のハミルトニアンで使われている $J_2$ とは表記が異なる）	前半では $[0,0.5]$ で走査、後半では 0.25 に固定
`CONSERVED_QUANTUMNUMBERS`、`Sz_total`	一重項（ $S_z=0$ ）と三重項（ $S_z=1$ ）セクターをそれぞれ分解する	`Sz`, 0 と 1
`NUMBER_EIGENVALUES`	`sparsediag` に要求する固有状態数	2（前半）または 5（後半）
`L`	鎖長	6, 8（前半）または 10, 12（後半）

格子

nnn chain lattice は、周期鎖に各サイトをその次近接サイトと結ぶ第2のボンド集合を追加したものです。

        J          J          J
   o---------o---------o---------o---(periodic)
   0         1         2         3
    \___________________________/
        J1 (=J2)   next-nearest-neighbour bonds

この内蔵格子は、ALPS格子ライブラリの他の部分とともに、これまでのチュートリアルで使用した単純な chain lattice の拡張として説明されています。

手法

ED-02 と同様に、単一項ギャップと三重項ギャップを得るために、sparsediag を用いて $S_z=0$ と $S_z=1$ のセクターをそれぞれ独立に対角化します。ここで最大のセクター（ $L=12$ 、 $S_z=0$ 、次元924）は Lanczos アルゴリズムにとって取るに足らないものです。ギャップの交差が見えるように十分密に $J_2$ を走査し、これによって臨界点の有限サイズ推定値を特定します。その後、この交差点に近い $J_2/J_1=0.25$ に固定して再度対角化を行い、CFT の塔を分解するためにより多くの固有状態を要求します。

そこでまず、基底状態と第一励起状態のエネルギー、および一重項（ $S_z = 0$ ）セクターと三重項（ $S_z=1$ ）セクターのギャップをプロットしてみましょう。

コマンドラインを用いる方法

parm5a は、 $L=6$ と $L=8$ の両方について、 $S_z=0$ と $S_z=1$ の両方のセクターで $J_2/J_1\in[0,0.5]$ を走査します。

MODEL="spin"
LATTICE="nnn chain lattice"
CONSERVED_QUANTUMNUMBERS="Sz"
local_S=1/2
J=1
NUMBER_EIGENVALUES=2
Sz_total=0
{ L=6; J1=0.0 }
{ L=6; J1=0.1 }
{ L=6; J1=0.2 }
{ L=6; J1=0.3 }
{ L=6; J1=0.4 }
{ L=6; J1=0.5 }
{ L=8; J1=0.0 }
{ L=8; J1=0.1 }
{ L=8; J1=0.2 }
{ L=8; J1=0.3 }
{ L=8; J1=0.4 }
{ L=8; J1=0.5 }
Sz_total=1
{ L=6; J1=0.0 }
{ L=6; J1=0.1 }
{ L=6; J1=0.2 }
{ L=6; J1=0.3 }
{ L=6; J1=0.4 }
{ L=6; J1=0.5 }
{ L=8; J1=0.0 }
{ L=8; J1=0.1 }
{ L=8; J1=0.2 }
{ L=8; J1=0.3 }
{ L=8; J1=0.4 }
{ L=8; J1=0.5 }

parameter2xml parm5a
sparsediag --write-xml parm5a.in.xml

Python

走査とデータ解析を自動化するために、スクリプト tutorial5a.py を使用します。いつも通りのインポートから始めます。

import pyalps
import pyalps.plot
from pyalps.dict_intersect import dict_intersect
import numpy as np
import matplotlib.pyplot as plt
import copy
import math

ここでも $S_z$ 量子数を使いますが、今回は異なるセクター $(S_z=0,1)$ でシミュレーションを実行します。システムサイズ $L=6,8$ で計算します。というのも、探している効果はすでに非常に小さいシステムサイズで現れるからです。

prefix = 'alps-nnn-heisenberg'
parms = []
for L in [6,8]:
    for Szt in [0,1]:
        for J1 in np.linspace(0,0.5,6):
            parms.append({
                'LATTICE'              : "nnn chain lattice",
                'MODEL'                : "spin",
                'local_S'              : 0.5,
                'J'                    : 1,
                'NUMBER_EIGENVALUES'   : 2,
                'CONSERVED_QUANTUMNUMBERS' : 'Sz',
                'Sz_total'             : Szt,
                'J1'                   : J1,
                'L'                    : L
            })

input_file = pyalps.writeInputFiles(prefix,parms)
res = pyalps.runApplication('sparsediag', input_file)
# res = pyalps.runApplication('sparsediag', input_file, MPI=4)
data = pyalps.loadEigenstateMeasurements(pyalps.getResultFiles(prefix=prefix))

この場合のデータ解析は、これまでのものより少し込み入っています。特に、階層的データセットという機能に大きく依存します。物理を理解するには、実際には基底状態と第一励起状態だけを見れば十分なので、ギャップの計算があまりに紛らわしいと感じても、あまり気にする必要はありません。

最初のステップとして、与えられた J1、L、Sz_total の組に対するすべてのエネルギーをまとめてソートします。まず、パラメータ——J1、L、Sz_total——でグループ化します。grouped についてのループの各要素は、したがって異なる運動量に対応するデータセットのリストを含むことになります。これらをまとめた後、dict_intersect 関数を使って結果のデータセットのプロパティを求めます。この関数は辞書のリストを受け取り、それらすべてに共通する部分だけを返します。numpy の argsort 関数を使って y をソートするインデックスのリストを取得し、これによって x も対応してソートできますが、おそらくこれは必要ないでしょう。

grouped = pyalps.groupSets(pyalps.flatten(data), ['J1', 'L', 'Sz_total'])
nd = []
for group in grouped:
    ally = []
    allx = []
    for q in group:
        ally += list(q.y)
        allx += list(q.x)
    r = pyalps.DataSet()
    sel = np.argsort(ally)
    r.y = np.array(ally)[sel]
    r.x = np.array(allx)[sel]
    r.props = dict_intersect([q.props for q in group])
    nd.append( r )
data = nd

次に、 $S_z=1$ セクターに現れる状態を $S_z=0$ セクターから取り除く必要があります。J1、L でグループ化することで、各グループが2つの異なる Sz_total セクターのスペクトルを含むようにします。そして関数 subtract_spectrum を使います。これは、第一引数として渡されたデータセットから、第二引数にも含まれる要素を取り除くものです。オプション引数として、最大相対差を指定できます。

grouped = pyalps.groupSets(pyalps.flatten(data), ['J1', 'L'])
nd = []
for group in grouped:
    if group[0].props['Sz_total'] == 0:
        s0 = group[0]
        s1 = group[1]
    else:
        s0 = group[1]
        s1 = group[0]
    s0 = pyalps.subtract_spectrum(s0, s1)
    nd.append(s0)
    nd.append(s1)
data = nd

ここで、基底状態（‘gs’）または第一励起状態（‘fe’）のエネルギーのみを含む新しいデータセットのリスト（sector_E）を作成します。この情報はプロパティ which に格納します。これにより後で collectXY 関数を使い、各 L について基底状態と第一励起状態のエネルギーを結合の関数としてプロットできるようになります。

sector_E = []
grouped = pyalps.groupSets(pyalps.flatten(data), ['Sz_total', 'J1', 'L'])
for group in grouped:
    allE = []
    for q in group:
        allE += list(q.y)
    allE = np.sort(allE)
    
    d = pyalps.DataSet()
    d.props = dict_intersect([q.props for q in group])
    d.x = np.array([0])
    d.y = np.array([allE[0]])
    d.props['which'] = 'gs'
    sector_E.append(d)
    
    d2 = copy.deepcopy(d)
    d2.y = np.array([allE[1]])
    d2.props['which'] = 'fe'
    sector_E.append(d2)

sector_energies = pyalps.collectXY(sector_E, 'J1', 'Energy', ['Sz_total', 'which', 'L'])
plt.figure()
pyalps.plot.plot(sector_energies)
plt.xlabel('$J_1/J$')
plt.ylabel('$E_0$')
plt.legend(prop={'size':8})

最後のステップとして、一重項ギャップと三重項ギャップを計算します。これらはそれぞれ、系の最低状態のエネルギーと、a) 一重項（ $S_z=0$ ）セクターの第一励起状態、b) 三重項（ $S_z=1$ ）セクターの最低状態、とのエネルギー差として定義されます。

grouped = pyalps.groupSets( pyalps.groupSets(pyalps.flatten(data), ['J1', 'L']), ['Sz_total'])

gaps = []
for J1g in grouped:
    totalmin = 1000
    for q in pyalps.flatten(J1g):
        totalmin = min(totalmin, np.min(q.y))
    
    for Szg in J1g:
        allE = []
        for q in Szg:
            allE += list(q.y)
        allE = np.sort(allE)
        d = pyalps.DataSet()
        d.props = pyalps.dict_intersect([q.props for q in Szg])
        d.props['observable'] = 'gap'
        print totalmin,d.props['Sz_total']
        if d.props['Sz_total'] == 0:
            d.y = np.array([allE[1]-totalmin])
        else:
            d.y = np.array([allE[0]-totalmin])
        d.x = np.array([0])
        d.props['line'] = '.-'
        gaps.append(d)

gaps = pyalps.collectXY(gaps, 'J1', 'gap', ['Sz_total', 'L'])

plt.figure()
pyalps.plot.plot(gaps)
plt.xlabel('$J_1/J$')
plt.ylabel('$\Delta$')
plt.legend(prop={'size':8})

plt.show()

出力データ

各セクターを独立に対角化すると、 $J_2/J_1$ （上述の通り ALPS パラメータでは J1 と表記）の関数として、次のような一重項ギャップと三重項ギャップが得られます。

$L$	$J_2/J_1$	一重項ギャップ $/J_1$	三重項ギャップ $/J_1$
6	0.0	1.3028	0.6847
6	0.1	1.0300	0.6688
6	0.2	0.7621	0.6455
6	0.25	0.6302	0.6302
6	0.3	0.5000	0.6119
6	0.5	0.0000	0.5000
8	0.0	0.9515	0.5227
8	0.1	0.7551	0.5039
8	0.2	0.5583	0.4800
8	0.25	0.4601	0.4663
8	0.3	0.3624	0.4516
8	0.5	0.0000	0.4045

どちらのシステムサイズについても、一重項ギャップ（ $J_2=0$ ではギャップレスで最低励起が $L\to\infty$ でのみゼロに近づくため、当初はより大きい）は $J_2/J_1=0.2$ から $0.3$ の間で三重項ギャップを下回り、ほぼちょうど $J_2/J_1\approx0.25$ で交差します——このような小さいサイズであってもすでに、受け入れられている熱力学極限の値 $J_2/J_1\approx0.2411$ に近い値です。この準位交差は、二量体化相への転移の有限サイズにおける兆候です。また、一重項ギャップが Majumdar-Ghosh 点 $J_2/J_1=0.5$ でちょうどゼロになることにも注目してください。これはその点での厳密な二重縮退基底状態と整合しています。

次近接結合を持つハイゼンベルク鎖：CFT の帰属

フラストレートした J1-J2 鎖の臨界点に調節することで、有限サイズ補正を大幅に減らすことができます。結合が異なるにもかかわらず、この模型は同じ連続臨界場の理論を持つことが示せるため、この極限でスケーリング次元を抽出できます。この点についての詳しい議論は、参考文献 I. Affleck, D. Gepner, H.J. Schulz and T. Ziman, J. Phys. A: Math. Gen. 22, 511 (1989) をご覧ください。

上で得られたスペクトルと比較してみてください。期待されるスケーリング次元との対応がはるかに見やすく、システムサイズが大きくなるにつれてはるかに速く収束することがわかるでしょう。

コマンドラインを用いる方法

parm5b は、上で見つけた交差点に近い $J_2/J_1=0.25$ に固定し、2つのシステムサイズについて $S_z=0$ セクターで5個の固有状態を要求します。

MODEL="spin"
LATTICE="nnn chain lattice"
CONSERVED_QUANTUMNUMBERS="Sz"
local_S=1/2
J=1
J1=0.25
NUMBER_EIGENVALUES=5
Sz_total=0
{ L=10 }
{ L=12 }

parameter2xml parm5b
sparsediag --write-xml parm5b.in.xml

Python

スクリプト tutorial5b.py で使われている新しいパラメータは次の通りです。

parms_ = {
    'LATTICE'              : "nnn chain lattice",
    'MODEL'                : "spin",
    'local_S'              : 0.5,
    'J'                    : 1,
    'J1'                   : 0.25,
    'NUMBER_EIGENVALUES'   : 5,
    'CONSERVED_QUANTUMNUMBERS' : 'Sz',
    'Sz_total' : 0
}
prefix = 'nnn-heisenberg'
parms = []
for L in [10,12]:
    parms_.update({'L':L})
    parms.append(copy.deepcopy(parms_))

スクリプトの残りの部分は、ED-04 のハイゼンベルク鎖 CFT 解析と全く同じように進みます：入力ファイルを書き出して実行し、固有状態の測定値を読み込み、 $S_z=0$ セクターの最低2状態間のギャップで再スケールし、期待されるスケーリング次元を重ねて表示します。

input_file = pyalps.writeInputFiles(prefix,parms)
res = pyalps.runApplication('sparsediag', input_file)
data = pyalps.loadEigenstateMeasurements(pyalps.getResultFiles(prefix=prefix))

E0 = {}
E1 = {}
for Lsets in data:
    L = pyalps.flatten(Lsets)[0].props['L']
    allE = []
    for q in pyalps.flatten(Lsets):
        allE += list(q.y)
    allE = np.sort(allE)
    E0[L] = allE[0]
    E1[L] = allE[1]

for q in pyalps.flatten(data):
    L = q.props['L']
    q.y = (q.y-E0[L])/(E1[L]-E0[L]) * (1./2.)
spectrum = pyalps.collectXY(data, 'TOTAL_MOMENTUM', 'Energy', foreach=['L'])

for SD in [0.5, 1, 1.5, 2]:
    d = pyalps.DataSet()
    d.x = np.array([0,4])
    d.y = SD+0*d.x
    spectrum += [d]

pyalps.plot.plot(spectrum)
plt.legend(prop={'size':8})
plt.xlabel("$k$")
plt.ylabel("$E_0$")
plt.xlim(-0.02, math.pi+0.02)
plt.show()

出力データ

$J_2/J_1=0.25$ において独立に対角化されたスペクトルを上と同じ方法で再スケールすると、次のようになります。

$L$	準位1	準位2	準位3	準位4, 5
10	0	0.500（定義）	0.947	1.401
12	0	0.500（定義）	0.960	1.432

これを、ED-04 における通常のハイゼンベルク鎖（ $J_2=0$ ）の同じサイズでの値と比較してみましょう：1に収束するはずの準位について、 $L=10$ で0.880、 $L=12$ で0.857でした。 $J_2/J_1=0.25$ では、同じ準位はすでに0.947と0.960に達しています——目標値により近く、かつ $L$ が大きくなるにつれて正しい方向に動いています。これはまさに、ED-04 で対数補正を引き起こしていた辺境演算子が調節によって取り除かれた場合に期待される通りです。

まとめ

フラストレートした $J_1$ – $J_2$ 鎖の一重項ギャップと三重項ギャップを比較することで、有限サイズの準位交差により、 $L=6,8$ の鎖だけを使って、二量体化転移を受け入れられている値 $J_2/J_1\approx0.2411$ から数パーセント以内に特定できます。さらに、その点付近での対角化により、ED-04 で同定されたのと同じ $c=1$ CFT 演算子内容がここではるかにきれいに現れることがわかります。これは、 $J_2$ を調節することで、フラストレーションのない鎖における緩やかな対数的有限サイズ補正を引き起こす辺境演算子が取り除かれるためです。

問題

一重項ギャップと三重項ギャップの表から、隣接するデータ点の間で線形補間することにより、 $L=6$ と $L=8$ について交差点 $J_2/J_1$ をそれぞれ推定してください。 $L$ が大きくなるにつれて、この推定値は $0.2411$ に近づいていきますか。
なぜ三重項ギャップではなく一重項ギャップが、Majumdar-Ghosh 点 $J_2/J_1=0.5$ でちょうどゼロになるのでしょうか。
$J_2/J_1=0.25$ における $L=10,12$ の再スケールされた準位3の値を、ED-04 の $J_2=0$ における対応する値と比較してください。どちらが CFT の予言する1により近いですか。どちらが $L$ に対してより速く収束しますか。

ED-04 Criticality ED-06 Full Diagonalization