hp140128 「京」以外 HPCI 一般利用 HPCI other than K General Use
動的密度行列繰り込み群法による三角格子ハバード模型の研究
Dynamical DMRG study of triangular lattice Hubbard model
曽田繁利1), 白川知功1,2), 遠山貴巳1,3), 柚木清司1)
Shigetoshi Sota1), Tomonori Shirakawa1,2), Takami Tohyama1,3), and Seiji Yunoki1) 1)理化学研究所, 2)SISSA, 3)東京理科大学
1)RIKEN, 2)SISSA, 3)Tokyo Univ. of Sci.
要旨 本研究では、三角格子ハバード模型のスピン励起状態の解析を動的密度行列繰り込み群法により 行った。三角格子ハバード模型のような二次元強相関電子系に対する動的密度行列繰り込み群法 の適用は非常に巨大な計算コストを要するため、必要な大規模並列プログラムの開発も同時に進 めた。本研究では、特に量子スピン液体状態におけるスピンギャップの有無を明らかにすること を目的とした。しかしながら、スピンギャップの有無を明らかにするためには、サイズ依存性等 の議論などのため、さらに巨大な計算資源を必要とすることが明らかになった。そのため、ここ では二次元強相関量子系に対する動的密度行列繰り込み群法の適用方法と本研究課題で得られ た結果について報告する。 キーワード：ハバード模型、スピン・フラストレート系、強相関系、三角格子模型、 密度行列繰り込み群法 Abstract
We have studied spin excitations on the half-filled triangular lattice Hubbard model by using dynamical DMRG method. The dynamical DMRG method on two-dimensional strongly correlated electron systems requires large computational resources. Thus we also developed massively parallel dynamical DMRG code. Our original main purpose was to examine if there is spin excitation gap in a possible quantum spin liquid phase. However, in order to obtain the conclusive numerical results, we had to systematically explore the size dependence. In this report, we show the algorithm of the dynamical DMRG for two-dimensional strongly correlated quantum systems and the results for the spin excitation spectra.
Keywords: Hubbard model, spin frustration, strongly correlated system, triangular lattice, DMRG © 2019 Research Organization for Information Science and Technology All rights reserved.
Received: 30 October 2018 Accepted: 6 June 2019 Available online: 10 June 2019
1. 研究の背景と目的 電子密度が半充填の三角格子を有する有機物質や無機物質が数多く発見されている。特に、三 角格子を有する有機物質については、そのモット転移点よりもわずかに絶縁体側に入ったところ に位置していると考えられる。その領域での基底状態がどのように特徴付けられるのかが、今日 の物性基礎研究の最重要課題のひとつとなっている。これまでの密度行列繰り込み群法[1]による 基底状態計算により、同一サイトの電子間クーロン斥力 U に対する相図を作成した[2]。本研究で は、この相図をもとに各相に対して動的密度行列繰り込み群法による計算を行い、各相のスピン 励起構造を明らかにすることを目的とする。 2. 計算モデル 本研究では、強相関電子系の代表的な量子格子模型である電子の飛び移り t と同一サイトの電 子間のクーロン斥力 U から構成されるハバード模型を三角格子に適用する。ハミルトニアンは、 次のように与えられる。
†
, , , , , ,h.c.
i j i i i j iH
t
c c
U
n n
ここで、c
i,(
c
i†,)
は i 番目のサイトのスピンの消滅(生成)演算子、ni, は i 番目のサイトのスピ ンの数演算子である。本研究では6×6=36 サイトの三角格子ハバード模型に対して計算を行っ た。図1 にその模式図を示す。この三角格子ハバード模型は、いわゆるスピン・フラストレーシ ョンが存在する系である。また、クーロン斥力 U が無限大の極限に対応する三角格子ハイゼンベ ルク模型の基底状態は隣り合うサイトのスピンが120°ずつ傾いた秩序状態(120°構造)を示すこ とが知られている。これまでの研究により、この120°構造と金属状態の間にスピン液体状態の 存在を明らかにしている[2]。そこで、本研究では密度行列繰り込み群法を動的な物理量の計算に 拡張した動的密度行列繰り込み群法により、動的スピン相関関数S(q,)の計算を行った。 図1 本研究における 6×6 サイト三角格子ハバー ド模型の模式図。青丸は各格子点を示し、数字は各 格子点の番号であり、密度行列繰り込み群法で一次 元的に導入される際の番号に対応する。赤線は密度 行列繰り込み群法のアルゴリズム上の一次元的な 経路を示し、青線は三角格子を形成するために導入 した密度行列繰り込み群法上の長距離相互作用に 対応する。また、本研究ではb 軸方向のみ周期的境 界条件を採用した。3. 並列計算の方法と効果（性能） 本研究では、「京」の利用に向けて開発された 2 次元強相関系に適用可能な密度行列繰り込み 群法(2D-DMRG)により、九州大学 FX10 上において、その計算を実行した。密度行列繰り込み群 法では、計算の目的に対応したターゲット状態の計算が最も計算コストを要する。本研究の動的 密度行列繰り込み群法では、基底状態と計算する動的な物理量に対応する励起状態がターゲット 状態となる。このターゲット状態は、本研究では基底状態に対してはランチョス法、励起状態に 対しては直交多項式展開法で計算される。ただし、ここで与えられるハミルトニアンは、密度行 列繰り込み群法によりターゲット状態に対して最適化された基底で表現される。また、密度行列 繰り込み群法で扱われるハミルトニアンの次元は、計算精度に対応して決定される密度行列繰り 込み群法の打ち切り次数 m に対して O(m2)の次元で与えられる。このターゲット状態の計算では、 ハミルトニアン(行列)とランチョス法、または直交多項式展開法による計算の各次数に対応した ベクトルの積が繰り返し現れる。この行列・ベクトル積が密度行列繰り込み群法を実行する上で 計算の中心となるため、この部分が大規模並列計算において如何に効率的に計算できるかが開発 のポイントとなる。そこで、本研究で用いられた2D-DMRG では、まずこの行列・ベクトル積に ついて、全系のハミルトニアンを密度行列繰り込み群法のアルゴリズムに対応した二つのブロッ クの直積で表現できるよう変換する。その結果、行列・ベクトル積は行列・行列積に変換される。 これにより、メモリバンド幅がボトルネックとなっていた本計算の単体性能の向上が図られた。 次に、この部分の大規模並列化として、この行列・行列積に対して、ハミルトニアンを行列の 縦方向と横方向、さらに元はベクトルであった行列の縦方向に対して 3 つの MPI のコミュニケ ータのグループを作成する。そして、それぞれのコミュニケータ内の通信だけで必要な並列計算 を実行するようにした。これにより、行列・行列積で必要とされる全対全通信は、ノード数を N として各グループに所属するノード数 N1/3の全対全通信となるため、利用するノード数の増加に 対して並列化効率の向上が図られる。例えば「京」を利用する場合、Tofu インターコネクトに対 して、この3 つのコミュニケータのグループを x 軸、y 軸、z 軸方向に振り分けることで、非常に 効率的な大規模並列計算が実現される。同じTofu インターコネクトを持つ九州大学 FX10 におい ても同様に非常に効率的な大規模並列計算が実現される。さらに高効率な大規模並列計算を実現 するため、2D-DMRG では実空間並列法[3]を採用した。この手法は、最適な基底を得るために実 行されるいわゆる有限系アルゴリズムに対する並列化手法であり、完全な基底を持つサイトを系 の中でスイープさせることで、ターゲット状態に対するハミルトニアンの表現を最適化する。特 に、本研究のような2 次元強相関系の密度行列繰り込み群法による計算では目的とする格子形状 に依存して導入される長距離相互作用に対する計算の収束を得るため、この有限系アルゴリズム によるプロセスが重要であり、計算時間のほとんどはこのプロセスで消費される。実空間並列法 [3]により、有限系アルゴリズムは、系の基底を部分ごとにほぼ完全な並列計算で最適化すること ができ、効率的な大規模並列計算が可能となる。
4. 研究成果 本研究では、三角格子ハバード模型のスピン励起状態解析に必要な動的物理量を、動的密度行 列繰り込み群法により計算する。特に、次のような動的スピン相関関数S(q,)を計算する。 0 1 1 ( , ) Im 0 ( ) ( ) 0 S S S E H i
q q q ここで、0
は基底状態、E0は基底状態のエネルギー、は微小な正の数、S(q)は運動量空間にお けるスピン相関関数であり、シリンダー状の境界条件を採用した本研究では、 , ,2
( )
sin(
)
(
)
(
1)
y y iq l x x l l l a bS
q l e
n
n
L
L
q
と与えられる。ここで、La、および Lbは図1 における a 軸、および b 軸方向のサイト数、lxおよ び lyは l 番目のサイトの x および y 座標を示す。動的密度行列繰り込み群法は、これまで一次元 強相関量子系の励起ダイナミクスの研究に用いられ、数多くの研究成果が報告されている。その 一方、動的密度行列繰り込み群法を本研究で対象とする三角格子ハバード模型のような二次元強 相関へ適用した例はほとんどない。その理由は、まず量子揺らぎを精密に取り扱うために必要な 密度行列繰り込み群法の打ち切り次数 m が非常に巨大になることが挙げられる。この事情は、エ ンタングルメント・エントロピーの面積則から理解される。この m について、密度行列繰り込み 群法の計算量はO(m3)と与えられるため、非常に大きな m による密度行列繰り込み群法の実行は 非常に大きな計算資源を必要とする。さらに、動的密度行列繰り込み群法では計算する動的な物 理量に対応した励起状態の計算も行う。この励起状態計算は、波動関数変換による高速化が期待 される基底状態と比較して100 倍以上の計算量を必要とする。したがって、効率的な大規模並列 動的密度行列繰り込み群法のプログラム開発が必要である。さらに、本研究を通じて確認された 注意点として、正しく収束した基底状態を出発点とした励起状態の計算が必要である点が挙げら れる。そのために、まず基底状態のみをタ ーゲット状態とした密度行列繰り込み群 法を実行する。そして、この基底状態計算 の収束が得られた後、計算する動的スピン 相関関数に対応する励起状態を含む 0 1 0 , S( ) 0 , S( ) 0 E H i
q q の3 状態をマルチターゲットとして計算を 実行する。本研究では、上述のような大規 模並列プログラム開発、および手法の確立 を進めつつ、目的となる動的スピン相関関 数の計算を行った。 -6 -4 -2 0 2 4 6 -6 -4 -2 0 2 4 6 qy qx M K 図2 6×6 サイト三角格子（シリンダー形状の境界条 件）において、運動量空間で独立な点を示したもの。 実線はブリルアンゾーン境界を示す。本研究では、これまでの研究で得られた電子間クーロン斥力 U に対する基底状態相図[2]をも とに、基底状態が金属、スピン液体、120°ネール秩序状態の 3 つに対応する U=6t, 8.5t, 12t につ いて、図1 に示した 6×6 サイト三角格子ハバード模型において動的スピン相関関数の計算を行 った。この6×6 サイト三角格子において運動量空間で独立な点を図 2 に示す。本研究課題では、 利用可能な計算資源(九州大学 FX10)の限界から図 2 の赤色の点で示した K 点、及び M 点近傍の 運動量に限定して計算を行った。 密度行列繰り込み群法の打ち切り次数 m としては、計算資源の制限から最大 m=4,000 として計 算を行った。図3 に本研究により得られた動的スピン相関関数 S(q,)を示す。U/t=12 の結果は、 スピン S=1/2 三角格子 Heisenberg 模型でよく記述できると考えられる Ba3CoSb2O9の中性子非弾 性散乱実験で観測されている[4,5]のスピン励起構造の特徴を再現している。特に/t=1 より高エ ネルギー側へ続く連続的なピーク構造は、これまでのスピン波理論では説明が困難であり、マグ ノン励起とは異なる量子揺らぎに起因する励起状態と考えられている。この連続状態の詳細を明 らかにすることは今後の重要な課題であり、現在、三角格子スピンS=1/2 反強磁性 Heisenberg 模 型において厳密対角化法、および動的密度行列繰り込み群法による解析を実行中である。U/t=8.5、 および U/t=6 の結果も同様の高エネルギー側へ続く連続状態が確認される。しかしながら、スピ ン液体状態におけるスピンギャップの有無を議論するためには、サイズ依存性の議論が必要であ る。本研究では利用できる計算資源の限界から 6×6 サイトの系に限定して計算を行わざるを得 なかった。今後の課題として、より大きな計算資源を用いてそのサイズ依存性を明らかにし、ス ピン液体状態におけるギャップの有無を議論したい。 0.0 0.5 1.0 1.5 2.0 2.5 3.0 0 5 10 15 20 25 30 35 40 S( q , ) (a rb . un it) / t U/t=12 U/t=8.5 U/t=6 (a) 0.0 0.5 1.0 1.5 2.0 2.5 3.0 0 20 40 60 80 100 120 140 (b) U/t=12 U/t=8.5 U/t=6 / t S( q , ) (a rb. un it) 図3 動的密度行列繰り込み群法により得られた動的スピン相関関数 S(q,ω) (a)および(b)はそれぞれ図 2 で示した M 点近傍および K 点近傍での計算結果を示す。
5. まとめと今後の課題 本研究では、動的密度行列繰り込み群法を三角格子ハバード模型に適用し、スピン液体状態と 考えられる領域に対するスピンギャップの有無を明らかにすることを目的に、動的スピン相関関 数の計算を行った。三角格子ハバード模型に対する動的密度行列繰り込み群法の適用はこれまで ほとんど前例がないこと、また非常に巨大な計算コストを必要とすると考えられるため、その手 法の確立と大規模並列計算を行うためのプログラム開発を進めつつ、その計算を行った。本研究 で得られた結果は、スピンS=1/2三角格子Heisenberg模型でよく記述できると考えられる Ba3CoSb2O9の中性子非弾性散乱実験で観測されている結果の特徴を再現するなど結果の妥当性 が認められる。しかしながら、本研究の目的であった量子スピン液体状態におけるギャップの有 無を明らかにするためには、そのサイズ依存性を明らかにするなど更なる計算と解析が必要であ り、これは今後の課題である。 参考文献
[1] S. R. White, Phys. Rev. Lett. 69, 2863 (1992); Phys. Rev. B 48, 10345 (1993).
[2] T. Shirakawa, T. Tohyama, J. Kokalj, S. Sota, and S. Yunoki, Phys. Rev. B 96, 205130 (2017). [3] E. M. Stoudenmire and S. R. White, Phys. Rev. B 87, 155137 (2013).
[4] J. Ma, Y. Kamiya, Tao Hong, H. B. Cao, G. Ehlers, W. Tian, C. D. Batista, Z. L. Dun, H. D. Zhou, and M. Matsuda, Phys. Rev. Lett. 116, 087201 (2016).
[5] S. Ito, N. Kurita, H. Tanaka, S. Ohira-Kawamura, K. Nakajima, S. Itoh, K. Kuwahara, and K. Kakurai, Nat. Communi. 8, 235 (2017).
