ED-04 Criticality

このチュートリアルでは、臨界スピン鎖を扱い、それを共形場理論の言葉で記述することとの関係を見ていきます。連続量子相転移では、スピン鎖はスケール不変になり、その長距離物理は $(1+1)$ 次元の共形場理論（CFT）によって記述されます。有限鎖が本当の意味で転移点にあることは決してありませんが、それでも転移の普遍的な痕跡を受け継いでいます：低エネルギーの有限サイズスペクトルは、CFT の初等演算子とその子孫演算子のスケーリング次元 $\Delta$ でラベル付けされた塔状の構造に組織化され、準位間隔は $1/L$ として消えていきます。厳密対角化はこれを直接見るための理想的な道具です。数十サイト程度の鎖について、運動量も含めた完全な低エネルギースペクトルが得られ、それだけですでに普遍的な演算子内容を分解するのに十分だからです。

イジング鎖

最初に考える模型は、次のハミルトニアンで与えられる臨界イジング鎖です。

H=J_{z} \sum_{\langle i,j \rangle} S^i_z S^j_z + \Gamma \sum_i S^i_x

ここで最初の和は最近接対について取られます。 $\Gamma$ は横磁場と呼ばれ、系は $\Gamma/J=\frac{1}{2}$ で臨界になります。 $\Gamma=0$ のとき、基底状態は $J\gt 0$ では反強磁性的、 $J \lt 0$ では強磁性的です。この系は厳密に解けます（P. Pfeuty, Annals of Physics 57, 79 (1970)）。

上の式で、 $\Delta$ はその場のスケーリング次元を表します。スケーリング場は群として現れます：最も低いものは初等場（primary field）と呼ばれ、スケーリング次元 $\Delta + m$ （ $m \in \lbrace 1, 2, 3, ... \rbrace$ ）を持つ無限個の子孫場を伴います。

イジング模型の厳密解（Pfeuty の論文の式 (3.7)）では、長距離相関は次のように減衰することがわかっています。

\langle S^i_z S^{i+n}_z \rangle \sim n^{-2\times 1/8}

\langle S^i_y S^{i+n}_y \rangle \sim n^{-2\times(1+1/8)}

\langle S^i_x S^{i+n}_x \rangle \sim n^{-2\times 1}

さらに、恒等演算子のスケーリング次元は 0 であると期待されます。

したがって、イジング模型の CFT には 0、1/8、1、1+1/8 というスケーリング次元が現れると期待されます。これを見るために、 $E \rightarrow \frac{E-E_0}{(E_1-E_0)8}$ に従ってスペクトルのすべてのエネルギーを再スケールします。これにより、最も低い2つの状態が期待されるスケーリング次元の位置にちょうど来るようになり、残りのスペクトルがこれと整合しているかどうかを確認できます。

パラメータ

パラメータ	意味	値
`LATTICE`	周期鎖	`chain lattice`
`MODEL`	量子スピン模型	`spin`
`local_S`	各サイトのスピン量子数	1/2
`Jxy`	面内（ $S^xS^x+S^yS^y$ ）結合	0
`Jz`	イジング（ $S^zS^z$ ）結合	$-1$
`Gamma`	横磁場	0.5（臨界点、$\Gamma/
`NUMBER_EIGENVALUES`	`sparsediag` に要求する低エネルギー固有状態の数	5
`L`	鎖長	10, 12

$\Gamma\neq 0$ になると（横磁場がスピンを反転させるため） $S_z$ は保存されないことに注意してください。したがって、これまでのチュートリアルとは異なり CONSERVED_QUANTUMNUMBERS の行はなく、sparsediag は全運動量のみで分解された鎖全体を対角化します。

格子

臨界横磁場イジング鎖は、これまでのチュートリアルと同じ ALPS格子ライブラリの周期的な chain lattice を用い、ここではイジング結合 $J_z$ とオンサイトの横磁場 $\Gamma$ を持ちます。

    Jz    Jz    Jz          Jz
o------o------o------ ... ------o     (each o also feels an on-site field Gamma S^x)
0      1      2                L-1
|______________________________________|
                Jz   (bond L-1 -- 0, periodic)

手法

$\Gamma=0.5$ では鎖は保存される $S_z$ を持たないため、sparsediag の Lanczos 反復は完全な $2^L$ 次元のヒルベルト空間全体で実行されます（ $L=10,12$ に対してそれぞれ $1024$ 、 $4096$ 状態のみ）。NUMBER_EIGENVALUES=5 を要求すると、初等場とその最初の子孫場を見るのに必要な最低5状態が返されます。完全なスペクトルではなく少数の低エネルギー状態のみが必要なため、対称性セクターを使って問題規模を縮小していなくても、スパース Lanczos 対角化は全対角化よりはるかに安価です。

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

パラメータファイル parm_ising は2種類のシステムサイズを設定します。

LATTICE    = "chain lattice"
MODEL      = "spin"
local_S    = 0.5
Jxy        = 0
Jz         = -1
Gamma      = 0.5
NUMBER_EIGENVALUES = 5

{L=10}
{L=12}

parameter2xml parm_ising
sparsediag --write-xml parm_ising.in.xml

Python 版

イジングシミュレーション用のパラメータファイルはこちらにあります。シミュレーションは Python スクリプトで実行されます。

まずいくつかのモジュールをインポートします。

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

続いて、2つのシステムサイズについてパラメータを設定しましょう。縦磁場 $h$ ではなく横磁場 $\Gamma$ を使うように注意してください。

parms = []
for L in [10,12]:
    parms.append({
        'LATTICE'    : "chain lattice",
        'MODEL'      : "spin",
        'local_S'    : 0.5,
        'Jxy'        : 0,
        'Jz'         : -1,
        'Gamma'      : 0.5,
        'NUMBER_EIGENVALUES' : 5,
        'L'          : L
    })

ご覧の通り、2つのシステムサイズをシミュレーションします。では入力ファイルを設定してシミュレーションを実行しましょう。

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

最後から2番目の行のコメントを外し、ジョブ数と mpirun 実行ファイル名（デフォルトは mpirun）を適切に調整すれば、ALPS に複数の CPU を同時に使わせることができます。

まず、各 L の値について最低状態と第一励起状態を抽出し、辞書にまとめます。

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]

上記のコードが機能するのは、ALPS がシミュレーションごとにグループ化されたリストとしてエネルギーを読み込むことがわかっているからです——つまり data はリストのリストであり、最上位に異なるシミュレーション、その下に異なる運動量が並びます。ここで、上で与えた式に従ってエネルギーを再スケールし、データを運動量の関数として集めます。

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

比較のために、初等場とその最初の数個の子孫場も表示しましょう。

for SD in [0.125, 1, 1+0.125, 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()

出力データ

同じハミルトニアンを独立に対角化すると、2つの鎖長について次のような再スケール後の低エネルギー準位（運動量について和を取り、ソートしたもの）が得られます。

$L$	準位1	準位2	準位3	準位4
10	0（定義）	0.125（定義）	0.994	1.107, 1.107
12	0（定義）	0.125（定義）	0.996	1.112, 1.112

構成上、最初の2つの再スケール準位はちょうど0と0.125に位置します。3番目と4番目の準位は対角化による真の予言であり、このように小さいサイズであってもすでに期待される CFT 値 1 と 1.125 から数パーセント以内に収まっており、 $L$ が大きくなるにつれてさらに収束していきます——これはイジング CFT の演算子内容 $\{0, \tfrac18, 1, 1+\tfrac18, \ldots\}$ の直接的な数値的証拠です。

ハイゼンベルク鎖

次に、より複雑な例として、次式で記述されるスピン1/2自由度の反強磁性ハイゼンベルク鎖を考えましょう。

H = J\sum_{\langle i,j \rangle} \mathbf{S}^i \cdot \mathbf{S}^j , \qquad J=1 .

この模型の臨界理論は中心電荷 $c=1$ を持ち、初等場は 0、0.5、1 です。イジング模型とは異なり、有限サイズ補正は対数的にしか消えないため、われわれが到達できるシステムサイズではかなり顕著に現れます。

パラメータ

パラメータ	意味	値
`LATTICE`	周期鎖	`chain lattice`
`MODEL`	量子スピン模型	`spin`
`local_S`	各サイトのスピン量子数	1/2
`Jxy`、`Jz`	面内結合とイジング結合（ここでは等方的）	1, 1
`NUMBER_EIGENVALUES`	各運動量セクターについて要求する固有状態数	5
`CONSERVED_QUANTUMNUMBERS`、`Sz_total`	必要な CFT の塔全体を含む $S_z=0$ セクターに制限する	`Sz`, 0
`L`	鎖長	10, 12

格子

上と同じ周期的な chain lattice（ALPS格子ライブラリ参照）ですが、ここでは等方的なハイゼンベルク結合を持ち、横磁場はありません。

    J     J     J           J
o-------o-------o--- ... ---o
0       1       2          L-1
|_________________________________|
                J   (bond L-1 -- 0, periodic)

手法

イジング点とは異なり、ハイゼンベルク鎖は $S_z$ を保存するため、sparsediag を $S_z=0$ セクターに制限できます——このセクターにはすでに基底状態と、CFT 内容を読み取るのに十分な低エネルギーの塔が含まれており、同時にセクターの次元を（ $L=12$ の場合） $2^{12}=4096$ から $\binom{12}{6}=924$ に減らすことができます。

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

パラメータファイルは2種類のシステムサイズを設定します。

LATTICE    = "chain lattice"
MODEL      = "spin"
local_S    = 0.5
Jxy        = 1
Jz         = 1
NUMBER_EIGENVALUES = 5
CONSERVED_QUANTUMNUMBERS = 'Sz'
Sz_total   = 0

{L=10}
{L=12}

parameter2xml parm_heisenberg
sparsediag --write-xml parm_heisenberg.in.xml

Python 版

この Python ファイルは上記の説明と同様の構成なので、ここでは詳しく説明しません。主な違いは、系の U(1) 対称性を利用できることです。また、関連するすべての状態を含む $S_z = 0$ セクターについてのみシミュレーションを実行します。heisenberg.py で使われている再スケーリングでは、ここでの最も低い初等場が $\Delta_1=1/2$ であるため、第一励起準位を 0.125 ではなく 0.5（スピン演算子のスケーリング次元）に固定します。

スペクトルを見て、さまざまなスケーリング場を同定し、システムサイズが大きくなるにつれてそれらが正しい値にどのように近づいていくかを調べてみてください。これがかなり難しいことに気づくでしょう。この系の詳しい議論については、I. Affleck, D. Gepner, H.J. Schulz and T. Ziman, J. Phys. A: Math. Gen. 22, 511 (1989) および次のチュートリアルを参照してください。

出力データ

独立に対角化を行い、 $(E-E_0)/\Delta_{S_z=1}\times\tfrac12$ （ $\Delta_{S_z=1}=E_0(S_z{=}1)-E_0(S_z{=}0)$ は三重項ギャップ）で再スケールした $S_z=0$ セクターの準位は次のようになります。

$L$	三重項ギャップ $\Delta/J$	準位2	準位3	準位4, 5
10	0.4232	0.500（定義）	0.880	1.149, 1.149
12	0.3558	0.500（定義）	0.857	1.149, 1.149
14	0.3071	0.500（定義）	0.839	1.148, 1.148
16	0.2702	0.500（定義）	0.825	1.145, 1.145

スケーリング次元1の子孫場に収束すると期待される3番目の準位は、イジングの準位がそうであったように単調に近づくのではなく、 $L=10$ の0.880から $L=16$ の0.825へと、むしろゆっくりと目標値から遠ざかっていきます。この鈍く非単調な振る舞いは、上で述べた辺境的無関係演算子の実際の症状であり、その対数的に抑制された補正は到達可能な $L$ の範囲でサイズ依存性を支配し、天文学的に大きなサイズになって初めて薄れていきます。ED-05 では、近くのフラストレーションのある点に調整することでこの辺境演算子を取り除き、速い教科書通りの $1/L$ 収束を回復する方法を示します。

まとめ

厳密に解ける横磁場イジング鎖も、収束がより困難なハイゼンベルク鎖も、いずれも有限サイズの厳密対角化スペクトルだけから、それぞれの場の理論が予言する普遍的な CFT 演算子内容（イジング： $c=1/2$ 、次元 $\{0,\tfrac18,1,\ldots\}$ ；ハイゼンベルク： $c=1$ 、次元 $\{0,\tfrac12,1,\ldots\}$ ）を確認します——ハイゼンベルクの例はまた、辺境的無関係演算子がその収束を数値的に紛らわしいものにし得ることも示しています。

問題

イジング鎖について、なぜ再スケールされたエネルギーは3番目の準位から縮退した対として現れるのでしょうか（例えば1.107付近の2つの状態）。（ヒント：対応する子孫演算子の運動量量子数について考えてみてください。）
臨界点から離れた $\Gamma$ の値（例えば $\Gamma=0.2$ や $\Gamma=1.0$ ）についてイジング CFT 解析を繰り返してください。 $L$ を大きくしても、再スケールされたスペクトルは依然として CFT 的に見えるでしょうか。
ハイゼンベルク鎖について、3番目の再スケール準位は $L$ が10から16に増えるにつれて、CFT の目標値に落ち着いていくのではなく遠ざかっていきます。可能であれば、さらにいくつかのサイズ（例えば $L=18,20$ ）を計算してみてください——この漂移は最終的に反転するでしょうか。また、手に入るシステムサイズが少数しかない状況で、「非常に緩やかな単調収束」と「本質的に非単調な振る舞い」をどのように区別すればよいでしょうか。

ED-03 Spectra ED-05 ED Phase Transition