CeB6

まずは、SOC結合を考慮しないScalar relativistic計算を行ってみる。計算用ディレクトリ（私ならnonrel_k12x12x12とでも名付ける）を作成し、ディレクトリ内でfeditを実行して=.inファイルを用意する。XC汎関数はPBE-GGAで、k点は12*12*12程度で良いだろう。fploを実行する。結晶の対称性が良いので、すぐに計算が終わる。化学的な電子配置から、Ce 4f, Ce 5d, B 2pがフェルミ準位近傍で重要な役割を果たすと予想される。これらのPDOSファイルを探してプロットしてみよう。（ヒント：head -n 1 +dos.sort001.nl*のようなコマンドで、各ファイルの1行目を表示させると、どのファイルにどの軌道のPDOSが書かれているか一目でわかる。）また、バンドの軌道ウェイト計算も行い、各軌道の重みもプロットしてみよう。重み付きバンドは、+bweightsファイルに出力される。軌道(site, n, l, m)ごとに分解されている。ちょっと自由度が多すぎて理解しにくいので、mに対して和をとって(site, n, l)ごとの重みに再構成するのが良い。この重みの足し合わせには、faddweiが役立つ（詳細はFPLOマニュアルを見よ）。今回は、Ce 4f, Ce 5d, B 2pごとの重みを計算したいので、以下のような=.addweiファイルを用意する。これを直下に置き、faddweiを実行すると、重みが計算される。なお、faddweiはデフォルトで=.addweiを読み込む。別のファイル名で保存している場合は、実行時に引数として指定する。

weightinfile  +bweights

weightoutfile +bweights_Ce4f_Ce5d_B2p

name Ce4f

  atom Ce sites 1 orbitals 4f

name Ce5d

  atom Ce sites 1 orbitals 5d

name B2p

  atom B sites 2 orbitals 2p

生成した新bweightsファイルをxfbpやgnuplotなどでプロットしてみよう。以下には、重み付きバンドとPDOSをプロットした結果を示す。まず、DOSを見てみると、化学的な予想通り、フェルミ準位近傍~10 eVのエネルギー領域では、Ce 4f, 5d, B 2p軌道が主要な寄与であることが分かる。特に、フェルミ準位では、Ce 4f軌道が飛び抜けて大きな状態を持っている。4f電子は約0.5 eVの狭いバンドを持つ一方で、Ce 5dおよびB 2pは比較的ブロードなバンド分散を示している。この結果は、重い電子系の電子構造の特徴をよく捉えている。

続いて、SOCも考慮したfully relativistic計算を行う。新たなディレクトリ（私ならrel_k12x12x12などと名付ける）を作成し、nonrel_k12x12x12ディレクトリ内から=.inおよび=.densファイルをコピーする。feditを実行して=.in中の相対論設定をfully relativisticに変更する。また、再度SCF計算を行うので、Band plotsはOFFにしておくことをお薦めする。なお、電子密度はscalar relativistic計算で収束した値を用いることで、早く安全に収束させることができる。fploを実行してSCF計算。fully relativistic計算にしたことでスピノル自由度が増加し、scalar relativistic計算時よりも1ループが遅くなる。しかし、結晶の対称性が良いので、数分で収束するはずだ。収束したら、Band plotsをONにして再度fploを実行し、バンド分散を計算する。解析のため、weightsもONにしておく。バンド計算が終わったら、scalar relativistic (SR) 計算とfully relativistic (FR) 計算のバンド分散を同時プロットして比較してみる。XFBPを使うなら、 $ xfbp ../nonrel_k12x12x12/+band ./+band。Ce 4f軌道へのSOCの効果を見たいので、-1 eV < E - E_F < 1 eVほどの窓で表示してみよう。E = E_F + 0.5 eV周辺の準位がSOCによってエネルギー上昇しているように見える（下図）。

このようなSOCの影響を詳しく確認するために、バンド重みをプロットしよう。fully relativistic計算でbweightsをオンにすると、+bweightsおよび+bweightslmsが生成される。前者は全角運動量量子数jおよびその磁気量子数mjの基底への射影で、後者は（近似的に）軌道角運動量l、その磁気量子数m、スピンsに射影した場合の重みが出力されている。状況に合わせて解析する量を変えれば良い。ここでは、SOCによるCe 4f軌道の分裂を見たいので、j, mj基底の解析を行う。SOCによって、j = 5/2の6重項とj = 7/2の8重項への分裂が予測される。+bweightsファイルには(site, n, l, j, mj)分解されたウェイトが出力されている。解析を簡単にするため、再びfaddweiを用いる。入力ファイル=.addweiは以下の通り。

weightinfile  +bweights

weightoutfile +bweights_Ce4f

name Ce4f_j5/2

  atom Ce sites 1 orbitals 4f5/2

name Ce4f_j7/2

  atom Ce sites 1 orbitals 4f7/2

出力された+bweights_Ce4fファイルをXFBPやgnuplotでプロットすると、以下のような図が得られる。各準位が分散を持つため完全には分離してはいないが、j = 5/2と7/2の間に分裂が生じていることが見て取れる。このSOC分裂は、見積り方にも依存するが、数百meV程度である。

上述のfully relativistic GGAでは、フェルミ準位で4f電子が大きなDOSを示しており、4f電子が遍歴的である様子が結論付けられる。一方、磁性・輸送・dHvAなどの測定では、CeB6においてCe4f電子は非常に局在的であることが示唆されている。したがって、この4f電子の局在性を正しく説明するためには、プレーンのDFT（LDA/GGA）を超えた計算が必要である。そこで、ここではDFT+U法を適用して、平均場の範疇で4f軌道に局所Coulomb相互作用を考慮する。

一般に、球対称Coulombポテンシャルは、Slater-CondonパラメーターF^kを用いて指定される。f軌道（l = 3）では、k = 0, 2, 4, 6の4自由度がある。k = 0は単極子成分、k > 0は多極子成分である。F^0は軌道平均した直接相互作用に相当し、F^k (k > 0)は軌道平均した交換相互作用を規定する。FPLOでは、このSlater-CondonパラメーターF^kを直接的に指定することができる。DFT+U法では、ふつうパラメーターとしてF^kを扱う。注目したい物理に即してパラメーターを選ぶべきである。例えば、磁気揺らぎや超伝導などフェルミ面の形状が重要となる物理を議論する場合には、dHvA効果やARPESなどで観測されたフェルミ面を再現するようにF^kを選ぶ。Heisenberg模型の結合パラメーターJ_ijなどの高エネルギーの量を議論するには、Curie-Weiss温度などの同程度のエネルギースケールの量を再現するように選ぶ。いずれにせよ、Ce 4fのバンド幅が数eVのオーダーであるから、DFT+Uで追加考慮する4f電子相関も数eVが限界であろう。

ここでは、簡単のためF^k (k > 0)をゼロと置く。非常に単純化した近似であるが、GGA+Uにおける相関軌道の振る舞いを見ることはできる。まずはF^0 = 1 eVを導入してGGA+U計算を回してみよう。ここでは作業ディレクトリをrel_u01ev_k12x12x12と名付ける。プレーンGGAの電子密度を引き継げるので、=.densをディレクトリ内にコピーしておこう。feditのサブメニューから(L)SDA+Uを選択することで、DFT+Uの設定ができる。これが収束したら、ディクレクトリごとコピーして、F^0を2 eVに増加させて計算を繰り返す。これを手動でやるのは大変なので、以下のようなスクリプトを使い自動化するのが良いだろう。F^0 = 1 eVはすでに収束していて、F^0 = 12eVまで1 eV刻みで計算するスクリプトである。

#! /bin/sh

prev_u="01"

for u_val in $(seq 1 1 12); do

    u_dir=$(printf "%02d" "$u_val")

    prev_dir="rel_u${prev_u}ev_k12x12x12"

    new_dir="rel_u${u_dir}ev_k12x12x12"

    cp -r "${prev_dir}" "${new_dir}"

    cd "${new_dir}" || exit 1

    sed -e "s/_F0_/${u_val}/" ../=.in.tmp > ./=.in

    fplo > fplo.out 2>&1

    cd ../ || exit 1

    prev_u="${u_dir}"

done

収束した電子密度から軌道重み付きのバンド図を描写する。j = 5/2の6重項とj = 7/2の8重項を観察したいので、以前使った=.addweiファイルを使ってfaddweiしておこう。j = 5/2成分を赤、j = 7/2を青で色付けした結果を以下に示す。まず、U = 0（GGAに相当）では、j = 5/2状態が部分占有されている。これは、CeB6中のCeイオンがおおよそCe3+、すなわち4f^1の価電子配置を持つことに対応している。

GGA+Uでは、このような4f状態に対して、局所クーロン相互作用Uを平均場的に導入する。その結果、局所4f状態密度における占有的な成分と非占有的な成分の間にエネルギー分裂が生じる。直感的には、占有4f状態は低エネルギー側へ、非占有4f成分は高エネルギー側へ移動する。今回の計算では、j = 5/2成分の一部が占有的な4f状態に対応し、j = 5/2の残りやj = 7/2成分が主に非占有4f状態として振る舞っていると解釈できる。実際、Uを大きくすると、占有4f成分と非占有4f成分の分離が大きくなる。

ただし、この分裂幅がそのままUに等しくなるわけではない。たとえば、単純な1軌道Hubbard模型を平均場近似で扱うと、占有状態と非占有状態の分裂幅はおおよそUで与えられる。しかし、実際のGGA+U計算では、Ce 4f軌道はCe 5d軌道やB 2p軌道と混成しており、バンド固有状態は純粋な4f状態ではない。また、電子密度、フェルミ準位、混成の大きさも自己無撞着に変化する。そのため、バンド構造上で観測される4f成分のエネルギーシフトや分裂幅は、入力したUの値よりもかなり小さく見えることがある。

さらに、Uを8 eV程度まで大きくしても、フェルミ準位近傍に明瞭なMott絶縁体的バンドギャップが生成されない。この点は、一見すると意外に見えるかもしれない。しかしCeB6では、フェルミ準位近傍にCe 4f軌道だけでなく、Ce 5dやB 2p軌道に由来する遍歴的な伝導バンドも存在する。今回のGGA+Uでクーロン相互作用を導入しているのはCe 4f軌道であり、伝導電子に対応するバンドはフェルミ準位を横切り続ける。また、Ce 4f軌道と伝導バンドとの混成も残っているため、フェルミ準位近傍のバンドには4f重みが残り得る。

このように、GGA+UによってCe 4f状態の占有・非占有分離が強調される一方で、CeB6が単純な4f Mott絶縁体にはならない。CeB6の低エネルギー電子構造は、局在的な4f電子と遍歴的な伝導電子が混成した金属状態として理解する必要がある。

続いて、フェルミ面を計算・プロットしてみよう。比較のため、CeをLaで置き換えた場合の結果（fully relativistic GGA）も(e)に示す。LaはCe3+から4f電子を全て取り除いた元素に相当する。(a)から(d)にかけてUを増加させると、(e) LaB6のフェルミ面に近づいてゆくことがわかる。LaB6は、4f電子が完全に局在した極限ともみなすことができるので、これは直感通りの結果とも言える。では、実際に観測されたフェルミ面と比較してみよう。T_Q = 3.4 K以下では反強四極子秩序が、T_N = 2.3 K以下では反強磁性秩序が存在するので、10 K程度の常磁性データと比較するのが望ましい。たとえば、ARPESの論文 https://doi.org/10.1103/PhysRevB.92.104420 でFig. 2 (a) にkz = 0カットが示されている。興味深いことに、我々のDFT計算の結果の中で、最もARPESのフェルミ面を再現するのは、(e) LaB6の結果である！これは、主要なフェルミ面がCe 5dおよびB 2p由来の伝導電子によって形成され、Ce 4f電子はかなり局在的に振る舞うことを示している。したがって、GGA+UでCeB6のフェルミ面を再現するには、Ce 4f電子がフェルミ準位から遠く離れるほど大きなUを入れるか、LaB6の電子構造を用いるしかない。しかし、こうすると、CeB6の低エネルギー準粒子や質量増大を完全に記述することはできない。これは静的DFT+Uの限界であって、動的平均場近似（DMFT）の導入などが必要になる。