密度行列繰り込み群法による二次元三角格子ハバード模型のスピン液体状態の研究Study of Spin Liquid States in Two-Dimensional Triangular Hubbard Model by Density-Matrix Renormalization Group

(1)

hp150112 「京」一般利用 K General Use

密度行列繰り込み群法による二次元三角格子ハバード模型の

スピン液体状態の研究

Study of Spin Liquid States in Two-Dimensional Triangular Hubbard Model

by Density-Matrix Renormalization Group

遠山貴巳 Takami Tohyama

東京理科大学 Tokyo University of Science

要旨半充填二次元三角格子ハバード模型の金属・絶縁体転移近傍では、三角格子特有のフラストレーション効果によりスピン液体状態が出現していると考えられている。一方向に周期境界条件、もう一方向に自由境界条件を課した系に対する密度行列繰り込み群法 (DMRG) のスピン相関関数の計算からはスピン液体を示唆する結果が得られていた。今回、運動量空間における磁気構造を明らかにするため、どちらの方向にも周期境界条件を課した系に DMRG を適用した。予想に反してスピン液体と思われる領域において磁気秩序が現れた。サイズ効果や打ち切り数依存性などの DMRG 計算の条件を詳細に検討した結果、系の形状や一次元化を施す手続きに結果が依存することがわかり、二次元三角格子ハバード模型特有の数値計算上の困難さが明らかとなった。キーワード：三角格子、ハバード模型、密度行列繰り込み群法、フラストレーション、スピン液体 Abstract

It is expected that there is a spin-liquid phase in the insulating region near metal-insulator transition in the two-dimensional triangular Hubbard model at half filling. It has been previously reported by density-matrix renormalization group (DMRG) calculations that there will be a spin-liquid behavior evidenced by spin-spin correlation for a system with periodic boundary condition in one direction and open boundary condition in the other direction. In order to clarify magnetic properties in the momentum space, we perform DMRG calculations for a system with periodic boundary condition in both directions. In contrast to an expected spin-liquid behavior, we find a magnetic order. After detailed investigations in terms of system size, truncation number, and so on, we find that the geometry of the systems as well as the choice of one-dimensional chain strongly influences on the results, indicating numerical difficulties in the frustrated two-dimensional Hubbard systems.

Keywords: Triangular lattice, Hubbard model, Density-matrix renormalization group, Frustration, Spin liquid

(2)

1. 研究の背景と目的 電子密度が半充填（ハーフ・フィーリング）の二次元三角格子を有する有機物質や無機物質が数多く発見されている。これらの物質は二次元三角格子上のハバード模型で記述できると考えられている。格子点上の電子間クーロン相互作用が小さいときは金属だが、クーロン相互作用を増加させると金属からモット絶縁体への転移を起こす。有機物質は、そのモット転移点よりもわずかに絶縁体側に入ったところに位置している。その領域での基底状態がどのように特徴付けられるのかが物性基礎研究の最重要課題のひとつとなっており、現在、理論・実験により様々なタイプの基底状態が提案されている。特に、三角格子特有のフラストレーション効果によりどのような磁気状態が出現しているのかが焦点であり、磁気長距離秩序の有無、スピン液体状態（磁気秩序は無いが電荷ギャップ有）でのスピンギャップの有無、最低励起はスピン一重項か三重項か等々、興味深い話題がつきない [1]。これらの課題を解決するには、大きな格子点数に対して正確な数値計算を行うことで基底状態、励起状態の情報を得ることが必要である。しかし、フラストレーションのある強相関電子系は、数値的厳密手法での取り扱いが最も難しい問題とされている。例えば、強相関電子系の手法として最も有力な計算手法に量子モンテカルロ法が挙げられるが、フラストレーションのある系では負符号問題が発生するため精度の保障は見込めない。これに対し、密度行列繰り込み群法(DMRG) は一次元系に対して有効であるとされてきたが、もし大規模計算が可能ならば二次元系に対しても適用可能である。二次元三角格子ハバード模型の金属・絶縁体転移近傍の磁気的性質に対しても DMRG による研究は何らかの情報を提供できる可能性がある。実際、我々は以前の課題で DMRG を用いて二次元三角格子ハバード模型の基底状態を調べ、一方向に周期境界条件、もう一方向に自由境界条件 (シリンダー境界条件) を課した 6x6=36 格子点の系に対するスピン相関関数の計算からスピン液体を示唆する結果を得ている [2]。本研究では、二次元三角格子ハバード模型のスピン液体状態の特徴を明らかにすることを目的とした。特に、以前のシリンダー境界条件とは異なり、どちらの方向も周期境界条件とした系を用いて、スピン液体状態での運動量空間での磁気構造を明らかにすることを中心課題として取り扱う。 2. 計算モデル 二次元三角格子上のハバード模型は以下のように与えられる。 † i i i i i i

H

t

c c

_ _

U

n n

   

 







ここで、 † i c はサイト

i

、スピン



の電子の生成演算子、t は最近接格子点間のホッピング、

U

は

(3)

や DMRG の計算 [2] から、

U

~ 7.5

t

より小さな

U

では金属、大きな

U

では絶縁体となることが指摘されている。さらに、_7.5_t_U _9.7_tの絶縁体領域では 120°磁気構造をもつ反強磁性秩序ではなく、スピン秩序のないスピン液体状態の兆候が得られている [2]。この計算結果は、

y

軸方向は周期境界条件だが

x

軸方向は自由境界条件であるシリンダー境界条件のもとでの6 6 格子で得られた。本研究では、運動量空間における磁気構造を明らかにして、スピン液体状態の特徴に迫ることを目的としたのでどちらの方向にも周期境界条件を施し DMRG 計算を実行した。また、サイズ効果も調べるため必要に応じて3 6 格子点や64格子点などの計算も行った。運動空間における磁気構造は、静的スピン構造因子

 

₀ z z ₀ S q  S Sq q を計算した。ここで、0 は基底状態、 z Sqは運動量qを持つスピン演算子のz成分である。 3. 並列計算の方法と効果 本研究では、京コンピュータの利用に向けて開発された 2 次元強相関系に適用可能な DMRG により計算を実施した。DMRG の計算は、要求する計算精度に対応して決定される打ち切り次数をmとして、 2 m のオーダーの次元をもつハミルトニアン行列に対して実行される。基底状態に ついては、ランチョス法を用いて計算した。ランチョス法では、ハミルトニアンとベクトルの積が現れる。この行列・ベクトル積が最も計算コストを要する部分であり、この大規模並列化が開発のポイントとなる。全系のハミルトニアンを二つのブロックの直積で表現できるようすると、そこに現れる行列・ベクトル積は行列・行列積に変換される。さらに行列・行列積ではハミルトニアンを行列の縦方向、横方向、元はベクトルであった行列の縦方向に対して 3 つの MPI のコミュニケータのグループを作成し、それぞれのコミュニケータ内の通信で大規模並列化を可能としている。これにより、行列・行列積で必要とする通信はそれぞれのグループ内で実行され、ノード数を_N として各グループに所属するノード数_N1 3_{の通信となるため、並列化効率の向上が図ら} れる。また、特に京コンピュータを利用する場合、その Tofu インターコネクトに対して、この 3 つのコミュニケータのグループを

x

軸、

y

軸、

z

軸方向に振り分けることで、非常に効率的な大規模並列計算を実現している。さらに高効率な大規模並列計算を実現するため、実空間並列法 [4] を採用した。この実空間並列化法は有限系アルゴリズムに対する並列化手法である。有限系アルゴリズムでは、完全な基底を持つサイトを系の中でスイープさせることで基底の最適化が行われる。特に、二次元ハバード模型の計算では、一筆書きで一次元形状にマッピングしたことによる仮想的な長距離相互作用のため、この有限系アルゴリズムによる最適化が重要である。以上の高度化を行った DMRG のプロダクトランでは、京コンピュータ 4800 ノードで平均実行効率 20%程度であった。

(4)

4. 研究成果 図 1 のような

x

方向、

y

方向とも周期境界条件を施した6 6 格子を採用する。ここで、格子点の番号は、DMRG を適用するため一筆書きした際の順番を示す。まず、金属・絶縁体転移と 120°スピン構造が生じる

U

の値を、エネルギーの

U

依存性から求めた。その結果、それぞれ

U

~ 5

t

、

U

~ 10

t

となった。前者の値は以前のシリンダー条件の結果 [2] (U ~ 7.55t) とは値が少し異なるが、m の増加により値が増加する ことが予想されるので以前の結果とは矛盾しない。次に、周期境界条件の特長を生かして運動量空間での静的スピン構造因子

S( )

q

の計算を実行した。図 2 は、いくつ かの m に対する q(4



3, 0),

( , 0)





q

での

S( )

q

の

U

依存性である。 (4



3, 0)  q は 120°磁気構造に対応する。予想通り、

U



10 t

では、120°磁気構造の強度が強くなる。一方、

U



10 t

では

_q



_{( ,0)}



の磁気構造を示す結果となった。

12 U



t

と

U



7 t

のときの第一ブリリュアン・ゾーンでの

S( )

q

を図 3 に示す。

12 U



t

では 120°磁気構造で期待されるようにq(4



3, 0)とその対称点で

_{S( )}

_q

は最大となる。一方、

U



7 t

では

q



( ,0)



だけで最大値を示している。そのため、方向性を持った磁気構造が実現している。このような結果はこれまで報告され ていなかったので、まず m 依存性を調 べたが、図 2 に示すように顕著な違いは得られなかった。また、サイズ効果の可能性を調べるため、3 6 格子点や 図 1 周期境界条件を課した6 6 32格子。番号は、DMRG を適用するため一筆書きした際の順番を示す。 図 2 q(4



3,0), ( ,0)



での_{S( )}_q の

U

依存性。いくつかの

m

での結果を示す。 (a) (b) 図 3 図 1 に示された6 6 格子に対するS( )q の第一ブリリュアン・ゾーンでの等高線図。(a) U12t、(b) 7 U  t。青色ほど値が小さく、赤色ほど値が大きい。第一ブリリュアン・ゾーンの右端はq(4 3,0)に対応し、そこでの値とその近傍の_q_{( ,0)} での値は図 2 に示されている。

(5)

の計算も行った。

_q



_{( ,0)}



での

S( )

q

が大きくなる場合もあったが、そのときは m の増大ととも に、その強度は弱くなっていった。したがって、6 6 格子での

_q



_{( ,0)}



の磁気構造の出現は、 6 6 格子特有の現象ではないかと予想された。実際、6 6 格子での一筆書きによる一次元鎖へのマッピングには任意性がある。その効果ではないかと疑い、

q



( ,0)



の磁気構造が一筆書きの違いによって生き残るか調べた。図 4 にはいくつかの場合のエネルギーと

_{S( )}

_q

のいくつかのmに対する結果を示す。ここで、“firstly vertical” は図1 の番号付けである。“firstly transvers” は図 1 の 1-12-･･･-25-35-36-26-･･･のように横に最初に一筆書きを始めた場合である。この2 つの場合は一筆書きに顕著な方向依存性があり、それが磁気構造に敏感に反映する可能性がある。“circular” は 1-2-3-4-5-6-7-18-19-30-31-32-33-･･･のように外周から内側に入り込んでいく一筆書きである。“random” はなるべく方向性を排除するように一筆書きしたものである (1-12-2-3-10-4-5-8-6-･･･)。エネルギーは “firstly vertical” が最低であるが、それだけではそのマッピングが最適とは言えない。なぜなら、上で述べたように

( ,0)





q

での大きな

S( )

q

の値は明らかに異常だからである。マッピングを変えると

q



( ,0)



での

S( )

q

は大きく変化する。方向性が小さいと期待される “random” では確かにその値は小さくなっており、q(4



3, 0)での値に近い。スピン液体状態での

S( )

q

には

q

依存性はないと期待されることと矛盾しない。以上の結果より、図 1 で採用した一筆書きでは

_q



_{( ,0)}



に対応するスピン配置を好むようなマッピングを採用していることが判明した。したがって、図 1 に示した境界条件での

q



( ,0)



の磁気構造の出現は、計算条件の設定不備による artifact の可能性が高いと考えている。以上のことを確認するため、割り当てられた京コンピュータでの計算資源を使い切ってしまい、本来の目的を達成できなかった。周期境界条件の下での運動量依存性の計算の困難さが明らかになっただけとは言え、周期境界条件の下での半充填二次元ハバード模型の DMRG 計算は京コンピュータなどでの大規模計算でのみ実行できるものであり、今後の計算を進める上での重要な情報となったと考えている。 図 4 周期境界条件を持つ6 6 格子に対する

U



7 t

での格子点当たりのエネルギー ( 青色菱印 ) とq(4



3, 0) ( , 0)  q での_{S( )}_q (それぞれ黒四角、赤丸)。いくつかの一筆書きの場合と、いくつかの

m

での結果を示す。

(6)

5. まとめと今後の課題 二次元三角格子ハバード模型のスピン液体状態の特徴を明らかにすることを目的として、スピン液体状態での運動量空間での磁気構造を明らかにしようと試みた。周期境界条件を施した6 6 格子の系を中心として、そのエネルギー、静的スピン構造因子を DMRG で用いて計算した。予想に反してスピン液体と思われる領域において磁気構造が現れた。サイズ依存性や精度を決定す る m 依存性、一筆書きによる一次元系へのマッピング依存性などの詳細な検討の末、その結果は 計算条件の設定の不備による artifact の可能性が高いことが判明した。以前用いたシリンダー境界条件では物理量の正確な運動量依存性を調べるのは難しい。特に、フラストレーションが強くこれまでもなかなか議論が収束してこなかった二次元三角格子ハバード模型の正確な基底状態の性質を議論するには運動量依存性は重要である。今後、DMRG のプログラムの改良を進めて、周期境界条件の下でのより高い精度での二次元三角格子ハバード模型の計算ができるように努力したい。 参考文献

[1] 例えばレビューとして P. A. Lee, Science 321, 1306 (2008); L. Balents, Nature 464, 199 (2010). [2] T. Shirakawa, T. Tohyama, J. Kokalj, S. Sota, and S. Yunoki, to be published.

[3] J. Kokalj and R. H. McKenzie, Phys. Rev. Lett. 110, 206402 (2013). [4] E. M. Stoudenmire and S. R. White, Phys. Rev. B 87, 155137 (2013).