第一原理バンド計算からの強束縛模型化

第一原理バンド計算によって得られたバンド構造を、そのまま理論計算で用い ることは難しい。特に電子相関効果を議論する場合や、本質となるバンド構造を抽 出するような場合は、何らかの形で模型化を行わなくてはならない。本研究では 強束縛模型とハバード模型を扱うが、第一原理バンド計算と強束縛模型への橋渡 しの役目をする手段として最局在ワニエ関数がある。最局在ワニエ関数はMarzari とVanderbiltが提唱したもので[45, 46]、本論文では提唱した論文から説明する。

ワニエ関数はブロッホ関数をフーリエ変換したものである。ワニエ関数|Rn⟩は

|Rn⟩= V 8π³

∫

dke^−ik·R|ψnk⟩ (2.90) と表し、Rを中心とした波動関数が得られる。Vは単位格子の体積、|ψnk⟩はn番 目のバンドに対応するブロッホ関数である。

|ψnk⟩=e^ik·r|u_nk⟩ (2.91)

2.3. 第一原理バンド計算からの強束縛模型化 33

ワニエ関数による位置の平均と位置の2乗平均は、

⟨Rn|r|0m⟩=i V 8π³

∫

dke^ik·R⟨u_nk|∇k|u_mk⟩ (2.92) を用いて、

¯r_n =i V 8π³

∫

dk⟨u_nk|∇k|u_nk⟩ (2.93)

⟨r²⟩n= V 8π³

∫

dk||∇k|u_nk⟩||² (2.94) となる。

今までの話ではワニエ関数は一意に定まっていたが、一般的に、ワニエ関数は ユニタリ行列Uによる任意性がある。そのため、

|u_nk⟩ →∑

m

U_mn^(k)|u_mk⟩ (2.95)

と今までの式を置き換える。このユニタリ行列の任意性に対して、ワニエ関数が 最も局在しているという条件を課すことで「最局在ワニエ関数」を得る。そのた め、最局在ワニエ関数を表す指標は分散Ωをできるだけ小さくすることになる。

Ω =∑

n

[⟨r²⟩n− ¯r²_n] (2.96) 以下の議論は全て、分散を最小化するための手順を探すことになる。まず、Ωを2 つの項に分離する。

Ω = ΩI+Ω˜ (2.97)

この2項はそれぞれ

ΩI = ∑

n



⟨r²⟩n−∑

Rm

|⟨Rm|r|0n⟩|²



 (2.98)

Ω =˜ ∑

n

∑

Rm,0n

|⟨Rm|r|0n⟩|² (2.99)

これらの式から、ΩIは、任意のユニタリ変換に対して不変量になることがわかる。

∑

n⟨r²⟩nと、∑

n

∑

Rm|⟨Rm|r|0n⟩|² = ∑

n⟨0n|r∑

Rm|Rm⟩⟨Rm|r|0n⟩ = ∑

n⟨0n|A|0n⟩はト

レース和となっているためである。簡略化で∑

Rm· · ·をAと置いた。よって、分散 を最小にするためには、Ω˜ を最小にすることになる。Ω˜ をさらに2つに分離する。

Ω = Ω˜ OD+ ΩD (2.100) これらはそれぞれ

ΩOD = ∑

m,n

∑

R

|⟨Rm|r|0n⟩|² (2.101) ΩD = ∑

n

∑

R,0

|⟨Rn|r|0n⟩|² (2.102) となる。

今までの議論は波数空間が連続的な場合だったが、実際の数値計算では離散的 な数値計算に焼き直さなければならない。そのため、波数空間において、ある波 数kに対する最隣接点へのベクトルをbを導入する。そして、

∑

b

w_bb_αb_β = δαβ (2.103)

を満たすように重みw_bを定義する。単純立方格子ならば、wb =3/Z|b|²となり、Z は最隣接数である。微分は

∇f (k) = ∑

b

w_bb[ f (k+b)− f (k)] (2.104)

|∇f (k)|² = ∑

b

w_b[ f (k+b)− f (k)]² (2.105) となる。ここから、

⟨Rn|r|0m⟩=i1 N

∑

k,b

e^−ik·R⟨u_nk|w_bb(|u_nk+b⟩ − |u_nk⟩) (2.106)

これらの式を用いることで、位置の平均値と2乗平均は以下の式になる。

¯r_n = i N

∑

k,b

wbb[⟨unk|unk+b⟩ −1] (2.107)

⟨r²⟩n = 1 N

∑

k,b

w_b[2−2Re⟨u_nk|u_nk+b⟩] (2.108)

2.3. 第一原理バンド計算からの強束縛模型化 35

Ω˜ を求めるために、新しい行列Mを導入する。

M^(k,b)_mn = ⟨u_mk|u_nk+b⟩ (2.109)

bは微少量であるためM^(k_nn^,b)は1に近い値になっていることが予想される。b離れ た波動関数は位相がe^−ib·r異なる。そこで、bで展開すると以下の式となる。

M_nn^(k,b)=1+ixb+ 1

2yb²+O(b³) (2.110) ここで、x,yは適当な変数である。ここから、

M_nn^(k,b)−1 = ixb+O(b²) (2.111)

|M_nn^(k,b)|² = 1+x²b²+yb²+O(b³) (2.112)

となることから、

ixb = i Im lnM_nn^(k,b)+O(b²) (2.113)

−yb² = 1− |M_nn^(k,b)|²+x²b²+O(b³) (2.114) が成り立つ。これらを位置の平均値に代入すると、

¯r_n = i N

∑

k,b

wbb[M^(k_nn^,b)−1] (2.115)

= −1 N

∑

k,b

w_bbIm lnM_nn^(k,b) (2.116) となり、2乗平均は

⟨r²⟩n = 1 N

∑

k,b

w_b[2−2ReM_nn^(k,b)]

= 1 N

∑

k,b

wb[−yb²]

= 1 N

∑

k,b

w_b[

1− |M_nn^(k,b)|²+ x²b²]

= 1 N

∑

k,b

w_b[

1− |M_nn^(k,b)|²+[Im lnM_nn^(k,b)]²]

(2.117)

と求まる。また、M行列を用いてΩそれぞれの値は ΩI = 1

N

∑

k,b

w_b



J−∑

mn

|M_mn^(k,b)|²



 (2.118)

ΩOD = 1 N

∑

k,b

wb

∑

m,n

|M_mn^(k,b)|² (2.119) ΩD = 1

N

∑

k,b

w_b∑

n

(−Im InM_nn^(k,b)−b· ¯r_n)² (2.120) Jはバンドの本数を表す。最局在するための問題をMを最小化するということに 言い換えている。

実際の計算手順としては、まずブロッホ関数に対するユニタリ行列Uを単位行 列から微小変化させる。

U^(k)_mn =δmn+dW_mn^(k) (2.121) ここでdW は無限小のアンチエルミート行列で、dW^† = −dW を満たす。これに より、

|u_nk⟩ → |u_nk⟩+∑

m

dW^(k)mn|u_mk⟩ (2.122) と変換されることになる。ここで、dΩ/dW_mn^(k)について考える。簡略化として、任 意の関数Fに対して

(dF dW )

nm

= dF

dW_mn (2.123)

とすると、任意の行列Bに対して

dtr[dW B]

dW = B (2.124)

d Re tr[dW B]

dW = A[B] (2.125)

d Im tr[dW B]

dW = S [B] (2.126)

が成り立つ。ここで、A[B] =(B−B^†)/2、S [B] =(B+B^†)/2iである。この式から、

ΩI、ΩOD、ΩDについてのdWによる微小変化を考えると、ΩI,OD= ΩI + ΩODに関 しては

dΩI,OD= 4 N

∑

k,b

w_bRe tr[dW^(k)R^(k_mn^,b)] (2.127)

2.3. 第一原理バンド計算からの強束縛模型化 37

ここで、

R^(k_mn^,b)= M^(k_mn^,b)M_nn^(k,b)∗ (2.128) である。また、ΩDに関しては

q^(k,b)_n = Im lnM^(k,b)_nn +k·¯r_n (2.129)

R˜^(k_mn^,b) = M_mn^(k,b)

M_nn^(k,b) (2.130)

T_mn^(k,b) = R˜^(k,b)_mn q^(k,b)_n (2.131)

を用いて

dΩD= −4 N

∑

k,b

wbIm tr[dW^(k)T^(k,b)] (2.132)

となる。ここから、

G^(k)= dΩ

dW^(k) = 4∑

b

w_b(A[R^(k,b)]−S [T^(k,b)]) (2.133)

となる。最小化するための手順として、dW = ϵGとなるように定めれば、

dΩ =−ϵ∑

k,mn

|G^(k)_mn|² (2.134)

となり、常にΩが減り続けることになる。また、最小ではない鞍点で計算が収束 してしまうことを避けるために、通常は波動関数のプロジェクションを指定する。

プロジェクション関数|g_n⟩を初期のブロッホ関数に射影させることによって、実 際の物質の波動関数に近づけることができる。最終的に、強束縛模型のハミルト ニアンを得るには、最局在ワニエ関数を構築するときのユニタリ行列で、エネル ギーが対角のみにある行列H˜_mn(k)= ϵmkδmnを回転させることになる。

H(k)= U^(k)†H(k)U˜ ^(k) (2.135) これをフーリエ変換させて実空間に直したものがホッピングパラメータになる。

H(R)=∑

k

e^−ik·RH(k) (2.136)

これまでの計算は、最局在ワニエ関数を求めたいバンドが絶縁体のように他の バンドとギャップが開いている場合のみ使える。しかし、例えばfcc格子を持つ銅 のd軌道のみを最局在ワニエ関数を用いて抽出したい場合は、d軌道と同じエネル ギーにいるs軌道が邪魔をしてくる。このような場合でも最局在ワニエ関数を求 めるように拡張を行う。まず、抜き出したい軌道を持つバンドが必ず入っている ようなエネルギー範囲を選ぶ。これをアウターウィンドウと呼ぶ。このとき、ユ ニタリ変換されたブロッホ関数は

|ψnk⟩= ∑

m∈N_win^(k)

U_mn^dis(k)|ψmk⟩ (2.137)

となる。N_win^(k)は波数kのアウターウィンドウ内にあるバンドの数を表す。ユニタリ 行列が正方行列でないために、ΩIは不変量ではなくなってしまう。そのため、ま ずΩIを一意に決定付けることになる。ΩIを決定するために、自己無撞着の手法を 用いる。i番目の繰り返しによって得たΩIをΩ⁽ⁱ⁾_I とする。ラグランジュの未定乗 数法から、

δΩI

δu⁽ⁱ⁾_mk^∗ +∑

n

Λ⁽ⁱ⁾_nmk δ

δu⁽ⁱ⁾_mk^∗[⟨u⁽ⁱ⁾_mk|u⁽ⁱ⁾_nk⟩ −δmn]=0 (2.138) ここで、Λは定数の行列である。Ω⁽ⁱ⁾_I は

Ω⁽ⁱ⁾_I = 1 N_{k p}

Nk p

∑

k=1

ω⁽ⁱ⁾_I (k) (2.139)

ω⁽ⁱ⁾_I (k) = ∑

b

wb

∑

m



1−∑

n

|⟨u⁽ⁱ⁾_mk|u⁽ⁱ_k+⁻¹⁾_b⟩|²



 (2.140)

となり、Nk pは波数の数を表す。i番目の情報をi−1番目のk+bの情報から更新 する。この式を未定乗数法に代入すると、



∑

b

w_bP⁽ⁱ⁻¹⁾_k+b



|u⁽ⁱ⁾_mk⟩=∑

n

Λ⁽ⁱ⁾_nmk|u⁽ⁱ⁾_nk⟩ (2.141) ここでP⁽ⁱ⁾_k = ^N₂^{k p} ∑

n|u⁽ⁱ⁾_nk⟩⟨u⁽ⁱ⁾_nk|である。Λに対角化した値を用いれば固有値問題に帰 着できるため、



∑

b

w_bP⁽ⁱ⁻¹⁾_k+b



|u⁽ⁱ⁾_mk⟩=λ⁽ⁱ⁾_mk|u⁽ⁱ⁾_nk⟩ (2.142)

2.3. 第一原理バンド計算からの強束縛模型化 39

結局ω⁽ⁱ⁾_I は

ω⁽ⁱ⁾_I (k)= N

∑

b

wb−∑

m

λ⁽ⁱ⁾_mk (2.143)

より求まる。これを自己無撞着に解くことで、ΩIが決定される。ΩIが決定した 後は以前と同様にΩD,ΩODを求めることで、最局在ワニエ関数を求めることがで きる。

本博士論文では、pwscf、またはwien2kパッケージを用いて第一原理バンド計 算を行い、wannier90パッケージ[47]を用いて最局在ワニエ関数を得た。wien2k から最局在ワニエ関数を得る際にはwien2kパッケージから得られた波動関数を wannier90で読み込める形に変換するために、wien2wannierパッケージを使用した [48]。

2.3.1 強束縛模型と多軌道ハバード模型

複数軌道の電子系のハミルトニアンの運動エネルギー項H_Kは場の演算子を用 いて

H_K =∑

σ

∑

µν

∫

drψ^µ†_σ(r)K(r)ψ^ν_σ(r) (2.144) と表される。ここで、K(r)は一粒子の運動エネルギー演算子を表し、場の演算子 はそれぞれ

ψ^µ†_σ = ∑

i

ϕ^µ∗(r−R_i)c^µ†_iσ (2.145) ψ^µ_σ = ∑

i

ϕ^µ(r−R_i)c^µ_iσ (2.146) ϕ^µ(r−Ri)はiサイトを中心に持つµ軌道の波動関数を表し、以下ϕ^µ_i と表す。c^µ_iσ,c^µ†_iσ はサイトiにスピンσを持つ電子の生成消滅演算子である。生成消滅演算子はフー リエ変換して波数表示

c^µ†_iσ = ∑

k

c^µ†_kσe^−ik·r (2.147)

c^µ_iσ = ∑

k

c^µ_kσe^ik·r (2.148)

で表す。第二量子化することで、運動エネルギー項は H_K =∑

i j

∑

µν

∑

σ

t_{i j}^µνc^µ†_iσc^ν_jσ (2.149) となる。ここで、t^µν_{i j} はν軌道の jサイトからµ軌道のiサイトへの電子の飛び移り やすさを表し、ホッピングパラメータと呼ばれ、

t^µν_{i j} =⟨ϕ^µ_i|K|ϕ^ν_j⟩ (2.150) となる。運動エネルギー項のみのこのハミルトニアンを強束縛模型という。強束 縛模型におけるエネルギーバンドは強束縛模型のハミルトニアンを対角化するこ とで得られる。また、強束縛模型におけるバンドn番目の波数kの群速度vⁿ_kは

vⁿ_k = 1

ℏ∇kϵkⁿ (2.151)

で定義されるが、ヘルマン・ファインマンの定理から、

∇kϵkⁿ = ∇k⟨ψⁿk|H_K|ψⁿk⟩

= ∇k(⟨ψⁿk|ψⁿk⟩)+⟨ψⁿk|∇kH_K|ψⁿk⟩

= ⟨ψⁿk|∇kH_K|ψⁿk⟩ (2.152) より、微分した強束縛模型のハミルトニアンを波動関数で挟めば得ることができる。

強束縛模型に電子間相互作用の演算子V(r,r^′)を加えることを考える。これも場 の演算子より

HI = 1 2

∑

σσ^′

∑

µνµ^′ν^′

∫ ∫

drdr^′ψ^µ†_σ(r)ψ^ν†_σ′(r^′)V(r,r^′)ψ^ν_σ^′′(r^′)ψ^µ_σ^′(r) (2.153)

と表される。式(2.145)(2.146)を代入し、

HI = 1 2

∑

i ji^′j^′

∑

µνµ^′ν^′

V_{i ji}^µνµ′j^′′ν′

∑

σσ^′

c^µ†_iσc^ν†_jσ′c^ν_j^′′σ^′c^µ_i′^′σ (2.154)

となる。ここで、

V_{i ji}^µνµ′j^′′ν′ =

∫ ∫

drdr^′ϕ^µ∗(r−Ri)ϕ^ν∗(r^′−Rj)V(r,r^′)ϕ^ν′(r−Rj^′)ϕ^µ′(r−Ri^′) (2.155)

2.3. 第一原理バンド計算からの強束縛模型化 41

となる。本論文において、電子間相互作用の起源としてクーロン相互作用を考え る。多体電子系におけるクーロン相互作用は遮蔽効果を受けるため、二個の自由 電子間クーロン相互作用の距離減衰よりも急激な減衰をする。これは自由電子よ りも相対的に短距離の電子間相互作用を考慮することが現実の物質を模型化する 上で適切であることを示している。

多体電子の遮蔽効果がある場合は、同じサイト同士の電子間相互作用のみを考 えることができる(i= j= i^′ = j^′)。この場合H_Iは以下の4項に書きなおすことが できる。

H_I =∑

i

[∑

µ

U_µn_iµ↑n_ıµ↓+∑

µ>ν

∑

σσ^′

U_µν^′ n_iµσn_iνσ′

−∑

µ,ν

J_µνS_iµ·S_iν+∑

µ,ν

J_µν^′ c^µ†_i↑c^µ†_i↓c^ν_i↓c^ν_i↑ ]

(2.156) それぞれ、U は同軌道間で働くオンサイト相互作用、U^′は他軌道間で働くオン サイト相互作用、Jはフント則によるエネルギー、J^′は電子が相互に移動するこ とによるペアホッピングエネルギーを表す(図2.6)。Siµはスピン演算子で、Siµ =

1 2

∑

s1,s2c^µ†_is

1σs1s2c^µ_is

2で、σs1s2はパウリ行列である。運動エネルギー項と、同サイト内 相互作用項のみを残した相互作用項を合わせた模型を多軌道ハバード模型という。

H = H_K+H_I =∑

i j

∑

µν

∑

σ

t^µν_{i j}c^µ†_iσc^ν_jσ

+∑

i

[∑

µ

U_µniµ↑nıµ↓+∑

µ>ν

∑

σσ^′

U_µν^′ niµσniνσ^′

−∑

µ,ν

J_µνS_iµ·S_iν+∑

µ,ν

J_µν^′ c^µ†_i↑c^µ†_i↓c^ν_i↓c^ν_i↑ ]

(2.157) また、多軌道ハバード模型のハミルトニアンをスピンによって符号が変化する部 分をSˆ、符号が変化しない部分をCˆというように置き換えると、

HI = 1 4

∑

i

∑

µ1µ2µ3ν4

∑

σ1σ2σ3σ4

(−S_µ1µ2µ3µ4σ1σ2+C_µ1µ2µ3µ4)

× δσ1σ4δσ2σ3c^µ_iσ^1†

1c^µ_iσ^2†

2c^µ_iσ³

3c^µ_iσ⁴

4 (2.158)

を得る[49, 50]。このσは上向きスピンでは1、下向きスピンでは-1をとる。ここ

で、Sˆ とCˆはそれぞれ、

Sˆ =











U U^′

J J^′











,Cˆ =











U

−U^′+2J 2U^′−J

J^′











(µ1 =µ2 =µ3= µ4) (µ1 =µ3 ,µ2= µ4) (µ1 =µ4 ,µ2= µ3) (µ1 =µ2 ,µ3= µ4)

(2.159)

である。Sˆ とCˆはそれぞれ4階のテンソル量だが、例えば2軌道模型の場合を考 えると、行列

Sˆ =











S₁₁₁₁ S₁₁₂₂ S₁₁₁₂ S₁₁₂₁ S₂₂₁₁ S₂₂₂₂ S₂₂₁₂ S₂₂₂₁ S1211 S1222 S1212 S1221

S₂₁₁₁ S₂₁₂₂ S₂₁₁₂ S₂₁₂₁











(2.160)

のように定義することで、以後述べる感受率との積を行列の積として表現するこ とができる。本論文においては、強束縛模型、または多軌道ハバード模型を用い、

ホッピングパラメータに必要な波動関数に最局在ワニエ関数を用いる。

2.4 多体計算

非従来型超伝導や、3d電子系など、電子相関効果が重要な役割を果たす系では、

多軌道ハバード模型のような多体相互作用のハミルトニアンを考慮しなくてはな らない。一般的に1体のハミルトニアンをH₀、2体のハミルトニアンをH₁し、

H = H₀+H₁での摂動論を考え、第二量子化されたハミルトニアンを用いて、本 論文ではH₀ =H_K、H1 = H_Iとする。本節では簡単のため、1軌道モデルで考える が、後の揺らぎ交換近似において多軌道に拡張する。大分配関数は以下となる。

Z= Tr[ e^−βH]

= e^{−βΩ(V,T,µ)} (2.161) Ωは系全体の熱力学ポテンシャル、Vは体積、T は温度、µは化学ポテンシャルを 表す。分配関数のボルツマン因子におけるHを非摂動項と摂動項に分けると、

Z =e^−βΩ =Tr[

e^−βH0U(β)]

= e^−βΩ0⟨U(β)⟩0 (2.162)

ドキュメント内 電気通信大学大学院電気通信学研究科 (ページ 52-66)