《寄稿論文》

物性研究・電子版 Vol. 5, No. 2, 052101 （2016年5月号）

密度行列繰り込み群無限系法初期推定に関する研究

横浜国立大学大学院工学府物理情報工学専攻 府川慎之介1

目次
1

はじめに

1

2

無限系法初期推定

2

2.1 無限系法繰り込み群・・・・・・・・・・・・・・・・・・・・・・・・・・・・・2
2.2 無限系法初期推定・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・4
2.3 初期推定符号調整の改良・・・・・・・・・・・・・・・・・・・・・・・・・・・6
2.4 近似推定による初期推定状態の改良・・・・・・・・・・・・・・・・・・・・・・10
2.5 基底対応付けの改良・・・・・・・・・・・・・・・・・・・・・・・・・・・・・12
2.6 高磁場領域計算時の初期推定・・・・・・・・・・・・・・・・・・・・・・・・・14345

基底・低励起状態計算
3.1

𝑆 = 1/2ハイゼンベルグ鎖・・・・・・・・・・・・・・・・・・・・・・・・・・・17

3.2

𝑆 = 1ハイゼンベルグ鎖・・・・・・・・・・・・・・・・・・・・・・・・・・・・18

3.3

𝑆 = 1/2ジグザグ鎖・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・20

高磁場領域計算23

4.1

𝑆 = 1/2ハイゼンベルグ鎖・・・・・・・・・・・・・・・・・・・・・・・・・・・23

4.2

𝑆 = 1ハイゼンベルグ鎖・・・・・・・・・・・・・・・・・・・・・・・・・・・・25
26

まとめ

参考文献

1

17

E-mail:fukawa-shinnosuke-nm@ynu.jp

《寄稿論文》

物性研究・電子版 Vol. 5, No. 2, 052101 （2016年5月号）

1 はじめにWhite によって 1992 年に考案された密度行列繰り込み群法(density matrix renormalizationgroup method,DMRG 法)[1]は、現在では大規模スピン系計算の主流となっている。厳密対角化ではメモリや計算時間の問題でせいぜい数十サイトの計算が限度であるが、DMRG 法では数百サイ
トの計算が可能であり、熱力学的極限の性質を十分な精度で論じることができるようになった。
DMRG 法では保持基底数によって計算精度が決まる。系が大きくなるほど精度を出すために必要な保持基底数が増え計算時間が増加する。
DMRG 計算の高速化としては、行列計算の並列化[2]、ハミルトニアンの粗行列性の利用[3]、基底変換による固有状態の初期推定[4]などがある。
DMRG 法ではランチョス法などの反復法で固有状態を求めるが、あらかじめ固有状態に近い状態を用意することで反復回数を減らし計算を高速化できる。DMRG 無限系法と呼ばれる計算では繰り込み群としての性質を利用した初期推定方法[5,6]が提案されているが、一般的な系への適用が不可能であることや、初期推定計算における符号調整が不完全であるなどの問題がある。
そこで本論文では初期推定計算の改良を行い、𝑆 = 1/2 ,1反強磁性ハイゼンベルグ鎖
𝐿−1

𝐻𝐻𝑒𝑖𝑠𝑒𝑛𝑏𝑒𝑟𝑔 = ∑ 𝑺𝑖 ∙ 𝑺𝑖+1

(1)

𝑖=1

及び𝑆 = 1/2ジグザグ鎖
𝐿−1

𝐿−2

𝐻𝑧𝑖𝑔𝑧𝑎𝑔 = 𝐽1 ∑ 𝑺𝑖 ∙ 𝑺𝑖+1 + 𝐽2 ∑ 𝑺𝑖 ∙ 𝑺𝑖+2
𝑖=1

(2)

𝑖=1

の基底・低励起状態の計算に適用した。ここで𝑺𝑖 はサイト𝑖のスピン演算子、𝐿はサイト数、𝐽1, 𝐽2は最隣接,次隣接相互作用定数である。またこれまで行われていなかった高磁場領域計算時の初期推定法を考案し、𝑆 = 1/2 ,1ハイゼンベルグ鎖に対し適用した。
本改良初期推定による計算の高速化は劇的である。例えば𝑆 = 1/2ハイゼンベルグ鎖の基底状態
計算では、初期ベクトルをランダムに設定した場合のランチョス法の反復回数は 100 回を超えるが、本論文の改良初期推定を行うと平均 10 回程度になり、DMRG 無限系法全体の計算で約 7 倍高速化した。また初期推定の精度は有限サイズ効果と関係があり、系が大きくなるほど精度が上がり計算時間が減る。
本論文の構成は次のようになる；２章では DMRG 法と初期推定方法の概要を解説したのち、初期推定の改良と高磁場領域での推定について論じた。３・４章では(1)(2)式のモデルへ適用し初期推定が確かに改良されることを示したほか、現れる現象から初期推定計算の性質を調査した。

1

《寄稿論文》

物性研究・電子版 Vol. 5, No. 2, 052101 （2016年5月号）

2 無限系法初期推定
2.1 DMRG 無限系法
DMRG 法には𝐿を次々に大きくする無限系法と、𝐿は一定のまま精度を高めてゆく有限系法の２種類がある。今回行った初期推定は無限系法での計算であるため、はじめにその概要を説明する。
DMRG 法では系全体(superblock)を system と environment の２つに分ける。system およびenvironment は共にブロックとスピンからなる。ランチョス法などによりターゲットとなる系全体のある固有状態|𝜓⟩を計算したのち、system に対する密度行列𝜌𝑖𝑖 ′ を計算する。
|𝜓⟩ = ∑ 𝜓𝑖𝑗 |𝑖⟩𝑠 |𝑗⟩𝑒𝑖,𝑗𝜌𝑖𝑖 ′ = ∑ 𝜓𝑖𝑗 𝜓𝑖 ′ 𝑗

(3)
(4)

𝑗

ここで|𝑖⟩𝑠 , |𝑗⟩𝑒 は system、environment の基底であり、𝜌𝑖𝑖 ′ は密度行列の𝑖𝑖 ′成分である。密度行列の対角化を行い、固有値の大きい𝑚個の固有状態を新たな基底とする基底変換を system に施す(繰り込み)。こうしてサイト数が１つ大きいブロックを計算し、このブロックから再び system とenvironment を構成し上記の計算を繰り返すことで状態数を増やさず次々と計算してゆく方法がDMRG 法である。特に DMRG 無限系法では environment は system を反転させたものを用いる。初めに４サイトでの対角化を行い、上に書いた手順に従ってブロックに繰り込み、系全体のサイト数を各ステップで 2 ずつ増やしながら目的のサイト数に達するまで計算を続ける。
DMRG 法における計算誤差は密度行列固有値𝜆𝑖 から見積もることができる。(4)式から𝜆𝑖 の総和
は∑𝑖 𝜆𝑖 = ∑𝑖 𝜌𝑖𝑖 = 1となるが、𝜆𝑖 を大きさで並び替えたとき切り捨てられる基底の固有値の総和𝑒𝑡𝑟𝑚

𝑒𝑡𝑟 = 1 − ∑ 𝜆𝑖
𝑖=1がおおよその計算精度となる。ここで𝑚は保持基底数である。

図１

DMRG 無限系法の概要

2(5)

《寄稿論文》

物性研究・電子版 Vol. 5, No. 2, 052101 （2016年5月号）

なお図１の系は開いた境界条件(OBC)の場合のものである。ブロックとスピンの設定によっては周期境界条件(PBC)の計算も可能だが、OBC に比べ𝑒𝑡𝑟 が 1,2 桁大きくなることが報告されている[1]。今回の計算では PBC で行う特別の理由はないため、すべての計算を OBC で行った。
多くのスピン系の計算ではハミルトニアンの性質から計算量を減らすことができる。サイト𝑖の𝑧
スピン演算子の z 成分𝑆𝑖𝑧 に対し、スピンの z 成分の総和𝑆𝑡𝑜𝑡𝑎𝑙は𝐿𝑧
𝑆𝑡𝑜𝑡𝑎𝑙= ∑ 𝑆𝑖𝑧

(6)

𝑖=1

𝑧
𝑧
で与えられるが、𝑆𝑡𝑜𝑡𝑎𝑙とハミルトニアン𝐻が交換するとき、𝐻は𝑆𝑡𝑜𝑡𝑎𝑙
が同じ部分空間でのみ値を𝑧持つ行列となり、対角化計算は𝑆𝑡𝑜𝑡𝑎𝑙
一定の部分空間でのみ行えばよい。この時計算されるエネル𝑧
ギー固有状態は一定の𝑆𝑡𝑜𝑡𝑎𝑙を持つ同時固有関数になっている。またそのような状態から密度行列𝑧を作ると、(4)式からわかるように system 部分のスピン z 成分の総和𝑆𝑠𝑦𝑠𝑡𝑒𝑚
𝐿 ⁄2𝑧
𝑆𝑠𝑦𝑠𝑡𝑒𝑚(7)= ∑ 𝑆𝑖𝑧
𝑖=1𝑧が同じ部分空間でのみ値を持つ行列となり、固有状態も一定の𝑆𝑠𝑦𝑠𝑡𝑒𝑚を持つ。この状態はDMRGの𝑧𝑧(𝐿 + 2)は𝑆𝑠𝑦𝑠𝑡𝑒𝑚
次ステップのブロック基底であり、ブロックのスピンz成分の総和𝑆𝑏𝑙𝑜𝑐𝑘(𝐿)と対応している。
ハミルトニアンがこの性質を持っているとき、容易に磁化曲線を計算することができる。磁場ℎによるエネルギー項(ゼーマン項)；𝑧
𝐻𝑍𝑒𝑒𝑚𝑎𝑛 = −ℎ ∑ 𝑆𝑖𝑧 = −ℎ𝑆𝑡𝑜𝑡𝑎𝑙

(8)

𝑖
𝑧
を考えると、元のハミルトニアン𝐻が𝑆𝑡𝑜𝑡𝑎𝑙と可換であれば𝐻 + 𝐻𝑍𝑒𝑒𝑚𝑎𝑛 も可換であり、𝐻𝑍𝑒𝑒𝑚𝑎𝑛 を加えても各固有状態は不変であるが、固有値のみが磁場によって以下の様に変化する；𝑧𝑧𝑧
𝐸0 (𝑆𝑡𝑜𝑡𝑎𝑙, ℎ) = 𝐸0 (𝑆𝑡𝑜𝑡𝑎𝑙, ℎ = 0) − ℎ𝑆𝑡𝑜𝑡𝑎𝑙

(9)

𝑧
𝑧
部分空間の最低エネルギーである。磁化曲線は各磁場に対して最低のエ
𝐸0 (𝑆𝑡𝑜𝑡𝑎𝑙, ℎ = 0)は各𝑆𝑡𝑜𝑡𝑎𝑙𝑧𝑧
ネルギーを持つ状態の磁化の値を示したものであり、𝐸0 (𝑆𝑡𝑜𝑡𝑎𝑙, ℎ) = 𝐸0 (𝑆𝑡𝑜𝑡𝑎𝑙+ 1, ℎ)を満たすℎで𝑧𝑧𝑧から𝑆𝑡𝑜𝑡𝑎𝑙+ 1に変化するようプロットすればよい。つまり各𝑆𝑡𝑜𝑡𝑎𝑙
部分空間での最低固磁化が𝑆𝑡𝑜𝑡𝑎𝑙

有値が分かっていれば、わざわざゼーマン項を加えた計算を行わずとも磁化曲線を得られるのである。
𝑧
DMRG 法で𝑆𝑡𝑜𝑡𝑎𝑙が大きい値での最低固有値を計算するとき１つ問題が生じる。厳密対角化で𝑧は計算したいサイト数𝐿𝑙𝑎𝑠𝑡 及び𝑆𝑡𝑜𝑡𝑎𝑙(𝐿𝑙𝑎𝑠𝑡 )での計算を初めから行うことができる。しかし DMRG

無限系法ではサイト数は徐々に増やしてゆくため、0 や 1 などの小さな値でない限り計算を進め𝑧
るなかで𝑆𝑡𝑜𝑡𝑎𝑙の値を増加させてゆく必要がある。この論文では DMRG 無限系法の途中での𝑧𝑆𝑡𝑜𝑡𝑎𝑙(𝐿)は次の様に設定した；𝑧𝑧(𝐿) = ceil[𝑆𝑡𝑜𝑡𝑎𝑙(𝐿𝑙𝑎𝑠𝑡 )𝑆𝑡𝑜𝑡𝑎𝑙

𝐿
𝐿𝑙𝑎𝑠𝑡

]

(10)

ここで ceil は切り上げである。このような設定にした理由は、異なるサイト数でも磁化𝑀 =𝑧⁄𝐿が同じ値ならば同じ性質を持ったエネルギー固有状態となり、DMRG 計算の途中での保𝑆𝑡𝑜𝑡𝑎𝑙

持基底を決定するのに妥当な状態と考えられるためである。
3

《寄稿論文》

物性研究・電子版 Vol. 5, No. 2, 052101 （2016年5月号）

2.2

無限系法初期推定

DMRG 無限系法では小さな系の固有状態から決定した保持基底を用いて大きな系の状態を精度良く求めている。これを繰り込み群の観点から見ると、系が大きくなるにつれパラメータ空間
上で熱力学的極限の固定点へ近づいていると解釈される。つまり系が大きくなると繰り込み変換𝑧で残される基底はほとんど同じ𝑆𝑏𝑙𝑜𝑐𝑘と密度行列固有値を持つようになり、またこれらの基底を用いて計算される全系の固有状態の各成分も前のステップでの固有状態の各成分と近い値になる。
この考えを元にして DMRG 無限系法の初期推定が考案された[5]。
無限系法初期推定を行うには、まず基底の対応付けが必要になる。先に述べたように、ハミル𝑧𝑧
トニアンが𝑆𝑡𝑜𝑡𝑎𝑙と交換するならば、密度行列固有状態(=保持基底)も一定の𝑆𝑏𝑙𝑜𝑐𝑘を持つ。これら𝑧𝑧を各𝑆𝑏𝑙𝑜𝑐𝑘で固有値の大きさ順に並べ替えたとき、𝑛番目に大きい基底を|𝑆𝑏𝑙𝑜𝑐𝑘, 𝑛⟩𝐿 と表す。この基𝑧𝑧底に対し、同じ𝑆𝑏𝑙𝑜𝑐𝑘, 𝑛でサイト数が𝐿′ のときのステップでの基底|𝑆𝑏𝑙𝑜𝑐𝑘, 𝑛⟩𝐿′ を対応付ける。そして𝐿′での全系の固有状態|𝜓⟩𝐿′ =∑

𝜓𝑖1𝑖2𝑖3𝑖4 |𝑖1⟩𝐿′ |𝑖2⟩|𝑖3⟩|𝑖4⟩𝐿′(11)

𝑖1,𝑖2,𝑖3,𝑖4

より、この𝜓𝑖1𝑖2𝑖3𝑖4を用いて初期推定状態|𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩𝐿|𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩𝐿 =∑𝜓𝑖1𝑖2𝑖3𝑖4 |𝑖1⟩𝐿 |𝑖2⟩|𝑖3⟩|𝑖4⟩𝐿(12)

𝑖1,𝑖2,𝑖3,𝑖4を設定する。ここで|𝑖1⟩𝐿(𝐿′ ) , |𝑖4⟩𝐿(𝐿′ ) はそれぞれサイト数𝐿 (𝐿′ )での system 及び environment のブ
ロックの基底であり、|𝑖1⟩𝐿′ と|𝑖1⟩𝐿 、|𝑖4⟩𝐿′ と|𝑖4⟩𝐿 は対応している。また|𝑖2⟩, |𝑖3⟩は system 及び𝑧environment のスピンの基底である。なお𝑆𝑏𝑙𝑜𝑐𝑘となるブロックの基底の数が一致しないことがあるが、サイト数𝐿′ での基底が余るならサイト数𝐿での基底の数に合わせて𝑛の大きいものを打ち切り、逆に足りないならば対応付けできない初期推定状態の成分は 0 と設定する。また基底の対応𝑧付けには𝑆𝑏𝑙𝑜𝑐𝑘の一致が必要なため、半整数スピン鎖の場合はブロック部分のサイト数の偶奇性が𝑧一致する𝐿′ = 𝐿 − 4に設定するが、整数スピンではこの偶奇性によらず𝑆𝑏𝑙𝑜𝑐𝑘を一致させることができるので𝐿′ = 𝐿 − 2とできる。
以上のようにして初期推定状態の元を作ることができるが、このままでは初期推定はできない。
なぜならブロックの基底は密度行列の固有状態であり、符号の任意性が残っているためである。
これにより初期推定状態の各成分の絶対値は真の固有状態と一致しているものの、符号がランダムに変わってしまっている。Ref. 5 では１次元反強磁性鎖で成り立つ Marshall’s sign rule を利用した符号調整を行い、𝑆 = 2ハイゼンベルグ鎖の基底状態と第一励起状態の初期推定を行っている。しかし Marshall’s sign rule は特定の場合でのみ成り立つ法則であり、例えばフラストレーションがある場合には厳密には成立しない。また Marshall’s sign rule の表式からスピンが半整数の場合には虚数が現れてしまうのでそのままでは利用できない。Ref. 6 ではハミルトニアンや上昇演算子の行列要素から符号調整を行っている。この方法はスピンやモデルによらず適用できるが、調整がうまくいかないことが多い。しかし本論文で提案する改良方法はある程度の符号調整がなされているほうが好ましく、この方法と共に用いることでほぼ完璧に調整ができるようになる。そこで Ref. 6 の演算子からの符号調整方法について説明しておく。

4

《寄稿論文》

物性研究・電子版 Vol. 5, No. 2, 052101 （2016年5月号）

図２ ブロック部分のハミルトニアン𝐻𝑏 (左)とブロック右端スピンの上昇演算子𝑆𝑟+ (右)。ハミルトニアンは𝑧𝑧
𝑆𝑏𝑙𝑜𝑐𝑘一定の部分空間でのみ値を持ち、上昇演算子は𝑆𝑏𝑙𝑜𝑐𝑘が 1 だけ大きい基底との要素のみ値を持つ。

𝑧𝑧
ハミルトニアン𝐻が𝑆𝑡𝑜𝑡𝑎𝑙と交換するとき、ブロック部分の基底も一定の𝑆𝑏𝑙𝑜𝑐𝑘を持つ。この基底
𝑧を𝑆𝑏𝑙𝑜𝑐𝑘
ごとに密度行列の大きさ順に並べ替えると、ブロック部分のハミルトニアン𝐻𝑏 とブロック𝑧
右端スピンの上昇演算子𝑆𝑟+ は図２のようになる。𝐻𝑏 は𝑆𝑏𝑙𝑜𝑐𝑘が同じ部分空間のみでのみ値を持ち𝑧𝑧(図２左)、またブロック右端スピンの上昇演算子𝑆𝑟+は𝑆𝑏𝑙𝑜𝑐𝑘と𝑆𝑏𝑙𝑜𝑐𝑘+1 の行列要素のみ値をもつ(図𝑧𝑧２ 左 ) 。 こ の 性 質 を 利 用 し て 、 同 じ 𝑆𝑏𝑙𝑜𝑐𝑘
で は 、 最 大 固 有 値 を 持 つ 基 底 |𝑆𝑏𝑙𝑜𝑐𝑘, 1⟩ に 対 し𝑧𝑧𝑧𝑧⟨𝑆𝑏𝑙𝑜𝑐𝑘, 1|𝐻𝑏 |𝑆𝑏𝑙𝑜𝑐𝑘, 𝑛⟩ (𝑛 > 1)が正になるように|𝑆𝑏𝑙𝑜𝑐𝑘, 𝑛⟩の符号を設定し、また異なる𝑆𝑏𝑙𝑜𝑐𝑘
では、
𝑧
𝑧
𝑧
⟨𝑆𝑏𝑙𝑜𝑐𝑘, 1|𝑆𝑟+ |𝑆𝑏𝑙𝑜𝑐𝑘+ 1,1⟩が正になるよう|𝑆𝑏𝑙𝑜𝑐𝑘+ 1,1⟩の符号を調整する。図２において○で表されている行列要素が調整に利用されるものである。全基底の符号反転を除けばこの計算によって一意に符号を決定できる。
この調整の妥当性について原論文では詳しく触れられていないが、定性的には以下の様に説明できる；基底の密度行列固有値が大きいほど対応するブロックの演算子の行列要素は大きな値を持ち、計算される全系の固有状態に大きく影響する傾向があるだろう。そのような固有状態に大きく影響する行列要素の符号を固定することで符号調整を行うことができる。
しかしこの計算がうまくいく場合は少ない。図３に𝑚 = 150での𝑆 = 1/2ハイゼンベルグ鎖の基𝑧𝑧
底状態(𝑆𝑡𝑜𝑡𝑎𝑙= 0最低エネルギー状態)と𝑆 = 1ハイゼンベルグ鎖の第一励起状態(𝑆𝑡𝑜𝑡𝑎𝑙= 1)の計算例を示す。横軸はスピン鎖の長さ𝐿、縦軸は初期推定状態のエネルギー期待値𝐸𝑖𝑔 と計算された基底エネルギー𝐸0 の相対誤差の常用対数log10(|(𝐸0 − 𝐸𝑖𝑔 )⁄𝐸0 |)である。白丸は上記の手順で得られた初期推定状態の相対誤差であり、青線は符号調整が限界まで行われた場合の初期推定状態の相対誤差である。この限界符号調整状態|𝜓𝑚𝑎𝑥 ⟩は初期推定状態と各成分の絶対値は同じであるが、符号を基底状態の成分のものと同じになるよう変えたものである。つまり真の基底状態|𝜓0 ⟩ =∑𝑖 𝜓𝑖 |𝑒𝑖 ⟩と初期推定状態|𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩ = ∑𝑖 𝜓𝑖′ |𝑒𝑖 ⟩に対し、|𝜓𝑚𝑎𝑥 ⟩ = ∑𝑖 𝜓𝑖𝑚𝑎𝑥 |𝑒𝑖 ⟩は|𝜓𝑖𝑚𝑎𝑥 | = |𝜓𝑖′ |𝑠𝑔𝑛(𝜓𝑖𝑚𝑎𝑥 ) = 𝑠𝑔𝑛(𝜓𝑖 )(13)を満たす。符号の任意性を考えたとき、もし完璧に符号調整ができているならば⟨𝜓0 |𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩~1であり、ほぼすべての𝑖に対し𝜓𝑖 𝜓𝑖′ ≥ 0が期待される。ゆえにこの|𝜓𝑚𝑎𝑥 ⟩の相対誤差と比較することでどれだけ符号調整ができていたか見積もることができる。
5

《寄稿論文》

物性研究・電子版 Vol. 5, No. 2, 052101 （2016年5月号）

図３ 𝑆 = 1/2反強磁性ハイゼンベルグ鎖の基底状態計算における演算子からの符号調整での初期推定誤差。
白丸は実際の初期推定状態のエネルギーの相対誤差であり、青線は限界符号調整状態のエネルギーの相対誤差である。どちらのグラフも白丸と青線は一致しておらず、符号調整は失敗している。

グラフからこの初期推定方法の不完全さは明らかである。𝑆 = 1/2の場合(図３左)は約半分程の計算では初期推定がなされているが、そうでない点も多々あり値がばらついている。このとき各演算子の符号と計算される基底状態の間に単純な関係は存在していないと考えられる。また𝑆 = 1の場合(図３右)は相対誤差の値は悪いながらもほぼ一定の値になっている。このとき演算子の符号と基底状態の間には比較的単純な関係が成り立っているものの、単に演算子成分を正にするような関係ではないと思われる。もちろんモデルや計算したい状態によってはこのアルゴリズムでも完璧に符号調整ができることもあるかもしれないが、毎回符号の関係を調べるのは手間である
し簡単な関係とも限らない。ゆえに符号法則によらず適用可能な符号調整方法が求められる。

2.3

初期推定符号調整の改良

一般的に適用可能な符号調整として、エネルギー期待値を利用する方法を考えた。初期推定状(𝑗)

態 |𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩ = ∑𝑖 𝜓𝑖′ |𝑒𝑖 ⟩ 及 び そ の 第 𝑗 成 分 を 反 転 さ せ た 状 態 |𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩ = ∑𝑖 𝜓𝑖′ |𝑒𝑖 ⟩ − 2𝜓𝑗′ |𝑒𝑗 ⟩ =(𝑗)|𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩ − 2𝜓𝑗′ |𝑒𝑗 ⟩と、真の基底状態|𝜓0 ⟩との内積⟨𝜓0|𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩, ⟨𝜓0 |𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩を考えると、符号が正し(𝑗)いものの方が 1 に近くなる。例えば、第𝑗成分の符号が|𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩の方が正しいならば(𝑗)|⟨𝜓0|𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩| < |⟨𝜓0 |𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩| < 1(14)となる。ここで|𝜓0 ⟩が基底状態であることを考えると、これらのエネルギー期待値の間に
(𝑗)

(𝑗)

⟨𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 |𝐻|𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩ > ⟨ 𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 |𝐻|𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩(15)という関係が成り立つと期待される。もちろん厳密な式ではないが、変分法的に考えても(15)式の6

《寄稿論文》

物性研究・電子版 Vol. 5, No. 2, 052101 （2016年5月号）(𝑗)

関係が成り立つとき、|𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩はより基底状態に近いだろう。
この考えを元にした改良方法の具体的な手順は次のようになる；(𝑗)

１ 初期推定状態|𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩と第𝑗成分が反転した状態|𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩のエネルギー期待値の差𝑑𝐸 =
(𝑗)

(𝑗)

⟨ 𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 |𝐻|𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩ − ⟨𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 |𝐻|𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩を計算
(𝑗)

２ 𝑑𝐸が負なら|𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩を新たな初期推定状態に採用
３ この計算を全成分に対し繰り返す
４ １～３の反転操作を 1 ループとし、ループ間で状態が変わらなくなるまで繰り返す
ここでいくつか注意事項がある。まずエネルギー差𝑑𝐸について、この値を毎回内積計算で求め
(𝑗)

(𝑗)

(𝑗)

る必要はない。|𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩ = |𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩ − 2𝜓𝑗′ |𝑒𝑗 ⟩となることに注意すれば、⟨ 𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 |𝐻|𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩ =⟨𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 |𝐻|𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩ − 4𝜓𝑗′ ⟨𝑒𝑗 |𝐻|𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩ + 4(𝜓𝑗′ )2⟨𝑒𝑗 |𝐻|𝑒𝑗 ⟩より
(𝑗)

(𝑗)

𝑑𝐸 = ⟨ 𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 |𝐻|𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩ − ⟨𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 |𝐻|𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩= −4𝜓𝑗′ ⟨𝑒𝑗 |𝐻|𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩ + 4(𝜓𝑗′ )2⟨𝑒𝑗 |𝐻|𝑒𝑗 ⟩(16)となる。あらかじめ𝐻|𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩を計算しておけば⟨𝑒𝑗 |𝐻|𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩はその第𝑗成分を呼び出すだけで済む。そして(16)式の𝑑𝐸が負になった時のみ反転を実行し、𝐻|𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩も対応したものに変化させる。
(𝑗)

反転後の𝐻|𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩は(𝑗)

𝐻|𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩ = 𝐻|𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩ − 2𝜓𝑗′ 𝐻|𝑒𝑗 ⟩(17)であるので、ハミルトニアン行列の第𝑗列目のみを呼び出し、それを2𝜓𝑗′ 倍したものを𝐻|𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩から引けばよい。これらのことに注意すれば、1 ループあたりの計算量はハミルトニアン行列と状態ベクトルの積の計算１回分以下となる。
またループの終了条件も重要になる。厳密に状態が変わらなくなるまでループ計算を繰り返していては計算時間がかかりすぎてしまう。初期推定状態はあくまで前ステップからの推定であり、数値計算誤差や有限サイズ効果により厳密には真の基底状態に一致しない。ゆえにある程度符号調整を行い、ループ間であまり状態が変わらなくなれば処理を終えてよい。具体的には以下のど
ちらかの条件を満たした時ループを終了するよう設定した；
1 − ⟨𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 (𝑙)|𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 (𝑙 − 1)⟩ < 𝑒𝑠(18)|(〈𝐸𝑖𝑛𝑖𝑡𝑖𝑎𝑙 (𝑙)〉 − 〈𝐸𝑖𝑛𝑖𝑡𝑖𝑎𝑙 (𝑙 − 1)〉)/〈𝐸𝑖𝑛𝑖𝑡𝑖𝑎𝑙 (𝑙)〉| < 𝑒𝑠(19)

ここで𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 (𝑙)は𝑙回目のループ終了時の初期推定状態であり、〈𝐸𝑖𝑛𝑖𝑡𝑖𝑎𝑙 (𝑙)〉はそのエネルギー期待値、また𝑒𝑠 は符号調整の設定誤差である。この𝑒𝑠 を適切な値に設定することでほとんど計算コス
トをかけずに符号調整を行うことができる。(18)式の条件はループ間で状態が変わらない場合の式であるが、(19)式はあまり基底状態にかかわらない基底まで保持されている場合を考えての条件式である。密度行列固有値の大きな基底は対応する初期推定状態の成分も大きく、符号反転によってエネルギーも大きく変化する。しかし固有値の小さな基底に対応する成分は反転エネルギ
ー差が小さくループを繰り返してもなかなか符号が決定されないため、(18)式の条件に引っかかりづらい。そのような場合でもエネルギー期待値はほとんど同じ値になっており、状態だけでなく(19)式のようにエネルギー期待値についても条件を課すことで効率的に計算できる。

7

《寄稿論文》

物性研究・電子版 Vol. 5, No. 2, 052101 （2016年5月号）

図４

𝑆 = 1/2ハイゼンベルグ鎖の基底状態計算における１成分毎の期待値調整。期待値調整のみでは初期推定は不完全だが(左)、演算子からの調整と併用することで完璧に調整される(右)。

図３で見たように演算子からの調整は不完全であるが、期待値からの調整も単体ではうまくいかない。図４は𝑚=150 での𝑆 = 1/2ハイゼンベルグ鎖の基底状態の初期推定計算における演算子
からの調整の適用・非適用時の比較である。これからわかるように、期待値からの符号調整のみ
でも初期推定ができていない点が散見される(図４左)。これは各成分を反転させ状態を変えてい
く中で、基底状態から離れているがどの１成分を反転させてもそれ以上エネルギーが下がらない
極小点に陥ってしまったためと思われる。一方、両方の調整を共に適用した場合(図４右)は、すべての点で限界符号調整状態の誤差と一致しており完璧に符号調整ができていることがわかる。演
算子調整のみ、あるいは期待値調整のみでの計算で初期推定ができていない点は一致してはいない。つまり両者の方法を同時に適用することでお互いに調整ができない場合を補い合い、すべての点で完璧に符号調整ができたといえる。なお誤差が大きくなっている点がいくつか見られるが、符号調整とは別の原因によるものであり、これについては 2.5 で詳しく述べる。
このように演算子からの調整と期待値からの調整を併用すれば完璧に符号調整ができるが、極まれに失敗してしまうこともある。事実、我々が行った計算のうち 4.2 の𝑆 = 1ハイゼンベルグ鎖の高磁場領域の計算において 1,2 点そのような場合が存在した。もちろんほとんどの点では失敗していないため、1,2 点符号調整ができておらずとも全体の計算時間にはほとんど影響しない。しかしこの改良を施しても初期推定を失敗することがあるという実証であり、他のモデルでは頻繁に起こることもあるかもしれない。そこでさらに安定的に符号調整ができる方法を考えた。
上記の改良方法は初期推定状態のある１成分のみに着目しエネルギー期待値から符号調整を行っている。しかし問題なのは初期推定状態成分の符号ではなくブロックの基底の符号であるため、エネルギー期待値からの符号調査もブロックの基底毎に行う方が趣旨に沿っている。ただし、そ8《寄稿論文》

物性研究・電子版 Vol. 5, No. 2, 052101 （2016年5月号）のような計算を行うには１度に初期推定状態の複数の成分を反転させてエネルギー差を調べる必要がある。その場合のアルゴリズムについて考えよう。あるブロック基底|𝑗⟩𝑏 を反転させたとき符号が反転する全系の基底の集合を𝐴𝑖 をする。つまり、|𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩ = ∑𝑖 𝜓𝑖′ |𝑒𝑖 ⟩に対し|𝑗⟩𝑏 を反転させた
(𝑗)

𝑏
⟩は
状態|𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙(𝑗)𝑏⟩ = |𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩ − ∑ 2𝜓𝑖′ |𝑒𝑖 ⟩|𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙

(20)

𝑖∈𝐴𝑗
(𝑗)

𝑏
⟩のエネルギー期待値はとなる。|𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙

(𝑗)

(𝑗)

⟨ 𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 |𝐻|𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩ = ⟨𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 |𝐻|𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩ − 4 ∑ 𝜓𝑖′ ⟨𝑒𝑖 |𝐻|𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩ + 4 ∑ 𝜓𝑖′ ⟨𝑒𝑖 |𝐻|𝑒𝑘 ⟩𝑖∈𝐴𝑗𝑖,𝑘∈𝐴𝑗= ⟨𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 |𝐻|𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩ − 4 ∑ 𝜓𝑖′ ⟨𝑒𝑖 |𝐻|𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩ + 4 ∑ 𝜓𝑖′ ⟨𝑒𝑖 |𝐻|𝑒𝑖 ⟩ + 8 ∑ 𝜓𝑖′ ⟨𝑒𝑖 |𝐻|𝑒𝑘 ⟩
𝑖∈𝐴𝑗

𝑖∈𝐴𝑗

(21)

𝑖<𝑘∈𝐴𝑗となる。１行目から２行目にはハミルトニアンのエルミート性を利用した。上式の２段目４項にハミルトニアンの非対角成分が現れるが、ハミルトニアンの疎行列性を利用することでこの計算量を減らすことができる[3]。これよりエネルギー差𝑑𝐸は
𝑑𝐸 = −4 ∑ 𝜓𝑖′ ⟨𝑒𝑖 |𝐻|𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩ + 4 ∑ 𝜓𝑖′ ⟨𝑒𝑖 |𝐻|𝑒𝑖 ⟩ + 8 ∑ 𝜓𝑖′ ⟨𝑒𝑖 |𝐻|𝑒𝑘 ⟩
𝑖∈𝐴𝑗

𝑖∈𝐴𝑗

𝑖<𝑘∈𝐴𝑗

(𝑗)

(22)

(𝑗)

𝑏
⟩を新たな初期推定状態に採用し、𝐻|𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩の計算を行う；となる。𝑑𝐸 < 0のとき|𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙(𝑗)

𝐻|𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩ = 𝐻|𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩ − ∑ 2𝜓𝑖′ 𝐻|𝑒𝑖 ⟩
𝑖∈𝐴𝑗

(23)

以上の計算をすべてのブロック基底に対して行い、それを１ループとして(17)(18)式の条件に引っかかるまで繰り返せばよい。
１成分毎の計算と比べると１度に多くの成分を変化させており、また反転の自由度が全系の状態数から保持基底数𝑚に減少しているため、極小点に陥りづらくなると期待される。図５は𝑚=150
での𝑆 = 1/2ハイゼンベルグ鎖の基底状態計算において、演算子からの調整を行わずにこちらの方法のみで初期推定を行った場合のグラフである。１成分毎の場合とは異なり、演算子からの調整を行わずとも完璧に符号調整がされている。このようにブロック基底毎に反転調査を行うことで高い精度で符号調整ができるが、演算子からの符号調整を共に適用することでさらに安定するだろう。またブロック基底毎の符号調整は１成分毎の場合より計算量が多いため、これ単体で調整が可能な場合でも演算子からの符号調整をしておくことで計算時間の短縮につながる。完璧に符号調整ができる３つの初期推定計算
Ⅰ演算子からの調整＋１成分毎の調整

Ⅱ

ブロック基底毎の調整

Ⅲ演算子からの調整＋ブロック基底毎の調整での計算時間はおおよそⅡ > Ⅲ ≥ Ⅰとなった。Ⅰの調整で初期推定に失敗することがよくある場
合ではⅡ, Ⅲの方が速くなることもあるかもしれないが、今回行った計算ではほとんどがⅠで調整できたため、特別の断りがない限り初期推定とはⅠの計算を指していると注意されたい。

9

《寄稿論文》

物性研究・電子版 Vol. 5, No. 2, 052101 （2016年5月号）

図５

𝑆 = 1/2ハイゼンベルグ鎖の基底状態計算におけるブロック基底毎の調整による初期推定誤差。演算子からの調整を行わずとも完璧に符号調整がされる。

2.4

近似推定による初期推定状態の改良

ここまで初期推定の符号問題を取り扱ってきたが、ここでは符号調整を行った上で初期推定状
態成分の絶対値を変化させ、より基底状態に近い状態を作ることを考える。3 章で示すが、初期推定誤差は有限サイズ効果と同じ𝐿依存性を持ち𝐿 → ∞(1⁄𝐿 → 0)で0になる。つまり真の基底状
態|𝜓0 ⟩𝐿 と初期推定状態|𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩𝐿 の差|∆𝜓1 ⟩𝐿 は𝐿−1の関数で表され、𝐿−1が十分小さければ線形近似が成り立つと考えられる；|∆𝜓1 ⟩𝐿 = |𝜓0 ⟩𝐿 − |𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩𝐿 ≅|𝜓0 ⟩𝐿′ − |𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩𝐿′
𝐿′ −1 − 𝐿′′ −1(𝐿−1 − 𝐿′−1)=|∆𝜓1 ⟩𝐿′∆𝐿′ −1∆𝐿−1(24)

ここで|𝜓0 ⟩𝐿′ は|𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩𝐿 の元になるサイト数𝐿′での基底状態であり、∆𝐿−1は初期推定元の状態のサ−1−1

イト数との逆数の差𝐿−1 − 𝐿′ 、∆𝐿′ は𝐿′での推定元状態のサイト数𝐿′′ との逆数の差𝐿′ある。この式から改良初期推定状態|𝜓𝑖𝑚𝑝1 ⟩𝐿|𝜓𝑖𝑚𝑝1 ⟩𝐿 = |𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩𝐿 +|∆𝜓1 ⟩𝐿′∆𝐿′ −1

∆𝐿−1

−1

− 𝐿′′

−1

で

(25)

を考えると、|𝜓0 ⟩𝐿 のテイラー展開と∆𝐿−1の１次まで一致する。つまり|𝜓𝑖𝑚𝑝1 ⟩𝐿 は∆𝐿−1の１次近似であり、この観点では|𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩𝐿 は０次近似であるといえる。また|𝜓0 ⟩𝐿 と|𝜓𝑖𝑚𝑝1 ⟩ の差|∆𝜓2 ⟩𝐿𝐿|∆𝜓2 ⟩𝐿 = |𝜓0 ⟩𝐿 − |𝜓𝑖𝑚𝑝1 ⟩𝐿(26)を考え、これの線形近似からさらに改良した|𝜓𝑖𝑚𝑝2 ⟩𝐿|𝜓𝑖𝑚𝑝2 ⟩ = |𝜓𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ⟩𝐿 +𝐿|∆𝜓1 ⟩𝐿′
′ −1

∆𝐿
10

∆𝐿−1 +|∆𝜓2 ⟩𝐿′∆𝐿′ −1∆𝐿−1(27)

《寄稿論文》

物性研究・電子版 Vol. 5, No. 2, 052101 （2016年5月号）

を作ることもできる。このとき|𝜓0 ⟩𝐿 と|𝜓𝑖𝑚𝑝2 ⟩𝐿 は𝐿が十分大きければ∆𝐿−1の２次まで一致する。同様の手順で次々に高次の改良初期推定状態を作ることができるが、ある程度のオーダーになると数値計算誤差の影響で推定精度は上がらなくなるだろう。
図６は𝑚=150 での𝑆 = 1/2ハイゼンベルグ鎖の基底状態計算における改良初期推定誤差である。
非改良状態（０次推定）に比べ改良推定（１次推定）は明らかに誤差が小さくなった。切り捨て
誤差は10−10程であったが、同じオーダー程の初期推定ができている。多少値がばらついているのは何らかの原因で非改良推定状態や基底状態に誤差が出ると改良初期推定状態もそれを引き継いで誤差が生じるためである。このばらつきの原因やその改良については 2.5 で述べる。１次推定の時点で値がばらつき𝐿−1にあまり依存しなくなっているため、２次推定の計算では誤差はさらに大きくばらつきあまり計算量は減らなかった。系によっては０次推定の時点で値がばらつくこともあるだろう。その場合は１次推定を行っても改良されないと思われる。

図６ 𝑆 = 1/2ハイゼンベルグ鎖の基底状態計算における近似改良初期推定誤差。非改良初期推定（白丸）に対し、１次推定（四角）は誤差が減っているが、２次推定（三角）は大きくばらついている。

11

《寄稿論文》

物性研究・電子版 Vol. 5, No. 2, 052101 （2016年5月号）

2.5

基底対応付けの改良

図４,５で符号調整は完璧にされているのにもかかわらず、誤差が大きくなる点がいくつかある。
この異常性の原因は単なる数値計算誤差ではない。もし数値計算誤差により異常な基底状態や密度行列固有状態が計算されたとすれば、その状態から初期推定を行う点でも誤差は大きくなるはずである。しかし実際には異常点は独立しており、その前後で初期推定誤差の𝐿依存性はなめらかに変化している。詳しく調べてみると初期推定での基底の対応付けが問題であるとわかった。2.2
で説明したように初期推定は保持基底及びエネルギー固有状態が繰り込みの固定点へ向かうことを利用しており、密度行列の各固有値も一定値に収束してゆく。しかし収束値や収束速度はそれぞれ異なり、徐々に系を大きくしてゆく中で固有値の大きさが入れ替わることがある。すると初期推定での基底の対応付けは密度行列固有値の大きさ順で行うため、本来対応していない基底同士が対応付けられてしまいこのような異常性が現れたのである。
この異常性を解消するには基底の対応付けを修正する必要がある。そこで対応付ける基底間で𝑧
内積計算を行うことを考えた。保持基底|𝑆𝑏𝑙𝑜𝑐𝑘, 𝑛⟩𝐿 は密度行列固有状態であり１ステップ前のブロック基底|𝑖1′ ⟩𝐿−2と system のスピン基底|𝑖2′ ⟩を用いて(𝑖1,𝐿)

𝑧|𝑆𝑏𝑙𝑜𝑐𝑘, 𝑛⟩𝐿 = |𝑖1⟩𝐿 = ∑ 𝐶𝑖1′ 𝑖2′ |𝑖1′ ⟩𝐿−2 |𝑖2′ ⟩(28)

𝑖1′ ,𝑖2′𝑧と表現される。同じく対応付け候補の基底|𝑆𝑏𝑙𝑜𝑐𝑘, 𝑛⟩𝐿′ も(𝑖1,𝐿′ )𝑧|𝑆𝑏𝑙𝑜𝑐𝑘, 𝑛⟩𝐿′ = |𝑖1⟩𝐿′ = ∑ 𝐶𝑖1′ 𝑖2′ |𝑖1′ ⟩𝐿′ −2|𝑖2′ ⟩(29)

𝑖1′ ,𝑖2′′(𝑖1,𝐿)(𝑖1,𝐿 )と表される。そして各成分の絶対値|𝐶𝑖1′ 𝑖2′ |, |𝐶𝑖1′ 𝑖2′ |をその大きさ順に並べ替え、その内積(𝑖1,𝐿)𝑧𝑧
𝑃(|𝑆𝑏𝑙𝑜𝑐𝑘, 𝑛⟩𝐿 , |𝑆𝑏𝑙𝑜𝑐𝑘, 𝑛⟩𝐿′ ) = ∑ |𝐶(𝑖)(𝑖1,𝐿′ )| |𝐶(𝑖)|(30)

𝑖(𝑖1,𝐿)(𝑖1,𝐿)
𝑧を計算した。ここで 𝐶(𝑖)は 𝐶𝑖1′𝑖2′ のうち絶対値が𝑖 番目に大きい成分である。 もし|𝑆𝑏𝑙𝑜𝑐𝑘, 𝑛⟩𝐿 と
(𝑖1,𝐿)(𝑖1,𝐿′)

𝑧|𝑆𝑏𝑙𝑜𝑐𝑘, 𝑛⟩𝐿′ が対応付けるべき基底同士なら、𝐶𝑖1′𝑖2′ , 𝐶𝑖1′𝑖2′ の絶対値の分布は一致するはずなので、(30)式の内積をとれば非常に 1 に近くなると予想される。実際に計算してみたところ、絶対値順
内積の値は対応付けるべき基底間では 0.9999 程、そうでないときには 0.8~0.9 程になり、対応付
けるべき基底かどうかを判定することができた。
以上の計算を基にした基底対応付けの改良手順は次のようになる；
１ 非改良時と同じく保持基底を固有値の大きさ順に対応付け
𝑧
𝑧
𝑧
𝑧
𝑧
２ |𝑆𝑏𝑙𝑜𝑐𝑘, 𝑛⟩𝐿 に対し|𝑆𝑏𝑙𝑜𝑐𝑘, 𝑛⟩𝐿′ 及び|𝑆𝑏𝑙𝑜𝑐𝑘, 𝑛 + 1⟩𝐿′ との絶対値順内積𝑃(|𝑆𝑏𝑙𝑜𝑐𝑘, 𝑛⟩𝐿 , |𝑆𝑏𝑙𝑜𝑐𝑘, 𝑛⟩𝐿′ ),𝑧𝑧
𝑃(|𝑆𝑏𝑙𝑜𝑐𝑘, 𝑛⟩𝐿 , |𝑆𝑏𝑙𝑜𝑐𝑘, 𝑛 + 1⟩𝐿′ )を計算
𝑧𝑧𝑧𝑧
３ もし𝑃(|𝑆𝑏𝑙𝑜𝑐𝑘, 𝑛⟩𝐿 , |𝑆𝑏𝑙𝑜𝑐𝑘, 𝑛 + 1⟩𝐿′ )の方が１に近ければ|𝑆𝑏𝑙𝑜𝑐𝑘, 𝑛⟩𝐿 は|𝑆𝑏𝑙𝑜𝑐𝑘, 𝑛 + 1⟩𝐿′ と対応付け、𝑧𝑧|𝑆𝑏𝑙𝑜𝑐𝑘, 𝑛 + 1⟩𝐿 は|𝑆𝑏𝑙𝑜𝑐𝑘, 𝑛⟩𝐿′ と対応付けるように修正

４ 以上の修正を𝑛 = 1から順に実行𝑧𝑧
この計算で対応付けの修正候補基底は|𝑆𝑏𝑙𝑜𝑐𝑘, 𝑛 + 1⟩𝐿′ のみだが、系によっては|𝑆𝑏𝑙𝑜𝑐𝑘, 𝑛 + 2⟩𝐿′ と対応付ける必要が出てくることも考えられる。しかし今回行った計算では上記の改良対応付けによってすべての場合で修正できた。
12

《寄稿論文》

物性研究・電子版 Vol. 5, No. 2, 052101 （2016年5月号）

図７

𝑆 = 1/2ハイゼンベルグ鎖の基底状態計算における対応付け改良初期推定誤差。

すべての計算で誤差のばらつきが改善されている。

図７は基底対応付け改良を施した場合の𝑚=150 での𝑆 = 1/2ハイゼンベルグ鎖の基底状態計算
における初期推定誤差である。図６がこの改良を行わない場合であるのでそちらも参照されたい。
０次近似では誤差が大きくなる点が全て消え、確かに改良することができた。また０次近似で完璧に初期推定できているため１次近似でもばらつきが抑えられ、その結果２次近似でさらに誤差を減らすことができた。しかし𝐿~170で誤差が大きくなり近似推定があまりうまくいっていない点が存在する。これは系を大きくしていったときの状態の変化が原因であると思われる。誤差が𝑧大きくなる点では基底状態計算での𝑆𝑡𝑜𝑡𝑎𝑙
一定の状態数が変化していて、また初期推定時の(30)式の内積計算で突出して１に近い値を示す対応基底が存在しない基底が現れた。密度行列固有値の大きさの入れ替えが保持されていた基底と切り捨てられていた基底の間で起こったり、あるいは系を大きくしていくとエネルギー固有状態の𝐿依存性が変化したりする場合、初期推定状態と真の状態の差を𝐿−1 の線形近似で表すことはできないため、誤差が大きくなったと考えられる。このような異常性は DMRG 法の特性と系の性質に由来するため、単なる初期推定計算の工夫では改善は難しいだろう。
なお密度行列固有値の入れ替わりは𝑆 = 1/2ハイゼンベルグ鎖の基底状態では頻繁に起こるが、第一励起や𝑆 = 1の場合にはあまり起こらず上記の修正を行ってもほとんど改良されない。ゆえに以下の計算で基底対応付けの改良を行っているのは 3.1 の基底状態計算のみである。

13

《寄稿論文》

物性研究・電子版 Vol. 5, No. 2, 052101 （2016年5月号）

2.6

高磁場領域計算時の初期推定𝑧2.1 で述べたように DMRG 法によって高磁場領域の基底状態、つまり大きな𝑆𝑡𝑜𝑡𝑎𝑙
での最低エ𝑧
ネルギー状態を計算したい場合には、(10)式の様にサイト数を増やしてゆく中で𝑆𝑡𝑜𝑡𝑎𝑙も大きくし𝑧ていく必要があるが、そのような場合に初期推定は行えないだろうか。計算途中で𝑆𝑡𝑜𝑡𝑎𝑙が変わるということは保持基底を決めるためのターゲット状態が変化するということである。そうすると
2.2 の冒頭で述べた固定点への漸近という考え方が成り立たなくなるように思われる。また初期𝑧𝑧
推定状態は推定元の状態と同じ𝑆𝑡𝑜𝑡𝑎𝑙となるため、そもそも𝑆𝑡𝑜𝑡𝑎𝑙が変化した初期状態を作ることが𝑧できない。しかし熱力学的極限で考えると、磁化曲線が滑らかになるなら𝑆𝑡𝑜𝑡𝑎𝑙
が１や２違っても𝑧⁄𝐿が一致するならば同じ性質を持った状態が実現するはず状態は同じであり、また磁化𝑀 = 𝑆𝑡𝑜𝑡𝑎𝑙である。そこで高磁場領域での初期推定に利用できる状態が存在するのか調べてみた。
初期推定状態は初期推定元の状態を各成分の絶対値はそのままに符号を適切に変化させたものである。高磁場領域の初期推定計算で問題なのは基底の対応付けの仕方であり、初期推定ができるならば推定元と推定対象の基底状態の成分の絶対値の分布は一致するはずである。そこで(30)式の保持基底の内積計算と同様に、DMRG 計算で得られる各基底状態に対する絶対値順内積計算
𝑧(𝐿𝑙𝑎𝑠𝑡 ) = 10の𝑆 = 1/2ハイゼンベルグ鎖での計算結果であをしてみた。表１は𝐿𝑙𝑎𝑠𝑡 = 200、𝑆𝑡𝑜𝑡𝑎𝑙𝑧𝑧𝑧𝑧(𝐿 = 22) = 2, 𝑆𝑡𝑜𝑡𝑎𝑙(𝐿 = 42) = 3, ⋯ , 𝑆𝑡𝑜𝑡𝑎𝑙(𝐿 = 182) =る。この計算では𝑆𝑡𝑜𝑡𝑎𝑙(𝐿 = 4) = 1から𝑆𝑡𝑜𝑡𝑎𝑙𝑧𝑧
10と 20 サイト大きくなる毎に𝑆𝑡𝑜𝑡𝑎𝑙
を 1 ずつ上げてゆく。𝐿 = 178は𝑆𝑡𝑜𝑡𝑎𝑙= 9となってから 9 ス𝑧
テップ目の計算であり、同じく𝑆𝑡𝑜𝑡𝑎𝑙= 9の𝐿 = 174から初期推定をしている。このとき表１の178－174の値（青文字）はほぼ 1 になっており、初期推定ができるならばこの計算で 1 に近い値が𝑧得られることが分かる。𝑆𝑡𝑜𝑡𝑎𝑙= 10となった直後の𝐿 = 182の行を見ると、通常の初期推定で利用する 4 サイト小さい𝐿 = 178の基底状態との値は178－174の時ほどは 1 に近くない。いろいろな状態との内積計算を行ってみたところ、𝐿 = 142との値が 1 に近く（赤文字）、この状態から初期推定ができる可能性があることが分かった。これらの状態をより詳しく調べると、中心２スピン𝑧𝑧の対応はそのままに、ブロック基底の対応が|𝑆𝑏𝑙𝑜𝑐𝑘, 𝑛⟩𝐿=142 ↔ |𝑆𝑏𝑙𝑜𝑐𝑘+ 1, 𝑛⟩𝐿=182 となっていた。

これは全系の状態はほとんど同じであり、ブロック基底も類似しているが、サイト数が増えた分
𝑧𝑧
𝑆𝑏𝑙𝑜𝑐𝑘の値が一様に大きくなったのであろう。また𝑆𝑡𝑜𝑡𝑎𝑙= 8となってから 9 ステップ目の計算である𝐿 = 158からは初期推定ができそうにないのは、基底の違いが原因であると思われる。
𝐿 = 178

𝐿 = 174

𝐿 = 162

𝐿 = 158

𝑧(𝑆𝑡𝑜𝑡𝑎𝑙= 9 第 1step)𝐿 = 178𝑧(𝑆𝑡𝑜𝑡𝑎𝑙= 9 第 9step)𝐿 = 182𝑧(𝑆𝑡𝑜𝑡𝑎𝑙= 10 第 1step)𝐿 = 142𝑧(𝑆𝑡𝑜𝑡𝑎𝑙= 8 第 1step)1

0.999990.907560.98046

0.93389

0.925130.92534

0.98027

0.92293

0.99927

表１ 様々な状態間の絶対値順内積計算。初期推定ができる178－174間の値はほとんど 1 にな𝑧𝑧っている（青文字）。𝐿 = 182は𝑆𝑡𝑜𝑡𝑎𝑙= 10になってから 1 ステップ目の計算であるが、𝑆𝑡𝑜𝑡𝑎𝑙=8になってから 1 ステップ目の計算である𝐿 = 148との値が他のものに比べ 1 に近い（赤文字）。

14《寄稿論文》

物性研究・電子版 Vol. 5, No. 2, 052101 （2016年5月号）

𝑧
𝐿 = 182,142はどちらも𝑆𝑡𝑜𝑡𝑎𝑙
増加直後であり、この時のブロック基底は１ステップ前の異なる𝑧
𝑆𝑡𝑜𝑡𝑎𝑙での基底状態から決められているため、このステップの基底状態を計算するのに最適化され𝑧ていない。対して𝐿 = 158のブロック基底は同じ𝑆𝑡𝑜𝑡𝑎𝑙の基底状態から計算されたものであり、この基底状態を計算するのに適切なものになっている。このように高磁場領域の初期推定では推定対象の状態がどのような基底で計算されるかも重要になる。
𝑧𝑧以上より、高磁場領域の初期推定では𝑆𝑡𝑜𝑡𝑎𝑙− 2 の状態の中で𝑆𝑡𝑜𝑡𝑎𝑙
増加後ステップ数が同じもの
𝑧𝑧から基底の対応付けを|𝑆𝑏𝑙𝑜𝑐𝑘+ 1, 𝑛⟩𝐿 ↔ |𝑆𝑏𝑙𝑜𝑐𝑘, 𝑛⟩𝐿′ として初期推定状態を作ればよいと思われる。
𝑧
もちろん𝑆𝑡𝑜𝑡𝑎𝑙
が変化せず通常の初期推定が適用できる部分はそちらを使えばよい。図８は𝐿𝑙𝑎𝑠𝑡 =𝑧(𝐿𝑙𝑎𝑠𝑡 ) = 10の𝑆 = 1/2ハイゼンベルグ鎖の初期推定誤差である。保持基底数𝑚 = 200で200、𝑆𝑡𝑜𝑡𝑎𝑙𝑧
切り捨て誤差は10−11であった。20 サイト大きくなる毎（DMRG10 ステップ毎）に𝑆𝑡𝑜𝑡𝑎𝑙を1ず𝑧𝑧
つ上げてゆくが、𝑆𝑡𝑜𝑡𝑎𝑙
増加後の 3 ステップは推定元との𝑆𝑡𝑜𝑡𝑎𝑙の違いや基底の違いから通常の初𝑧期推定ができない。ゆえにこれらに対しては𝑆𝑡𝑜𝑡𝑎𝑙− 2の状態からの推定法を適用した。グラフか𝑧
らわかるように、𝑆𝑡𝑜𝑡𝑎𝑙が変化する点においてもそうでない点ほどではないが初期推定ができてい𝑧𝑧
る。同一の𝑆𝑡𝑜𝑡𝑎𝑙での推定は𝐿 − 4の状態から計算しており、異なる𝑆𝑡𝑜𝑡𝑎𝑙での推定は𝐿 − 40の基底状

態から初期推定状態を求めている。有限サイズ効果等を考えれば推定誤差は前者の方が小さくなるのは理解できよう。

図８

𝑆 = 1/2ハイゼンベルグ鎖の高磁場領域の初期推定誤差。𝑚 = 200、切り捨て誤差は10−11 だった。白丸は通常の初期推定時、バツは高磁場初期推定時の誤差、青線は符号調整が限界まで行わ𝑧
れた場合の誤差である。𝑆𝑡𝑜𝑡𝑎𝑙が増加した点においてもある程度初期推定ができている。

15

《寄稿論文》

物性研究・電子版 Vol. 5, No. 2, 052101 （2016年5月号）

このようにして高磁場領域でも確かに初期推定を行うことができた。なお、図８の計算では𝐿𝑙𝑎𝑠𝑡𝑧𝑧(𝐿𝑙𝑎𝑠𝑡 )の値により常に𝑆𝑡𝑜𝑡𝑎𝑙と𝑆𝑡𝑜𝑡𝑎𝑙− 2の状態から初期推定ができたが、磁化曲線を求める場合は𝑧𝑧
すべての𝑆𝑡𝑜𝑡𝑎𝑙
での最低エネルギー固有値が必要であり、必ずしも𝑆𝑡𝑜𝑡𝑎𝑙− 2の状態から推定できる𝑧とは限らない。スピンが半整数の場合、𝑆𝑡𝑜𝑡𝑎𝑙− 2の状態のブロックのサイト数が異なっていると𝑧𝑧
基底の対応付けができないこともある。その場合には𝑆𝑡𝑜𝑡𝑎𝑙− 4など𝑆𝑡𝑜𝑡𝑎𝑙− 2𝑛 (𝑛 = 1,2, ⋯ )で表される状態のなかでブロックのサイト数の偶奇が一致するものから初期推定状態を作ればよい。
𝑧
なお、この方法はあくまで高磁場領域計算時の無限系法を高速化する方法であり、途中で𝑆𝑡𝑜𝑡𝑎𝑙

を変化させている以上、無限系法のみではその結果は信頼できないことに注意されたい。きちんと性質等を議論するには DMRG 有限系法[1]による精度の向上が不可欠であるが、そちらの計算においても別の初期推定方法[4]が考案されており、大きな成功を収めている。

16

《寄稿論文》

物性研究・電子版 Vol. 5, No. 2, 052101 （2016年5月号）

3 基底・低励起状態計算
3.1

𝑺 = 𝟏/𝟐ハイゼンベルグ鎖𝑧𝑧
𝑆 = 1/2 ハイゼンベルグ鎖(1)式での基底状態(𝑆𝑡𝑜𝑡𝑎𝑙= 0)と第一励起状態(𝑆𝑡𝑜𝑡𝑎𝑙= 1)の計算を行った。ここでは𝐿 − 4での状態から初期推定を行った。図９は𝑚 = 200, 𝐿 = 400での計算結果であり、上図は横軸を線形軸、下図は対数軸にとったものである。(5)式の切り捨て誤差𝑒𝑡𝑟 の最大値はどちらもおおよそ10−10であった。基底状態も第一励起状態も初期推定誤差（白丸）と最大推定誤差（青線）は一致しており完璧に符号調整ができている。この系はギャップレスで基底状態も第一励起状態も Tomonaga-Luttinger (TL)状態だが、それを反映してか推定誤差の𝐿依存性はどちらも同じであった。またグラフは横軸を対数軸にすると線形になり（下図）、初期推定誤差は𝐿のべき乗で表されることがわかる。𝑆 = 1/2 反強磁性ハイゼンベルグ鎖のエネルギーやその他の物理量の多くは𝐿に対しべき乗型の依存性を示すことが知られているが[7][8]、初期推定誤差の依存性もそれと一貫したものになっている。

図９ 𝑆 = 1/2 ハイゼンベルグ鎖の基底状態（左図）と第一励起状態（右図）での初期推定誤差。白丸は０次、
四角は１次の初期推定状態のエネルギーの相対誤差であり、青線は限界符号調整状態の誤差である。

17

《寄稿論文》

物性研究・電子版 Vol. 5, No. 2, 052101 （2016年5月号）

また１次近似での初期推定誤差（□）はほぼすべての点で減少しており、ばらつきはあるがこちらもべき乗的な依存性を示した。𝐿~200の時点で切り捨て誤差と同じオーダーになり、さらに大きくなると対角化計算にほとんどコストがかからなくなった。また 2.4 でも述べたが、１次推定の時点ですでにばらつきがあり誤差を𝐿−1の関数で表すことができないため、グラフにはしていないが２次推定を行ってもあまり改良されなかった。

3.2

𝑺 = 𝟏ハイゼンベルグ鎖𝑧
𝑆 = 1 反強磁性ハイゼンベルグ鎖の場合、初期推定誤差の𝐿依存性は各𝑆𝑡𝑜𝑡𝑎𝑙
で大きく異なる。図𝑧𝑧𝑧１０は基底状態(𝑆𝑡𝑜𝑡𝑎𝑙= 0)と第一励起状態(𝑆𝑡𝑜𝑡𝑎𝑙= 1)、図１１は第二励起状態(𝑆𝑡𝑜𝑡𝑎𝑙= 2)の初期推定誤差である。𝑚 = 150, 𝐿 = 200、切り捨て誤差はそれぞれ10−10 , 10−12 , 10−10であった。𝐿 − 2の状態から推定を行ったが、すべての場合で完璧に符号調整されている。図１０左の基底状態について、初期推定誤差は特殊な振る舞いをしている。𝑆 = 1/2も含め他の状態ではなめらかに変化しているのに対し、この状態では少なくとも２つの kink が存在し、𝐿~100で収束しかけるがその後再び減少してゆく。このような異常性は縮退が原因と思われる。𝑆 = 1 反強磁性ハイゼンベルグ
鎖の基底状態は Valence Bond Solid(VBS)状態であり、OBC では疑似四重縮退を持つ。これは端のスピンの自由度に由来するものであり、バルク極限ではどれも VBS 状態になるため、本来𝑧𝑧𝑧
𝑆𝑡𝑜𝑡𝑎𝑙= 1の振る舞いは一致するはずである。しかし四重縮退のうち𝑆𝑡𝑜𝑡𝑎𝑙= 0に２つ、= 0と𝑆𝑡𝑜𝑡𝑎𝑙𝑧𝑧
𝑆𝑡𝑜𝑡𝑎𝑙= ±1にそれぞれ１つ準位が存在するため、𝑆𝑡𝑜𝑡𝑎𝑙= 0の計算では𝐿が大きくなると縮退している２つの状態が混ざってしまい、このような異常性が現れたと考えられる。
第一励起状態（図１０右）の場合、推定誤差の log は𝐿~150まで線型的に減少したのち一定値をとる。これは VBS 状態が𝐿に対し指数的な依存性で特徴付けられること[9]と一致しており、一定

図１０ 𝑆 = 1ハイゼンベルグ鎖の基底状態（左）と第一励起状態（右）計算での初期推定誤差。𝐿 − 2の状態から推定を行った。基底状態では異常な振る舞いを見せるが、第一励起状態では誤差は指数的に減少してゆく。

18

《寄稿論文》

物性研究・電子版 Vol. 5, No. 2, 052101 （2016年5月号）

図１１

𝑆 = 1ハイゼンベルグ鎖の第二励起状態計算での初期推定誤差。𝐿~70に kink があり（左）、対数軸にするとその前後で𝐿依存性が異なることが分かる（右）。

値になるのは切り捨て誤差や数値計算誤差などによる推定精度の限界に達したためと考えられる。
また第二励起状態（図１１左）では一見滑らかだが、𝐿~70で kink が存在しその前後で依存性が変化している。横軸を対数にとると（図１１右）その振る舞いから推定誤差の𝐿依存性は𝐿 ≲ 70
では指数的、𝐿 ≳ 70ではべき乗的になると予想される。この状態は連続準位の底にあたり、第一励起状態とのエネルギー差は𝐿が大きい時には𝐿−2 に依存することが分かっている[9]。また𝐿が小さい時には上で述べた指数的な依存性が残っているだろう。つまり、各𝐿でのエネルギー依存性と初期推定誤差の依存性は見事に一致しており、その振る舞いを説明付けることができる。
図１２は𝐿 − 4の状態から推定を行った時の初期推定誤差である。比較のため𝐿 − 2の状態からの𝑧
推定誤差も載せた。𝑆𝑡𝑜𝑡𝑎𝑙= 0では図１０と同様に異常な振る舞いを示し、改良推定ができなかったため除いた。非改良時（０次推定）では𝐿 − 4と𝐿 − 2の推定で値が異なっている。有限サイズ効𝑧
果を考えれば誤差は𝐿 − 4の方が大きくなりそうだが、𝑆𝑡𝑜𝑡𝑎𝑙= 1では逆に小さくなっている。詳しく調べると基底状態|𝜓0 ⟩𝐿 の各成分は振動しながら固定点に漸近していた。このとき𝐿と𝐿 − 2は逆位相、𝐿と𝐿 − 4は同位相になっており、𝐿 − 4では振動による誤差が消え推定誤差が小さくなった𝑧と思われる。𝑆𝑡𝑜𝑡𝑎𝑙= 2では誤差はスピン鎖が短いときは𝐿 − 4の方が小さいが、長くなると𝐿 − 2の𝑧方が小さくなる。また𝐿 − 4の誤差は常にべき乗的な依存性を示した。つまり𝑆𝑡𝑜𝑡𝑎𝑙= 2においても振動が存在するが、振動による誤差は指数的であり振動を除いた誤差はべき乗的と考えられる。
サイト数𝐿が小さいときは振動による指数的な誤差が大きく、それを除去できる𝐿 − 4の方が良い
推定ができ、𝐿が大きい時には有限サイズ効果が小さくなる𝐿 − 2の方が良くなるのであろう。
𝑧𝑧
𝐿 − 4の１次近似による改良推定誤差は特に𝑆𝑡𝑜𝑡𝑎𝑙= 2で小さくなった。つまり𝑆𝑡𝑜𝑡𝑎𝑙= 2の推定誤𝑧𝑧
差は∆𝐿−1の依存性が強いが、𝑆𝑡𝑜𝑡𝑎𝑙= 1は高次の項の影響が強いと思われる。𝑆𝑡𝑜𝑡𝑎𝑙= 1の有限サイ

ズ効果が指数的であることを考えれば妥当であり、これゆえ２次の推定を行ってもさらに高次の項の影響であまり改善されなかった。なお𝐿 − 2での近似推定では上記の振動性により線形近似が成り立たずまったく改良されなかった。
19

《寄稿論文》

物性研究・電子版 Vol. 5, No. 2, 052101 （2016年5月号）

𝑧
図１２ 𝑆 = 1ハイゼンベルグ鎖の𝐿 − 4の状態からの初期推定誤差。𝐿 − 2からの推定と比べ𝑆𝑡𝑜𝑡𝑎𝑙= 1（左図）𝑧ではすべての𝐿で誤差が減るが、𝑆𝑡𝑜𝑡𝑎𝑙= 2（右図）では𝐿が大きい時は悪化する。１次近似を行えばこちらもすべての𝐿で誤差が減る。

3.3

𝑺 = 𝟏/𝟐ジグザグ鎖𝑆 = 1/2ジグザグ鎖(2)式の基底状態計算での初期推定を行った。この系の相図は調べられており[10]、図１３左のようになる。𝐽2⁄𝐽1 = 0では系は１本のハイゼンベルグ鎖と一致するが、
𝐽2⁄𝐽1 ~0.27付近で TL 相から dimer 相へ転移する。また𝐽2⁄𝐽1 = 0.5のとき系は Majumdar-Ghoshmodel と呼ばれ、基底状態は図１３右のように１つ飛ばしの最隣接ボンドでシングレットを組む
厳密な Dimer 状態になる。さらに𝐽2⁄𝐽1を大きくすると並進対称性の破れが現れ[11]、𝐽2⁄𝐽1 → ∞の極限では２本の独立なハイゼンベルグ鎖となる。今回、𝐽2⁄𝐽1 = 0.2, 0.4, 1.0, 10.0および𝐽2⁄𝐽1 =0.5の５点において初期推定を行った。𝐽2⁄𝐽1 = 0.5以外の４つの保持基底数𝑚は全て 200、サイト
数は𝐿 = 400、切り捨て誤差はそれぞれ10−10 , 10−11 , 10−9 , 10−6であった。𝐽2⁄𝐽1 = 0.5を除いているのは DMRG 計算と初期推定が特殊な振る舞いを示すためであり、その原因も含め後に説明する。

図１３ 𝑆 = 1/2ジグザグ鎖の基底状態の相図と Majumdar-Ghosh model の基底状態（右図）。右図の実線は最隣接ボンド、点線は次隣接ボンドである。赤丸で示したペアがシングレットを組む。

20

《寄稿論文》

物性研究・電子版 Vol. 5, No. 2, 052101 （2016年5月号）

図１４に上記の４つの計算における初期推定誤差を示す。まず𝐽2⁄𝐽1 = 0.2(左上)について、系はハイゼンベルグ鎖と同じ TL 状態であり、初期推定誤差も同じ振る舞いを示す（初期推定誤差はべき乗的な依存性を持つ）。𝐽2⁄𝐽1 = 0.4(右上)では系は dimer 状態でありギャップを持つ。このとき𝑆 = 1 ハイゼンベルグ鎖と同じくエネルギーの𝐿依存性は指数的になる[11]が、初期推定誤差も同様に指数的に振る舞う。𝐽2⁄𝐽1 = 1.0(左下)も同じく dimer 状態であり、限界精度まで指数的に減少している。しかし𝐽2⁄𝐽1 = 0.4よりがたつきが大きい。これは並進対称性の破れが原因と思われる。𝑆 = 1/2ジグザグ鎖では𝐽2⁄𝐽1が大きくなると並進対称性が破れスピン相関などで数十サイ
トの周期の大きな振動が現れる[11]。DMRG 無限系法では 2 サイトずつ増やしてゆくが、このよ
うな大きな周期性により計算される状態が少しずつ変わってくる。これによる初期推定誤差ががたつきとなって表れたと考えられる。𝐽2⁄𝐽1 = 10.0となると相図上では dimer だが、実際にはほとんど独立な２本のハイゼンベルグ鎖と見なせるだろう。この時𝐿 − 4からの初期推定(右下)ではほとんど初期推定ができていない。スピンが半整数のとき初期推定状態はブロック部分のサイト

図１４

𝑆 = 1/2ジグザグ鎖の初期推定誤差。それぞれ

𝐽2⁄𝐽1 = 0.2, 0.4, 1.0, 10.0の基底状態計算時の値である。

21

《寄稿論文》

物性研究・電子版 Vol. 5, No. 2, 052101 （2016年5月号）

数の偶奇性が一致している必要があるが、２本鎖であることを考えると𝐿 − 4からの推定ではスピン鎖１本あたりの偶奇性は異なってしまっている。そこで図１５のように𝐿 − 8からの推定を行ったところ、計算誤差や並進対称性の破れが原因と思われるばらつきがあるものの、ある程度初期推定をすることができた。
最後に𝐽2⁄𝐽1 = 0.5の場合について言及しておく。この場合の初期推定誤差はゼロ、つまり初期推定状態が真の基底状態と厳密に一致した。上で述べたようにこの系の基底状態|𝜓𝑑 ⟩は図１３右のようにサイト2𝑖 − 1と2𝑖がシングレットを組む；|𝜓𝑑 ⟩ = 2−𝐿⁄2 (|+⟩1 |−⟩2 − |−⟩1 |+⟩2 ) ⊗ ⋯ ⊗ (|+⟩𝐿−1 |−⟩𝐿 − |−⟩𝐿−1 |+⟩𝐿 )(31)

このとき DMRG 計算では保持基底数𝑚 ≥ 2で厳密な状態を計算できる[11]。(31)式の様に基底状
態は２サイトごとのシングレットの直積で表されるため、ブロックの基底として|𝑏⟩+ = 2−𝐿⁄4 (|+⟩1 |−⟩2 − |−⟩1 |+⟩2 ) ⊗ ⋯ ⊗ (|+⟩𝐿/2−3 |−⟩𝐿/2−2 − |−⟩𝐿/2−3 |+⟩𝐿/2−2 )⨂|+⟩𝐿⁄2−1|𝑏⟩− = 2−𝐿⁄4 (|+⟩1 |−⟩2 − |−⟩1 |+⟩2 ) ⊗ ⋯ ⊗ (|+⟩𝐿/2−3 |−⟩𝐿/2−2 − |−⟩𝐿/2−3 |+⟩𝐿/2−2 )⨂|−⟩𝐿⁄2−1(32)の２つを保持していれば次ステップでの厳密な基底状態を表現できるためである。そしてこの直積性により𝐿 − 4と𝐿の基底状態はこれらのブロック基底と中心 2 スピン基底の全く同じ結合で表されるので、初期推定状態が基底状態と一致するのである。

図１５𝐽2⁄𝐽1 = 10.0におけるジグザグ鎖基底状態計算の初期推定誤差。𝐿 − 8の基底状態から初期推定を行った。

22

《寄稿論文》

物性研究・電子版 Vol. 5, No. 2, 052101 （2016年5月号）

4 高磁場領域計算
4.1

𝑺 = 𝟏/𝟐ハイゼンベルグ鎖𝑧𝑧
この系の計算では𝑆𝑡𝑜𝑡𝑎𝑙が増加してから３ステップは 2.6 の𝑆𝑡𝑜𝑡𝑎𝑙
が異なる場合の、それ以外は𝑧𝑧
2.2 の𝑆𝑡𝑜𝑡𝑎𝑙が同じ場合の初期推定を行った。2.6 で説明したように𝑆𝑡𝑜𝑡𝑎𝑙
増加後１，２ステップ目
𝑧𝑧は推定元状態と𝑆𝑡𝑜𝑡𝑎𝑙が異なり、３ステップ目は基底を決める状態の𝑆𝑡𝑜𝑡𝑎𝑙が異なる。図１６は𝐿 =𝑧𝑧400, 𝑆𝑡𝑜𝑡𝑎𝑙= 50,100,150の計算結果である。保持基底数は𝑆𝑡𝑜𝑡𝑎𝑙= 50では𝑚 = 200(左上)と50(右上)𝑧の２通り、𝑆𝑡𝑜𝑡𝑎𝑙= 100,150ではそれぞれ𝑚 = 200(左下),150(右下)である。切り捨て誤差は順に𝑧10−11, 10−8, 10−11,10−11であった。また𝑆𝑡𝑜𝑡𝑎𝑙= 20の計算が図８なのでそちらも参照されたい。
𝑧𝑧
𝑆𝑡𝑜𝑡𝑎𝑙= 50, 𝑚 = 50以外の結果について、図８のように𝑆𝑡𝑜𝑡𝑎𝑙
が小さく初期推定元の状態が何十サ
𝑧
イトも小さいと推定誤差は明らかに大きいが、𝑆𝑡𝑜𝑡𝑎𝑙= 50(𝐿 − 16から推定)や100(𝐿 − 8から推定)𝑧になるとあまり変わらなくなる。しかし𝑆𝑡𝑜𝑡𝑎𝑙= 150になるとむしろ悪化する。これはターゲット𝑧
状態の𝑆𝑡𝑜𝑡𝑎𝑙が頻繁に変化し、適切な基底があまり残されていないためと思われる。

図１６

𝑆 = 1⁄2ハイゼンベルグ鎖の高磁場領域計算の初期推定誤差。

23

《寄稿論文》

物性研究・電子版 Vol. 5, No. 2, 052101 （2016年5月号）

𝑧
𝑧
𝑆𝑡𝑜𝑡𝑎𝑙= 50, 𝑚 = 50の振る舞いは奇妙なものとなった。𝑆𝑡𝑜𝑡𝑎𝑙
増加時の推定誤差は指数的に減少し、
𝑧
𝐿~100程で𝑆𝑡𝑜𝑡𝑎𝑙
一定時の初期推定より誤差が小さくなる。有限サイズ効果を考えれば𝐿が近い状

態からの推定の方が誤差は小さくなるはずだが、実は初期推定誤差は他の要因で振る舞いが変化することが分かった。図１７は𝑆 = 1/2 ハイゼンベルグ鎖の基底状態計算をあえて小さな保持基
底数𝑚で行った時の初期推定誤差であるが、このときも誤差減少が起こっている。𝑚 = 100(白丸)では誤差の減少はだんだん緩やかになるが、𝑚 = 20(青三角)では𝐿 ≳ 50から指数的に減少している。𝑚が小さい時に DMRG 計算がどうなるか考えてみよう。保持基底は密度行列固有値の大きなものから採用するので、𝑚が小さいとどの基底も大きな固有値を持つことになる。このとき密度行列(4)式から分かるように、これらの基底を用いて計算される全系の基底状態はどの成分も大きな値を持つだろう。もし有限サイズ効果による変化より十分大きな値ならば、初期推定状態は真の基底状態と非常に近くなるため、𝑚が大きい時より誤差が減少したと考えられる。これと同じ𝑧𝑧ことが𝑆𝑡𝑜𝑡𝑎𝑙が変化する部分の初期推定でも起こったのだろう。𝑆𝑡𝑜𝑡𝑎𝑙が増加し DMRG 法のターゲット状態が変わると、その状態を計算するのに適切でない基底が保持されることになり、実質的な保持基底数が小さくなる。ゆえに高磁場推定部分のみが図１７と同じように減少していったと
思われる。

図１７

𝑆 = 1/2 ハイゼンベルグ鎖の基底状態計算における初期推定誤差。

保持基底数𝑚が小さいほど誤差が小さくなっている。

24

《寄稿論文》

物性研究・電子版 Vol. 5, No. 2, 052101 （2016年5月号）

4.2

𝑺 = 𝟏ハイゼンベルグ鎖𝑧
𝑆 = 1ハイゼンベルグ鎖の高磁場領域初期推定計算は、𝑆𝑡𝑜𝑡𝑎𝑙が大きい時のばらつきや𝑚が小さい𝑧⁄𝐿 = 0.5で時の誤差減少など𝑆 = 1/2とほぼ同じ振る舞いを示した。この項では唯一異なった𝑆𝑡𝑜𝑡𝑎𝑙𝑧の計算について述べる。図１８は𝐿 = 200, 𝑚 = 150, 𝑆𝑡𝑜𝑡𝑎𝑙= 100での初期推定誤差である。この時切り捨て誤差は10−9であった。左図は 2.2 の演算子からの調整と 2.3 前半の１成分毎の期待値からの調整で初期推定計算を行ったものである。このとき𝐿 = 130において符号調整に失敗していた。
この計算を詳しく調べると、演算子からの調整のみでも８割方の点で符号調整ができており、そうでない点も初期推定誤差は小さくなっていて基底状態に近い状態が得られていたが、𝐿 = 130では全く調整ができていなかった。また１成分毎の調整のみを適用してみると半数以上の点で符号調整ができなかった。ある程度調整できていて基底状態に近い状態になっていれば１成分毎の調整でも正しい初期推定状態を得ることができるため、演算子からの調整が全くうまくいかなかっ
た𝐿 = 130では期待値からの調整を行っても符号は間違ったままだったのであろう。そこで演算子からの調整と 2.3 の後半で述べたブロック基底毎の調整で初期推定を行ったところ、すべての点で符号調整ができた（右図）。もし基底毎の調整でなければ符号調整ができない系であれば計算の高速化になると思われるが、今回は符号調整に失敗した点が１つだけだったため、計算量の多い
こちらの方法では数%~10 数%程計算時間が増えてしまった。一般的に利用するなら計算量を減らす改良が必要であろう。

図１８

𝑆 = 1ハイゼンベルグ鎖の高磁場領域計算での初期推定。１成分毎の調整（左図）では矢印で示した𝐿 = 130において符号調整を失敗しているが、ブロック基底毎の調整（右図）ではすべての点で調整がされている。

25

《寄稿論文》

物性研究・電子版 Vol. 5, No. 2, 052101 （2016年5月号）

5 まとめ
本論文では無限系法初期推定の改良と各種モデルへの適用を行った。２章ではこれまでの初期推定法の不完全性を示したのち、エネルギー期待値による符号調整や初期推定状態の近似推定、および密度行列固有値の入れ替えを考慮した基底対応付けなどの改良を行った。またこれまで行われていなかった高磁場領域の初期推定について考察した。３・４章では実際の系への適用を通じ、確かに初期推定計算が改良されることを示したほか、ここで現れる現象から初期推定計算の性質について調査した。推定誤差は系のサイズ依存性と関連しているほか、保持基底数𝑚が小さいときは DMRG の性質に由来する誤差減少が起こることが分かった。𝑆 = 1ハイゼンベルグ鎖の高磁場領域の計算では改良初期推定でも失敗することがあり、2.3 後半のブロック基底毎の調整方法を用いることで完璧な符号調整ができた。
このようにさまざまな改良を行うことで、DMRG 無限系法でも有限系法での初期推定計算[4]と同等まで計算時間を減らせる。無限系法初期推定の原理的に系が大きくなるほど有限サイズ効
果が抑えられ対角化の計算量が劇的に少なくなる。特に近似推定による改良の効果は大きいが、これを行うには基底の符号調整が完璧でなければならない。しかし１成分毎の符号調整では失敗することがある一方、ブロック基底毎では完璧に調整できるが計算時間が増えてしまう。また
DMRG 法の特性や系の性質など、初期推定とはあまり関係のない部分が原因で誤差が出てしまうこともある。更なる計算時間の短縮には、初期推定計算だけでなく DMRG 法の計算にも踏み込んだ改良が必要である。

参考文献[1] S.R. White, Phys. Rev. Lett. 69, 2863 (1992); Phys. Rev. B 48, 10345 (1993).
[2] 山田進, 五十嵐亮, 奥村雅彦, 今村俊幸, 町田昌彦, 応用数理 20(2), 132-147 (2010).
[3] S. Qin and J. Lou, Commun. Theor. Phys. 36, 635 (2001).
[4] S. R. White, Phys. Rev. Lett. 77, 3633 (1996).
[5] U. Schollwö ck, Phys. Rev. B 58, 8194 (1998); 59, 3917(E) (1999).
[6] L. Sun, J. Zhang, S. Qin, and Y. Lie, Phys. Rev. B 65, 132412 (2002).
[7] K. Nomura, Phys. Rev. B 48, 16814 (1993).
[8] H. W. J. Blö te, J. L. Cardy, and M. P. Nightingale, Phys. Rev. Lett. 56, 742 (1986).
[9] S. Qin, Y. Liu, and L. Yu, Phys. Rev. B 55, 2721 (1997).
[10] K. Okunishi, T. Tonegawa, J. Phys. Soc. Jpn. 72, pp. 479-482 (2003)[11] S. R. White, I. Affleck, Phys. Rev. B 54, 9862 (1996).

26