《講義ノート》

物性研究・電子版 Vol. 7, No. 2, 072205 （2018年11月号）

密度行列繰り込み群の原理と応用
東北大学 理学研究科
柴田 尚和

密度行列繰り込み群は、1 次元量子多体系の基底状態を精密に求める方法として 1992 年に

S. White によって考案された数値的解析法である [1, 2]。小さな系の基底状態を出発点として、密度行列の対角化によって重要な基底だけを残しながら系を大きくするこの計算法は、既存の方法では難しかった電子系やフラストレート量子スピン系の取り扱いを可能にし、解析的に解けない量子多体問題を数値的に解く新たな道を切り開いた。その最も重要な特徴は基底数を制限しながら系の拡張を繰り返すことで行列積型の波動関数を構成し、その最適化を変分原
理を満たすように行うところにある。その優れたアイディアは後に高次元古典系 [3, 4]、有限温度量子系の熱力学量 [5, 6, 7]、2 次元量子系の基底状態 [8, 9]、時間発展 [10] の計算に応用され、テンソルネットワーク [11, 12] へと発展した。このテキストは、こうした多くの研究分野で使われている密度行列繰り込み群の原理と応用 [13] について解説したものである。

密度行列繰り込み群の原理

1多電子系の波動関数を表現するためには、一般に任意の電子配置を表現できる基底が必要であり、その数は電子数や系のサイズが大きくなると指数関数的に増大する。例えば、格子点上の各原子あたりに 1 つの電子が局在する S = 1/2 量子スピン系を考えた場合、独立なスピン
配置の数 N は電子数 Ne の増加に対して N = 2Ne のように増加し、30 個程度のスピンの系を考えてもその基底の数 (異なるスピン配置の数) は 10 億を超えてしまう。このような基底数の急激な増大は系のハミルトニアンの次元を増大させ、その数値的対角化を困難にすること
で多体問題の直接的解法に実質上の制限を与える。この一般的な数値的解法の限界を超えた解析を行うためには発散的に増大する多体基底数の問題を解決する必要があり、その最良の解決策の一つが密度行列を用いた基底選択である。

1.1

密度行列の特徴

この基底選択の方法について考える前に、密度行列から量子多体系のどのような特徴を知ることができるか確認してみよう。波動関数はそこに含まれる変数によって一意に値が定まる関数であり、その絶対値の 2 乗は個々の電子の配置や量子化されたスピンの向き等を指定する
波動関数の変数が示す状態の実現確率を表す。多体系の場合、電子やスピンは複数存在するた

《講義ノート》

物性研究・電子版 Vol. 7, No. 2, 072205 （2018年11月号）

め、結果的に電子やスピンの相互の関係によって変化する確率を波動関数は表すことになる。
例えば、反強磁性相互作用が働いている量子スピン系の低エネルギー状態を考えれば、隣り合
うスピンが互いに反対向きのとき、波動関数の振幅は大きくなり、逆に同じ向きのときにはその振幅は小さくなっているはずである。まず、このことを具体的に数式を使って確認してみよう。
いま、系の中に含まれる 2 つの S = 1/2 スピンに注目し、その状態を σ1 = {↑1 , ↓1 }, σ2 ={↑2 , ↓2 } で表す。この 2 つ以外のスピンの状態はその独立なスピン配列の一つ一つに番号を付け、(例えば ↓ は 0、↑ は 1 に対応させる変換によってスピン配列を２進数に置き換えて) 整数

j で表す。j は複数のスピンの状態を表すので、j で指定される数はその異なるスピン配列の数だけ存在する。系全体の任意のスピン配置は σ1 , σ2 , j を使って表すことができ、多体の波動
関数 |Ψ⟩ は、それらの重ね合わせによって

|Ψ⟩ =

∑

Ψσ1 ,σ2 ,j |σ1 ⟩|σ2 ⟩|j⟩(1)

σ1 ,σ2 ,jと表現することができる。ここで Ψσ1 ,σ2 ,j は線形結合の係数である。この波動関数において、σ1 が ↑、σ2 が ↓ の状態である量子力学的な実現確率 P↑1 ,↓2 は、
P↑1 ,↓2 = ⟨Ψ|δσ1 ,↑1 δσ2 ,↓2 |Ψ⟩ =∑Ψ∗↑1 ,↓2 ,j Ψ↑1 ,↓2 ,j =j∑|Ψ↑1 ,↓2 ,j |2(2)|Ψ↑1 ,↑2 ,j |2(3)jであり、σ1 が ↑、σ2 も ↑ の状態である確率 P↑1 ,↑2 は、

P↑1 ,↑2 = ⟨Ψ|δσ1 ,↑1 δσ2 ,↑2 |Ψ⟩ =∑Ψ∗↑1 ,↑2 ,j Ψ↑1 ,↑2 ,j =j∑jとなる。このような確率は相関関数 ⟨Ψ|S1z S2z |Ψ⟩ やエネルギーの期待値の計算においても現れ、例えば、反強磁性交換相互作用の大きさを J > 0 とする σ1 , σ2 のスピン間相互作用の対角成分 JS1z S2z によって生じるエネルギーの期待値 E1,2 はE1,2= ⟨Ψ|JS1z S2z |Ψ⟩

(4)

=

J∑
(|Ψ↑1 ,↑2 ,j |2 − |Ψ↑1 ,↓2 ,j |2 − |Ψ↓1 ,↑2 ,j |2 + |Ψ↓1 ,↓2 ,j |2 )
4 j

(5)

=

J
{(P↑1 ,↑2 + P↓1 ,↓2 ) − (P↑1 ,↓2 + P↓1 ,↑2 )}4(6)となり、エネルギー E1,2 が負になることと、個々のスピン配置の実現確率に差が生じること(P↑1 ,↓2 + P↓1 ,↑2 ) > (P↑1 ,↑2 + P↓1 ,↓2 ) が関係していること、言い換えれば、ハミルトニアンの特徴を反映して低エネルギー状態におけるスピン配置の実現確率に差が現れることが分かる。
この実現確率の差に注目して、確率の低い状態を捨てることで基底数を減らし、系を拡張しても基底数が増えないようにして基底数の発散の問題を回避することが、密度行列を用いた基底選択の意図である。ここで重要な点は、エネルギーの低い状態にもともと含まれない確率の《講義ノート》

物性研究・電子版 Vol. 7, No. 2, 072205 （2018年11月号）

低い状態だけを捨てるため、基底数を減らしても近似的に正しい基底状態の波動関数とエネルギーが得られることである。
このことを踏まえて、局所的なエネルギーの期待値の計算式をもう一度確認すると、σ1 が ↑ でσ2 のスピンが ↓ の基底 | ↑1 ⟩| ↓2 ⟩|j⟩ と、σ1 のスピンが ↓ で σ2 のスピンが ↑ の基底 | ↓1 ⟩| ↑2 ⟩|j⟩は、σ1 と σ2 の交換に対するハミルトニアンの対称性によって同じ局所的なエネルギーを与え、基底状態への寄与に関してこの部分だけでは優劣がつかないことが分かる。これは、局所的なハミルトニアンにこの二つの基底をつなぐ行列要素が含まれないことから生じる縮退の存在を示しており、このような場合に一方の基底だけを選ぶと対称性を破る基底状態が得られることになる。しかしここで、通常の交換相互作用に含まれるようなスピンの交換項がハミルトニアンの中に存在すれば、この縮退は解け、二つの基底の線形結合によってより低いエネル
ギーの状態が実現する。実際にハミルトニアンに交換項が含まれると局所的なエネルギーの
期待値は、

E1,2

= ⟨Ψ|J

=

+
−
+

}
{
1
S1z S2z + (S1+ S2− + S1− S2+ ) |Ψ⟩
2

(7)

)
J ∑(
|Ψ↑1 ,↑2 ,j |2 + |Ψ↓1 ,↓2 ,j |2
4 j

)
J ∑(
|Ψ↑1 ,↓2 ,j |2 + |Ψ↓1 ,↑2 ,j |2
4 j)J ∑( ∗
Ψ↑1 ,↓2 ,j Ψ↓1 ,↑2 ,j + Ψ∗↓1 ,↑2 ,j Ψ↑1 ,↓2 ,j
2 j

(8)

であるから、式 (8) の第 3 項に注目して Ψ↑1 ,↓2 ,j と Ψ↓1 ,↑2 ,j の間の位相差が π のときにエネル
ギーを下げることが分かる。このことから、低エネルギー状態を表現するときには、Ψ↑1 ,↓2 ,jの基底 | ↑1 ⟩| ↓2 ⟩|j⟩ と Ψ↓1 ,↑2 ,j の基底 | ↓1 ⟩| ↑2 ⟩|j⟩ の二つを独立に残す代わりに、それらの線形結合である √12 [| ↑1 ⟩| ↓2 ⟩|j⟩ − | ↓1 ⟩| ↑2 ⟩|j⟩] という基底を一つ残せば、エネルギーの低い状態が得られる。
このようにしてエネルギーの期待値と波動関数の係数 Ψσ1 ,σ2 ,j との関係を確認すると、局
所的エネルギーは波動関数 Ψσ1 ,σ2 ,j から得られる密度行列

ρ(σ1 ,σ2 )(σ1′ ,σ2′ ) =∑Ψ∗σ1 ,σ2 ,j Ψσ1′ ,σ2′ ,j(9)jの各要素 ρ(σ1 ,σ2 )(σ1′ ,σ2′ ) の関係によって定められ、エネルギーの低い状態に含まれる基底の特徴は、基底状態の密度行列 ρ の中に織り込まれていることが分かる。

1.2

密度行列を用いた基底選択と系の拡張

上に述べた密度行列の特徴に着目して、基底状態への寄与の大きな基底を最も効果的に求める方法を考えよう。密度行列の対角項は波動関数の各成分の絶対値の 2 乗であり、それは各

《講義ノート》

物性研究・電子版 Vol. 7, No. 2, 072205 （2018年11月号）

図 1: 基底選択を行う領域 (ブロック) の設定
基底が表現する状態の実現確率であるから、全ての対角項の和をとれば、全ての確率の和とな
り、規格化された波動関数であればその値 (密度行列のトレース) は 1 になっている。ここで、密度行列の非対角成分を 0 にする正規直交変換を行ない、密度行列を対角化すると、その対角要素の和であるトレースは変化しない一方で、非対角成分によって繋がる基底の対角成分同士は一般に分裂するように変化して固有値に移行する。すなわち、非対角要素を少しずつ大きくしていくと、大きい対角成分はより大きく、小さい対角成分はより小さい値へと移行する。
したがって、対角化の後に得られる最も大きな固有値に対応する変換後の基底を求めれば、その固有値は密度行列の定義によって実現確率であるから、最も実現確率が高い基底ということになる。つまり密度行列を求めてそれを対角化し、その固有値が大きな基底から順番に選び出
せば、限られた基底数の範囲内で基底状態を表現するための最も良い基底の選択ができることになる。
そこで、この基底選択の方法を利用して大きな系の基底状態を求めることを考える。基底の選択を行い基底数を制限することができれば、系を拡張して新しくスピンや電子軌道を加えても系全体の基底数を一定数以下に保つことができる。しかし、拡張前の基底状態を使って選んだ基底が、系の拡張の後の基底状態でも重要な基底になっていなければ、基底数の制限は単に誤差を生み出すだけであるから、拡張前の密度行列と拡張後の密度行列はほぼ同じになっていてほしい。したがって、基底選択を行う適切な領域について次に考えなければならない。
密度行列は波動関数によって定まるため、ほぼ同じ密度行列を得るためには拡張の前と後で波動関数の特徴が変化しない領域を選択する必要がある。微分方程式の解である波動関数の位相や振幅は境界条件によって定まるので、境界条件が設定されている場所を起点として、そこからの相対的な位置が移動しない領域を選べば良さそうである。また、基底選択を行なう
領域に隣接する新しい自由度を付け加える場所は、境界条件の影響が現れにくい場所に設定
すると、系を大きくしたときに並進対称な無限系の特徴を再現する波動関数に近づくと考えられる。したがって、図 1 に示すように、基底選択を行う領域の一端は境界条件を設定する
系の端に固定し、もう一端はその系の端から最も離れた場所に設定して、そこに新しい自由度を加えて拡張していくのが良さそうである。具体的には、図 2 に示すように、はじめに数値対角化が可能な小さな開放端の系を考え、その基底状態を求めた後、系を半分に分けて、それぞれの領域で基底選択を行う。すると、その領域の一端は境界条件によって振幅が固定され、もう一端はその系の境界から最も離れた中央部分になり、そこは一様な状態に最も近い場所になる。したがって、系を半分に分けてできる両側のブロックで基底選択を行なった後、2 つの《講義ノート》

物性研究・電子版 Vol. 7, No. 2, 072205 （2018年11月号）

図 2: 無限系の波動関数を求める方法 (infinite method)
ブロックの間の系の中央部分に新しい自由度を加えると、拡張部分は十分大きな系になったときに並進対称な状態に近づくと考えられる。この基底数の制限と系の拡張を繰り返すことで、いくらでも系を大きくしていくことができ、最終的に拡張の前後で系の中央部分の状態が変化しなくなれば、端の影響を受けない無限系の性質が得られる。このようにして無限系の波動関数を得る方法が infinite method である。
この方法の問題点は、基底選択の際に用いる波動関数と系の拡張後の波動関数がそれぞれサイズの異なる系の固有状態になっていることである。すなわち、同じ領域に対する密度行列を求めても、系の拡張の前後で波動関数は多少異なったものになるため、密度行列の固有値や固有ベクトルには差が現れてしまい、拡張前の系で最適化された基底は一般に拡張後の系では最適化されていない。そのため、infinite method において系の拡張ごとに密度行列の小さな固有値に対応する状態を捨てていくと、不完全な最適化により生じる誤差が累積してしまう。
この問題を解決するためには、系の拡張後に改めて基底選択をやり直し、再度最適化する必要がある。そのため、infinite method を用いて求めたいサイズに到達したら、図 3 に示すように、系の大きさを固定するために左右のブロックの一方のみを大きくし、もう一方は逆に小さくする。こうして２つのブロックの大きさを非対称にすることで、全体のサイズを固定したま
ま片側のブロックを最初の出発点まで小さくし、そこでブロックの基底をリセットする。ここから、今度はそれぞれ逆方向にブロックの大きさを変化させると、基底選択をやり直しながらブロックを作り直していくことができ、系の大きさを固定して基底選択をやり直すことが可能になる。重要な点は、常に同じ系の基底状態から密度行列を求めていることである。
このような操作を、左右のブロックの役割を交替させながら行い、それぞれの基底選択を最適化していくことで、最終的に各ブロックに残す状態数だけで定まる最もエネルギーの低い状態を得ることができる。最後に残る誤差は、各ブロックの拡張の際に部分的に基底を捨てたことで生じるもので、このときヒルベルト空間が狭くなっていることから、得られる最低エネル
ギーは厳密な基底状態のエネルギーより必ず高くなっている。したがって、基底状態のエネル
ギーの上限値がこの計算によって得られることになる。厳密な基底状態のエネルギーを下回ることはないので、残す状態数を増やしてエネルギーが下がれば、より真の基底状態に近づくことを意味し、その値が収束すれば真の値を確認したと判断できるだろう。

《講義ノート》

物性研究・電子版 Vol. 7, No. 2, 072205 （2018年11月号）

図 3: 系の大きさを固定して最適化する方法 (finite system method)1.3

密度行列を用いた基底選択の精度

実際にこのような基底の選択と系の拡張を行うためには、ハミルトニアンやそこに含まれる演算子の表現を基底変換に合わせて変更し、行列表現された演算子の次元を基底選択後の次元に合わせていかなければならない。密度行列の固有ベクトルは、対応する固有値 (実現確率) をもつ新たな基底への変換行列となるので、この固有ベクトルを使って各演算子の基底変換を行うことができる。いま、密度行列 ρ の α 番目に大きな固有値が wα であり、その規格化された固有ベクトルの i 番目の成分を vαi とすると、異なる固有値に属する固有ベクトルは直
交関係

∑

∗
vαi
vα′ i = δαα′

(10)

i

をみたす。したがって、vαi は変換前の基底 i から α 番目の固有値を与える状態への直交基底
B
変換行列を構成している。ブロックのハミルトニアン Hii′ は、この行列で基底変換すること

ができ
B
Hαα
′ =

∑

∗
HiiB′ vαi
vα′ i′(11)ii′のように密度行列の固有ベクトルの基底で表すことができる。この基底変換によってブロックのハミルトニアンの次元が変換前の基底 i の次元数から、残す固有値の数である α の最大値に変わる。例えば、α の最大値を m として、その数を i の基底数の半分に制限すれば、このブロックに新しく S = 1/2 のスピンを一つ加えても、付け加えるスピンには ↑ か ↓ の 2 通りの状態しかないので、基底の数は 2m になるだけであり、結局のところ、スピンを新たに加えてブロックを拡張しても、そのハミルトニアンの次元は変換前の基底 i の次元数 2m に戻るだけ
で次元の増加は生じない。したがって、この変換と拡張を繰り返せば、いくらでも系を大きくしていくことが可能になる。
このような演算子やハミルトニアンの基底変換を行う際に考慮すると良い点は、ハミルトニ
アンに対称性があり保存量子数が存在する場合、ハミルトニアンの固有状態や基底はその保《講義ノート》

物性研究・電子版 Vol. 7, No. 2, 072205 （2018年11月号）

図 4: 1 次元近藤格子模型 (J/t = 1, nc = 1) の場合の密度行列の固有値分布
存量子数によって分類可能になっていることである。例えばハイゼンベルグ模型では、全スピ
ンとその z 成分が保存するため、固有状態を構成する基底の全スピンの z 成分は全て同じ大きさになっている。この事情はブロックのハミルトニアンに対しても同様であり、ハミルトニ
アンの非対角行列要素は同じ量子数の基底の間にしか存在しない。したがって、全ての基底を
この保存量子数によって分類しておけば、ハミルトニアンや演算子の表現をコンパクトな形にまとめることができ、対称性によって 0 になることが保証されている行列要素をはじめから省いておくことができる。また、異なる量子数の状態がハミルトニアンによって混ざらないということは、その固有状態になっている波動関数の基底も異なる量子数の間で混ざることはな
く、密度行列の対角化も、個々の量子数毎に別々に行うことができる。このような工夫は数値計算の精度や速さを大幅に向上させるため、実際に計算を行う際には大変重要になる。
以上のように、基底の選択と演算子の変換を行うことで系を大きくしていくことができる
が、このような基底数の制限を行うことで生じる誤差についてどのように評価すればよいか
ここで考えてみよう。基底数の制限によって波動関数に導入される誤差は、捨ててしまった基底を本来の基底状態がどの程度含んでいたか、すなわち、捨てた基底が表現する状態空間の中に基底状態が観測される確率によって生じると考えることができる。それは、もともと基底状態の中に含まれていたけれども捨ててしまった基底の実現確率であり、捨てた基底が持っていた密度行列の固有値の大きさそのものである。したがって、m 個の基底を残したときに生じる誤差 Em は、捨てた基底に対応する基底状態の密度行列の固有値の和

Em =

α∑
max
α=m+1

wα = 1 −

m
∑

wα

(12)

α=1

によって与えられ、それが切り捨てた実現確率となって誤差の原因になると考えることができる。実際にどの程度の誤差が生じるかは、密度行列の固有値分布によって決まることになる。
一例として、近藤格子模型 (J/t = 1) の基底状態の場合で密度行列の固有値分布を確認する

《講義ノート》

物性研究・電子版 Vol. 7, No. 2, 072205 （2018年11月号）と、図 4 に示すように、64 電子の場合でも wα は速やかに減衰し、100 個程度の状態数を残せ
は波動関数のノルムに生じる誤差は 10−6 程度に抑えられることが分かる。局在スピン 32 個と遍歴電子の軌道 32 個からなるこの系の本来の基底自由度の数は 232 × 432 = 296 ∼ 1029 であるから、極めて効率の良い基底選択ができていることが分かる。

1.4

波動関数の特異値分解と密度行列の関係

このようにして得られる波動関数の構造について、ここで改めて考えてみる。初期状態として系全体の基底が左右のブロックの基底 |l⟩ と |r⟩ の直積 |l⟩|r⟩ で表現される場合を考えると、基底状態の波動関数 |Ψ⟩ は、一般に

|Ψ⟩ =

∑

Ψlr |l⟩|r⟩

(13)

lr

のように表されることになる。係数行列 Ψlr は、波動関数に含まれる各基底の成分の大きさを行列形式で表したものであり、波動関数の規格化条件によって⟨Ψ|Ψ⟩ =∑

Ψ∗lr Ψlr = 1

(14)

lr

が成り立っている。ここで、この係数行列 Ψlr を特異値分解1 すると Ψlr =∑α ulα dαα vαr と表現することができ、波動関数を

|Ψ⟩ =

∑

Ψlr |l⟩|r⟩ =∑∑lrulα dαα vαr |l⟩|r⟩

(15)

α

lr

のように書き直すことができる。ここで ulα , vαr は特異値分解したときの直交変換行列であり、それぞれ∑u∗lα′ ulα= δα′ α(16)vα∗ ′ r vαr= δα′ α

(17)

l

∑
r

の関係をみたすため、|l⟩, |r⟩ から |uα ⟩, |vα ⟩ への基底変換を|uα ⟩ =∑

ulα |l⟩(18)vαr |r⟩(19)l|vα ⟩ =∑rによって定義すると

|Ψ⟩ =

∑
lr

Ψlr |l⟩|r⟩ =∑dαα |uα ⟩|vα ⟩

(20)

α

1 特異値分解は非対称行列を含む任意の行列 A に対して、2 種類の直交変換行列 U, V と特異値を対角項とする対角行列 D を用いて A = U DV の形に表現することであり、A が対称行列であれば A の対角化 (固有値分解) と等価になる。

《講義ノート》

物性研究・電子版 Vol. 7, No. 2, 072205 （2018年11月号）となる。これは |uα ⟩, |vα ⟩ で表現される状態が dαα の大きさで波動関数に含まれていることを表しており、波動関数を左右のブロックの基底を使って行列表現し、その行列を特異値分解す
ると、波動関数を構成する各ブロックの主要な成分を特異値 dαα の大きさを基準にして選び出せることを示している。すなわち、dαα の大きなものに対応する左右のブロックの基底 |uα ⟩と |vα ⟩ を利用すれば、波動関数の主要部分は表現できるのである。
では、この特異値 dαα と基底変換ベクトル |uα ⟩ と |vα ⟩ は、波動関数の密度行列とどのよRうな関係になっているのだろう。いま、左右のブロックの密度行列 ρLl,l′ , ρr,r′ をそれぞれの基

底 l, r を使って表すと

ρL
l,l′

∑

=

Ψ∗lr Ψl′ r(21)

Ψ∗lr Ψlr′

(22)

r

ρR
r,r ′

∑

=

l

であるから、特異値分解した表式 Ψlr =∑ ∗r vα′ r vαr = δα′ α であることを用いて

ρL
l,l′

=

∑

∑

α ulα dαα vαr を代入すると、∑∗
l ulα′ ulα = δα′ α 、Ψ∗lr Ψl′ r

r

=

∑∑
r

=

∑

∗
u∗lα d∗αα vαrul′ α′ dα′ α′ vα′ rαα′u∗lα d∗αα ul′ α′ dα′ α′ δα′ α

αα′

=

∑

u∗lα |dαα |2 ul′ α

(23)

α

ρR
r,r′

=

∑

Ψ∗lr Ψlr′

l

=

∑∑
l

=

∑

∗
u∗lα d∗αα vαr
ulα′ dα′ α′ vα′ r′αα′∗
δα′ α d∗αα vαrdα′ α′ vα′ r′

αα′

=

∑

∗
vαr
|dαα |2 vαr′

(24)

α

と書くことができる。この式は、密度行列 ρL , ρR を対角化するときの表式 ρL = U † DU , ρR =V † DV と同じであり、そのことから、密度行列の固有値 (対角行列 D の対角要素) が特異値
dαα の絶対値の 2 乗になり、その固有ベクトル (V, U の行ベクトル) が特異値分解したときの基底変換行列になっていることがわかる。つまり、密度行列を求めて対角化することで、特異値分解をせずにその特異値と変換行列が得られていたのである。

《講義ノート》

物性研究・電子版 Vol. 7, No. 2, 072205 （2018年11月号）

1.5

行列積型変分波動関数

波動関数と基底変換行列の関係が明らかになったので、この変換行列を使って系を拡張した
ときに得られる波動関数の表式について最後に確認してみよう。まず、出発点となる 4 スピンの波動関数 |Ψ(4)⟩ を i 番目のスピン自由度 σi (例えば S = 1/2 量子スピンであれば σi = {↑, ↓})を使って|Ψ(4)⟩ =∑Ψ(4)(σ1 ,σ2 ),(σ3 ,σ4 ) |σ1 ⟩|σ2 ⟩|σ3 ⟩|σ4 ⟩(25)

σ1 ,σ2 ,σ3 ,σ4のように表し、左右のブロックの基底変換 (σ1 , σ2 ) → α2 と (σ3 , σ4 ) → α2′ を、それぞれ変換
行列 u(4)σσ21 ,α2 と v(4)σα3′ ,σ4 を使って定義する。特異値行列を d(4)α2 とすると、特異値分解し
2

た表式から

Ψ(4)(σ1 ,σ2 ),(σ3 ,σ4 )∑=u(4)σσ21 ,α2 d(4)α2 v(4)σα32 ,σ4

(26)

α2

と書けることが分かる。6 スピンの波動関数 |Ψ(6)⟩ は 4 スピンのときの左側ブロックの基底

α2 と新しく加えた 3 番目のスピン自由度 σ3 、および 4 スピンのときの右側ブロックの基底 α2′と新しく加えた 4 番目のスピン自由度 σ4 を使って Ψ(6)(α2 ,σ3 ),(σ4 ,α′2 ) で表現できるが、これをやはり特異値分解して左右のブロックの基底変換 (α2 , σ3 ) → α3 と (α2′ , σ4 ) → α3 を変換行列

u(6)σα32 ,α3 と v(6)σα43 ,α′ を使って表し、特異値行列 d(6)α3 を用いると2

Ψ(6)(α2 ,σ3 ),(σ4 ,α′2 )∑=u(6)σα32 ,α3 d(6)α3 v(6)σα43 ,α′

(27)

2

α3

を得る。ここで、左右のブロックの α2 と α2′ の基底を元の σ1 , σ2 および σ5 , σ6 を使って表すと、変換行列は u(4)σσ21 ,α2 と v(4)σα5′ ,σ6 であったから、2

Ψ(6)(σ1 ,σ2 ,σ3 ),(σ4 ,σ5 ,σ6 )=∑u(4)σσ21 ,α2 u(6)σα32 ,α3 d(6)α3 v(6)σα43 ,α′ v(4)σα5′ ,σ62α2 ,α′2 ,α3

(28)

2

と書き直すことができる。そして、8 スピンの波動関数 |Ψ(8)⟩ は 6 スピンのときの左側ブロックの基底 α3 と新しく加えた 4 番目のスピン自由度 σ4 、および 6 スピンのときの右側ブロックの基底 α3′ と新しく加えた 5 番目のスピン自由度 σ5 を使って Ψ(8)(α3 ,σ4 ),(σ5 ,α′3 ) で表現でき、やはり特異値分解して左右のブロックの基底変換 (α3 , σ4 ) → α4 と (α3′ , σ5 ) → α4 の変換行列

u(8)σα43 ,α4 と v(8)σα54 ,α′ と、特異値行列 d(8)α4 を用いて
3

Ψ(8)(α3 ,σ4 ),(σ5 ,α′3 )∑=u(8)σα43 ,α4 d(8)α4 v(8)σα54 ,α′

(29)

3

α4

と書き表すことができる。ここで、また、左右のブロックの α3 と α3′ の基底を元の α2 , σ3 お
よび σ6 , α2′ を使って表すと、変換行列は u(6)σα32 ,α3 と v(6)σα6′ ,α′ であるから

Ψ(8)(α2 ,σ3 ,σ4 ),(σ5 ,σ6 ,α′2 )

=

∑

3

2

u(6)σα32 ,α3 u(8)σα43 ,α4 d(8)α4 v(8)σα54 ,α′ v(6)σα6′ ,α′α3 ,α′3 ,α4

3

3

2

(30)

《講義ノート》

物性研究・電子版 Vol. 7, No. 2, 072205 （2018年11月号）となり、さらに、左右のブロックの α2 と α2′ の基底を元の σ1 , σ2 および σ5 , σ6 を使って表すと、変換行列 u(4)σσ21 ,α2 と v(4)σα7′ ,σ8 を使って
2

Ψ(8)(σ1 ,σ2 ,σ3 ,σ4 ),(σ5 ,σ6 ,σ7 ,σ8 ) =∑u(4)σσ21 ,α2 u(6)σα32 ,α3 u(8)σα43 ,α4 d(8)α4 v(8)σα54 ,α′ v(6)σα6′ ,α′ v(4)σα7′ ,σ83α2 ,α′2 ,α3 ,α′3 ,α4

3

2

(31)

2

のように表すことができる。
このようにして 2n スピンの系の波動関数を変換行列を使って表現すると

Ψ(2n)(σ1 ,σ2 ,··· ,σn ),(σn+1 ,···σ2n−1 ,σ2n ) =
∑

σ

σ

σ

2n−1
u(4)σ1 ,α2 u(6)σα32 ,α3· · · u(2n)σαnn−1 ,αn d(2n)αn v(2n)αn+1· · · v(6)α2n−2′′ v(4)α′ ,σ
n ,α
3 ,α
2n

α2 ,α′2 ,α3 ,α′3 ,··· ,αn

n−1

2

2

(32)
となる。こうしてみると、多体の波動関数が u(2m) や v(2m) で表される変換行列の積で表現されていることが分かる。このような波動関数は行列積型波動関数 [11, 12] と呼ばれるものであり、密度行列繰り込み群によって得られる波動関数は結局のところ、有限の次元数で打ち切った行列の積の形で波動関数を近似的に表現したものであることが分かる。基底変換の最適化はこの行列積型の変分波動関数の最適化と同じであり、密度行列繰り込み群の変分法としての側面が理解できる。

有限温度量子系への拡張

2

前章では、密度行列繰り込み群による量子 1 次元系の基底状態の計算法を説明したが、ここでは、密度行列繰り込み群の応用例として、有限温度における熱力学量を求める計算法を解説する。第１章で述べた基底状態の計算では、ハミルトニアンの行列次元を制限しながら系を実空間方向に拡張し、最終的に大きな系の基底状態を求めた。しかしながら、行列次元の制限は基底状態に寄与しない状態の排除によって行われているため、そこで得られたハミルトニアン
を用いてエネルギーの高い励起状態を正しく得ることは難しい。このことは、熱的に励起された状態が多数現れる有限温度の計算が困難であることを意味するが、実は量子転送行列を利用
すると密度行列繰り込み群を用いて無限系の有限温度の計算を正確に行うことが可能になる。

2.1

量子転送行列法

はじめに、この量子転送行列法について説明する。熱力学的性質を理論的に明らかにするためには、熱平衡状態における各固有状態の重みを決める分配関数が必要になる。大抵の熱力学

《講義ノート》

物性研究・電子版 Vol. 7, No. 2, 072205 （2018年11月号）

量はこの分配関数の微分によって求めることができるが、分配関数 Z は系のハミルトニアン

H を指数部に含む指数関数のトレース}{= Tr e−βHZ(33)であるから、相互作用が働く多体系の場合には、ハミルトニアンの基底数が粒子数に対して指数関数的に増大するため、極めて大きな行列の指数関数を求めなければならない。そのため、一般に全ての固有値が得られていない限り、まともに計算するのは困難となる。そこで、ハミルトニアンを 2n や 2n + 1 を空間座標とする局所ハミルトニアン h2n,2n+1 と h2n+1,2n+2 に分解し、[h2n,2n+1 , h2n′ ,2n′ +1 ] = [h2n+1,2n+2 , h2n′ +1,2n′ +2 ] = 0
を満たすように Hodd =∑L/2−1n=0

h2n+1,2n+2 と Heven =∑L/2−1n=0(34)

h2n,2n+1 に分類して、さらに Hodd と Heven の非可換性から生じる補正項を微小量に抑えるように十分大きな分割数 NT
で逆温度 β = 1/(kB T ) を分割して、分配関数 Z をZ= Tr e−βH=[

lim Tr e−βHodd /NT e−βHeven /NT

NT →∞



∏

L/2−1

=

lim Tr 

NT →∞

(35)

]NT

(36)
∏

NT

L/2−1

e−βh2n+1,2n+2 /NTn=0e−βh2n,2n+1 /NT 

(37)

n=0

のように小さな行列−βhe(τj+1i,i+1
,τj )

/NT

≡

∑

∑

∑ ∑

′
′
σi,τj+1 σi+1,τj+1 σi,τσi+1,τ
jj′′|σi,τj+1 σi+1,τj+1 ⟩⟨σi σi+1 |e−βhi,i+1 /NT |σi′ σi+1⟩⟨σi,τ
σ′|j i+1,τj(38)の積の形で表す鈴木・トロッター分解 [14, 15, 16] を利用する。ここで τi は、逆温度 β を分割数 NT で等間隔に分割したことで生じる逆温度軸上の座標を表し、その間隔を τi+1 − τi =β/NT = β0 と表す。この分割によって、分配関数は図 5 に示す対角化が可能なサイズの局所
ハミルトニアン hi,i+1 を指数部にもつ指数関数 e−β0 h2n,2n+1 および e−β0 h2n+1,2n の積のトレースとして表現され、ハミルトニアンが並進対称であれば、この指数関数の積も並進対称であり、heven = h2n,2n+1 , hodd = h2n+1,2n として逆温度 β = 1/(kB T ) の方向に NT 分割された

1 列分の指数関数の積
T

=

lim

M →∞

N∏
T −1

−βh

/N

−βh

/N

even
T
T
odd
(e(τ2j+2,τ2j+1 ) e(τ2j+1 ,τ2j ) )j=0

を転送行列とみなして分配関数 Z をこの転送行列 T の積のトレースとして(39)

《講義ノート》

物性研究・電子版 Vol. 7, No. 2, 072205 （2018年11月号）

図 5: (a) 鈴木・トロッター分解された分配関数、(b) 式 (38) に対応する行列{}
Z = Tr T L/2(40)のように表すことができる。この転送行列 T は量子転送行列 [17, 18, 19] と呼ばれるものであり、Z の値はこの量子転送行列の固有値を L/2 乗したものの和で表され、L が十分大きい場
合、最大固有値からの寄与は最大固有値以外からの寄与と比較して十分大きくなり、L が無限大になる極限では転送行列の最大固有値 λmax だけで

Z=

∑

L/2

λi

= λL/2max(L → ∞)(41)iのように定まることになる。このことは、無限系の有限温度の分配関数が量子転送行列の最大固有値だけで決まることを意味し、例えば自由エネルギー F がF= − ln Z/(Lβ)(42)= − ln λmax /(2β)(43)と表現できるように、有限温度における熱力学量の計算は、量子転送行列の最大固有値を求める計算に置き換えられることになる。
量子転送行列は逆温度 β = 1/(kB T ) の方向に NT 分割されている指数関数の行列 e−βheven /NTおよび e−βhodd /NT の積で表現されている。そのため NT に比例する分割点においてこの hodd
や heven の中に含まれる局所自由度 σi,τj などが現れ、NT の増加に対して量子転送行列の次元は指数関数的に増大する。この問題はハミルトニアンを実空間方向に拡張したときに、その行列次元が電子数の増加に伴い指数関数的に増大するときの問題と同じであり、直接量子転送行
列を対角化して分配関数を求めようとすると、行列次元の制約によって NT を大きくすることができず、β を大きくしたときに β/NT を十分小さくできない。β/NT を十分小さくしなけれ《講義ノート》

物性研究・電子版 Vol. 7, No. 2, 072205 （2018年11月号）

ば鈴木・トロッター分割で無視した非可換項の影響が無視できなくなるため、分配関数の信頼性が失われてしまう。そこで、基底状態の波動関数を求めるときと同じように密度行列を使って基底の選択を行い、量子転送行列の NT を少しずつ増やして最終的に NT の大きな量子転送行列の最大固有値を得るのが有限温度密度行列繰り込み群 [5, 6, 7] である。

2.2

密度行列繰り込み群の適用

分配関数の計算に必要なものは量子転送行列の最大固有値であるから、最大固有値を与える
固有ベクトルの密度行列を求め、その対角化によって重要な基底だけを選び出せばよい。ただし、チェッカーボード型に鈴木・トロッター分解した分配関数の量子転送行列は、左右の反転対称性がなく、その行列はエルミートではない。そのため、その固有ベクトルは左側から作用させる VjL′ と右側から作用させる VjR で∑

VjL Tjj ′= λmax VjL′

(44)

Tjj ′ VjR′

= λmax VjR

(45)

j

∑
j′

のように異なったものになり、左右のベクトルの規格化条件は、左右のベクトルの内積 ⟨VL |VR ⟩ =∑ L R
j Vj Vj = 1 によって定義される。これは固有ベクトルが左右で複素共役の関係になるエル∗ミート行列の場合 vjL = vjR と異なるものであり、この非エルミート性により密度行列

ρA
kk′ =

∑

VklL VkR′ l

(46)

l

も非対称になり、その対角化によって得られる固有ベクトルも左右で異なるものになる。したがって、チェッカーボード型に鈴木・トロッター分解した量子転送行列の基底ベクトルは左右で異なったものになり、その規格直交性も左右の基底ベクトルの内積によって定義される。このように有限温度の計算では非対称性の取り扱いが基底状態を求めるときと異なるが、基本的には局所自由度を加えてハミルトニアンを実空間方向に拡張するときと同じように量子転送行列を逆温度軸方向に拡張して分配関数を求めることができる。
量子転送行列を用いた計算法では、鈴木・トロッター分解によって十分小さな逆温度 β0 で定義される式 (38) の行列を逆温度軸方向に追加していくため、分配関数全体の β は逆温度軸方向の要素数を N としたとき、β = N β0 のように N に比例するように増加する。そのため、分配関数の温度は量子転送行列の拡張とともに低下し、高温から低温までの温度変化を量子転送行列の拡張の過程で得ることができる。このメリットは、比熱のような温度に対する変化率が必要な場合や、物理量の温度変化を計算する場合に生かすことができる。
局所的な相関関数や密度のような物理量の期待値は転送行列の最大固有値を与える固有ベ
クトルを使って直接計算することができ、また 3 サイト以内の演算子によって定義される物理《講義ノート》

物性研究・電子版 Vol. 7, No. 2, 072205 （2018年11月号）

量や相関関数であれば転送行列と最大固有値の固有ベクトルを用いて⟨S0z S1z ⟩ ≡zzTr{e−βH S0z S1z }/Z = ⟨V L |S0,τ
T S1,τ
|V R ⟩/λmax
0
0

(47)

⟨S0z S2z ⟩ ≡zz
T S2,τ
|V R ⟩/λmax
Tr{e−βH S0z S2z }/Z = ⟨V L |S0,τ
0
0

(48)

のように求めることができる。したがって、エネルギーに寄与する相関関数を求め、内部エネ
ルギーを直接求めて、その温度微分から比熱を計算することもできる。比熱の計算は自由エネ
ルギー

F (T ) = − ln λmax /(2β)(49)の 2 階微分C(T ) = −T

∂2
F
∂T 2

(50)

から求めることもできるが、内部エネルギーの微分から求めた方が微分の階数が少ないため、数値的に安定している。また、帯磁率や電荷感受率についても自由エネルギーの 2 階微分

χs
χc

∂2
F
∂hz 2
∂2
= − 2F
∂µ
= −

(51)
(52)

から求めることもできるが、微小な外場を加えた状況で磁化や密度の期待値⟨S0z ⟩ ≡z
Tr{e−βH S0z }/Z = ⟨V L |S0,τ
|V R ⟩
0(53)⟨n0 ⟩ ≡

Tr{e−βH n0 }/Z = ⟨V L |n0,τ0 |V R ⟩(54)を求め、外場がない状況からの差分

χs

=

χc

=

lim ⟨S0z ⟩/hz

(55)

lim (⟨n0 ⟩µ+∆µ − ⟨n0 ⟩µ )/∆µ

(56)

hz →0

∆µ→0

によって求めることもできる。例えば、スピンの反転対称性がある場合のように外場がない状
況での値が対称性により 0 になる場合には、一回の計算で帯磁率や感受率を得ることができるため、分配関数の微分から求めるより効率が良い。図 6 に 1 次元近藤格子模型のスピン液体相にホールをドープしたときの電荷感受率の温度依存性を示すが、伝導電子のバンド幅 (4t)の 1/100 の温度まで計算できることが分かる。
また、量子転送行列の最大固有値を与える固有ベクトルからは、虚時間軸に対応する逆温度
軸方向の相関関数

χAB (τ2j ) ≡ Tr{e−βH Ai (τ2j ) Bi (0)}/Z= ⟨V L |Ai (τ2j ) Bi (0)|V R ⟩(57)

《講義ノート》

物性研究・電子版 Vol. 7, No. 2, 072205 （2018年11月号）

図 6: 1 次元近藤格子模型の絶縁相 (スピン液体相) にホールをドープしたときの電荷感受率の温度依存性 [20]。縦軸上の黒丸は通常の密度行列繰り込み群で求めたドープ量 δ が 0.3 のときの基底状態の電荷感受率。
を求めることができる。この相関関数のフーリエ変換

β ∑ iωn τ2jeχAB (τ2j )M j=0MχAB (iωn ) =(58)

によって松原振動数 ωn で表した相関関数が得られるので、その相関関数を解析接続すること
で実時間の相関関数に変換することができ [21]、得られる動的相関関数の ω および温度依存性から、準粒子状態密度やエネルギーの吸収、放出に関わるスペクトルとその温度変化を知ることができる。解析接続の方法としてはパデ近似や最大エントロピー法といった計算法があり、特に解析接続の精度が確保できる低エネルギー領域でその力を発揮する。図 7 に最大エントロピー法で求めた 1 電子状態密度の温度依存性を示すが、近藤絶縁体と呼ばれるスピン液体
相にホールをドープしたときの低エネルギー領域の特徴的な変化が確認できる。低温における高エネルギー領域の動的相関関数の波数や周波数依存性については、有限系における基底状態からの励起の時間発展を計算し、境界条件の影響を受けない範囲での短時間の相関関数の時間依存性から求めることができる。このとき必要になる時間発展についても密度行列繰
り込み群を応用して計算することができ、大きな系での解析が可能になっている [10]。

《講義ノート》

物性研究・電子版 Vol. 7, No. 2, 072205 （2018年11月号）

図 7: 1 次元近藤格子模型の絶縁相 (スピン液体相) にホールをドープしたときの 1 電子状態密度の温度依存性 [21]3

量子 2 次元系への応用
この章では、密度行列繰り込み群の量子 2 次元系への適用について考える。密度行列繰り込

み群の本質は、基底の変換と選択を繰り返して系を拡張するところにあり、その基底変換と選択の特徴は最終的に得られる行列積型の波動関数に反映される。したがって、2 次元系において精度を確保するためには、この一連の基底変換の流れに応じて形成される波動関数の最適化をしなければならない。そのため、この行列積型の波動関数に適した 2 次元の自由度の 1 次
元化をはじめに考えたいが、行列積型の波動関数の一般的な特徴は等比数列のような相関関数の単調な減衰であるから、1 次元化した後の自由度の間での相関が短距離型になる自然な一次元化をしなければならない。このような 1 次元化は一般には簡単ではないが磁場中の 2 次
元電子系では比較的容易に実現できる。
半導体の界面やグラフェンなどの原子層に閉じ込められた 2 次元電子系に垂直磁場を加えると、電子の運動はローレンツ力を受けることでサイクロトロン運動に移行する。その閉じた電子軌道に課せられる量子化条件は古典的サイクロトロン半径を量子化し、その結果として
運動エネルギーが離散化され、ランダウ準位が形成される。このときのランダウ準位の指数をN とし、ベクトルポテンシャルの表現としてランダウゲージを用いると、1 電子状態は][][(x − X)2
x−X(59)
ΨN X = CN exp iky y −
H
N
2ℓ2
ℓ
のように、N 次のエルミート多項式 HN (x)、規格化因子 CN 、磁気長 ℓ = (h̄c/eH)1/2 、サイクロトロン運動する電子の中心座標 X によって表現される。ここで、ky は y 方向の運動を表す波数であり X = ky ℓ2 の関係式によって中心座標 X と結び付いている。このことから、y方向の長さ Ly の周期境界条件により y 方向の波数が ky = 2πn/Ly のように離散化されると、

《講義ノート》

物性研究・電子版 Vol. 7, No. 2, 072205 （2018年11月号）

図 8: 最低ランダウ準位の占有率 ν = 1/3 の分数量子ホール状態の電子対相関関数

X = ky ℓ2 の関係によって中心座標も(n = 0, 1, 2, 3, · · · , M − 1)

Xn = 2πℓ2 n/Ly(60)のように離散化され、連続自由度の系が離散化された格子自由度の系に変換される。ここで、M はランダウ準位の縮退度 (1 電子軌道の数) で 2πM ℓ2 = Lx Ly によって系の面積 Lx × Lyと関係している。したがって、全ての 1 電子状態はランダウ準位の指数 N と同一ランダウ準位内の離散化された量子数 Xn によって一意に指定でき、電子が占有しているランダウ準位の指数を固定すれば、全ての 1 電子状態は離散化された有限個の中心座標 Xn の値だけで指定できるようになる。
この量子数 Xn は、まさに 1 次元的に並べられた格子系の電子軌道の中心座標と等価に扱うことができ、この量子数を用いてハミルトニアンを表現すれば短距離の相関が強く現れる 1次元系への変換が可能になる。電子間の相互作用は、離散化された電子の X 座標の指数 n を
使って

H=S

∑
n

c†n cn +1 ∑∑∑∑
An1 n2 n3 n4 c†n1 c†n2 cn3 cn4
2 n n n n
1

2

3

(61)

4

のように表現することができる。ここで、一体のエネルギー S は、x,y 方向に周期境界条件を課すことで生じる Lx × Ly の格子によって定義される古典的なウィグナー結晶のエネルギー[22] であり、c†n は、中心座標が Xn の軌道の電子の生成演算子、An1 n2 n3 n4 は 2 電子間のクー
ロン相互作用の行列要素

An1 n2 n3 n4

1 ∑ ′
2πe2
δn1 −n4 ,qy Ly /2π
Lx Ly q
ϵq
[ 2 2
]
[
]2
q ℓ
q x Lx
LN (q 2 ℓ2 /2) exp −− i(n1 − n3 )2M

= δn′ 1 +n2 ,n3 +n4(62)になっている [22]。ここで、LN (x) はランダウ準位に依存するラゲール多項式であり、N がランダウ準位の指数になる。δn′ 1 ,n2 は n1 = n2 (mod M ) のときに 1 になりそれ以外の場合は0 となる。

《講義ノート》

物性研究・電子版 Vol. 7, No. 2, 072205 （2018年11月号）

図 9: N = 2 ランダウ準位の電荷密度波状態の電子対相関関数 [9]このようにして一次元化した磁場中 2 次元電子系の基底状態を、通常の密度行列繰り込み群によって求めると、図 8 に示すように、最低ランダウ準位における占有率 ν = 1/3 での分数量子ホール液体 [23] や、図 9 に示す、高次のランダウ準位におけるウィグナー結晶、バブル
相、ストライプ相 [9, 24, 25] といった多彩な秩序状態の出現を確認することができる。

局所エネルギー変換を用いた外場応答の計算

4

密度行列繰り込み群を用いて基底状態を求めると、行列積型の波動関数が得られることを第1 章で示した。このことは、長距離の相関関数は全て漸近的に指数関数でスケールされ、臨界的なベキ乗則に従う長距離揺らぎは最終的に失われることを意味する。そのため、そのような近似法に依存した特徴が現れる前に系の拡張を止め、行列次元に対する外挿を行うことで臨界的な長距離揺らぎが現れるかどうか有限の大きさの系で確認しなければならない。しかしながら一方で、有限の長さに閉じ込められた系では量子化に伴うエネルギー準位の離散化によって磁化や電荷密度の外場応答に境界条件に依存した構造が現れ、特異性のない連続的な応答が得られる無限系の特徴を再現することが難しくなる。ここでは、このような境界条件に依存するエネルギー準位の離散化が作り出す有限サイズ効果を抑制し、無限系の外場応答を効率的に求める方法について紹介する。
有限系の境界条件や離散的な保存量子数によって生じる離散準位のサイズ効果は、局所エネ
ルギスケール変換の一つである SSD (Sine Square Deformation)[26, 27, 28, 29] を適用することで極めて効果的に抑制することができる。SSD と呼ばれるエネルギースケール変換は、系を構成するスピンや軌道に端から順番に j = 1, 2, 3, · · · , L と番号を付けて、サイン 2 乗の関
数形{
f (j) = sin2}1π(j − )
L
2(63)

でハミルトニアンのスケールを j に応じて少しずつ滑らかに変化させる変換である。この関数
は、系の中央 j = (1 + L)/2 でオリジナルのエネルギーの大きさを保ち、輪を作るように両端

《講義ノート》

物性研究・電子版 Vol. 7, No. 2, 072205 （2018年11月号）-0.5780

S =1/2 Spin Ladder
J⊥/J = 1.0

E/L-0.578100.00001

0.000020.00003

1/L3

図 10: SSD を適用したときのエネルギー密度のサイズ依存性を仮想的に繋ぐボンドの中心 j = L + 1/2 でエネルギーを 0 にする。周期境界条件をみたすハ
ミルトニアンにこの関数をかけると L 番目と 1 番目を繋ぐボンドが切れて、開放端条件のハ
ミルトニアンに変換されるところがポイントである。この SSD を適用したハミルトニアンを用いると、一様な基底状態が得られないように思えるが、実は並進対称な周期境界条件の基底状態を正確に得ることができる。
実際にこのことを自由電子の系

H = −tL−1∑(c†j+1 cj + c†j cj+1 ) − µjL∑c†j cj(64)

jで確認してみよう。j 番目と j + 1 番目を繋ぐ遷移要素 t には f (j + 1/2) のスケール因子がかかることに注意すると、SSD を適用することでハミルトニアンは

HSSD

= −t

L−1
∑

sin2

j

= −t

L−1
∑
j

1
2

}
{
L
{π }
∑
π
1
j (c†j+1 cj + c†j cj+1 ) − µsin2(j − ) c†j cj
L
L
2
j

(65)

{
}
L
∑
1
2π
1
1 − cos(j − ) c†j cj
2
L
2
j

(66)

{
}
2π
1 − cosj (c†j+1 cj + c†j cj+1 )
L
−µ

= −t

L−1
∑
j

1
2

{
}
2π
1 2π
1 − (ei L j + e−i L j ) (c†j+1 cj + c†j cj+1 )
2
−µ

{
}
L
∑
1
2π
1
1
1 2π
1 − (e L (j− 2 ) + e− L (j− 2 ) ) c†j cj
2
2
j

(67)

《講義ノート》

物性研究・電子版 Vol. 7, No. 2, 072205 （2018年11月号）となる。このハミルトニアンのフーリエ変換は

1 ∑ ik(j−1/2)
√
e
cj
L j
1 ∑ −ik(j−1/2)
√
e
ck
L k
L

ck

=

cj

=

(68)
(69)

と δ = 2π
L の定義を用いて

HSSD

=

1∑
(−2t cos k − µ)c†k ck
2
k
{
(
)
}
δ
1 ∑ −iδ/2−2t cos k ++e− µ c†k+δ ck
2
2
k
{
(
)
}
δ
1 ∑ iδ/2−2t cos k −e− µ c†k−δ ck
+
2
2

(70)

k

となり、SSD によって生じた式 (70) の第 2 項と第 3 項は個々の電子の波数を隣の波数に移す
遷移要素を生成することが分かる。ここで、本来の並進対称性が残る第 1 項のハミルトニアン
の基底状態

|G⟩ =

∏

ck |0⟩(71)|k|≤kFを考え、その状態に第 2 項を作用させると、|k| ≤ kF の状態には電子が完全に詰まっている
ため、パウリ原理によってほとんどの遷移要素は消失し、唯一フェルミエネルギー直下の電子がフェルミ面直上の状態に遷移する要素だけが残る。そのため、基底状態において電子が−kF から kF までの状態を占有しているとすると、第 2 項から kF → kF + δ の遷移要素が、第 3 項から −kF → −kF − δ の遷移要素がそれぞれ一つずつ現れることになる。しかし、この 2 つの項にはそれぞれ −2t cos (kF + δ/2) − µ と −2t cos (−kF − δ/2) − µ の係数が付いているため、化学ポテンシャル µ を µ = −2t cos(kF + δ/2) に選ぶことで、この二つ要素を同時に消すことができる。結局、HSSD を第一項の基底状態 |G⟩ に作用させても第一項の係数∑1|k|≤kF (−2t cos k − µ) だけが固有値として現れるだけで、並進対称な第一項の基底状態 |G⟩
2がそのまま HSSD の固有状態になる。このように、開放端境界条件で計算しているにもかかわらず周期境界条件を適用した解が得られる SSD の計算法は、開放端境界条件で計算を行う密度行列繰り込み群を利用して、境界条件の影響が少ない並進対称な基底状態を得るために利用することができる。
開放端の境界条件でも、並進対称な波動関数が得られるということは、開放端の存在によって生じる基底エネルギーの有限サイズ効果も抑制されていると考えられる。実際に SSD を適用して基底状態のサイト当たりのエネルギーを計算すると、図 10 に示すように、通常現れる

1/L2 に比例する開放端条件に起因する有限サイズの補正が、1/L3 の補正量に変化して、サイズ L の増加とともに速やかに無限系のエネルギー密度に収束していることが確認できる。このようなサイズ依存性の変化は、1/L2 の寄与を生み出す端の効果が抑制されるだけでなく、

《講義ノート》

物性研究・電子版 Vol. 7, No. 2, 072205 （2018年11月号）

0.54

L = 50µ = 0 ~ −0.22

0.5

n(x)
0.46
SSD
OBC

0.42
10

20

30

40

x
図 11: SSD を適用したとき (実線) と通常の開放端のとき (破線) の電子密度分布端から最も離れた系の中央付近のエネルギーが周期境界条件と同様に最適化されていることを示し、端のエネルギースケールが 0 に漸近することによって、端の効果が系の中央付近の状態に影響し難くなっていることを反映している。実際にここで化学ポテンシャルを連続的に変化させると、図 11 に示すように、系の中央付近の電子密度は連続的に変化し、有限サイズ効
果の影響が現れない無限系の結果をほぼ正確に再現することが確認できる。このように、エネ
ルギーや電子密度が化学ポテンシャルの連続的変化に応じて無限系と同じように系統的に変化することから、特定の化学ポテンシャルにおいて並進対称な基底状態が正確に再現される自由電子のときの状況が、一般のモデルや化学ポテンシャルの場合においても近似的に成り立つことを期待することができる。
そこで最後に、SSD を用いることで帯磁率や電荷感受率といった外場応答が正確に得られることを、自由電子の局所状態密度を求めることで確認してみよう。外場応答は基底状態における量子相転移や臨界的特性、励起ギャップの有無を実験的に確認する際に利用する重要な物
理量であるが、有限系において磁場や電場に対する応答係数を求めようとすると、保存量子数の存在によって、誘起される磁化や電子密度の変化が離散的に生じ、通常はサイズスケーリン
グを行わないと外場の大きさと誘起される変化量との間の比例係数を正しく求めることはできない。この離散的変化の原因は系に与えた量子数と境界条件によって外場の大きさと無関係に量子数の密度が固定されてしまう境界条件の影響であるため、SSD によってこの境界条件の縛りを解いてしまえば、正しい応答係数が有限系で計算できると考えられる。

SSD を適用すると、系の端と中央にはエネルギーに対する重みの差が生じるため、中央付近の状態を外場に対して最もエネルギーが下がるように調整し、全電子数を整数値にするための余分な電子密度は端の部分に押し付けるという中央部分の最適化が自動的に行われる。この最適化によって、外場に対する中央領域の密度の連続的な変化が可能になるが、このときの変化の大きさは、中央付近の状態密度の大きさによって定まるはずである。そこで最後に、系の中央付近の局所状態密度を自由電子の場合で計算して、厳密な無限系での状態密度と比較《講義ノート》

物性研究・電子版 Vol. 7, No. 2, 072205 （2018年11月号）

0

Exact
SSD

-0.2

t =1
L = 50

-0.4

E
-0.6
-0.8
-1

0

0.1

0.2

0.3

0.4

ρ(E)

図 12: SSD を適用した 50 サイトの系 (破線) と無限系 (実線) の局所状態密度 ρ(E)してその信頼性を確認してみよう。
図 12 に SSD を適用して計算した系の中央での局所状態密度を示すが、得られた有限サイズの系の結果は、ほぼ正確に無限系の結果を再現している。通常の有限系で見られる離散的な効果が消失しているとともに、わずか数十サイトの系で定量的にも無限系の結果を再現しており、電場や磁場といった外場に対する応答係数が確かに正しく得られることが理解できる。

参考文献[1] S. R. White: Phys. Rev. Lett. 69 2863 (1992)[2] S. R. White: Phys. Rev. B 48 10345 (1993)[3] T. Nishino and K. Okunishi: J. Phys. Soc. Jpn. 65 891 (1996)[4] T. Nishino and K. Okunishi: J. Phys. Soc. Jpn. 66 3040 (1997)[5] R. J. Bursill, T. Xiang and G. A. Gehring: J. Phys.: Condns. Matter 8 L583 (1996)[6] X. Wang and T. Xiang: Phys. Rev B 56 5061 (1997)[7] N. Shibata: J. Phys. Soc. Jpn. 66 2221 (1997)[8] S. R. White and D. J. Scalapino: Phys. Rev. B 55 6504 (1997)[9] N. Shibata and D. Yoshioka: Phys. Rev. Lett. 86 5755 (2001)[10] G. Vidal: Phys. Rev. Lett. 93 040502 (2004)[11] S. Östlund and S. Rommer: Phys. Rev. Lett. 75 3537 (1995)

《講義ノート》

物性研究・電子版 Vol. 7, No. 2, 072205 （2018年11月号）[12] T. Nishino, Y. Hieida, K. Okunishi, N. Maeshima, Y. Akutsu, and A. Gendiar: Prog.
Theor. Phys. 105 409 (2001)[13] U. Schollwöck: Rev. Mod. Phys. 77 259 (2005)[14] H. F. Trotter: Proc. Amer. Math. Soc. 10 545 (1959)[15] M. Suzuki: Commun. Math. Phys. 51 183 (1976)[16] M. Suzuki: Phys. Rev. B 31 2957 (1985)[17] H. Betsuyaku: Phys. Rev. Lett. 53 629 (1984)[18] H. Betsuyaku: Prog. Theor. Phys. 73 319 (1985)[19] T. Koma: Prog. Theor. Phys. 81 783 (1989)[20] N. Shibata, B. Ammon, M. Troyer, M. Sigrist and K. Ueda: J. Phys. Soc. Jpn. 67 1086
(1998)[21] T. Mutou, N. Shibata and K. Ueda: Phys. Rev. Lett. 81 4939 (1998)[22] D. Yoshioka: Phys. Rev B 29 6833 (1984)[23] R. B. Laughlin: Phys. Rev. Lett. 50 1395 (1983)[24] A. A. Koulakov, M. M. Fogler, and B. I. Shklovskii: Phys. Rev. Lett. 76 499 (1996)[25] M. P. Lilly, K. B. Cooper, J. P. Eisenstein, L. N. Pfeiﬀer, and K. W. West: Phys. Rev.
Lett. 82 394 (1999)[26] A. Gendiar, R. Krcmar and T. Nishino: Prog. Theor. Phys. 122 953 (2009)[27] T. Hikihara and T. Nishino: Phys. Rev. B 83 060414 (2011)[28] H. Katsura: J. Phys. A 44 252001 (2011)[29] C. Hotta and N. Shibata: Phys. Rev. B 86 041108 (2012)