第18章分光スペクトルの第一原理計算

ここまでの17章で，我々は多電子問題を密度汎関数理論という形に整理し(第I〜III部)，固体のバンド構造として読み解き(第IV部)，実際に計算機の上で解く技術を整え(第V部)，原子核を動かし(第16章)，電流を流した(第17章)．しかし理論が正しいかどうかを決めるのは，最終的には実験である．そして電子状態を直接のぞき込む実験の代表が分光——物質に光を当て，出てくる電子や透過・反射する光を測る手法——である．

本章の主題は，計算した電子状態を「実験と直接比較できる量」に翻訳することである．翻訳すべき量は4つある．内殻電子を真空へ叩き出すX線光電子分光(XPS)，内殻電子を非占有状態へ持ち上げるX線吸収端近傍構造(XANES/NEXAFS)，価電子を真空へ叩き出す角度分解光電子分光(ARPES)，そして価電子を伝導帯へ持ち上げる紫外・可視・近赤外吸収(UV-Vis-NIR)である．それぞれについて，実験が測っているものは何か，計算では何を求めればよいか，そして基底状態の理論であるDFTでどこまで踏み込めるのかを，式のレベルで詰めていく．

本章はまた本書の最終章でもある．18.9節では，新しい二次元物質の構造を複数の分光量を突き合わせて確定していく実際のワークフローを見る．単一の観測量では構造は決まらず，すべての分光量を同時に説明できるモデルだけが生き残る——この方法論こそが，第一原理計算が実験と組む理由である．そして18.10節で，本書全体を振り返り，基底状態DFTの先に広がる地図を描いて筆を置く．

この章で学ぶこと

4種類の電子遷移(XPS・XANES・ARPES・UV-Vis-NIR)の統一的な整理
光電効果のエネルギー保存則と，仕事関数・化学ポテンシャルの正しい扱い
内殻束縛エネルギーの $\Delta$SCF 定式化:$E_B = E^{(0)}(N{-}1) - E^{(0)}(N) + \mu_0$ の完全導出
Janakの定理からの積分表示 $E_B = -\int_0^1\varepsilon_c(n)\dd n$ と Slater の遷移状態近似の誤差評価
コアホールを局在化させる技術(ペナルティ汎関数)と，荷電セルの静電発散を絶つ技術(厳密クーロンカットオフ)
化学シフトを「価電子電荷 + Madelungポテンシャル + 緩和」に分解する解釈の枠組み
Gunnarsson–Lundqvistの議論:各対称性クラスの最低状態にDFTが使える理由
双極子選択則 $\Delta l = \pm1$ の導出と，XANESが局所部分状態密度を映す仕組み
ARPESの運動学と，スーパーセルのバンドを解きほぐすアンフォールディング重み $W_{\bm{K}J}(\kk)$
独立粒子近似の誘電関数 $\varepsilon_2(\omega)$ の完全導出，Kubo–Greenwood形式，Kramers–Kronig関係
複数の分光量を突き合わせて構造を確定するワークフロー

本章でもHartree原子単位系($\hbar = m_e = e = 4\pi\varepsilon_0 = 1$，第1章)を基本に用いる．ただし分光の世界では慣習的に eV が使われるので，要所では eV 換算を併記する(1 Hartree $= 27.2114$ eV，1 Rydberg $= 13.6057$ eV)．長さは Bohr(1 Bohr $= 0.52918$ Å)を基本とし，実験値との比較では Å を用いる．電磁場を扱う18.8節では，$4\pi$ 因子の置き場所を明示するため Gauss単位系(cgs)の表式を併記する．光子エネルギーと波長の換算は，覚えておくと便利な関係

$$ \lambda\,[\text{Å}] \times E\,[\mathrm{eV}] = 12398 $$

を使う(可視光 $2$ eV は $6200$ Å，炭素K端の $285$ eV は $43.5$ Å，Al K$\alpha$ 線の $1486.7$ eV は $8.34$ Å)．この数値は18.6節と18.8節で双極子近似の妥当性を判定するときに効いてくる．

18.1 4種類の電子遷移と分光法

分光法の名前は無数にあるが，電子状態計算の立場から見ると，光が引き起こす一電子遷移は始状態がどこかと終状態がどこかの2つで分類できる．始状態は「内殻」か「価電子帯」か，終状態は「真空(試料の外)」か「物質内の非占有状態」かである．この $2\times2$ が本章の全体構成そのものになる．

図18.1 光が引き起こす4種類の電子遷移．始状態(内殻か価電子帯か)と終状態(真空か物質内の非占有状態か)の組合せで，主要な分光法が整理される．本章はこの図の4本の矢印を，左から順に第一原理計算の対象へ翻訳していく．

それぞれの矢印を，必要な光子エネルギー・測れる物理量・計算すべき量の3つで整理しておこう．

分光法	遷移	光子エネルギー	実験で測る量	計算で求める量
XPS(X線光電子分光)	内殻 → 真空	$10^2$–$10^4$ eV	内殻束縛エネルギー $E_B$，その化学シフト	コアホールを持つ終状態と基底状態の全エネルギー差
XANES / NEXAFS	内殻 → 非占有状態	$10^2$–$10^4$ eV	吸収端近傍の微細構造	コアホール存在下の非占有局所部分状態密度と双極子行列要素
UPS / ARPES	価電子帯 → 真空	$10$–$10^2$ eV	$E_B$ と面内波数 $k_\parallel$ の関係(バンド分散)	バンド構造 $\varepsilon_n(\kk)$，スーパーセルならアンフォールドしたスペクトル関数
UV-Vis-NIR 吸収	価電子帯 → 伝導帯	$0.1$–$6$ eV	吸収係数・反射率・屈折率	複素誘電関数 $\varepsilon(\omega)$(バンド間遷移行列要素の和)

この表の右端の列に注目してほしい．XPSだけが全エネルギーの差を要求しており，他の3つは一電子的な量(固有値と行列要素)で書ける．この違いは本質的である．XPSでは電子が1個減る——電子数が変わる——のに対し，他の3つでは系の電子数は変わらないか(XANES，UV-Vis)，あるいは変わっても価電子帯という「電子の海」から1個抜けるだけで系の場がほとんど変化しない(ARPES)からである．内殻から電子を抜くと，残った電子が空孔に向かって一斉に引き寄せられ，有効ポテンシャルが数十 eV も変化する．この巨大な緩和を扱うには，一電子固有値では足りず，全エネルギーの差を正面から計算しなければならない．18.2〜18.5節はこの一点をめぐる話である．

物理的意味:なぜ「4本の矢印」で整理するのか

実験装置のカタログには，XPS・ESCA・UPS・ARPES・NEXAFS・XANES・EXAFS・EELS・UV-Vis・楕円偏光解析…と多数の名前が並ぶ．しかし電子状態計算の立場から見れば，必要な情報は「どの準位から，どこへ，電子が動くのか」だけである．装置の名前ではなく遷移の始点と終点で分類すると，計算すべき量が一意に決まる．たとえば電子エネルギー損失分光(EELS)は，光の代わりに高速電子を使う点を除けばXANES/UV-Vis吸収と同じ遷移を見ており，計算上は同じ枠組みで扱える．逆に，同じ「光電子分光」という名前でも，内殻を見るXPSと価電子帯を見るARPESでは，計算の難しさの質がまったく異なる．

18.2 XPSの基礎:何が測れるのか

本節では，計算に入る前に実験の側を整理する．XPSが測っているのは何か，なぜ表面に敏感なのか，なぜ元素が識別できるのか，そしてなぜ同じ元素でも化学環境によってピークがずれるのか．これらを押さえておかないと，次節で導く式が何と比較されるべき量なのか分からなくなる．

18.2.1 光電効果とエネルギー保存則

物質に振動数 $\nu$ の光を当てると電子が飛び出す．これが光電効果であり，飛び出した電子(光電子)の運動エネルギーを測るのが光電子分光である．素朴には，光子のエネルギー $h\nu$ が電子の束縛を解くのに使われ，残りが運動エネルギーになる:

$$ \begin{equation} E_{\mathrm{kin}} = h\nu - E_B - \phi . \label{eq:18-ebind-exp} \end{equation} $$

ここで $E_B$ が束縛エネルギー(binding energy)，$\phi$ が仕事関数である．この式を $E_B$ について解いた

$$ \begin{equation} E_B = h\nu - E_{\mathrm{kin}} - \phi \label{eq:18-eb-def} \end{equation} $$

が，実験で $E_B$ を決める処方箋である．光子エネルギー $h\nu$ は既知(実験室系のAl K$\alpha$ 線なら $1486.7$ eV，Mg K$\alpha$ なら $1253.6$ eV，放射光なら可変)，$E_{\mathrm{kin}}$ は測定量，$\phi$ は装置の較正から決まる．

ところが，この $\phi$ をどう扱うかに落とし穴がある．「仕事関数」は，電子を系の中から真空中の静止状態まで持ち出すのに要する仕事であり，どこを基準に測るかで意味が変わる．金属では占有準位の上端がフェルミ準位で明確に決まるので，$\phi \equiv E_{\mathrm{vac}} - \mu$($\mu$ は化学ポテンシャル $=$ フェルミ準位)と定義できる．ところが絶縁体や分子では，化学ポテンシャルはバンドギャップの中のどこかにあり，不純物や欠陥や表面吸着で簡単に動いてしまう．しかも試料が絶縁体だと，光電子が抜けた分だけ試料が正に帯電し，以後の光電子が減速されてピークが見かけ上高束縛エネルギー側にずれる(チャージアップ)．

実験では，これを次の巧妙な仕掛けで回避する．試料と測定器(電子エネルギー分析器)を電気的に接続(オーミック接触)する．すると両者の間で電子が自由にやりとりでき，熱平衡では両者の化学ポテンシャルが一致する．以後，束縛エネルギーはこの共通の化学ポテンシャルを基準に測る．この状況を図18.2に描いた．

図18.2 XPS測定のエネルギーダイアグラム．試料と測定器はオーミック接触しており，化学ポテンシャル $\mu$ を共有する．光電子の全エネルギー(横の破線)はどこでも一定だが，真空準位が試料側と測定器側で異なる($\phi$ の差)ため，測定器が測る運動エネルギー $K$ は試料表面直上での運動エネルギーとは異なる．結果として，式\eqref{eq:18-eb-def}に現れるのは測定器の仕事関数だけになる．

図18.2から読み取れる重要な帰結を式にしておこう．試料の外に出た電子は，試料と測定器の間の接触電位差によって加速あるいは減速される．したがって「試料表面直上での運動エネルギー」と「測定器の中で測られる運動エネルギー $K_{\mathrm{spec}}$」は違う．ところが，光電子の全エネルギー(運動エネルギー $+$ その場のポテンシャルエネルギー)は保存する．測定器の中で電子が静止しているとしたときのポテンシャルエネルギーを $V_{\mathrm{spec}}$ と書くと，これは測定器の真空準位に等しく，共通の化学ポテンシャル $\mu$ から測って測定器の仕事関数 $\phi_{\mathrm{spec}}$ だけ上にある:

$$ \begin{equation} V_{\mathrm{spec}} = \mu + \phi_{\mathrm{spec}} . \label{eq:18-vspec} \end{equation} $$

ここが要点である．試料の仕事関数 $\phi_{\mathrm{samp}}$ はどこにも現れない．試料側の真空準位が高かろうが低かろうが，電子は接触電位差で加速・減速されて帳尻が合い，最終的に測定器の真空準位を基準にした運動エネルギーが測られる．だから\eqref{eq:18-eb-def}の $\phi$ は，正しくは $\phi_{\mathrm{spec}}$ である．装置の較正(標準試料のフェルミ端を測る)とは，この $\phi_{\mathrm{spec}}$ を決める作業にほかならない．

18.2.2 表面敏感性:非弾性平均自由行程

XPSの第二の特徴は，極端な表面敏感性である．X線自体は物質を $\mu$m のオーダーで透過する．しかし，深いところで生成した光電子は，表面に到達するまでにプラズモン励起や電子-正孔対生成といった非弾性散乱を受けてエネルギーを失い，本来のピーク位置には現れない(バックグラウンドになる)．したがって「元のエネルギーを保ったまま」出てこられるのは，表面からごく浅い領域で生成した光電子だけである．

この深さを支配するのが非弾性平均自由行程(inelastic mean free path, IMFP)$\lambda$ である．深さ $z$ で生成した光電子が非弾性散乱なしに表面に達する確率は $\exp(-z/\lambda\cos\theta)$($\theta$ は表面法線からの脱出角)であり，信号の $95\,\%$ は深さ $3\lambda\cos\theta$ までから来る．

$\lambda$ を運動エネルギーの関数としてプロットすると，驚くほど物質によらないユニバーサル曲線が現れる(図18.3)．$50$–$100$ eV 付近で最小値 $4$–$6$ Å をとり，両側で増大する U 字型である．この形の起源は次のように理解できる．エネルギー損失の主要なチャネルは価電子系の集団励起(プラズモン)であり，その断面積は損失関数 $\mathrm{Im}\bigl[-1/\varepsilon(q,\omega)\bigr]$ に比例する:

$$ \begin{equation} \frac{1}{\lambda} \;\propto\; \int \frac{\dd q}{q}\int \dd\omega\; \mathrm{Im}\!\left[\frac{-1}{\varepsilon(q,\omega)}\right] , \label{eq:18-imfp} \end{equation} $$

ここで $\varepsilon(q,\omega)$ は物質の誘電関数である(この量は18.8節で正面から計算する)．高エネルギー側では，電子が速すぎて相互作用する時間が短くなるため断面積が落ち，$\lambda$ は増大する($\lambda \propto E^{1/2}$ 程度)．低エネルギー側では，電子のエネルギーがプラズモン励起エネルギー($10$–$20$ eV)を下回るとその損失チャネルが閉じるため，やはり $\lambda$ が増大する．両者の間に最小値が現れる．誘電関数の大まかな形はどの固体でも似ている(価電子密度がそれほど違わない)ため，曲線が「ユニバーサル」になるのである．

非弾性平均自由行程のユニバーサル曲線（両対数）．光電子の運動エネルギー 2〜10000 eV に対して，50〜100 eV で最小値約5 Å をとるU字型の曲線で，Al Kα 線（1486.7 eV）の位置に赤い破線と赤丸を示す． — 図18.3 非弾性平均自由行程のユニバーサル曲線(両対数)．$50$–$100$ eV で最小値 $4$–$6$ Å をとり，これは原子層にして2〜3枚分にすぎない．実験室のAl K$\alpha$ 線($1486.7$ eV)では $\lambda \approx 15$–$25$ Å で，これでも情報深さは数 nm にとどまる．XPSが「表面の分析手法」と呼ばれる理由である．

物理的意味:表面敏感性は長所でもあり短所でもある

$\lambda \approx 5$–$25$ Å という値は，XPSが表面数原子層の情報しか持たないことを意味する．これは表面科学にとっては最大の武器であり，吸着分子・酸化被膜・界面の化学状態をピンポイントで調べられる．しかし計算と比較するときには注意が要る．バルクの単位胞で計算した内殻束縛エネルギーと，表面第一層の原子の値は数百 meV から $1$ eV 程度ずれる(表面原子は配位数が少なく，電荷分布も緩和も異なる)．したがって「実験値」と比較するときは，スラブ模型で表面近傍のサイトを別々に計算し，$\lambda$ による深さの重みを付けて足し合わせるのが正しい手続きである．18.9節の事例研究では，まさにこの区別が構造決定の決め手になる．

18.2.3 元素選択性と化学シフト

内殻準位のエネルギーは，原子核の電荷 $Z$ でほぼ決まる．$1s$ 準位なら，水素様原子の粗い見積もりで $E_B \approx Z_{\mathrm{eff}}^2/2$ Hartree であり，$Z$ とともに急速に深くなる．実際，$1s$ 束縛エネルギーは C: $285$ eV，N: $400$ eV，O: $533$ eV，F: $686$ eV と，元素ごとにはっきり分離している．したがってスペクトルにピークが立つ位置を見れば，そこにどんな元素があるかが直ちに分かる．これが元素選択性である(ただし水素とヘリウムは内殻を持たないので，この方法では検出できない)．さらにピーク面積は，その元素の原子数に(光イオン化断面積で重みを付けて)比例するので，定量分析もできる．

ここからが本章の主題である．同じ元素でも，化学環境が違うと束縛エネルギーが数 eV ずれる．これを化学シフトという．古典的な例として，トリフルオロ酢酸エチル $\mathrm{CF_3\text{-}CO\text{-}O\text{-}CH_2\text{-}CH_3}$ を挙げよう．この分子には炭素原子が4個あるが，化学環境がすべて異なる．実測のC $1s$ スペクトルは4本のピークに分かれ，束縛エネルギーの低い順に

炭素サイト	結合相手	C $1s$ 束縛エネルギー(概値)	相対シフト
$\mathrm{CH_3}$	C, H	$291.2$ eV	$0$(基準)
$\mathrm{CH_2}$	O, C, H	$293.4$ eV	$+2.2$ eV
$\mathrm{C{=}O}$	O, O, C	$295.7$ eV	$+4.5$ eV
$\mathrm{CF_3}$	F, F, F, C	$298.9$ eV	$+7.7$ eV

となる．順序は明快である．電気陰性度の高い原子(F, O)に囲まれた炭素ほど価電子を奪われて正に帯電し，その分だけ内殻電子が原子核に強く束縛される．$\mathrm{CF_3}$ の炭素はフッ素3個に電子を引かれてもっとも正に帯電しているので，束縛エネルギーがもっとも大きい．同様に，シリコンの $2p$ 準位は，単体シリコン($\mathrm{Si^0}$)で $99.4$ eV，二酸化ケイ素($\mathrm{Si^{4+}}$)で $103.4$ eV と $4$ eV もずれ，この差を使って酸化膜の厚さや界面の化学状態が定量される．

この定性的な描像を，あとで定量的な式に翻訳する(18.5節)．いま押さえておくべきなのは，化学シフトには2つの起源があるということである．

始状態効果(initial-state effect):光を当てる前の基底状態において，内殻軌道のエネルギーそのものが化学環境でずれている．原子上の価電子電荷が減れば内殻電子の遮蔽が弱まり，内殻準位は深くなる．また周囲のイオンが作る静電ポテンシャル(Madelungポテンシャル)も内殻準位を一様にずらす．
終状態効果(final-state effect):内殻に空孔ができた瞬間，周囲の価電子がその正電荷に引き寄せられて空孔を遮蔽する．この緩和によって終状態のエネルギーが下がり，束縛エネルギーは緩和エネルギーの分だけ小さくなる．遮蔽の効きやすさは環境に依存する(金属中では価電子が動きやすく緩和が大きい，絶縁体では小さい)ので，これも化学シフトに寄与する．

両者を分離することがXPSの解釈の要であり，そのためには終状態を実際に計算しなければならない．次節でその定式化を行う．

補足:内殻ピークが分裂するその他の要因

化学シフト以外にも，内殻ピークが複数本に分かれる要因が3つある．実験スペクトルを解釈するときは，これらを混同しないことが重要である．

スピン軌道相互作用．内殻では原子核近傍の強い電場のためスピン軌道相互作用が大きく，$l \ge 1$ の準位は全角運動量 $j = l \pm 1/2$ の2本に分裂する．強度比は縮重度 $2j+1$ の比で決まり，$j = l-1/2$ と $j = l+1/2$ について $$ I_{j=l-1/2} : I_{j=l+1/2} = 2l : 2(l+1) $$ すなわち $p$ 準位($l=1$)では $2p_{1/2} : 2p_{3/2} = 1:2$，$d$ 準位($l=2$)では $3d_{3/2}:3d_{5/2} = 2:3$ である．Si $2p$ の分裂幅は $0.6$ eV 程度で，高分解能測定では明確に分離できる．
磁気的交換分裂．内殻に空孔ができると，残った内殻電子は不対スピンを持つ．磁性体では価電子の正味スピンとの交換相互作用によってエネルギーが分裂する．$\mathrm{O_2}$ 分子($S=1$)の O $1s$ や，$\mathrm{MnF_2}$, $\mathrm{MnO}$ の Mn $3s$ が典型例である．閉殻の $\mathrm{N_2}$ では分裂しない．
化学ポテンシャルのシフト．18.2.1節で述べたとおり，半導体・絶縁体の $\mu$ はギャップ中にあり，ドーピング・欠陥・表面吸着・バンド曲がりで動く．$\mu$ を基準に測る束縛エネルギーは，その分だけ丸ごとシフトする．$p$ 形と $n$ 形の同じ半導体で内殻ピークがギャップ幅に近いだけずれることがあり，これは電子構造の変化ではなく基準の移動である．

18.3 内殻束縛エネルギーの理論

本節では，実験で測る $E_B$ を第一原理計算で求める量に翻訳する．結論を先に述べると，$E_B$ はコアホールを持つ終状態と基底状態の全エネルギー差である．この一見当たり前の主張を，図18.2のエネルギーダイアグラムから丁寧に導き，そのうえで金属と絶縁体で基準の取り方がどう違うかを明らかにする．最後に，第10章のJanakの定理を使って別ルートからも同じ結果に到達し，Slaterの遷移状態近似の誤差の次数を評価する．

18.3.1 エネルギー保存則から $\Delta$SCF へ

まず記号を定める．$N$ 電子系の基底状態(始状態)のエネルギーを $E_i(N)$，内殻軌道 $c$ に空孔を持つ $N-1$ 電子系(終状態)のエネルギーを $E_f(N-1)$ と書く．図18.2の全エネルギー保存則は次のようになる．

導出:測定量と全エネルギー差の関係

ステップ1:エネルギー保存則を書く．光を当てる前の系は「$N$ 電子の試料 $+$ 光子」であり，そのエネルギーは $E_i(N) + h\nu$ である．光電子放出後は「$N-1$ 電子の試料(内殻に空孔) $+$ 測定器の中を運動する電子1個」となり，そのエネルギーは $E_f(N-1) + V_{\mathrm{spec}} + K_{\mathrm{spec}}$ である．$V_{\mathrm{spec}}$ は測定器の中で静止した電子のポテンシャルエネルギー，$K_{\mathrm{spec}}$ は測定される運動エネルギーである．したがって

$$ \begin{equation} E_i(N) + h\nu = E_f(N-1) + V_{\mathrm{spec}} + K_{\mathrm{spec}} . \label{eq:18-econs} \end{equation} $$

ここで2つの近似をしている．第一に，光電子が飛び出すときに原子核が受ける反跳エネルギーを無視した．運動量保存から反跳エネルギーは $E_{\mathrm{rec}} = K_{\mathrm{spec}}\,m_e/M_{\mathrm{atom}}$ 程度であり，$1$ keV の光電子と炭素原子($M \approx 12\times1836\,m_e$)なら $0.045$ eV，目標精度 $0.1$–$0.5$ eV に対して無視できる(ただし軽元素・高エネルギーでは効いてくる)．第二に，原子核の位置は光電子放出の前後で変わらないとした(Franck–Condon的な固定核近似)．光電子が原子を離れる時間スケールは $10^{-16}$ 秒程度，格子振動の周期は $10^{-13}$ 秒程度なので，この分離は妥当である．

ステップ2:測定器の真空準位を化学ポテンシャルで表す．式\eqref{eq:18-vspec}を\eqref{eq:18-econs}に代入する:

\begin{align} E_i(N) + h\nu &= E_f(N-1) + \mu + \phi_{\mathrm{spec}} + K_{\mathrm{spec}} \label{eq:18-econs2}\\[2pt] \therefore\quad h\nu - K_{\mathrm{spec}} - \phi_{\mathrm{spec}} &= E_f(N-1) - E_i(N) + \mu . \label{eq:18-econs3} \end{align}

\eqref{eq:18-econs3}は\eqref{eq:18-econs2}を移項しただけである．左辺は実験の測定量そのもの，すなわち式\eqref{eq:18-eb-def}で定義した束縛エネルギー $E_B$ である．よって

$$ \begin{equation} E_B = E_f(N-1) - E_i(N) + \mu . \label{eq:18-eb-mu} \end{equation} $$

ステップ3:基準の任意性を吸収する．ここで問題になるのは，$E_f$, $E_i$, $\mu$ のすべてが「電位の基準をどこに取るか」に依存することである．試料と測定器を接触させると界面付近で電荷移動が起き，試料全体の静電ポテンシャルが一定量 $\Delta\mu$ だけシフトする．この基準のずれを陽に書き出そう．ある標準的な基準(たとえば計算コードが採用する電位のゼロ点)で測ったエネルギーに上付き $(0)$ を付けると，電子1個あたり $\Delta\mu$ だけポテンシャルがずれるので，電子数に比例して

\begin{align} E_i(N) &= E_i^{(0)}(N) + N\,\Delta\mu , \label{eq:18-shift-i}\\ E_f(N-1) &= E_f^{(0)}(N-1) + (N-1)\,\Delta\mu , \label{eq:18-shift-f}\\ \mu &= \mu_0 + \Delta\mu \label{eq:18-shift-mu} \end{align}

となる(化学ポテンシャルは1電子あたりの量なので $\Delta\mu$ が1回だけ加わる)．これらを\eqref{eq:18-eb-mu}に代入する:

\begin{align} E_B &= \Bigl[E_f^{(0)}(N-1) + (N-1)\Delta\mu\Bigr] - \Bigl[E_i^{(0)}(N) + N\Delta\mu\Bigr] + \bigl[\mu_0 + \Delta\mu\bigr] \label{eq:18-sub1}\\[2pt] &= E_f^{(0)}(N-1) - E_i^{(0)}(N) + \mu_0 + \bigl[(N-1) - N + 1\bigr]\Delta\mu \label{eq:18-sub2}\\[2pt] &= E_f^{(0)}(N-1) - E_i^{(0)}(N) + \mu_0 . \label{eq:18-sub3} \end{align}

\eqref{eq:18-sub2}では $\Delta\mu$ の係数を集め，\eqref{eq:18-sub3}ではその係数が $(N-1)-N+1 = 0$ で厳密に消えることを使った．∎

いま得られた結果が，内殻束縛エネルギー計算の出発点である．

$$ \begin{equation} \boxed{\;E_B = E_f^{(0)}(N-1) - E_i^{(0)}(N) + \mu_0\;} \label{eq:18-eb-general} \end{equation} $$

この式の意味を噛みしめておこう．第一に，$\Delta\mu$ が厳密に消えた．これは，試料の帯電やバンド曲がりによる一様なポテンシャルシフトが束縛エネルギーに影響しないことを意味する．逆にいえば，計算する側は電位のゼロ点を好きに選んでよい——ただし $E_f^{(0)}$, $E_i^{(0)}$, $\mu_0$ の3つを同じ基準で計算しなければならない．第二に，右辺は3つとも第一原理計算で求まる量である:基底状態の全エネルギー $E_i^{(0)}(N)$，コアホールを持つ $N-1$ 電子系の全エネルギー $E_f^{(0)}(N-1)$，そして基底状態の化学ポテンシャル $\mu_0$(金属ならフェルミ準位，半導体・絶縁体ならギャップ中の適切な位置)．第三に，この式は金属・半導体・絶縁体・分子のすべてに共通である．

定義:$\Delta$SCF法

始状態と終状態のそれぞれについて自己無撞着(SCF)計算を独立に行い，その全エネルギーの差からエネルギー差の観測量を求める方法を$\Delta$SCF法という(第10章でイオン化エネルギーの評価法として登場した)．式\eqref{eq:18-eb-general}はまさに $\Delta$SCF の処方であり，内殻束縛エネルギーへの適用が本章の主題である．$\Delta$SCF の利点は，終状態における電子の再配置(緩和)が変分的に完全に取り込まれる点にある．欠点は，終状態を安定に収束させる技術が必要な点であり，それが18.4節のテーマになる．

18.3.2 金属の場合:$\mu_0$ が消える

式\eqref{eq:18-eb-general}には $\mu_0$ が残っている．金属の場合，この項をJanakの定理で消去して，より扱いやすい形にできる．以下でそれを示す．

導出:金属における $E_B = E_f^{(0)}(N) - E_i^{(0)}(N)$

ステップ1:何を比べているのかを整理する．\eqref{eq:18-eb-general}に現れる $E_f^{(0)}(N-1)$ は，「コアホールがあり，かつ電子が1個少ない」系のエネルギーである．一方，計算上扱いやすいのは電荷中性な系，すなわち「コアホールがあり，抜けた電子はフェルミ準位に置かれている」系 $E_f^{(0)}(N)$ である(電子数が $N$ のままなのでセルが中性で，周期境界条件の下でも静電的な発散が起きない)．両者の差を評価しよう．

ステップ2:Janakの定理を使う．コアホールを持つ系(これも一つの well-defined な自己無撞着系である)において，フェルミ準位近傍の状態の占有数 $n$ を $1$ から $0$ へ連続的に変える経路を考える．第10章のJanakの定理により，占有数に関する全エネルギーの微係数はその軌道のKS固有値に等しい:

$$ \begin{equation} \frac{\partial E_f^{(0)}}{\partial n} = \varepsilon_{\mathrm{F}}(n) . \label{eq:18-janak-metal} \end{equation} $$

金属では，この $\varepsilon_{\mathrm{F}}(n)$ は化学ポテンシャル $\mu_0$ そのものである．しかも，マクロな金属($N \sim 10^{23}$)から電子を1個抜いてもフェルミ準位はまったく動かない($\delta\mu_0 \sim 1/N \to 0$)．すなわち $\varepsilon_{\mathrm{F}}(n) = \mu_0$ は $n$ によらない定数である．

ステップ3:積分する．微積分学の基本定理により

\begin{align} E_f^{(0)}(N-1) - E_f^{(0)}(N) &= -\int_0^1 \dd n\,\frac{\partial E_f^{(0)}}{\partial n} \label{eq:18-mint1}\\[2pt] &= -\int_0^1 \dd n\,\mu_0 \label{eq:18-mint2}\\[2pt] &= -\mu_0 . \label{eq:18-mint3} \end{align}

\eqref{eq:18-mint1}で符号が負なのは，$n$ を $1$ から $0$ へ下げる(電子を抜く)向きに積分しているからである．\eqref{eq:18-mint2}でJanakの定理と $\varepsilon_{\mathrm{F}} = \mu_0 = \text{const}$ を使い，\eqref{eq:18-mint3}で定数の積分を実行した．

ステップ4:一般式に代入する．\eqref{eq:18-mint3}より $E_f^{(0)}(N-1) = E_f^{(0)}(N) - \mu_0$ だから，これを\eqref{eq:18-eb-general}に入れると

\begin{align} E_B &= \bigl[E_f^{(0)}(N) - \mu_0\bigr] - E_i^{(0)}(N) + \mu_0 \label{eq:18-mfin1}\\[2pt] &= E_f^{(0)}(N) - E_i^{(0)}(N) . \label{eq:18-mfin2} \end{align}

$\mu_0$ が相殺した．∎

$$ \begin{equation} \text{金属:}\qquad E_B = E_f^{(0)}(N) - E_i^{(0)}(N) . \label{eq:18-eb-metal} \end{equation} $$

この式は計算上たいへん都合がよい．どちらの計算も電荷中性のセルで行えるので，周期境界条件下での静電エネルギーの発散(18.4.1節)を心配しなくてよい．しかも，コアホールが価電子によって遮蔽される様子——正電荷が1個できて，フェルミ準位付近の電子がそこへ流れ込む——が計算にそのまま入っている．金属におけるコアホールの「価電子遮蔽描像」は，この式によって理論的に正当化される．

物理的意味:なぜ絶縁体では金属の式が使えないのか

ステップ2で使ったのは「フェルミ準位に電子を置ける」「その準位のエネルギーは電子数によらない」という2つの性質であり，いずれも金属に固有である．有限のギャップを持つ系では，抜いた電子を置くべき「フェルミ準位の状態」が存在しない．無理に伝導帯の底に置けば，それは電子-正孔対を余分に励起した状態であり，$\mu_0$ の位置(ギャップ中のどこか)とは異なるエネルギーを支払うことになる．実際，有限ギャップ系に\eqref{eq:18-eb-metal}を適用すると，固体アンモニアの N $1s$ で $1$ eV 以上の過大評価が残ることが知られている．一方，金属に対して一般式\eqref{eq:18-eb-general}を使うのは正しいが，荷電セルの扱いが必要になるためセルサイズ収束が遅くなる．系に応じて式を使い分けるのが実務である．両者が同じ極限値に収束することは，たとえば一次元炭素鎖でコアホール間距離を $225$ Å まで伸ばした計算で，差 $0.07$ eV まで確認されている．

18.3.3 Janakの定理による別ルート:遷移状態近似

いま導いた $\Delta$SCF は2回のSCF計算を要する．実は，1回の計算で(近似的に)同じ答えを出す方法がある．これがSlaterの遷移状態近似である．導出はJanakの定理をもう一度，今度は内殻軌道の占有数について使う．第10章と重なる部分もあるが，内殻という「緩和が極端に大きい」場合にこそこの議論の価値が現れるので，ここで完結した形で再構成する．

いま考えるのは孤立系(分子)，あるいは真空準位を基準に取れる系である．内殻軌道 $c$ の占有数を $n$ とし，$n$ を $1$(基底状態)から $0$(コアホール状態)まで連続的に変えながら，そのつど自己無撞着に解く．分数占有のKS形式は第10章10.7.2節で導入した．

導出:$E_B$ の Janak 積分表示と Slater の遷移状態

ステップ1:Janakの定理を内殻軌道に適用する．第10章の定理により，他の占有数を固定して $n_c = n$ だけ動かすと

$$ \begin{equation} \frac{\partial E}{\partial n} = \varepsilon_c(n) \label{eq:18-janak-core} \end{equation} $$

である．ここで $\varepsilon_c(n)$ は「内殻軌道の占有数を $n$ に固定して自己無撞着に解いたときの，その軌道のKS固有値」であり，$n$ に強く依存する量である．

ステップ2:積分して有限のエネルギー差にする．束縛エネルギーは $E_B = E(n{=}0) - E(n{=}1)$ だから，微積分学の基本定理より

\begin{align} E_B &= E(n{=}0) - E(n{=}1) \label{eq:18-eint1}\\[2pt] &= -\bigl[E(n{=}1) - E(n{=}0)\bigr] \label{eq:18-eint2}\\[2pt] &= -\int_0^1 \dd n\,\frac{\partial E}{\partial n} \label{eq:18-eint3}\\[2pt] &= -\int_0^1 \dd n\,\varepsilon_c(n) . \label{eq:18-eb-integral} \end{align}

\eqref{eq:18-eint3}は $E(1)-E(0) = \int_0^1 (\partial E/\partial n)\dd n$ を使い，\eqref{eq:18-eb-integral}で\eqref{eq:18-janak-core}を代入した．束縛エネルギーは，内殻固有値を占有数について $0$ から $1$ まで平均した値の符号を変えたものである．

ステップ3:被積分関数を中点まわりに展開する．厳密な積分を実行するには $\varepsilon_c(n)$ を全 $n$ で計算しなければならず，これでは1回のSCFで済まない．そこで $\varepsilon_c(n)$ を $n=1/2$ のまわりでTaylor展開する:

$$ \begin{equation} \varepsilon_c(n) = \varepsilon_c\!\left(\tfrac12\right) + \left(n-\tfrac12\right)\varepsilon_c' + \frac{1}{2}\left(n-\tfrac12\right)^2\varepsilon_c'' + \frac{1}{6}\left(n-\tfrac12\right)^3\varepsilon_c''' + \frac{1}{24}\left(n-\tfrac12\right)^4\varepsilon_c''''+\cdots \label{eq:18-taylor} \end{equation} $$

ここで $\varepsilon_c', \varepsilon_c'', \dots$ はすべて $n=1/2$ での微係数である．この展開が有効なのは，$\varepsilon_c(n)$ が $[0,1]$ 上で滑らかだからである(実際，後述するように $\varepsilon_c(n)$ はほぼ直線になる)．

ステップ4:項別に積分する．$u = n - 1/2$ と置換すると $\dd n = \dd u$，積分区間は $u \in [-1/2, 1/2]$ となる．奇数次の項は対称区間で消えるので，

\begin{align} \int_{-1/2}^{1/2}\dd u &= 1, \label{eq:18-mom0}\\ \int_{-1/2}^{1/2} u\,\dd u &= \left[\frac{u^2}{2}\right]_{-1/2}^{1/2} = \frac{1}{8}-\frac{1}{8} = 0, \label{eq:18-mom1}\\ \int_{-1/2}^{1/2} u^2\,\dd u &= \left[\frac{u^3}{3}\right]_{-1/2}^{1/2} = \frac{1}{24}+\frac{1}{24} = \frac{1}{12}, \label{eq:18-mom2}\\ \int_{-1/2}^{1/2} u^3\,\dd u &= 0, \label{eq:18-mom3}\\ \int_{-1/2}^{1/2} u^4\,\dd u &= \left[\frac{u^5}{5}\right]_{-1/2}^{1/2} = \frac{2}{5}\cdot\frac{1}{32} = \frac{1}{80} . \label{eq:18-mom4} \end{align}

\eqref{eq:18-mom1}と\eqref{eq:18-mom3}がゼロになるのは被積分関数が奇関数だからである．これらを\eqref{eq:18-taylor}の項別積分に代入すると

\begin{align} \int_0^1 \varepsilon_c(n)\,\dd n &= \varepsilon_c\!\left(\tfrac12\right)\cdot 1 + \varepsilon_c'\cdot 0 + \frac{\varepsilon_c''}{2}\cdot\frac{1}{12} + \frac{\varepsilon_c'''}{6}\cdot 0 + \frac{\varepsilon_c''''}{24}\cdot\frac{1}{80}+\cdots \label{eq:18-mid1}\\[2pt] &= \varepsilon_c\!\left(\tfrac12\right) + \frac{1}{24}\varepsilon_c'' + \frac{1}{1920}\varepsilon_c''''+\cdots . \label{eq:18-mid} \end{align}

これは数値積分でいう中点則とその誤差項にほかならない．1階微分の寄与が完全に消えていることが要点である．

ステップ5:束縛エネルギーの表式を得る．\eqref{eq:18-eb-integral}に\eqref{eq:18-mid}を代入して

$$ \begin{equation} E_B = -\varepsilon_c\!\left(\tfrac12\right) - \frac{1}{24}\varepsilon_c'' - \frac{1}{1920}\varepsilon_c'''' - \cdots \;\simeq\; -\varepsilon_c\!\left(\tfrac12\right) . \label{eq:18-slater} \end{equation} $$

すなわち，内殻軌道の占有数を $1/2$ にして1回だけ自己無撞着に解き，その固有値の符号を変えれば束縛エネルギーが得られる．この半占有状態を遷移状態と呼ぶ．∎

$$ \begin{equation} E_B = -\int_0^1 \varepsilon_c(n)\,\dd n \;\simeq\; -\varepsilon_c\!\left(n=\tfrac12\right) \qquad(\text{Slaterの遷移状態近似}) \label{eq:18-slater-key} \end{equation} $$

他の評価法と誤差を比べておこう．いずれも\eqref{eq:18-taylor}を代入して積分値\eqref{eq:18-mid}と比較するだけである．

導出:Koopmans近似・台形則・中点則の誤差の比較

(a) Koopmans型($-\varepsilon_c(1)$)．\eqref{eq:18-taylor}に $n=1$($u = 1/2$)を代入すると

$$ \begin{equation} \varepsilon_c(1) = \varepsilon_c\!\left(\tfrac12\right) + \frac{1}{2}\varepsilon_c' + \frac{1}{8}\varepsilon_c'' + \frac{1}{48}\varepsilon_c''' + \cdots \label{eq:18-koop1} \end{equation} $$

である．これと厳密な積分値\eqref{eq:18-mid}との差を取る:

\begin{align} \varepsilon_c(1) - \int_0^1\varepsilon_c(n)\dd n &= \frac{1}{2}\varepsilon_c' + \left(\frac{1}{8}-\frac{1}{24}\right)\varepsilon_c'' + \frac{1}{48}\varepsilon_c''' + \cdots \label{eq:18-koop2}\\[2pt] &= \frac{1}{2}\varepsilon_c' + \frac{1}{12}\varepsilon_c'' + \frac{1}{48}\varepsilon_c''' + \cdots \label{eq:18-koop-err} \end{align}

\eqref{eq:18-koop2}では $1/8 - 1/24 = 3/24 - 1/24 = 2/24 = 1/12$ を計算した．1階微分に比例する誤差が残ることが致命的である．内殻の場合，$\varepsilon_c'$ は電子を抜くことによる有効ポテンシャルの深化を表す量で，その大きさは数十 eV に達する．だから $-\varepsilon_c(1)$(基底状態のKS固有値をそのまま読む)は $E_B$ を大幅に過小評価する．炭素の $1s$ で言えば，KS固有値は $-270$ eV 程度なのに実測の $E_B$ は $285$ eV である．

(b) 台形則($-\frac{1}{2}[\varepsilon_c(0)+\varepsilon_c(1)]$)．\eqref{eq:18-taylor}に $n=0$($u=-1/2$)を代入すると

$$ \begin{equation} \varepsilon_c(0) = \varepsilon_c\!\left(\tfrac12\right) - \frac{1}{2}\varepsilon_c' + \frac{1}{8}\varepsilon_c'' - \frac{1}{48}\varepsilon_c''' + \cdots \label{eq:18-trap1} \end{equation} $$

\eqref{eq:18-koop1}と\eqref{eq:18-trap1}を足して2で割ると，奇数次の項が相殺して

\begin{align} \frac{\varepsilon_c(0)+\varepsilon_c(1)}{2} &= \varepsilon_c\!\left(\tfrac12\right) + \frac{1}{8}\varepsilon_c'' + \frac{1}{384}\varepsilon_c'''' + \cdots \label{eq:18-trap2}\\[2pt] \frac{\varepsilon_c(0)+\varepsilon_c(1)}{2} - \int_0^1\varepsilon_c\,\dd n &= \left(\frac{1}{8}-\frac{1}{24}\right)\varepsilon_c'' + \cdots = \frac{1}{12}\varepsilon_c'' + \cdots \label{eq:18-trap-err} \end{align}

(4次の係数は $\frac{1}{24}\cdot\frac{1}{2}\bigl[(\frac12)^4+(-\frac12)^4\bigr] = \frac{1}{24}\cdot\frac{1}{16} = \frac{1}{384}$ である．)1階微分の誤差は消えたが，2階微分の誤差は中点則の $2$ 倍で，しかも符号が逆である．

(c) 中点則(遷移状態)．\eqref{eq:18-mid}より誤差は $\frac{1}{24}\varepsilon_c''$ のみ．三者のなかで最良である．∎

図18.4 Janakの積分表示\eqref{eq:18-eb-integral}の図形的意味．青い曲線が $\varepsilon_c(n)$，その下の面積が $\int_0^1\varepsilon_c\dd n = -E_B$ である．緑の破線は中点則の長方形で，曲線がほぼ直線であるためこの長方形の面積は曲線下の面積をきわめてよく近似する．基底状態の固有値 $\varepsilon_c(1)$(上の赤破線)をそのまま使うKoopmans型評価は，緩和エネルギー $R$ の分だけ束縛エネルギーを過小評価する．

18.3.4 始状態効果と終状態効果の分離

18.2.3節で「化学シフトには始状態効果と終状態効果がある」と述べた．式\eqref{eq:18-eb-integral}を使うと，この分離を厳密に定義できる．

やることは単純で，被積分関数に $\varepsilon_c(1)$ を足して引くだけである:

\begin{align} E_B &= -\int_0^1 \varepsilon_c(n)\,\dd n \label{eq:18-sep1}\\[2pt] &= -\int_0^1 \Bigl[\varepsilon_c(1) + \bigl(\varepsilon_c(n) - \varepsilon_c(1)\bigr)\Bigr]\dd n \label{eq:18-sep2}\\[2pt] &= -\varepsilon_c(1)\int_0^1\dd n - \int_0^1\bigl[\varepsilon_c(n)-\varepsilon_c(1)\bigr]\dd n \label{eq:18-sep3}\\[2pt] &= -\varepsilon_c(1) - \int_0^1\bigl[\varepsilon_c(n)-\varepsilon_c(1)\bigr]\dd n . \label{eq:18-sep4} \end{align}

\eqref{eq:18-sep2}は恒等的な足し引き，\eqref{eq:18-sep3}は積分の線形性($\varepsilon_c(1)$ は $n$ によらない定数なので外に出せる)，\eqref{eq:18-sep4}は $\int_0^1\dd n = 1$ である．そこで緩和エネルギーを

$$ \begin{equation} R \equiv \int_0^1 \bigl[\varepsilon_c(1) - \varepsilon_c(n)\bigr]\dd n \;\ge\; 0 \label{eq:18-relax} \end{equation} $$

と定義すると($n$ が減ると遮蔽が弱まり有効ポテンシャルが深くなるので $\varepsilon_c(n) \le \varepsilon_c(1)$，したがって $R \ge 0$)，

$$ \begin{equation} E_B = \underbrace{-\varepsilon_c(1)}_{\text{始状態効果}} \;+\; \underbrace{R}_{\text{終状態効果}} . \label{eq:18-eb-init-final} \end{equation} $$

という分解が得られる．$-\varepsilon_c(1)$ は基底状態の情報だけで決まる始状態項，$R$ はコアホールができたあとの電子の再配置だけから来る終状態項である．

物理的意味:緩和はなぜ「符号が正」なのか

式\eqref{eq:18-eb-init-final}を素朴に読むと「緩和すると束縛エネルギーが増える」と見えるが，これは符号の取り方の問題である．$R$ の定義\eqref{eq:18-relax}で $\varepsilon_c(n)$ が負の方向に深くなることが $R \gt 0$ に対応する．物理的には，終状態(コアホールのある系)が緩和によって安定化し，そのぶん $E_f$ が下がる．$E_B = E_f - E_i$ だから，$E_f$ が下がれば $E_B$ は小さくなる．そして実際，$-\varepsilon_c(1)$ を「緩和を許さない場合の束縛エネルギー」と読むと，緩和を許した本当の $E_B$ はそれより小さいはずである——ところが\eqref{eq:18-eb-init-final}では $+R$ で大きくなっている．矛盾しているだろうか．していない．鍵は，$-\varepsilon_c(1)$ は「凍結軌道近似での束縛エネルギー」ではないことである．KS固有値 $\varepsilon_c$ は Hartree–Fock の軌道エネルギーとは異なり，自己相互作用の扱いのために「その軌道の電子1個分のエネルギー」を表さない(第10章)．Janakの定理が教えるのは，$\varepsilon_c$ が微分量だということだけである．したがって\eqref{eq:18-eb-init-final}は「Koopmans型の値 $+$ 補正」という数学的な分解と読むのが正確で，その補正 $R$ が化学環境に依存する部分が終状態寄与の化学シフトである．実測の内殻シフトの $70$–$90\,\%$ は始状態項で説明でき，残りが終状態項というのが多くの系での経験である．

化学シフトはこの分解の差として書ける．同一元素の2つのサイト $A$, $B$ について

$$ \begin{equation} \Delta E_B \equiv E_B^{A} - E_B^{B} = \bigl[-\varepsilon_c^{A}(1) + \varepsilon_c^{B}(1)\bigr] + \bigl[R^{A} - R^{B}\bigr] . \label{eq:18-chemshift-split} \end{equation} $$

第1項が始状態シフト，第2項が終状態シフトである．18.5節では第1項をさらに「原子上の価電子電荷」と「周囲からのMadelungポテンシャル」に分解し，化学シフトの系統性を理解する枠組みを作る．その前に，終状態を実際に計算する技術を整えなければならない．

18.4 コアホール計算の技術

式\eqref{eq:18-eb-general}と\eqref{eq:18-eb-metal}は，原理的には「2回SCFを回して引き算するだけ」に見える．しかし実際には，終状態を計算しようとすると2つの技術的な壁にぶつかる．第一に，特定の原子の特定の内殻軌道に空孔を作り，それを維持したまま自己無撞着に解くにはどうするか．第二に，周期境界条件のもとで電荷を持つセルをどう扱うか．本節ではこの2つを順に片付ける．

18.4.1 周期境界条件と荷電セルの静電発散

まず問題の所在をはっきりさせる．周期系のHartreeエネルギーは，逆格子ベクトル $\bm{G}$ の和で書ける．密度のFourier成分を $\tilde\rho(\bm{G}) = \int_\Omega \rho(\rr)e^{-i\bm{G}\cdot\rr}\dd^3r$($\Omega$ はセル体積)と定義すると，

$$ \begin{equation} E_H = \frac{1}{2}\int\!\!\int \frac{\rho(\rr)\rho(\rr')}{\abs{\rr-\rr'}}\dd^3r\,\dd^3r' = \frac{1}{2\Omega}\sum_{\bm{G}}\frac{4\pi}{G^2}\abs{\tilde\rho(\bm{G})}^2 \label{eq:18-hartree-G} \end{equation} $$

である($4\pi/G^2$ はCoulomb核 $1/r$ のFourier変換，第12章・第14章)．ここで $\bm{G}=\bm{0}$ の項を取り出すと

$$ \begin{equation} E_H^{(\bm{G}=\bm{0})} = \frac{1}{2\Omega}\lim_{G\to0}\frac{4\pi}{G^2}\abs{\tilde\rho(\bm{0})}^2 , \qquad \tilde\rho(\bm{0}) = \int_\Omega\rho(\rr)\dd^3r = Q \label{eq:18-G0-term} \end{equation} $$

となり，セルの正味電荷 $Q$ がゼロでない限り発散する．物理的には当たり前で，同じ符号の電荷を無限に並べればCoulombエネルギーは無限大になる．コアホールを持つ終状態は $Q = +1$ なので，まさにこの状況にある．

慣用的な処方は補償背景電荷(jellium背景)の導入である．一様な電荷密度 $-Q/\Omega$ をセル全体に敷き，全体を中性にする．これは $\tilde\rho(\bm{0}) \to 0$ とすること，すなわち\eqref{eq:18-hartree-G}の $\bm{G}=\bm{0}$ 項を単純に落とすことに等しい．計算は走るようになる．しかし，こうして得られたエネルギーは「電荷 $Q$ の点欠陥が周期的に並び，一様背景で中和された系」のエネルギーであって，我々が欲しい「孤立した1個のコアホール」のエネルギーではない．

数学ノート:補償背景を使ったときの有限サイズ誤差(Makov–Payne展開)

電荷分布 $\Delta\rho$ が原点まわりの狭い領域に局在し，正味電荷 $Q$ と2次モーメント $q_r = \int r^2\Delta\rho(\rr)\dd^3r$ を持つとする．これを一辺 $L$ の立方格子に周期的に並べ，一様背景で中和した系の全エネルギー $E(L)$ は，孤立系の値 $E_\infty$ から次のようにずれる:

$$ \begin{equation} E(L) = E_\infty - \frac{\alpha_M Q^2}{2L} - \frac{2\pi Q\,q_r}{3L^3} + O(L^{-5}) . \label{eq:18-makov-payne} \end{equation} $$

第1項は，電荷 $Q$ の点電荷格子と補償背景からなる系(Wigner格子)のMadelungエネルギーである．$\alpha_M$ はMadelung定数で，単純立方格子では $\alpha_M = 2.8373$，面心立方では $2.8883$，体心立方では $2.8883$ 程度である．第2項は，電荷分布の広がり(単極子でない部分)と背景の相互作用から来る．

この補正の大きさを見積もってみよう．$Q = 1$，単純立方セルで $L = 20$ Bohr($\simeq 10.6$ Å)とすると

$$ \frac{\alpha_M Q^2}{2L} = \frac{2.8373}{40}\ \mathrm{Hartree} = 0.0709\ \mathrm{Hartree} = 1.93\ \mathrm{eV} . $$

目標精度は $0.1$–$0.5$ eV である．しかもこの誤差は $1/L$ でしか減らない．$0.1$ eV 以下にするには

$$ L \gt \frac{\alpha_M}{2\times(0.1/27.211)}\ \mathrm{Bohr} = 386\ \mathrm{Bohr} = 204\ \text{Å} $$

が必要で，原子数にして $10^6$ 個規模になる．式\eqref{eq:18-makov-payne}で補正を掛ければ改善するが，高次項の評価が難しく，$0.1$ eV の再現性は保証できない．絶対値を狙う以上，補償背景では勝負にならない．これが18.4.3節の厳密クーロンカットオフ法を必要とする理由である．

図18.5 (a) 周期境界条件の下では，コアホールは無限個の周期像とともに現れ，像どうしのCoulomb相互作用がエネルギーに混入する．(b) 有限ギャップ系では，一般式\eqref{eq:18-eb-general}を厳密クーロンカットオフとともに用いると内殻ホール間距離 $15$–$20$ Å で収束する．一方，金属用の式\eqref{eq:18-eb-metal}を有限ギャップ系に適用すると，どこまでセルを大きくしても $1$ eV 規模の系統誤差が残る．

18.4.2 コアホールを局在させる:ペナルティ汎関数法

次に，狙った原子の狙った内殻軌道に空孔を作る方法を考える．素朴には「内殻軌道の占有数を $0$ に固定すればよい」のだが，自己無撞着計算の途中で軌道は混ざり合うので，どの解が「内殻軌道」なのかを毎回同定しなければならない．しかも周期系では，内殻準位はきわめて狭いバンドを作り，複数の等価な原子があると空孔が分散してしまう(1個の原子に $1$ の空孔ではなく，$M$ 個の原子に $1/M$ ずつ，という非物理的な解に落ちる)．

これを一挙に解決するのがペナルティ汎関数法である．アイデアは単純で，「狙った内殻軌道を占有すると莫大なエネルギー損をする」ようにエネルギー汎関数を書き換えてしまう．

定義:ペナルティ汎関数

狙った原子の狙った内殻軌道の(規格化された)波動関数を $\ket{\chi_c}$ とする．これを使って射影演算子

$$ \begin{equation} \hat{P} = \ket{\chi_c}\,\Delta\,\bra{\chi_c} \label{eq:18-projector} \end{equation} $$

を作る($\Delta$ は正の大きな定数，実用値は $\Delta \sim 100$ Rydberg $= 1361$ eV)．全エネルギー汎関数を

$$ \begin{equation} E_f = E_{\mathrm{DFT}}\bigl[\{\psi\}\bigr] + E_{\mathrm{pen}}, \qquad E_{\mathrm{pen}} = \frac{1}{V_B}\int_B\!\dd^3k \sum_{\mu} f_{\mu}^{(\kk)} \braket{\psi_\mu^{(\kk)}|\hat{P}|\psi_\mu^{(\kk)}} \label{eq:18-penalty} \end{equation} $$

と拡張する($V_B$ はBrillouin域の体積，$f_\mu^{(\kk)}$ は占有数)．$E_{\mathrm{pen}}$ は「占有された状態が $\ket{\chi_c}$ 成分をどれだけ含むか」を $\Delta$ 倍して足し上げた量であり，常に非負である．

この汎関数を変分すると，Kohn–Sham方程式に射影演算子が1個加わるだけである．以下でそれを確かめる．

導出:ペナルティ項が生む修正Kohn–Sham方程式

ステップ1:変分問題を立てる．規格直交条件 $\braket{\psi_\mu|\psi_\nu} = \delta_{\mu\nu}$ のもとで $E_f$ を停留させる．Lagrangeの未定乗数 $\lambda_{\mu\nu}$ を導入し(第5章・第10章と同じ手続き)，

$$ \begin{equation} \Xi = E_f - \sum_{\mu\nu}\lambda_{\mu\nu}\bigl(\braket{\psi_\mu|\psi_\nu} - \delta_{\mu\nu}\bigr) \label{eq:18-pen-lag} \end{equation} $$

を $\psi_\mu^*$ について汎関数微分してゼロと置く($\kk$ の添字は煩雑なので省く)．

ステップ2:DFT部分の微分．第10章で示したとおり

$$ \begin{equation} \frac{\delta E_{\mathrm{DFT}}}{\delta\psi_\mu^*(\rr)} = f_\mu\,\hat{h}_{\mathrm{KS}}\,\psi_\mu(\rr), \qquad \hat{h}_{\mathrm{KS}} = -\tfrac12\nabla^2 + v_{\mathrm{eff}}[\rho](\rr) \label{eq:18-pen-dft} \end{equation} $$

である．有効ポテンシャル $v_{\mathrm{eff}}$ は密度 $\rho$ を通じて $\psi$ に依存するが，その依存性は $\delta E_{\mathrm{DFT}}/\delta\rho = v_{\mathrm{eff}}$ という関係にまとめられている(連鎖律)．

ステップ3:ペナルティ項の微分．$E_{\mathrm{pen}} = \sum_\mu f_\mu\braket{\psi_\mu|\hat{P}|\psi_\mu}$ は $\psi_\mu^*$ について1次なので，微分は自明である:

\begin{align} \frac{\delta E_{\mathrm{pen}}}{\delta\psi_\mu^*(\rr)} &= f_\mu\frac{\delta}{\delta\psi_\mu^*(\rr)}\int\!\!\int \psi_\mu^*(\rr_1)P(\rr_1,\rr_2)\psi_\mu(\rr_2)\dd^3r_1\dd^3r_2 \label{eq:18-pen-d1}\\[2pt] &= f_\mu\int P(\rr,\rr_2)\psi_\mu(\rr_2)\dd^3r_2 \;=\; f_\mu\,\hat{P}\psi_\mu(\rr) . \label{eq:18-pen-d2} \end{align}

ここで $P(\rr_1,\rr_2) = \chi_c(\rr_1)\Delta\chi_c^*(\rr_2)$ は射影演算子の核である．\eqref{eq:18-pen-d1}から\eqref{eq:18-pen-d2}へは $\delta\psi_\mu^*(\rr_1)/\delta\psi_\mu^*(\rr) = \delta(\rr_1-\rr)$ を使っただけである．$\hat{P}$ がエルミートなので，$\psi_\mu$ についての微分からも同じ項が(複素共役で)出る．

ステップ4:停留条件を書く．\eqref{eq:18-pen-dft}と\eqref{eq:18-pen-d2}を\eqref{eq:18-pen-lag}に入れて

$$ \begin{equation} f_\mu\bigl(\hat{h}_{\mathrm{KS}} + \hat{P}\bigr)\psi_\mu = \sum_\nu\lambda_{\mu\nu}\psi_\nu . \label{eq:18-pen-stat} \end{equation} $$

右辺の $\lambda$ 行列はエルミートなのでユニタリ変換で対角化でき(第5章と同じ議論:占有数が等しい軌道の部分空間内での回転は密度も $E_f$ も変えない)，正準形

$$ \begin{equation} \bigl(\hat{h}_{\mathrm{KS}} + \hat{P}\bigr)\ket{\psi_\mu} = \varepsilon_\mu\ket{\psi_\mu} \label{eq:18-pen-ks} \end{equation} $$

を得る．∎

式\eqref{eq:18-pen-ks}の効果を読み取ろう．$\ket{\chi_c}$ は $\hat{P}$ の固有値 $\Delta$ の固有ベクトルである($\hat{P}\ket{\chi_c} = \Delta\ket{\chi_c}$，$\braket{\chi_c|\chi_c}=1$ より)．もともとの内殻準位はきわめて局在しており $\ket{\chi_c}$ とほとんど一致するから，その固有値は

$$ \varepsilon_c \;\longrightarrow\; \varepsilon_c + \Delta $$

と押し上げられる．炭素 $1s$ なら $-270\ \mathrm{eV} + 1361\ \mathrm{eV} = +1091$ eV となり，フェルミ準位のはるか上に飛ぶ．したがってアウフバウ原理に従えばこの状態は占有されない．抜けた1電子は，金属なら\eqref{eq:18-eb-metal}に従ってフェルミ準位に置かれ，絶縁体なら\eqref{eq:18-eb-general}に従ってセルから取り去られる．こうしてコアホールが自己無撞着に生成・維持される．

物理的意味:ペナルティ項はエネルギーに寄与しない

$E_f = E_{\mathrm{DFT}} + E_{\mathrm{pen}}$ と書いたが，我々が\eqref{eq:18-eb-general}に代入すべきなのは物理的な全エネルギーであって，人工的な $E_{\mathrm{pen}}$ を含む量ではない．ところが心配は要らない．自己無撞着解に到達すると $E_{\mathrm{pen}}$ は自動的にほぼゼロになる．なぜなら，占有されたどの状態も $\ket{\chi_c}$ 成分をほとんど持たないからである($\braket{\chi_c|\psi_\mu}$ の残差は典型的に $10^{-3}$ 以下，$E_{\mathrm{pen}} \sim \Delta\times10^{-6}$ Hartree)．$\hat{P}$ は「解を選ぶための道具」であって，エネルギーを歪める項ではない．これは第10章で見た「Lagrange未定乗数は拘束を課すが，拘束が満たされた解ではエネルギーに寄与しない」という構造とまったく同じである．

もう一つの利点は局在性である．$\ket{\chi_c}$ は特定の原子に張り付いた関数なので，$\hat{P}$ はその原子の内殻だけを押し上げる．等価な原子が何個あっても，空孔は狙った1個の原子に完全に局在する．

補足:スピン軌道分裂を含めるには

$2p$ 以上の内殻では，スピン軌道相互作用による $j = l\pm1/2$ の分裂を再現したい(Si $2p$ で $0.6$ eV)．そのためには，相対論的擬ポテンシャル(第15章)と2成分スピノール形式を使い，射影関数 $\ket{\chi_c}$ を原子Dirac方程式の角度解で作る．$J = l+1/2$，$M = m+1/2$ に対して

$$ \begin{equation} \ket{\Phi_{l+1/2}^{\,m+1/2}} = \sqrt{\frac{l+m+1}{2l+1}}\,\ket{Y_l^{m}}\ket{\alpha} + \sqrt{\frac{l-m}{2l+1}}\,\ket{Y_l^{m+1}}\ket{\beta}, \label{eq:18-dirac-plus} \end{equation} $$

$J = l-1/2$，$M = m+1/2$ に対して

$$ \begin{equation} \ket{\Phi_{l-1/2}^{\,m+1/2}} = -\sqrt{\frac{l-m}{2l+1}}\,\ket{Y_l^{m}}\ket{\alpha} + \sqrt{\frac{l+m+1}{2l+1}}\,\ket{Y_l^{m+1}}\ket{\beta} \label{eq:18-dirac-minus} \end{equation} $$

である($\ket{\alpha},\ket{\beta}$ は上向き・下向きスピン，全体の位相は任意)．係数の2乗和が $\frac{l+m+1}{2l+1}+\frac{l-m}{2l+1} = \frac{2l+1}{2l+1} = 1$ となり規格化されていること，$\ket{\Phi_{l+1/2}}$ と $\ket{\Phi_{l-1/2}}$ が直交することは直接確かめられる．射影関数を $\ket{\chi_c} = \ket{R_{nl}\,\Phi_J^M}$($R_{nl}$ は内殻の動径関数)と取れば，$j$ ごとに別々のコアホール計算ができ，$2p_{1/2}$ と $2p_{3/2}$ の絶対束縛エネルギーが個別に求まる．

補足:擬ポテンシャルにコアホールを埋め込む方法との違い

より古典的な処方として，あらかじめ「$1s$ 電子を1個抜いた配置」で原子の擬ポテンシャルを生成し，狙った原子だけをその擬ポテンシャルで置き換える方法がある(コアホール擬ポテンシャル法)．実装が簡単で，XANESのスペクトル形状を求めるには十分実用的である．しかし絶対束縛エネルギーには使えない．理由は，始状態用と終状態用の擬ポテンシャルが異なる原子計算から作られており，エネルギーの原点が揃っていないからである．擬ポテンシャルには任意定数の自由度があり，それが $E_i^{(0)}$ と $E_f^{(0)}$ で違えば，その差である $E_B$ に直接効いてしまう．同種原子どうしの相対シフト(化学シフト)なら定数がキャンセルするので使えるが，絶対値は出ない．ペナルティ汎関数法では始状態と終状態でまったく同じ擬ポテンシャルを使うので，この問題が原理的に生じない．

18.4.3 厳密クーロンカットオフ法

残る問題は荷電セルである．18.4.1節で見たとおり，補償背景では $1/L$ の誤差が消えない．ここでは，Coulomb相互作用そのものを有限距離で切ることで，周期像を厳密に切り離す方法を導く．

まず密度を分解する．終状態の密度 $\rho_f$ を，始状態の密度 $\rho_i$(中性・周期的)とその差分に分ける:

$$ \begin{equation} \rho_f(\rr) = \rho_i(\rr) + \Delta\rho(\rr), \qquad \int\Delta\rho(\rr)\dd^3r = -1 \label{eq:18-rho-split} \end{equation} $$

(電子が1個減るので $\Delta\rho$ の積分は $-1$．電子密度を正の量として数えている)．Hartreeポテンシャルは線形汎関数なので，対応して

$$ \begin{equation} V_H[\rho_f] = \underbrace{V_H[\rho_i]}_{\text{周期成分 }V_H^{(\mathrm{P})}} + \underbrace{V_H[\Delta\rho]}_{\text{非周期成分 }V_H^{(\mathrm{NP})}} \label{eq:18-vh-split} \end{equation} $$

と分かれる．$V_H^{(\mathrm{P})}$ は中性密度から作られるので，通常どおり高速Fourier変換(FFT)で問題なく計算できる．厄介なのは正味電荷を持つ $\Delta\rho$ の側だけである．

ここで重要な物理的観察がある．$\Delta\rho$ は狙った原子のまわりに局在している．コアホールという点電荷ができると，価電子がそこへ流れ込んで遮蔽するので，遮蔽電荷を含めた $\Delta\rho$ 全体の広がりは有限である．シリコンの $2p$ コアホールの場合，$\Delta\rho$ の球面平均は半径 $7$ Å 程度でほぼゼロになる(遮蔽の主役は同じ原子上の $3p$ 電子である)．

導出:打ち切りCoulomb核のFourier変換

ステップ1:核を定義する．カットオフ半径 $R_c$ を導入し，

$$ \begin{equation} v_c(r) = \begin{cases} \dfrac{1}{r} & (r \le R_c)\\[4pt] 0 & (r \gt R_c)\end{cases} \label{eq:18-vc-def} \end{equation} $$

と定める．

ステップ2:角度積分を実行する．Fourier変換を球座標で書き，$\bm{G}$ を極軸に取る．$u = \cos\theta$ と置くと $\dd\Omega = 2\pi\,\dd u$ であり，

\begin{align} \int\dd\Omega\, e^{-i\bm{G}\cdot\rr} &= 2\pi\int_{-1}^{1} e^{-iGru}\,\dd u \label{eq:18-ang1}\\[2pt] &= 2\pi\left[\frac{e^{-iGru}}{-iGr}\right]_{u=-1}^{u=1} \label{eq:18-ang2}\\[2pt] &= 2\pi\cdot\frac{e^{-iGr}-e^{iGr}}{-iGr} = 2\pi\cdot\frac{-2i\sin(Gr)}{-iGr} \label{eq:18-ang3}\\[2pt] &= \frac{4\pi\sin(Gr)}{Gr} . \label{eq:18-ang4} \end{align}

\eqref{eq:18-ang2}は指数関数の積分，\eqref{eq:18-ang3}ではEulerの公式 $e^{i x}-e^{-i x} = 2i\sin x$ を使った．

ステップ3:動径積分を実行する．\eqref{eq:18-ang4}を使って

\begin{align} \tilde{v}_c(\bm{G}) &= \int_{r\le R_c} \frac{1}{r}\,e^{-i\bm{G}\cdot\rr}\,\dd^3r \;=\; \int_0^{R_c}\dd r\,r^2\cdot\frac{1}{r}\cdot\frac{4\pi\sin(Gr)}{Gr} \label{eq:18-rad1}\\[2pt] &= \frac{4\pi}{G}\int_0^{R_c}\sin(Gr)\,\dd r \label{eq:18-rad2}\\[2pt] &= \frac{4\pi}{G}\left[\frac{-\cos(Gr)}{G}\right]_0^{R_c} \label{eq:18-rad3}\\[2pt] &= \frac{4\pi}{G^2}\bigl(1-\cos(G R_c)\bigr) . \label{eq:18-vc-G} \end{align}

\eqref{eq:18-rad1}で $r^2\cdot\frac1r\cdot\frac1{Gr} = \frac{1}{G}$ と約分し，\eqref{eq:18-rad3}で $\sin$ の不定積分を実行した．

ステップ4:$\bm{G}\to\bm{0}$ の極限を調べる．ここが要点である．$\cos x \simeq 1 - x^2/2$ を使うと

\begin{align} \lim_{G\to0}\tilde{v}_c(\bm{G}) &= \lim_{G\to0}\frac{4\pi}{G^2}\left[1-\left(1-\frac{(GR_c)^2}{2}+O(G^4)\right)\right] \label{eq:18-lim1}\\[2pt] &= \lim_{G\to0}\frac{4\pi}{G^2}\cdot\frac{G^2R_c^2}{2} \;=\; 2\pi R_c^2 . \label{eq:18-lim2} \end{align}

有限である．通常のCoulomb核 $4\pi/G^2$ が $G\to0$ で発散するのに対し，打ち切り核は $G=0$ でも有限値 $2\pi R_c^2$ を持つ．したがって $\bm{G}=\bm{0}$ 項を落とす必要も，補償背景を入れる必要もない．∎

$$ \begin{equation} \tilde{v}_c(\bm{G}) = \frac{4\pi}{G^2}\bigl(1-\cos(G R_c)\bigr), \qquad \tilde{v}_c(\bm{0}) = 2\pi R_c^2 \label{eq:18-cutoff-key} \end{equation} $$

あとは，この核を使って非周期成分を

$$ \begin{equation} V_H^{(\mathrm{NP})}(\rr) = \frac{1}{\Omega}\sum_{\bm{G}}\widetilde{\Delta\rho}(\bm{G})\,\tilde{v}_c(\bm{G})\,e^{i\bm{G}\cdot\rr} \label{eq:18-vhnp} \end{equation} $$

と計算すればよい．FFTがそのまま使えるので，計算コストは通常のHartreeポテンシャルと変わらない．

ただし，この処方が厳密であるためには条件がある．それを確かめておこう．

導出:厳密性の条件 $R_c = 2R$ かつ $L \gt 4R$

$\Delta\rho$ が原子位置を中心とする半径 $R$ の球内に完全に含まれているとする．求めたいのは孤立系の静電エネルギー

$$ E_H^{\mathrm{iso}} = \frac{1}{2}\int\!\!\int \frac{\Delta\rho(\rr)\Delta\rho(\rr')}{\abs{\rr-\rr'}}\dd^3r\,\dd^3r' $$

である．一方，FFTで\eqref{eq:18-vhnp}を計算して得られるのは，周期像をすべて含んだ巡回畳み込み

$$ E_H^{\mathrm{FFT}} = \frac{1}{2}\sum_{\bm{T}}\int\!\!\int \Delta\rho(\rr)\,v_c(\rr-\rr'-\bm{T})\,\Delta\rho(\rr')\dd^3r\,\dd^3r' $$

である($\bm{T}$ は格子ベクトル)．両者が一致するために必要な条件は2つある．

条件1($\bm{T}=\bm{0}$ の項が正しいこと)．$\rr$ と $\rr'$ はともに半径 $R$ の球内にあるから，三角不等式より $\abs{\rr-\rr'} \le 2R$ である．この範囲で $v_c(\rr-\rr') = 1/\abs{\rr-\rr'}$ となるためには

$$ R_c \ge 2R $$

が必要十分である．

条件2($\bm{T}\ne\bm{0}$ の項が消えること)．格子ベクトルの最短長を $L$ とすると，やはり三角不等式より

$$ \abs{\rr-\rr'-\bm{T}} \ge \abs{\bm{T}} - \abs{\rr-\rr'} \ge L - 2R $$

である．これが $R_c$ より大きければ $v_c = 0$ となり像の寄与は厳密にゼロになる:

$$ L - 2R \gt R_c . $$

両条件を組み合わせる．$R_c$ を最小値 $2R$ に取れば，条件2は $L - 2R \gt 2R$，すなわち

$$ \begin{equation} R_c = 2R,\qquad L \gt 4R \label{eq:18-cutoff-cond} \end{equation} $$

となる．$\Delta\rho$ の広がりの4倍以上のセルを用意すれば，静電エネルギーは厳密に孤立系の値になる．∎

物理的意味:なぜ「補償背景」ではなく「核を切る」のか

補償背景法と厳密クーロンカットオフ法の決定的な違いは，電位の基準がずれるかどうかである．補償背景を入れると，セル全体に一様なポテンシャルオフセットが加わり，しかもその大きさは電荷 $Q$ とセル体積に依存する．したがって $Q=0$ の始状態計算と $Q=+1$ の終状態計算で電位の原点が食い違い，その差である $E_B$ に直接誤差として入ってしまう．18.3.1節の導出で「$E_i^{(0)}$, $E_f^{(0)}$, $\mu_0$ を同じ基準で計算せよ」と強調したのは，まさにこの点である．厳密クーロンカットオフ法では，$\bm{G}=\bm{0}$ 成分を人工的に操作しないので，始状態と終状態が最初から同じ電位基準の上に乗っている．

18.4.4 収束の確認

実務上，確認すべきパラメータは(基底関数・$\kk$ 点・実空間格子といった通常のものに加えて)ただ一つ，コアホール間距離である．周期境界条件を使う以上，コアホールは周期的に並ぶ．厳密クーロンカットオフは静電的な像相互作用を消してくれるが，電子構造そのものの像間干渉——隣のコアホールが作る歪んだ電荷分布が重なる効果——は残る．したがってセルを系統的に大きくして $E_B$ が収束することを確認しなければならない．

系	必要なコアホール間距離	おおよその原子数
c-BN(絶縁体)	$\simeq 15$ Å	数百
固体 $\mathrm{NH_3}$(分子性結晶)	$\simeq 20$ Å	数百
Si(半導体)	$\simeq 27$ Å	$\simeq 500$
グラフェン・TiN・TiC(金属・半金属)	金属式\eqref{eq:18-eb-metal}で $\simeq 10$ Å	$\simeq 64$

収束が半導体で遅く金属で速い理由は明快である．金属では自由電子が即座にコアホールを遮蔽するので $\Delta\rho$ の広がりが小さく，しかも中性セル計算\eqref{eq:18-eb-metal}が使えるので静電的な問題自体が生じない．半導体・絶縁体では遮蔽が誘電応答(有限の $\varepsilon_1$)を通じてしか起こらず，$\Delta\rho$ が長い尾を引く．$\varepsilon_1 \simeq 12$ のシリコンで $\Delta\rho$ の有効半径が $7$ Å，条件\eqref{eq:18-cutoff-cond}の $L \gt 4R = 28$ Å は，上表の $27$ Å とよく整合している．

18.5 ベンチマークと解釈

技術が揃ったので，実際にどこまで当たるのかを見る．そのうえで，得られた数値を「なぜその値なのか」に翻訳する解釈の枠組みを作る．

18.5.1 絶対束縛エネルギーの精度

内殻束縛エネルギーは絶対値が数百 eV の量である．そこに $0.1$–$0.5$ eV の精度で切り込めるかどうかが，この手法の価値を決める．代表的な固体について，式\eqref{eq:18-eb-general}(有限ギャップ系)および\eqref{eq:18-eb-metal}(金属・半金属)にペナルティ汎関数法と厳密クーロンカットオフ法を組み合わせて得られる結果を，実験値と並べてみよう．

系	準位	計算値 (eV)	実験値 (eV)	差 (eV)
c-BN	N $1s$	$398.87$	$398.1$	$+0.77$
固体 $\mathrm{NH_3}$	N $1s$	$398.92$	$399.0$	$-0.08$
ダイヤモンド	C $1s$	$286.50$	$285.6$	$+0.90$
Si	$2p_{1/2}$	$100.13$	$99.8$	$+0.33$
Si	$2p_{3/2}$	$99.40$	$99.2$	$+0.20$
グラフェン	C $1s$	$284.23$	$284.4$	$-0.17$
TiN	N $1s$	$396.43$	$397.1$	$-0.67$
TiC	C $1s$	$281.43$	$281.5$	$-0.07$

平均絶対誤差は約 $0.4$ eV，絶対値に対する平均相対誤差は $0.16\,\%$ である．気相分子($\mathrm{CO}$, $\mathrm{N_2}$, $\mathrm{N_2O}$, $\mathrm{SiH_4}$, $\mathrm{SiF_4}$ など23例)についても同様の計算を行うと，平均絶対誤差 $0.5$ eV，相対誤差 $0.22\,\%$ が得られる．$400$ eV の量を $0.5$ eV の精度で当てるということは，5桁目まで合っているということである．

とくに注目すべきは Si $2p_{1/2}$ と $2p_{3/2}$ の分裂で，計算値 $100.13 - 99.40 = 0.73$ eV に対し実験値 $99.8 - 99.2 = 0.6$ eV である．相対論的擬ポテンシャルと2成分スピノールによる射影(式\eqref{eq:18-dirac-plus},\eqref{eq:18-dirac-minus})が，スピン軌道分裂を定量的に再現していることが分かる．

例:GGAの限界が見えるとき

すべてがうまくいくわけではない．$\mathrm{O_2}$ 分子($S=1$)の O $1s$ 交換分裂(18.2.3節の数学ノート参照)を計算すると，GGA-PBE では $0.5$ eV となり，実験値 $1.1$ eV を大幅に過小評価する．コアホールが残した不対スピンと分子の三重項スピンとの交換相互作用を，近似交換相関汎関数が正しく記述できていないためである．絶対束縛エネルギーは $0.5$ eV で当たるのに，交換分裂は半分しか出ない——同じ計算の中でも，量によって精度が違うことに注意しなければならない．一般に，スピン分極や局在した $d$/$f$ 電子が関わる量では，汎関数の限界が露呈しやすい．

18.5.2 化学シフトの解釈:電荷ポテンシャル模型

計算値が実験と合ったとして，次に問われるのは「なぜその値なのか」である．ここでは化学シフトを3つの寄与に分解する，古典的だが今なお有用な枠組みを導く．

導出:電荷ポテンシャル模型

ステップ1:内殻電子が感じる静電ポテンシャルを2つに分ける．内殻電子は原子核のごく近く($1s$ なら半径 $1/Z$ 程度)にいる．その位置での静電ポテンシャルへの寄与は，(i) 同じ原子の価電子から，(ii) 他の原子から，の2つに分けられる(原子核自身の寄与は化学環境によらないので，基準に繰り込む)．

ステップ2:自分の価電子からの寄与．原子 $A$ の価電子が半径 $r_v$ の球殻上に一様に分布し，その数を $N_v$ とする．半径 $r_v$，全電荷 $-N_v$ の球殻が中心に作る静電ポテンシャルは，Gaussの法則より

$$ \phi_{\text{valence}} = \frac{-N_v}{r_v} $$

である(球殻の内側ではポテンシャルは一定で，殻上の値に等しい)．電荷 $-1$ の内殻電子が持つポテンシャルエネルギーは

$$ \begin{equation} U_{\text{onsite}} = (-1)\times\phi_{\text{valence}} = \frac{N_v}{r_v} \label{eq:18-pot-onsite} \end{equation} $$

で，正の量(価電子との反発)である．

ステップ3:価電子が $q_A$ だけ減ったときの変化．原子 $A$ が正味電荷 $+q_A$ を帯びる，すなわち価電子が $q_A$ 個減ると，$N_v \to N_v - q_A$ となり

$$ \begin{equation} \Delta U_{\text{onsite}} = -\frac{q_A}{r_v} \equiv -k\,q_A , \qquad k \equiv \frac{1}{r_v} \label{eq:18-pot-k} \end{equation} $$

だけポテンシャルエネルギーが下がる．内殻固有値も同じだけ下がるので，$E_B = -\varepsilon_c$ は $+k q_A$ だけ増える．$k$ の大きさは，炭素の価電子半径 $r_v \simeq 1.6$ Bohr から $k \simeq 0.63$ Hartree $\simeq 17$ eV(電子1個あたり)と見積もられる．

ステップ4:他の原子からの寄与(Madelungポテンシャル)．他の原子 $B$ が正味電荷 $q_B$ を持てば，$A$ の位置に静電ポテンシャル $q_B/R_{AB}$ を作る．内殻電子のポテンシャルエネルギーは $-q_B/R_{AB}$，したがって $E_B$ への寄与は $+q_B/R_{AB}$ である．すべての $B$ について足して

$$ \begin{equation} V_M(A) = \sum_{B\ne A}\frac{q_B}{R_{AB}} \label{eq:18-madelung} \end{equation} $$

をMadelungポテンシャルと呼ぶ．

ステップ5:終状態効果を加えて完成．式\eqref{eq:18-eb-init-final}の分解を思い出し，緩和エネルギー $R_A$ を引く．以上をまとめて

$$ \begin{equation} E_B(A) = E_B^{0} + k\,q_A + V_M(A) - R_A \label{eq:18-potmodel} \end{equation} $$

を得る．$E_B^0$ は基準となる中性原子の値である．∎

$$ \begin{equation} E_B(A) = E_B^{0} + \underbrace{k\,q_A}_{\text{自分の電荷}} + \underbrace{\sum_{B\ne A}\frac{q_B}{R_{AB}}}_{\text{Madelung}} - \underbrace{R_A}_{\text{緩和}} \label{eq:18-potmodel-key} \end{equation} $$

この式は，化学シフトについて実験家が経験的に知っていることをきれいに説明する．

酸化数との相関．形式酸化数が上がれば $q_A$ が増え，$E_B$ が上がる．図18.6の硫黄がその典型で，$\mathrm{S^{2-}}$(硫化物)から $\mathrm{S^{6+}}$(硫酸)まで $S\,2p_{3/2}$ が単調に $7$ eV 上昇する．
傾きが $k$ より小さい理由．ステップ3で $k \simeq 17$ eV/電子と見積もったのに，実測の傾きは酸化数1あたり $1$–$2$ eV にすぎない．理由は2つある．第一に，実際に移動する電荷は形式酸化数よりずっと小さい(共有結合性が高いほど $q_A \ll$ 酸化数)．第二に，$A$ が正に帯電すれば周囲 $B$ は負に帯電するので $V_M(A) \lt 0$ となり，$k q_A$ を部分的に打ち消す．このMadelung相殺は，イオン結晶で特に顕著で，$k q_A$ の半分以上を食うことがある．
金属と絶縁体の違い．$R_A$ は遮蔽の効きやすさに比例する．金属では価電子が自由に動くので $R_A$ が大きく($10$ eV 規模)，絶縁体では小さい．同じ元素が金属状態と酸化物状態にあるとき，観測される化学シフトには $\Delta R_A$ も含まれている．

硫黄の形式酸化数と S 2p3/2 束縛エネルギーの相関．S(2−)，S(0)，S(4+)，S(6+) の 4 点が酸化数とともに 161.5 eV から 168.5 eV へ上がり，赤い破線は傾き約 0.9 eV/酸化数の直線． — 図18.6 硫黄 $2p_{3/2}$ 束縛エネルギーと形式酸化数の相関(概略値)．式\eqref{eq:18-potmodel-key}の $k q_A$ 項が支配的で，酸化数が上がるほど束縛エネルギーが上がる．傾きが $k \simeq 17$ eV/電子よりはるかに小さいのは，実際の電荷移動が形式酸化数より小さいことと，Madelung項による相殺の両方による．

物理的意味:「有効電荷」で束縛エネルギーを整理する

式\eqref{eq:18-potmodel-key}の $q_A$ は「原子上の電荷」だが，量子力学的にはこれは一意に決まらない量である(電子密度を原子に切り分ける方法は無数にある)．それでも，Wannier関数や Mulliken/Bader 解析(いずれも本書では扱わない)で定義した有効電荷を使うと，$E_B$ の系統性は驚くほどよく整理される．とくに，始状態の有効電荷 $N_i$ と終状態(コアホールあり)の有効電荷 $N_f$ の和 $N_i + N_f$ を横軸に取ると，同一元素の異なる化合物における内殻束縛エネルギーがほぼ直線に乗ることが知られている．$N_i$ が始状態効果，$N_f$ が終状態(遮蔽)効果を表しているという読み方ができる．式\eqref{eq:18-eb-integral}が「$\varepsilon_c(n)$ を $n=0$ から $1$ まで平均せよ」と教えていたことを思い出せば，始状態($n=1$)と終状態($n=0$)の情報が同じ重みで入るのは自然である．

18.6 XANESと励起状態のDFT

図18.1の2本目の矢印に移る．内殻電子を真空まで叩き出すのではなく，物質内の非占有状態へ持ち上げる過程である．光子エネルギーを掃引しながら吸収係数を測ると，ある値($=$ 内殻から最低非占有状態への励起エネルギー)で吸収が急激に立ち上がる．これが吸収端である．端から $50$ eV 程度までの微細構造をXANES(X-ray absorption near edge structure)またはNEXAFS(near edge X-ray absorption fine structure)と呼び，さらに高エネルギー側の緩やかな振動をEXAFSと呼ぶ．

本節では，(i) 基底状態の理論であるDFTで励起状態を扱ってよいのかという根本問題に答え，(ii) 吸収スペクトルを与える公式を黄金律から導き，(iii) 双極子選択則 $\Delta l = \pm1$ を導出して，XANESが「局所の $l\pm1$ 部分状態密度」を映していることを示す．

図18.7 炭素K端を模したX線吸収スペクトルの模式図．吸収端の直前・直後に現れる鋭い共鳴(π*，σ*)が XANES/NEXAFS 領域，端から $50$ eV 以上高エネルギー側の緩やかな振動が EXAFS 領域である．XANES は非占有状態の局所的な電子構造(結合の種類・対称性・酸化状態)を，EXAFS は光電子波の後方散乱干渉を通じて原子間距離を反映する．

18.6.1 基底状態理論で励起を扱えるか:Gunnarsson–Lundqvistの議論

正面から問題を立てよう．Hohenberg–Kohnの定理(第9章)は基底状態について定式化されている．基底状態密度が外部ポテンシャルを決め，エネルギー汎関数の最小値が基底状態エネルギーを与える——それだけである．ところがXANESの終状態は，内殻に空孔があり電子が伝導帯にいる励起状態であって，系全体の基底状態ではない．この状態にDFTを適用してよい理由はあるのか．

答えは「ある．ただし条件付きで」である．鍵になるのは，次のごく一般的な変分原理である．

定理:対称性クラスごとの変分原理(Gunnarsson–Lundqvist)

ハミルトニアン $\hat{H}$ と可換な観測量 $\hat{O}$($[\hat{H},\hat{O}]=0$)があるとする．$\hat{O}$ の固有値 $\lambda^{(A)}$ を持つ $\hat{H}$ の固有状態全体をクラス $A$ と呼び，それ以外をクラス $B$ と呼ぶ．このとき，クラス $A$ に属する試行波動関数 $\Psi$ に対して

$$ \begin{equation} \frac{\braket{\Psi|\hat{H}|\Psi}}{\braket{\Psi|\Psi}} \;\ge\; E_{\min}^{(A)} \label{eq:18-gl-var} \end{equation} $$

が成り立つ．ここで $E_{\min}^{(A)}$ はクラス $A$ の最低エネルギーである．すなわち，各対称性クラスの最低状態は，それが系全体の基底状態でなくても，変分的に決定できる．

導出:クラスごとの変分原理と，射影ハミルトニアンによる無拘束化

ステップ1:展開する．$\hat{H}$ の固有状態は完全系をなす．$[\hat{H},\hat{O}]=0$ なので，固有状態は $\hat{O}$ の固有値でラベル付けでき，クラス $A$ の固有状態 $\{\Psi_i^{(A)}\}$ とクラス $B$ の固有状態 $\{\Psi_j^{(B)}\}$ に分けられる．クラス $A$ の任意の試行関数は

$$ \ket{\Psi} = \sum_i c_i\ket{\Psi_i^{(A)}} $$

と展開できる(クラス $B$ の成分は $\hat{O}$ の固有値が違うので入らない)．

ステップ2:期待値を計算する．$\hat{H}\ket{\Psi_i^{(A)}} = E_i^{(A)}\ket{\Psi_i^{(A)}}$ と固有状態の直交性から

\begin{align} \braket{\Psi|\hat{H}|\Psi} &= \sum_{i,i'} c_{i'}^*c_i E_i^{(A)}\braket{\Psi_{i'}^{(A)}|\Psi_i^{(A)}} = \sum_i \abs{c_i}^2 E_i^{(A)}, \label{eq:18-gl1}\\ \braket{\Psi|\Psi} &= \sum_i \abs{c_i}^2 . \label{eq:18-gl2} \end{align}

ステップ3:重み付き平均は最小値以上．$\abs{c_i}^2 \ge 0$ かつ $E_i^{(A)} \ge E_{\min}^{(A)}$ だから

\begin{align} \sum_i \abs{c_i}^2 E_i^{(A)} \;\ge\; \sum_i \abs{c_i}^2 E_{\min}^{(A)} = E_{\min}^{(A)}\sum_i\abs{c_i}^2 \label{eq:18-gl3} \end{align}

であり，両辺を $\braket{\Psi|\Psi} = \sum_i\abs{c_i}^2 \gt 0$ で割ると\eqref{eq:18-gl-var}を得る．等号はクラス $A$ の最低状態そのもの($c_i = \delta_{i,\min}$)のときに限る．∎

ステップ4(補足):射影ハミルトニアン．実際の計算では「クラス $A$ に属する試行関数だけを動かす」という拘束が扱いにくい．そこで射影演算子

$$ \begin{equation} \hat{A} = \sum_{i\in A}\ket{\Psi_i^{(A)}}\bra{\Psi_i^{(A)}}, \qquad \hat{H}_A \equiv \hat{A}^\dagger\hat{H}\hat{A} \label{eq:18-gl-proj} \end{equation} $$

を導入する．任意の(クラス $B$ 成分も含む)試行関数 $\ket{\Psi} = \sum_i c_i\ket{\Psi_i^{(A)}} + \sum_j d_j\ket{\Psi_j^{(B)}}$ について，$\hat{A}\ket{\Psi} = \sum_i c_i\ket{\Psi_i^{(A)}}$ だから

$$ \frac{\braket{\Psi|\hat{H}_A|\Psi}}{\braket{\Psi|\Psi}} = \frac{\sum_i\abs{c_i}^2E_i^{(A)}}{\sum_i\abs{c_i}^2 + \sum_j\abs{d_j}^2} $$

となる．エネルギーの原点を全固有値が負になるようにずらしておけば，分母が大きくなるほど比は大きく(ゼロに近く)なるので，最小化は $d_j = 0$ を選ぶ．すなわち無拘束の最小化がクラス $A$ の最低状態を与える．18.4.2節のペナルティ汎関数 $\hat{P}$ は，まさにこの射影を実装する道具にほかならない．

次に，この変分原理を密度汎関数の言葉に翻訳する．第9章で学んだLevyの制約付き探索(constrained search)を，クラス $A$ に制限して繰り返せばよい．

導出:クラス $A$ に制限した密度汎関数

ハミルトニアンを $\hat{H} = \hat{T} + \hat{V}_{ee} + \hat{V}_{\mathrm{ext}}$ と分ける．クラス $A$ の最低エネルギーは，定理\eqref{eq:18-gl-var}より

\begin{align} E_{\min}^{(A)} &= \min_{\Psi\in A}\braket{\Psi|\hat{T}+\hat{V}_{ee}+\hat{V}_{\mathrm{ext}}|\Psi} \label{eq:18-levyA1}\\[2pt] &= \min_{\rho}\;\Bigl[\;\min_{\Psi\in A,\;\Psi\to\rho}\bigl(\braket{\Psi|\hat{T}+\hat{V}_{ee}|\Psi} + \braket{\Psi|\hat{V}_{\mathrm{ext}}|\Psi}\bigr)\Bigr] \label{eq:18-levyA2}\\[2pt] &= \min_{\rho}\;\Bigl[\;\min_{\Psi\in A,\;\Psi\to\rho}\braket{\Psi|\hat{T}+\hat{V}_{ee}|\Psi} + \int v_{\mathrm{ext}}(\rr)\rho(\rr)\dd^3r\Bigr] \label{eq:18-levyA3}\\[2pt] &\equiv \min_{\rho}\Bigl[F_A[\rho] + \int v_{\mathrm{ext}}(\rr)\rho(\rr)\dd^3r\Bigr] . \label{eq:18-levyA4} \end{align}

\eqref{eq:18-levyA2}では最小化を2段階に分けた:まず密度 $\rho$ を固定して，その密度を与えるクラス $A$ の波動関数の中で最小を取り，次に $\rho$ について最小化する．全体としての最小値は変わらない(第9章9.4節と同じ論法)．\eqref{eq:18-levyA3}では，$\hat{V}_{\mathrm{ext}}$ が1体演算子なのでその期待値が $\int v_{\mathrm{ext}}\rho$ と密度だけで書け，内側の最小化から外に出せることを使った．\eqref{eq:18-levyA4}で

$$ \begin{equation} F_A[\rho] \equiv \min_{\Psi\in A,\;\Psi\to\rho}\braket{\Psi|\hat{T}+\hat{V}_{ee}|\Psi} \label{eq:18-FA} \end{equation} $$

と定義した．∎

$F_A[\rho]$ は $v_{\mathrm{ext}}$ に依存しないという意味で普遍的である．したがってクラス $A$ の中でHohenberg–Kohn–Kohn–Sham の枠組み全体をそのまま組み立てられる．ただし決定的な注意がある．

注意:$F_A$ は「もとの $F$」ではない

$F_A[\rho]$ はクラス $A$ の選び方に依存する．すなわち，基底状態の普遍汎関数 $F[\rho]$ とは異なる汎関数であり，当然その交換相関部分も異なる．厳密には，コアホールを持つ系にはコアホール系専用の交換相関汎関数を使うべきである．実際の計算では，基底状態用に構成されたLDA/GGA汎関数をそのまま流用する．これは制御されていない近似であり，その妥当性は最終的にベンチマークによってしか判定できない．18.5.1節で見た $0.4$–$0.5$ eV という精度は，この流用がうまく働いていることの経験的な証拠である．逆に，$\mathrm{O_2}$ の $1s$ 交換分裂が半分しか出ないのは，この流用が破綻する場合の例である．

18.6.2 コアホール状態への適用

あとは $\hat{O}$ を具体的に選ぶだけである．第6章の第二量子化を使う．着目する内殻軌道の生成・消滅演算子を $\hat{a}_c^\dagger, \hat{a}_c$，数演算子を $\hat{n}_c = \hat{a}_c^\dagger\hat{a}_c$ とし，

$$ \begin{equation} \hat{A} = 1 - \hat{n}_c \label{eq:18-corehole-proj} \end{equation} $$

と取る．これがきちんと射影演算子になっていることを確かめよう．

導出:$\hat{A} = 1-\hat{n}_c$ は射影演算子である

ステップ1:$\hat{n}_c$ が冪等であることを示す．フェルミオンの反交換関係(第6章)$\{\hat{a}_c,\hat{a}_c^\dagger\} = \hat{a}_c\hat{a}_c^\dagger + \hat{a}_c^\dagger\hat{a}_c = 1$ と $\hat{a}_c^\dagger\hat{a}_c^\dagger = 0$ を使う:

\begin{align} \hat{n}_c^2 &= \hat{a}_c^\dagger\hat{a}_c\,\hat{a}_c^\dagger\hat{a}_c \label{eq:18-nc1}\\[2pt] &= \hat{a}_c^\dagger\bigl(1 - \hat{a}_c^\dagger\hat{a}_c\bigr)\hat{a}_c \label{eq:18-nc2}\\[2pt] &= \hat{a}_c^\dagger\hat{a}_c - \hat{a}_c^\dagger\hat{a}_c^\dagger\hat{a}_c\hat{a}_c \label{eq:18-nc3}\\[2pt] &= \hat{n}_c - 0 = \hat{n}_c . \label{eq:18-nc4} \end{align}

\eqref{eq:18-nc2}では反交換関係から $\hat{a}_c\hat{a}_c^\dagger = 1 - \hat{a}_c^\dagger\hat{a}_c$ を代入し，\eqref{eq:18-nc4}では $\hat{a}_c^\dagger\hat{a}_c^\dagger = 0$(同じ状態を2回作れない，Pauli原理)を使った．

ステップ2:$\hat{A}$ の冪等性とエルミート性．

$$ \hat{A}^2 = (1-\hat{n}_c)^2 = 1 - 2\hat{n}_c + \hat{n}_c^2 = 1 - 2\hat{n}_c + \hat{n}_c = 1 - \hat{n}_c = \hat{A}, \qquad \hat{A}^\dagger = 1 - \hat{n}_c^\dagger = \hat{A} . $$

よって $\hat{A}$ は射影演算子である．

ステップ3:何を射影しているか．$\hat{n}_c$ の固有値は $0$ か $1$ である．したがって

$$ \hat{A}\ket{\Psi} = \ket{\Psi}\ \Longleftrightarrow\ \hat{n}_c\ket{\Psi} = 0 \ \Longleftrightarrow\ \text{内殻軌道 } c \text{ が空である} $$

である．クラス $A$ は「狙った内殻軌道が空である状態全体」であり，その中の最低エネルギー状態がコアホール第一励起状態にほかならない．∎

これで論理が閉じた．18.4.2節で導入したペナルティ汎関数法は，単なる「計算のトリック」ではなく，Gunnarsson–Lundqvistの変分原理に裏打ちされた手続きである．$\hat{P}$ による罰則は，変分空間をクラス $A$ に制限する射影 $\hat{A}$ を実装している．

注意:$[\hat{H},\hat{n}_c] = 0$ は厳密ではない

上の議論は $\hat{O} = \hat{n}_c$ が $\hat{H}$ と可換であることを前提にした．しかし多体ハミルトニアンの電子間相互作用は配置を混ぜるので，厳密には $[\hat{H},\hat{n}_c] \ne 0$ である．したがって「内殻が空である状態」は $\hat{H}$ の固有状態の集合として閉じていない．この近似が許されるのは，内殻準位が価電子帯から数百 eV も離れており，混成の行列要素がエネルギー差に比べて桁違いに小さいからである．エネルギー分母が大きいので2次摂動論的な混ざり込みは無視できる．逆に，浅い内殻($3d$ 遷移金属の $3p$ など)や強相関系では，この分離が悪くなり，多重項構造・サテライトピークとして現れる．それらは本節の枠組みの外にある．

さらに高い励起状態(内殻電子がLUMOではなくLUMO$+k$ に入った状態)は，第一励起状態のエネルギー $E_1$ にKS固有値差を足して近似する:

$$ \begin{equation} \Delta E_k \simeq E_1 + \bigl(\varepsilon^{\sigma}_{\mathrm{LUMO}+k} - \varepsilon^{\sigma}_{\mathrm{LUMO}}\bigr) . \label{eq:18-higher-exc} \end{equation} $$

ここで $\sigma$ はコアホールと同じスピンチャネルである．光子は電子のスピンを反転させないので，遷移は $\Delta S = 0$ を満たさなければならない(18.6.4節)．つまり，上向きスピンの内殻電子を抜いたなら，行き先も上向きスピンの非占有状態である．式\eqref{eq:18-higher-exc}は，コアホールによって歪んだ終状態のKS固有値を使う点が要点で，基底状態の固有値を使ってはならない．

18.6.3 終状態則とZ$+1$近似

いま「終状態のKS固有値を使え」と述べた．この処方には終状態則(final-state rule)という名前が付いている．

物理的意味:終状態則

X線吸収スペクトルの形状は，コアホールが存在する状態(すなわち終状態)の電子構造で決まる．逆にX線発光(内殻空孔に価電子が落ちてX線を出す過程)のスペクトル形状は，その始状態——やはりコアホールがある状態——で決まる．いずれも「コアホールがある側」の電子構造が効く，というのが終状態則の内容である．

理由は直観的である．吸収の終状態では，励起された電子はまだ原子の近くにいて，コアホールが作る引力ポテンシャルを感じている．このポテンシャルは電子を原子側に引き寄せるので，非占有状態の局所部分状態密度が吸収原子の上に大きく積み上がり，鋭い共鳴(白線ピーク)を作る．基底状態(コアホールなし)の非占有状態を使って計算すると，この引力が抜け落ち，端近傍のピーク強度が大幅に過小評価され，ピーク位置も数 eV ずれる．「バンド計算の空バンドDOSをそのままXANESと比べる」という素朴な処方が失敗する理由がこれである．

コアホールの効果を手早く見積もる古典的な近似がZ$+1$近似(等価内殻近似)である．

例:Z$+1$近似の考え方

原子番号 $Z$ の原子の $1s$ 軌道から電子を1個抜くと，価電子が感じる核電荷の遮蔽が1単位分だけ弱まる．価電子から見れば，これは「核電荷が $Z$ から $Z+1$ に増えた」のとほぼ同じである．したがって，$1s$ ホールを持つ炭素は窒素のように振る舞う．同様に，$1s$ ホールを持つ窒素は酸素のように，$2p$ ホールを持つケイ素はリンのように振る舞う．

この近似の使い道は2つある．第一に，コアホール計算をせずに終状態の性質を定性的に予想できる．たとえば $\mathrm{CO}$ の C $1s$ 励起の終状態は $\mathrm{NO}$ に似ているはずだ，という具合である．第二に，化学シフトを熱力学サイクルで見積もれる．内殻束縛エネルギーの差は，$Z$ の化合物と $Z+1$ の化合物の生成エネルギーの差に関係づけられる．

もちろん定量的には粗い近似で，$0.1$ eV 精度が要求される現代の解析には使えない．しかし「コアホールとは核電荷が1つ増えることだ」という直観は，18.5.2節の電荷ポテンシャル模型と併せて，スペクトルを読む際の羅針盤になる．

18.6.4 フェルミの黄金律と双極子選択則

ここからは吸収強度の話である．まず，時間に依存する摂動に対する遷移確率の一般公式(フェルミの黄金律)を用意する．これは18.8節の誘電関数でもそのまま使うので，丁寧に導いておく．

数学ノート:フェルミの黄金律の導出

時間に依存しないハミルトニアン $\hat{H}_0$($\hat{H}_0\ket{n} = E_n\ket{n}$)に，$t=0$ から振動する摂動

$$ \hat{H}'(t) = \hat{V}e^{-i\omega t} + \hat{V}^\dagger e^{i\omega t} $$

を加える(実の摂動をこの形に書けることは，$\cos\omega t$ を指数関数に分解すれば分かる)．

ステップ1:展開係数の方程式．$\ket{\psi(t)} = \sum_n c_n(t)e^{-iE_nt}\ket{n}$ を時間依存Schrödinger方程式 $i\partial_t\ket{\psi} = (\hat{H}_0+\hat{H}')\ket{\psi}$ に代入し，左から $\bra{m}e^{iE_mt}$ を掛けると

$$ i\dot{c}_m(t) = \sum_n \braket{m|\hat{H}'(t)|n}\,e^{i(E_m-E_n)t}\,c_n(t) $$

を得る($\hat{H}_0$ の項は指数因子の微分と相殺する)．これはまだ厳密である．

ステップ2:1次摂動．初期条件 $c_i(0)=1$，他はゼロ．右辺の $c_n$ を初期値で置き換える(1次のBorn近似):

$$ i\dot{c}_f(t) = \braket{f|\hat{V}|i}e^{i(E_f-E_i-\omega)t} + \braket{f|\hat{V}^\dagger|i}e^{i(E_f-E_i+\omega)t} . $$

第1項は吸収($E_f \simeq E_i+\omega$)，第2項は誘導放出($E_f\simeq E_i-\omega$)に対応する．吸収を考えるので第1項だけ残す．

ステップ3:時間積分．$\Omega \equiv E_f-E_i-\omega$ と置くと

$$ c_f(t) = -i\braket{f|\hat{V}|i}\int_0^t e^{i\Omega t'}\dd t' = -i\braket{f|\hat{V}|i}\,\frac{e^{i\Omega t}-1}{i\Omega} . $$

ステップ4:確率と，$t\to\infty$ の極限．

$$ \abs{c_f(t)}^2 = \abs{\braket{f|\hat{V}|i}}^2\frac{\abs{e^{i\Omega t}-1}^2}{\Omega^2} = \abs{\braket{f|\hat{V}|i}}^2\frac{4\sin^2(\Omega t/2)}{\Omega^2} . $$

($\abs{e^{i\theta}-1}^2 = (\cos\theta-1)^2+\sin^2\theta = 2-2\cos\theta = 4\sin^2(\theta/2)$ を使った．)最後の因子は $\Omega=0$ に幅 $\sim 1/t$ の鋭いピークを持ち，その面積は

$$ \int_{-\infty}^{\infty}\frac{4\sin^2(\Omega t/2)}{\Omega^2}\dd\Omega = 4\cdot\frac{t}{2}\int_{-\infty}^{\infty}\frac{\sin^2 x}{x^2}\dd x = 2t\cdot\pi = 2\pi t $$

である($x = \Omega t/2$ と置換，$\dd\Omega = 2\dd x/t$，そして標準積分 $\int_{-\infty}^\infty(\sin x/x)^2\dd x = \pi$ を使った)．したがって $t\to\infty$ で

$$ \frac{4\sin^2(\Omega t/2)}{\Omega^2} \;\longrightarrow\; 2\pi t\,\delta(\Omega) $$

となり，単位時間あたりの遷移確率

$$ \begin{equation} w_{i\to f} = \frac{\abs{c_f(t)}^2}{t} = 2\pi\abs{\braket{f|\hat{V}|i}}^2\delta(E_f-E_i-\omega) \label{eq:18-golden} \end{equation} $$

が得られる．これがフェルミの黄金律である．$\delta$ 関数はエネルギー保存を表す．

次に，X線に対する摂動 $\hat{V}$ を書き下す．電磁場中の電子のハミルトニアンは，最小結合 $\bm{p}\to\bm{p}+\bm{A}/c$(電子の電荷が $-1$ であることに注意)により

$$ \hat{H} = \frac{1}{2}\left(\hat{\bm{p}}+\frac{\bm{A}}{c}\right)^2 + V(\rr) = \hat{H}_0 + \frac{1}{2c}\bigl(\hat{\bm{p}}\cdot\bm{A}+\bm{A}\cdot\hat{\bm{p}}\bigr) + \frac{A^2}{2c^2} $$

である．Coulombゲージ $\nabla\cdot\bm{A}=0$ では $[\hat{\bm{p}},\bm{A}] = -i\nabla\cdot\bm{A} = 0$ なので $\hat{\bm{p}}\cdot\bm{A}=\bm{A}\cdot\hat{\bm{p}}$，また $A^2$ の項は光子2個の過程(散乱)に対応し1光子吸収には寄与しない．よって

$$ \begin{equation} \hat{H}' = \frac{1}{c}\bm{A}\cdot\hat{\bm{p}} . \label{eq:18-Ap} \end{equation} $$

平面波 $\bm{A}(\rr,t) = A_0\hat{\bm{e}}\cos(\bm{q}\cdot\rr-\omega t)$ を代入し，$\cos$ を指数関数に分解すると，吸収に対応する部分は

$$ \hat{V} = \frac{A_0}{2c}\,e^{i\bm{q}\cdot\rr}\,\hat{\bm{e}}\cdot\hat{\bm{p}} $$

である．

導出:双極子近似の妥当性と長さ形式への変換

ステップ1:$e^{i\bm{q}\cdot\rr}\simeq1$ とおけるか．行列要素 $\braket{f|e^{i\bm{q}\cdot\rr}\hat{\bm{e}}\cdot\hat{\bm{p}}|c}$ で，$\ket{c}$ は内殻軌道なので，積分に効くのは吸収原子のごく近傍だけである．そこでの $\abs{\bm{q}\cdot\rr}$ の大きさを見積もる．炭素K端 $\omega = 285$ eV では波長 $\lambda = 12398/285 = 43.5$ Å，波数 $q = 2\pi/\lambda = 0.144$ Å$^{-1}$ である．炭素 $1s$ 軌道の半径はおよそ $a_0/Z_{\mathrm{eff}} \simeq 0.53/5.7 = 0.09$ Å だから

$$ \abs{\bm{q}\cdot\rr}\sim 0.144\times0.09 = 0.013 \ll 1 . $$

したがって $e^{i\bm{q}\cdot\rr} = 1 + i\bm{q}\cdot\rr + \cdots \simeq 1$ としてよい．これが双極子近似である($10$ keV 級の硬X線でも，重元素の内殻はさらに小さいので概ね成り立つ．破れると四重極遷移が現れ，$3d$ 遷移金属のK端前縁ピークとして観測される)．

ステップ2:運動量形式から長さ形式へ．$\hat{H}_0 = \hat{p}^2/2+V(\rr)$ に対して，正準交換関係 $[\hat{p}_\alpha,\hat{r}_\beta] = -i\delta_{\alpha\beta}$ を使うと

\begin{align} [\hat{H}_0,\hat{\bm{r}}] &= \left[\frac{\hat{p}^2}{2},\hat{\bm{r}}\right] = \frac{1}{2}\Bigl(\hat{\bm{p}}[\hat{\bm{p}},\hat{\bm{r}}] + [\hat{\bm{p}},\hat{\bm{r}}]\hat{\bm{p}}\Bigr) \label{eq:18-comm1}\\[2pt] &= \frac{1}{2}\bigl(-i\hat{\bm{p}} - i\hat{\bm{p}}\bigr) = -i\hat{\bm{p}} . \label{eq:18-comm2} \end{align}

\eqref{eq:18-comm1}では $[\hat{A}\hat{B},\hat{C}] = \hat{A}[\hat{B},\hat{C}]+[\hat{A},\hat{C}]\hat{B}$ を使い，$V(\rr)$ は $\hat{\bm{r}}$ と可換なので落ちる．したがって $\hat{\bm{p}} = i[\hat{H}_0,\hat{\bm{r}}]$ であり，固有状態間の行列要素は

\begin{align} \braket{f|\hat{\bm{p}}|c} &= i\braket{f|\hat{H}_0\hat{\bm{r}}-\hat{\bm{r}}\hat{H}_0|c} \label{eq:18-len1}\\[2pt] &= i(E_f-E_c)\braket{f|\hat{\bm{r}}|c} = i\omega\braket{f|\hat{\bm{r}}|c} \label{eq:18-len2} \end{align}

となる(\eqref{eq:18-len1}で左右の固有状態に $\hat{H}_0$ を作用させ，\eqref{eq:18-len2}でエネルギー保存 $E_f-E_c=\omega$ を使った)．よって $\abs{\braket{f|\hat{\bm{e}}\cdot\hat{\bm{p}}|c}}^2 = \omega^2\abs{\braket{f|\hat{\bm{e}}\cdot\hat{\bm{r}}|c}}^2$ である．∎

これらを黄金律\eqref{eq:18-golden}に入れ，吸収されるエネルギーを入射エネルギー流束で割ると吸収断面積が得られる．入射平面波の時間平均Poyntingベクトルは $S = c\abs{E_0}^2/(8\pi) = \omega^2A_0^2/(8\pi c)$($\abs{E_0} = \omega A_0/c$)なので，

\begin{align} \sigma(\omega) &= \frac{\omega\,\sum_f w_{c\to f}}{S} = \frac{\omega\cdot 2\pi\dfrac{A_0^2}{4c^2}\sum_f\abs{\braket{f|\hat{\bm{e}}\cdot\hat{\bm{p}}|c}}^2\delta(\varepsilon_f-\varepsilon_c-\omega)} {\omega^2A_0^2/(8\pi c)} \label{eq:18-xsec-a}\\[2pt] &= \frac{4\pi^2}{c\,\omega}\sum_f\abs{\braket{f|\hat{\bm{e}}\cdot\hat{\bm{p}}|c}}^2\delta(\varepsilon_f-\varepsilon_c-\omega) \label{eq:18-xsec-b}\\[2pt] &= 4\pi^2\alpha\,\omega\sum_f\abs{\braket{f|\hat{\bm{e}}\cdot\hat{\bm{r}}|c}}^2\delta(\varepsilon_f-\varepsilon_c-\omega) . \label{eq:18-xas-sigma} \end{align}

\eqref{eq:18-xsec-b}では $A_0^2$ と $\omega$ を整理し，\eqref{eq:18-xas-sigma}では長さ形式への変換 $\abs{M_p}^2=\omega^2\abs{M_r}^2$ を使い，原子単位系での微細構造定数 $\alpha = 1/c = 1/137.036$ を導入した．

$$ \begin{equation} \sigma(\omega) = 4\pi^2\alpha\,\omega\sum_f\abs{\braket{f|\hat{\bm{e}}\cdot\hat{\bm{r}}|c}}^2\,\delta(\varepsilon_f-\varepsilon_c-\omega) \label{eq:18-xas-key} \end{equation} $$

次に，この行列要素がどの終状態に対してゼロでないかを調べる．それが選択則である．

導出:双極子選択則 $\Delta l=\pm1$，$\Delta m = 0,\pm1$，$\Delta s=0$

ステップ1:双極子演算子を球面調和関数で書く．球座標で $z = r\cos\theta$，$x\pm iy = r\sin\theta\,e^{\pm i\varphi}$ である．$Y_1^0 = \sqrt{3/4\pi}\cos\theta$，$Y_1^{\pm1} = \mp\sqrt{3/8\pi}\sin\theta\,e^{\pm i\varphi}$ を使えば

$$ z = r\sqrt{\frac{4\pi}{3}}\,Y_1^0, \qquad x\pm iy = \mp r\sqrt{\frac{8\pi}{3}}\,Y_1^{\pm1} $$

と書ける．したがって $\hat{\bm{e}}\cdot\rr$ は，$r$ と $Y_1^{\mu}$($\mu=0,\pm1$)の線形結合の積である．

ステップ2:行列要素を動径部分と角度部分に分ける．内殻状態を $\ket{c} = R_{n_cl_c}(r)Y_{l_c}^{m_c}$，終状態の吸収原子まわりの成分を $R_l(r;\varepsilon)Y_l^m$ とすると

$$ \braket{f|\hat{\bm{e}}\cdot\rr|c} = \underbrace{\int_0^\infty R_l^*(r)\,r\,R_{n_cl_c}(r)\,r^2\dd r}_{\text{動径積分}} \times \underbrace{C_\mu\int Y_l^{m*}\,Y_1^{\mu}\,Y_{l_c}^{m_c}\,\dd\Omega}_{\text{角度積分(Gaunt係数)}} $$

となる．選択則は角度積分から出る．

ステップ3:$\varphi$ 積分から磁気量子数の規則．$Y_l^m \propto e^{im\varphi}$ だから，$\varphi$ に関する積分は

$$ \int_0^{2\pi} e^{-im\varphi}e^{i\mu\varphi}e^{im_c\varphi}\dd\varphi = 2\pi\,\delta_{m,\,m_c+\mu} $$

である．$\mu = 0,\pm1$ なので $\Delta m \equiv m-m_c = 0,\pm1$．

ステップ4:パリティから $l$ の偶奇．空間反転 $\rr\to-\rr$($\theta\to\pi-\theta$, $\varphi\to\varphi+\pi$)のもとで $Y_l^m(-\hat{\bm{r}}) = (-1)^lY_l^m(\hat{\bm{r}})$ である(証明:$Y_l^m\propto P_l^m(\cos\theta)e^{im\varphi}$ で，$P_l^m(-x)=(-1)^{l+m}P_l^m(x)$，$e^{im(\varphi+\pi)}=(-1)^me^{im\varphi}$ だから，積は $(-1)^{l+2m}=(-1)^l$ 倍される)．すると角度積分の被積分関数は反転で $(-1)^{l+1+l_c}$ 倍される．積分領域(球面)は反転で不変だから，積分値 $I$ は

$$ I = (-1)^{l+1+l_c}I $$

を満たさなければならない．$(-1)^{l+1+l_c} = -1$ なら $I=0$ である．ゆえに $l+l_c+1$ が偶数，すなわち $l-l_c$ は奇数でなければならない．

ステップ5:三角不等式から $l$ の範囲．角運動量の合成則により，積 $Y_1^\mu Y_{l_c}^{m_c}$ は $\abs{l_c-1}\le L\le l_c+1$ を満たす $Y_L^M$ の線形結合に展開できる．球面調和関数の直交性から，$\int Y_l^{m*}Y_1^\mu Y_{l_c}^{m_c}\dd\Omega$ がゼロでないためには $l$ がこの範囲になければならない:

$$ \abs{l - l_c} \le 1 . $$

ステップ6:両者を組み合わせる．ステップ4より $l-l_c$ は奇数，ステップ5より $\abs{l-l_c}\le1$．両立する値は $l-l_c = \pm1$ のみである($l-l_c=0$ は偶数なので排除される)．

ステップ7:スピン．双極子演算子 $\hat{\bm{e}}\cdot\rr$ はスピン空間に作用しないので，スピン部分の内積は $\braket{\sigma_f|\sigma_c} = \delta_{\sigma_f\sigma_c}$ を与える．すなわち $\Delta s = 0$．∎

$$ \begin{equation} \Delta l = \pm1,\qquad \Delta m = 0,\pm1,\qquad \Delta s = 0 \qquad(\text{双極子選択則}) \label{eq:18-selection} \end{equation} $$

この規則がXANESの解釈を決定する．吸収端の名前と，そこで見える終状態の対称性の対応は次のとおりである．

吸収端	始状態	$l_c$	許される終状態	実際に主役となるもの
K	$1s$	0	$l=1$	$p$ 成分の非占有局所DOS
L$_1$	$2s$	0	$l=1$	$p$ 成分
L$_2$, L$_3$	$2p_{1/2}$, $2p_{3/2}$	1	$l=0$ または $l=2$	$d$ 成分(動径積分が $s$ より大きい)
M$_4$, M$_5$	$3d_{3/2}$, $3d_{5/2}$	2	$l=1$ または $l=3$	$f$ 成分(希土類の解析に使われる)

18.6.5 XANESは局所部分状態密度を映す

選択則が分かったので，断面積\eqref{eq:18-xas-key}をもう一段書き直そう．内殻軌道 $\ket{c}$ は吸収原子 $A$ のまわり半径 $0.1$ Å 程度に閉じ込められている．したがって行列要素は，終状態 $\ket{f}$ の原子 $A$ の近傍における振る舞いしか見ていない．そこで終状態を原子 $A$ を中心とする球面調和関数で展開する:

$$ \ket{f} \;\simeq\; \sum_{lm} a^{f}_{lm}\,R_l(r;\varepsilon_f)\,Y_l^m(\hat{\bm{r}}) \qquad (\rr \text{ が } A \text{ の近傍のとき}) . $$

K端($l_c=0$)の場合，選択則より $l=1$ の成分だけが効く．動径積分

$$ M_{\mathrm{rad}}(\varepsilon) = \int_0^\infty R_1^*(r;\varepsilon)\,r\,R_{1s}(r)\,r^2\dd r $$

はエネルギーとともにゆっくりしか変化しない(内殻軌道の広がりの中でしか積分が効かず，その領域では終状態の動径関数の形はエネルギーに鈍感である)．これを積分の外に出すと

\begin{align} \sigma(\omega) &\simeq 4\pi^2\alpha\,\omega\,\abs{M_{\mathrm{rad}}(\varepsilon_c+\omega)}^2 \sum_f\sum_{m}\abs{a^f_{1m}}^2\,\delta(\varepsilon_f-\varepsilon_c-\omega) \label{eq:18-pdos1}\\[2pt] &= 4\pi^2\alpha\,\omega\,\abs{M_{\mathrm{rad}}}^2\;\rho^{A}_{l=1}(\varepsilon_c+\omega) , \label{eq:18-xas-pdos} \end{align}

となる．ここで最後の等号で，原子 $A$ 上の $l$ 射影局所状態密度

$$ \begin{equation} \rho^A_{l}(\varepsilon) \equiv \sum_f\sum_m \abs{a^f_{lm}}^2\,\delta(\varepsilon-\varepsilon_f) \label{eq:18-pdos-def} \end{equation} $$

の定義(第13章の部分状態密度)を使った．

$$ \begin{equation} \sigma_{\mathrm{K}}(\omega)\;\propto\;\omega\,\abs{M_{\mathrm{rad}}}^2\,\rho^{A}_{p}\bigl(\varepsilon_c+\omega\bigr) \qquad(\text{コアホールが存在する系のDOSを使う}) \label{eq:18-xas-pdos-key} \end{equation} $$

物理的意味:XANESが「局所プローブ」である3つの理由

元素選択的:吸収端のエネルギーが元素ごとに大きく違うので，混合物でも狙った元素だけを見られる．
サイト選択的:同じ元素でも化学環境が違えば内殻準位が化学シフトするので(18.5節)，異なるサイトの寄与を分離できる．
対称性選択的:選択則\eqref{eq:18-selection}により，特定の角運動量成分だけが見える．さらに偏光を使えば $m$ 成分まで選べる．

この3つが揃うので，XANESは「非晶質・液体・埋もれた界面」といった，回折が使えない系の局所構造解析に威力を発揮する．

例:炭素K端の π* と σ*，そして偏光依存性

不飽和炭化水素の炭素K端スペクトルには，吸収端の手前に鋭いピークが現れる．これは $1s\to\pi^*$ 遷移である．$\pi^*$ 軌道は $\mathrm{C}$–$\mathrm{C}$ 結合軸に垂直な $p$ 成分($p_z$)からできており，選択則 $\Delta l=\pm1$ を満たす．一方，飽和炭化水素には $\pi^*$ がなく，$\sigma^*$ 遷移だけが数 eV 高いところに現れる．実際，代表的な分子の第一ピーク位置は次のようになる．

分子	第一ピークの性格	計算値 (eV)	実験値 (eV)
エチレン $\mathrm{C_2H_4}$	$1s\to\pi^*$(C=C)	$284.17$	$284.3$
アセチレン $\mathrm{C_2H_2}$	$1s\to\pi^*$(C≡C)	$285.56$	$285.8$
エタン $\mathrm{C_2H_6}$	$1s\to\sigma^*$(C–H/C–C)	$287.00$	$287.2$
c-BN(B K端)	$1s\to\pi^*$	$195.28$	$195.5$

いずれも絶対値で $0.2$–$0.3$ eV の一致である．飽和分子(エタン)の第一ピークが不飽和分子より高エネルギー側にあるのは，$\pi^*$ 軌道が $\sigma^*$ 軌道より低いという結合論の常識(第2章)がそのままスペクトルに現れたものである．

偏光依存性．放射光の直線偏光を使い，試料を傾けながら測ると，$\hat{\bm{e}}$ の向きによって $\pi^*$ ピークと $\sigma^*$ ピークの強度比が劇的に変わる．$\abs{\braket{f|\hat{\bm{e}}\cdot\rr|c}}^2$ は $\hat{\bm{e}}$ が終状態の $p$ 軌道の向きと平行なときに最大になるからである．金属表面に平らに寝た芳香環分子なら，$\hat{\bm{e}}$ が表面法線を向くとき $\pi^*$ が最大，面内を向くとき $\sigma^*$ が最大になる．この角度依存性から，吸着分子の傾き角が $\pm5^\circ$ 程度の精度で決定できる．NEXAFSが表面科学の標準ツールである理由である．

c-BNのボロンK端は，コアホール計算のセルサイズ依存性を見るのに好例である．8原子の単位胞では実験のスペクトル形状を再現できないが，64原子胞まで大きくすると実験とよく一致する．18.4.4節で見た「コアホールの遮蔽電荷が広がる距離」が，ここでもセルサイズ収束を支配している．多重散乱法(FEFFなど)で同じ計算を行う場合も，クラスターサイズを5原子から239原子まで系統的に増やす必要があり，事情は同じである．

補足:振動子強度と多電子効果

式\eqref{eq:18-xas-key}は一電子的な行列要素で書かれているが，原理的には多電子の遷移双極子行列要素

$$ f_{G,E} \;\propto\; \frac{1}{\omega_{GE}}\abs{\braket{\Phi_G|\hat{\bm{e}}\cdot\hat{\bm{D}}|\Phi_E}}^2, \qquad \hat{\bm{D}} = \sum_{i}\rr_i $$

を計算すべきである($\Phi_G$ は基底状態，$\Phi_E$ は励起状態の多体波動関数)．ここで重要なのは，自己無撞着に得られた励起状態のSlater行列式 $\Phi_E$ を基底状態の軌道で展開すると

$$ \ket{\Phi_E} = c_G\ket{\Phi_G} + \sum_{\mu\le N}\sum_{\nu\gt N}c_{\mu\nu}\,\hat{a}_\nu^\dagger\hat{a}_\mu\ket{\Phi_G} + (\text{2電子以上の励起配置}) $$

という形になることである．すなわち，コアホール状態を自己無撞着に解くという行為は，基底状態から見ると無限個の多電子励起配置の重ね合わせを暗黙に含んでいる．これがコアホールの高次遮蔽過程(shake-up/shake-off サテライト)を部分的に取り込む仕組みであり，$\Delta$SCF が単なる一電子近似より良い結果を与える理由の一つである．

18.7 ARPESとバンドアンフォールディング

図18.1の3本目の矢印に移る．今度は始状態が価電子帯である．価電子帯から真空へ電子を叩き出す実験が紫外光電子分光(UPS)であり，飛び出す方向(角度)まで分解して測るのが角度分解光電子分光(ARPES)である．ARPESは，バンド構造 $\varepsilon_n(\kk)$ を直接「写真に撮る」ことのできるほとんど唯一の実験手法であり，第一原理バンド計算にとって最も直接的な検証手段である．

本節では，(i) 何を測ると $\kk$ が決まるのかという運動学，(ii) スーパーセル計算をすると計算結果が実験と直接比較できなくなるという問題，(iii) その解決策であるバンドアンフォールディングの定式化，を順に扱う．

18.7.1 ARPESの運動学

ARPESの測定量は2つ，光電子の運動エネルギー $E_{\mathrm{kin}}$ と脱出方向の極角 $\theta$(表面法線から測る)である．ここから束縛エネルギーと波数を決める．

図18.8 ARPES測定のジオメトリ．試料表面に光を当て，極角 $\theta$ の方向に飛び出す光電子の運動エネルギーを分析器で測る．表面平行方向には並進対称性が保たれているので運動量の面内成分が保存し，$E_{\mathrm{kin}}$ と $\theta$ から結晶運動量の面内成分 $k_\parallel$ が決まる．面直成分は表面で保存しないため，別の仮定が必要になる．

エネルギー．18.2.1節と同じエネルギー保存則から，束縛エネルギーは

$$ \begin{equation} E_B = h\nu - E_{\mathrm{kin}} - \phi \label{eq:18-arpes-EB} \end{equation} $$

で与えられる．価電子帯の場合，$E_B$ はフェルミ準位から測って $0$ から $10$ eV 程度の量である．

面内波数．結晶表面は，表面に平行な方向には周期性を保っている(表面再構成があってもその周期で)．したがって表面平行方向の運動量は，逆格子ベクトル $\bm{G}_\parallel$ を除いて保存する．真空中の光電子の運動量の大きさは $p = \sqrt{2mE_{\mathrm{kin}}}$，その面内成分は $p\sin\theta$ だから

$$ \begin{equation} \hbar k_\parallel = \sqrt{2mE_{\mathrm{kin}}}\,\sin\theta \qquad\Longleftrightarrow\qquad k_\parallel\,[\text{Å}^{-1}] = 0.5123\sqrt{E_{\mathrm{kin}}\,[\mathrm{eV}]}\,\sin\theta \label{eq:18-arpes-kpar} \end{equation} $$

である(数値係数は $\sqrt{2m_e\cdot1\,\mathrm{eV}}/\hbar$ をÅ$^{-1}$ 単位で書いたもの)．この式が，ARPESが「バンドを直接測る」実験である理由のすべてである．$\theta$ を掃引しながらエネルギー分布曲線を取れば，$(k_\parallel, E_B)$ 平面上の強度分布——すなわちバンド分散——が直接得られる．

例:測れる波数の範囲

典型的なブリルアン域の大きさは $\pi/a \simeq 1$ Å$^{-1}$($a=3$ Å)である．式\eqref{eq:18-arpes-kpar}で $\theta = 90^\circ$(表面すれすれ)としても $k_\parallel^{\max} = 0.5123\sqrt{E_{\mathrm{kin}}}$ なので，第一ブリルアン域の端($k_\parallel \simeq 1$ Å$^{-1}$)に到達するには $E_{\mathrm{kin}} \gtrsim 4$ eV が必要である．実際には $\theta \lesssim 60^\circ$ しか使えないので，$E_{\mathrm{kin}} \gtrsim 5$ eV，すなわち $h\nu \gtrsim 15$ eV 程度が必要になる．He I 共鳴線($21.2$ eV)がARPESの標準光源である理由がこれである．逆に，光子エネルギーを上げすぎると図18.3のユニバーサル曲線により $\lambda$ が最小になり，バルクの情報が取れなくなる．$20$–$100$ eV という光子エネルギー領域は，この2つの要求のバランスで決まっている．

面直波数の不定性．問題は $k_\perp$ である．表面は面直方向の並進対称性を破っているので，$k_\perp$ は保存しない．それでも $k_\perp$ を知りたい場合には，結晶内部での終状態を自由電子で近似する(自由電子終状態近似)．

導出:内部ポテンシャルによる $k_\perp$ の決定

ステップ1:結晶内部の終状態を仮定する．光電子は結晶内部では十分高いエネルギーにあるので，その分散を自由電子的とする．ただしエネルギーの原点は真空準位より $V_0$ だけ下にある(結晶の「内部ポテンシャル」，伝導帯の底の深さに対応する):

$$ E = \frac{\hbar^2 k^2}{2m} - V_0 = \frac{\hbar^2(k_\parallel^2 + k_\perp^2)}{2m} - V_0 . $$

ステップ2:真空中でのエネルギーを書く．真空中に出た電子については

$$ E_{\mathrm{kin}} = \frac{\hbar^2(k_\parallel^2 + k_{\perp,\mathrm{out}}^2)}{2m} . $$

ステップ3:保存量を突き合わせる．エネルギー $E = E_{\mathrm{kin}}$ と面内波数 $k_\parallel$ の両方が保存するので，上の2式を等しいと置くと $k_\parallel$ の項が消えて

$$ \frac{\hbar^2k_\perp^2}{2m} - V_0 = \frac{\hbar^2 k_{\perp,\mathrm{out}}^2}{2m} \qquad\Longrightarrow\qquad \frac{\hbar^2k_\perp^2}{2m} = \frac{\hbar^2 k_{\perp,\mathrm{out}}^2}{2m} + V_0 . $$

ステップ4:測定量で書く．$\hbar^2k_{\perp,\mathrm{out}}^2/2m = E_{\mathrm{kin}} - \hbar^2k_\parallel^2/2m = E_{\mathrm{kin}}(1-\sin^2\theta) = E_{\mathrm{kin}}\cos^2\theta$ だから

$$ \begin{equation} k_\perp = \frac{1}{\hbar}\sqrt{2m\bigl(E_{\mathrm{kin}}\cos^2\theta + V_0\bigr)} . \label{eq:18-arpes-kperp} \end{equation} $$

∎

ここで $V_0$ は実験から決めるパラメータである．光子エネルギーを掃引して同じ $k_\parallel$ でのバンド位置の周期性を見ると，$k_\perp$ がブリルアン域を1周する周期が読み取れ，そこから $V_0$($5$–$20$ eV 程度)がフィットされる．この点はARPESの解釈における不可避の曖昧さであり，二次元物質(面直方向に分散のない系)がARPES研究で好まれる理由の一つでもある．

注意:ARPESの強度はスペクトル関数そのものではない

ARPESの測定強度は，模式的に

$$ I(\kk,\omega) \;\propto\; \abs{M_{fi}(\kk,\omega)}^2\,A(\kk,\omega)\,f(\omega) $$

と書ける．$A(\kk,\omega)$ が一電子除去のスペクトル関数(第13章)，$f(\omega)$ がFermi分布(占有状態しか見えない)，$M_{fi}$ が遷移行列要素である．$M_{fi}$ は光子エネルギー・偏光・実験ジオメトリに依存し，あるバンドが「消える」ことすらある(行列要素効果)．したがって「ピークが見えない=バンドがない」ではない．逆に，偏光を変えて強度変化を追うことで，バンドの軌道対称性を同定できる(18.6.4節の偏光依存性と同じ原理である)．第一原理計算と比較する際は，まず $A(\kk,\omega)$ の比較から始め，必要に応じて $\abs{M_{fi}}^2$ を計算するのが実務的な順序である．

18.7.2 スーパーセル計算のジレンマ:バンドの折り畳み

ここからが本節の主題である．現実の試料は完全結晶ではない．欠陥がある，不純物が入っている，合金である，表面に分子が吸着している，長周期の再構成が起きている——こうした系を第一原理で扱うにはスーパーセル(基本単位胞を $N$ 倍に拡張したセル)を使うしかない．ところがスーパーセルを使うと，バンド構造が実験と直接比較できなくなる．

理由は単純である．実空間のセルを $N$ 倍にすると，逆空間のブリルアン域は $1/N$ に縮む．もとの基本セルのブリルアン域に入っていたバンドは，縮んだブリルアン域に折り畳まれて詰め込まれる．1本だったバンドが $N$ 本に増え，分散は寸断される．図18.9(a),(b)にこの様子を示した．

バンドの折り畳みとアンフォールディング．(a) 基本セルの cos 型バンド，(b) 2倍セルで上下 2 本に折り畳まれたバンド，(c) 基本セルのバンドと重みの大きい濃い丸，ゴーストバンドを示す薄い赤の小丸． — 図18.9 一次元強束縛鎖を例にしたバンドの折り畳みとアンフォールディング．(a) 基本セルの正しいバンド分散．(b) セルを2倍にすると，ブリルアン域が半分になり，1本のバンドが2本に折り畳まれる．この図をARPESと比べても一致しない．(c) スペクトル重み $W_{\bm{K}J}(\kk)$ を計算して基本セルのブリルアン域に戻すと，もとの分散が復元される．薄い赤の小円は「重みがほぼゼロのゴーストバンド」で，完全結晶では厳密にゼロ，摂動が入ると有限の重みを持つようになる．

ここで重要な認識がある．スーパーセルが単に基本セルを繰り返しただけ(摂動なし)なら，折り畳まれたバンドの情報量は基本セルのそれと完全に同じである．失われたものは何もなく，単に「表示のしかた」が変わっただけである．したがって，適切な変換を施せば元に戻せるはずである．そして摂動がある場合には，その変換が摂動の効果を定量的に可視化する道具になる．これがアンフォールディングである．

18.7.3 アンフォールディングの定式化

設定を明確にする．参照する基本セルの格子ベクトルを $\bm{a}_1,\bm{a}_2,\bm{a}_3$，スーパーセルの格子ベクトルを $\bm{A}_i = N_i\bm{a}_i$ とし，$N = N_1N_2N_3$ とおく．スーパーセル1個の中には $N$ 個の基本並進

$$ \bm{R}_{\bm{n}} = n_1\bm{a}_1+n_2\bm{a}_2+n_3\bm{a}_3, \qquad 0\le n_i \lt N_i $$

が含まれる．基本セルの逆格子ベクトルを $\bm{G}_p$，スーパーセルのそれを $\bm{G}_s$ と書く．$\bm{A}_i = N_i\bm{a}_i$ より $\bm{B}_i = \bm{b}_i/N_i$ だから，スーパーセルの逆格子は基本セルの逆格子より密であり，$\{\bm{G}_p\}\subset\{\bm{G}_s\}$ という包含関係が成り立つ．

スーパーセルのKohn–Sham固有状態を，スーパーセル波数 $\bm{K}$(縮んだブリルアン域内)とバンド指標 $J$ でラベルし，平面波で展開する:

$$ \begin{equation} \ket{\Psi_{\bm{K}J}} = \sum_{\bm{G}_s} C_{\bm{K}J}(\bm{G}_s)\,\ket{\bm{K}+\bm{G}_s}, \qquad \ket{\bm{q}} \equiv \frac{1}{\sqrt{\Omega}}e^{i\bm{q}\cdot\rr}, \qquad \sum_{\bm{G}_s}\abs{C_{\bm{K}J}(\bm{G}_s)}^2 = 1 . \label{eq:18-sc-expand} \end{equation} $$

問いはこうである．この状態は，基本セルのBloch波数 $\kk$ の性格をどれだけ持っているか．答えを与えるには，「基本セルのBloch波数 $\kk$ を持つ」ということを演算子で表現しなければならない．それが並進群の既約表現への射影である．

定義:基本並進群の既約表現への射影演算子

並進演算子を $(\hat{T}_{\bm{R}}\psi)(\rr) = \psi(\rr-\bm{R})$ と定める．基本セル波数 $\kk$ に属する成分を取り出す演算子を

$$ \begin{equation} \hat{P}_{\kk} = \frac{1}{N}\sum_{\bm{n}} e^{i\kk\cdot\bm{R}_{\bm{n}}}\,\hat{T}_{\bm{R}_{\bm{n}}} \label{eq:18-projector-k} \end{equation} $$

と定義する(和はスーパーセル内の $N$ 個の基本並進すべてにわたる)．これは有限巡回群 $\{\hat{T}_{\bm{R}_{\bm{n}}}\}$ の，指標 $e^{-i\kk\cdot\bm{R}}$ を持つ既約表現への標準的な射影演算子である．

導出:$\hat{P}_{\kk}$ の作用と射影演算子であることの確認

ステップ1:平面波への作用．平面波は並進演算子の固有関数である:

$$ \hat{T}_{\bm{R}}\,e^{i\bm{q}\cdot\rr} = e^{i\bm{q}\cdot(\rr-\bm{R})} = e^{-i\bm{q}\cdot\bm{R}}e^{i\bm{q}\cdot\rr} . $$

したがって

$$ \begin{equation} \hat{P}_{\kk}\ket{\bm{q}} = \left[\frac{1}{N}\sum_{\bm{n}}e^{i(\kk-\bm{q})\cdot\bm{R}_{\bm{n}}}\right]\ket{\bm{q}} . \label{eq:18-proj-pw} \end{equation} $$

ステップ2:括弧の中の幾何級数を評価する．$\bm{q}$ が\eqref{eq:18-sc-expand}に現れる形 $\bm{q} = \bm{K}+\bm{G}_s$ であり，かつ $\kk$ も $\bm{K}$ と同じ折り畳み類に属する($\kk = \bm{K}+\bm{G}_s'$)場合を考える．すると $\kk-\bm{q} = \bm{G}_s'-\bm{G}_s$ はスーパーセル逆格子ベクトルであり，$\bm{G}_s'-\bm{G}_s = \sum_i m_i\bm{B}_i$ と書ける．$\bm{B}_i\cdot\bm{a}_j = (2\pi/N_i)\delta_{ij}$ だから

\begin{align} \frac{1}{N}\sum_{\bm{n}}e^{i(\bm{G}_s'-\bm{G}_s)\cdot\bm{R}_{\bm{n}}} &= \prod_{i=1}^{3}\left[\frac{1}{N_i}\sum_{n_i=0}^{N_i-1}\exp\!\left(\frac{2\pi i\,m_i n_i}{N_i}\right)\right] \label{eq:18-geo1}\\[2pt] &= \prod_{i=1}^{3}\delta_{m_i \equiv 0 \ (\mathrm{mod}\ N_i)} . \label{eq:18-geo2} \end{align}

\eqref{eq:18-geo1}で3方向の和が積に分解し，\eqref{eq:18-geo2}では等比級数の公式

$$ \frac{1}{N_i}\sum_{n=0}^{N_i-1}z^n = \frac{1}{N_i}\cdot\frac{z^{N_i}-1}{z-1} = 0 \quad(z = e^{2\pi i m_i/N_i}\ne1), \qquad = 1\quad(z=1) $$

を使った($z^{N_i} = e^{2\pi i m_i} = 1$ なので分子は常にゼロ，$z\ne1$ なら和はゼロ)．$m_i$ が $N_i$ の倍数であることは $\sum_i m_i\bm{B}_i = \sum_i (m_i/N_i)\bm{b}_i$ が基本セルの逆格子ベクトルであることと同値である．したがって

$$ \begin{equation} \hat{P}_{\kk}\ket{\bm{q}} = \begin{cases} \ket{\bm{q}} & (\bm{q}-\kk \in \{\bm{G}_p\})\\ 0 & (\text{それ以外}) \end{cases} \label{eq:18-proj-action} \end{equation} $$

すなわち $\hat{P}_{\kk}$ は「波数が $\kk$ と基本逆格子ベクトルの分だけ違う平面波」を通し，それ以外を消す．

ステップ3:射影演算子であることの確認．\eqref{eq:18-proj-action}より，$\hat{P}_{\kk}$ は正規直交基底 $\{\ket{\bm{q}}\}$ において固有値 $0$ または $1$ の対角行列である．対角で実固有値なのでエルミート，固有値が $0,1$ なので $\hat{P}_{\kk}^2=\hat{P}_{\kk}$．よって射影演算子である．

ステップ4:完全性．各平面波 $\ket{\bm{q}}$ は，$\bm{q}$ を基本セルのブリルアン域に還元した波数 $\kk$ のちょうど1つのクラスに属する．したがって，$\bm{K}$ と同じ折り畳み類に属する $N$ 個の $\kk$ について足すと

$$ \begin{equation} \sum_{\kk} \hat{P}_{\kk} = \hat{1} \label{eq:18-proj-complete} \end{equation} $$

となる．∎

これで準備は整った．スペクトル重みを定義する．

定義:アンフォールディングのスペクトル重み

$$ \begin{equation} W_{\bm{K}J}(\kk) \equiv \braket{\Psi_{\bm{K}J}|\hat{P}_{\kk}|\Psi_{\bm{K}J}} = \sum_{\bm{G}_p}\abs{\braket{\kk+\bm{G}_p|\Psi_{\bm{K}J}}}^2 \label{eq:18-weight} \end{equation} $$

第2の等号は\eqref{eq:18-proj-action}と\eqref{eq:18-sc-expand}から直ちに従う:$\hat{P}_{\kk}$ が通す平面波はちょうど $\{\kk+\bm{G}_p\}$ の集合だからである．

$$ \begin{equation} W_{\bm{K}J}(\kk) = \sum_{\bm{G}_p}\abs{\braket{\kk+\bm{G}_p|\Psi_{\bm{K}J}}}^2, \qquad \sum_{\kk}W_{\bm{K}J}(\kk) = 1 \label{eq:18-weight-key} \end{equation} $$

和則 $\sum_{\kk}W = 1$ は完全性\eqref{eq:18-proj-complete}と規格化 $\braket{\Psi|\Psi}=1$ から出る．この和則は実装のチェックとして有用である(重みの合計が $1$ にならなければ，どこかが間違っている)．

最後に，アンフォールドしたスペクトル関数を

$$ \begin{equation} A(\kk,\omega) = \sum_{J} W_{\bm{K}J}(\kk)\,\delta\bigl(\omega - \varepsilon_{\bm{K}J}\bigr) \label{eq:18-unfolded-A} \end{equation} $$

と定義する($\bm{K}$ は $\kk$ をスーパーセルのブリルアン域に還元した波数)．第13章で導入したGreen関数 $\hat{G}(z) = \sum_{\bm{K}J}\ket{\bm{K}J}\bra{\bm{K}J}/(z-\varepsilon_{\bm{K}J})$ を使えば，$\hat{A}(\omega) = -\pi^{-1}\mathrm{Im}\,\hat{G}(\omega+i0^+)$ の対角要素を $\hat{P}_{\kk}$ で重み付けした量にほかならない．これがARPES強度と直接比較すべき理論量である．

物理的意味:「ぼやけたバンド」は何を表しているか

摂動がない場合．スーパーセルが基本セルの単純な繰り返しなら，各固有状態は純粋な基本セルBloch状態である．したがってある1つの $\kk$ について $W = 1$，他はすべて $0$ となり，折り畳まれたバンドは完全にほどけて元の分散を復元する(図18.9(c)の大きな丸)．

摂動がある場合．欠陥・不純物・合金の乱れ・吸着子が入ると，基本セルの並進対称性が破れる．すると固有状態は複数の $\kk$ 成分の混合になり，$W$ が複数の $\kk$ に分散する．$\kk$ 方向への広がりは，摂動によるBloch状態の寿命(あるいは平均自由行程 $\ell$)と $\Delta k \sim 1/\ell$ の関係で結ばれる．ARPESで観測される「ぼやけたバンド」は，まさにこの $W$ の分散を見ている．したがって，アンフォールドしたスペクトル関数の $\kk$ 方向の幅を実験の幅と比べれば，乱れの強さを定量的に議論できる．

ゴーストバンド．図18.9(c)の小さな赤い丸は，完全結晶では重みが厳密にゼロだが，摂動が入ると有限の重みを獲得するバンドである．実験でこれが弱く見えれば，それは対称性の破れの直接的な証拠になる．折り畳まれたバンド図をそのまま眺めていては，「本物のバンド」と「折り畳みの人工物」の区別さえつかない．アンフォールディングは，この区別を定量的に付ける手続きである．

補足:非直交局在基底(LCAO)でのスペクトル重み

平面波基底での式\eqref{eq:18-weight-key}は透明だが，局在基底(第14章のLCPAO)を使う実装では書き換えが必要である．スーパーセルのKohn–Sham状態を

$$ \ket{\Psi_{\bm{K}J}} = \sum_{\bm{R}}\sum_{M} e^{i\bm{K}\cdot\bm{R}}\,c^{\bm{K}J}_{M}\ket{\phi_{M\bm{R}}} $$

($M$ は原子軌道の指標，$\bm{R}$ はスーパーセル格子ベクトル，$\ket{\phi_{M\bm{R}}}$ は位置 $\bm{\tau}_M+\bm{R}$ に置かれた原子軌道)と書くと，同じ射影 $\hat{P}_{\kk}$ の期待値は

$$ \begin{equation} W_{\bm{K}J}(\kk) = \sum_{M,N}\sum_{\bm{R}} c^{\bm{K}J\,*}_{M}\,c^{\bm{K}J}_{N}\; e^{i\kk\cdot(\bm{\tau}_N+\bm{R}-\bm{\tau}_M)}\; S_{M\bm{0},\,N\bm{R}} \label{eq:18-weight-lcao} \end{equation} $$

という形になる．$S_{M\bm{0},N\bm{R}} = \braket{\phi_{M\bm{0}}|\phi_{N\bm{R}}}$ は重なり行列である．3つの構造的な要点を押さえておきたい．

位相因子が本質．$\kk$ 依存性はすべて $e^{i\kk\cdot(\bm{\tau}_N+\bm{R}-\bm{\tau}_M)}$ に入っている．この因子こそが「基本セルのBloch対称性への射影」を実装しており，これを落とすと単なる占有数の分解になってしまい，アンフォールディングにならない．
重なり行列が必要．基底が非直交なので，$S$ を落とすと和則 $\sum_{\kk}W=1$ が破れる．
参照セルのBloch関数を知る必要はない．式\eqref{eq:18-weight-lcao}に現れるのはスーパーセルの計算結果($c$, $S$)と参照セルの幾何(どの並進で写り合うか)だけである．参照セルの電子状態を別途計算する必要はない．だからこそ，「基本セルが存在しない系」(合金，非整合構造)にも形式的に適用できる．
軌道分解ができる．和を特定の軌道 $M$(たとえば表面第一層のZr $d$ 軌道)に制限すれば，その軌道の寄与だけを取り出した重み $W^{\kk}_{\bm{K}JM}$ が得られる．「このバンドは表面の $d$ 軌道由来である」という主張を，図の上で定量的に示せる．

例:スラブ計算における表面状態の同定

表面を扱うには，真空層を挟んだスラブ模型を使う．スラブは面直方向にバルクの単位胞を何層も積んだスーパーセルだから，面直方向のバンドは激しく折り畳まれる($n$ 層なら $n$ 本)．この折り畳まれたバンド図から「どれが表面状態でどれがバルク由来か」を読み取るのは，目視ではほとんど不可能である．

アンフォールディングを使えば話は明快になる．バルクの単位胞を参照セルとしてスラブのバンドをアンフォールドすると，バルク由来の状態はバルクバンドの位置に重みを集中させる．一方，表面状態はバルクバンドギャップの中に，バルクバンドとは無関係な位置に現れる．さらに軌道分解した重み($W^{\kk}_{\bm{K}JM}$，$M$ を表面第一層の軌道に制限)を円の大きさで描けば，その状態がどれだけ表面に局在しているかが一目で分かる．ZrB$_2$(0001) 表面ではこの手続きによって，2本の表面状態(Zr $d$ 軌道由来)と，量子閉じ込めによるスラブ固有の状態が明確に区別された．この解析は18.9節の事例研究の出発点になる．

18.8 誘電関数とKubo–Greenwoodの式

図18.1の最後の矢印である．価電子帯から伝導帯への遷移は，電子を試料の外に出さないので，測定は「光がどれだけ吸収・反射・屈折されるか」を通じて行われる．これらのマクロな光学量をすべて統括するのが複素誘電関数 $\varepsilon(\omega)$ である．本節では，$\varepsilon(\omega)$ の虚部を独立粒子近似(RPA)で完全に導出し，Kubo–Greenwood形式との対応を付け，Kramers–Kronig関係で実部を復元し，最後に近似の限界と，その先の理論への展望を述べる．

18.8.1 複素誘電関数とマクロ光学量

定義:複素誘電関数

物質に振動電場 $\bm{E}(\omega)$ をかけたとき，電気変位 $\bm{D}$ が線形に応答するとして

$$ \begin{equation} D_\alpha(\omega) = \sum_\beta \varepsilon_{\alpha\beta}(\omega)\,E_\beta(\omega), \qquad \varepsilon(\omega) = \varepsilon_1(\omega) + i\,\varepsilon_2(\omega) \label{eq:18-eps-def} \end{equation} $$

で定義されるテンソル $\varepsilon_{\alpha\beta}$ を複素誘電関数と呼ぶ．等方的な系(立方晶や多結晶)ではスカラー $\varepsilon(\omega)$ になる．実部 $\varepsilon_1$ は分極(エネルギーを蓄える成分)，虚部 $\varepsilon_2$ は散逸(吸収)を表す．以下ではGauss単位系を用いる．

誘電関数は，周波数帯によって異なる自由度が担う．$\mathrm{GHz}$ 帯では分子の配向分極，$\mathrm{THz}$ 帯ではイオンの変位(光学フォノン)，$\mathrm{PHz}$ 帯(可視光〜紫外)では電子の分極である．第一原理電子状態計算が担当するのはこの最後の電子分極であり，その周波数は $0.1$–$10$ eV に対応する．

18.8.2 $\varepsilon_2(\omega)$ の導出

これから，フェルミの黄金律\eqref{eq:18-golden}から出発して $\varepsilon_2(\omega)$ の表式を組み立てる．方針は「ミクロに計算した吸収パワーを，マクロな誘電関数で書いた吸収パワーと等置する」ことである．

導出:独立粒子近似の $\varepsilon_2(\omega)$

ステップ1:摂動と双極子近似．18.6.4節と同様，摂動は $\hat{H}' = (1/c)\bm{A}\cdot\hat{\bm{p}}$ である．入射光を

$$ \bm{A}(\rr,t) = A_0\hat{\bm{e}}\cos(\bm{q}\cdot\rr-\omega t) $$

とする．可視光では $\lambda \simeq 5000$ Å，$q = 2\pi/\lambda \simeq 1.3\times10^{-3}$ Å$^{-1}$ であり，ブリルアン域の大きさ $\sim 1$ Å$^{-1}$ の $10^{-3}$ にすぎない．したがって $e^{i\bm{q}\cdot\rr}\simeq1$(双極子近似)とすると同時に，遷移は $\kk$ を変えない「垂直遷移」とみなせる．吸収に対応する摂動演算子は

$$ \hat{V} = \frac{A_0}{2c}\,\hat{\bm{e}}\cdot\hat{\bm{p}} . $$

ステップ2:黄金律を適用する．占有バンド $v$ から非占有バンド $c$ への同じ $\kk$ での遷移確率は，\eqref{eq:18-golden}より

$$ \begin{equation} w_{v\kk\to c\kk} = 2\pi\left(\frac{A_0}{2c}\right)^2 \abs{\braket{c\kk|\hat{\bm{e}}\cdot\hat{\bm{p}}|v\kk}}^2\, \delta\bigl(\varepsilon_{c\kk}-\varepsilon_{v\kk}-\omega\bigr) . \label{eq:18-eps-golden} \end{equation} $$

ステップ3:単位体積あたりの吸収パワー(ミクロ)．1回の遷移でエネルギー $\omega$ が吸収される．すべての $\kk$，すべてのバンド対について足し，スピン縮重度 $2$ を掛け，体積 $\Omega$ で割ると

\begin{align} P_{\mathrm{micro}} &= \frac{2\omega}{\Omega}\sum_{\kk}\sum_{v,c} w_{v\kk\to c\kk} \label{eq:18-pmicro1}\\[2pt] &= \frac{2\omega}{\Omega}\cdot 2\pi\frac{A_0^2}{4c^2} \sum_{\kk,v,c}\abs{\braket{c\kk|\hat{\bm{e}}\cdot\hat{\bm{p}}|v\kk}}^2\delta(\varepsilon_{c\kk}-\varepsilon_{v\kk}-\omega) \label{eq:18-pmicro2}\\[2pt] &= \frac{\pi\omega A_0^2}{c^2\,\Omega} \sum_{\kk,v,c}\abs{\braket{c\kk|\hat{\bm{e}}\cdot\hat{\bm{p}}|v\kk}}^2\delta(\varepsilon_{c\kk}-\varepsilon_{v\kk}-\omega) . \label{eq:18-pmicro} \end{align}

\eqref{eq:18-pmicro}では $2\cdot2\pi/4 = \pi$ と整理した．

ステップ4:単位体積あたりの吸収パワー(マクロ)．マクロには，電場が分極電流にする仕事が吸収パワーである．$\bm{E}(t) = \mathrm{Re}\bigl[\bm{E}_0e^{-i\omega t}\bigr]$，分極 $\bm{P} = \chi\bm{E}$，$\chi = (\varepsilon-1)/4\pi$，分極電流 $\bm{j} = \partial\bm{P}/\partial t$ である．$\bm{j}(t) = \mathrm{Re}\bigl[-i\omega\chi\bm{E}_0e^{-i\omega t}\bigr]$ だから，2つの実部の積の時間平均の公式 $\braket{\mathrm{Re}[a e^{-i\omega t}]\,\mathrm{Re}[b e^{-i\omega t}]} = \frac12\mathrm{Re}[a^*b]$ を使って

\begin{align} P_{\mathrm{macro}} = \braket{\bm{j}\cdot\bm{E}} &= \frac{1}{2}\mathrm{Re}\Bigl[\bigl(-i\omega\chi\bm{E}_0\bigr)^*\cdot\bm{E}_0\Bigr] \label{eq:18-pmacro1}\\[2pt] &= \frac{1}{2}\mathrm{Re}\bigl[\,i\omega\chi^*\bigr]\abs{\bm{E}_0}^2 \label{eq:18-pmacro2}\\[2pt] &= \frac{\omega}{2}\,\mathrm{Re}\!\left[\,i\cdot\frac{\varepsilon_1-i\varepsilon_2-1}{4\pi}\right]\abs{\bm{E}_0}^2 \label{eq:18-pmacro3}\\[2pt] &= \frac{\omega}{8\pi}\,\varepsilon_2(\omega)\,\abs{\bm{E}_0}^2 . \label{eq:18-pmacro} \end{align}

\eqref{eq:18-pmacro3}から\eqref{eq:18-pmacro}へは，$i(\varepsilon_1-1-i\varepsilon_2) = i(\varepsilon_1-1)+\varepsilon_2$ の実部が $\varepsilon_2$ であることを使った．散逸は $\varepsilon_2$ だけが担うという重要な事実がここで確認された．

ステップ5:$A_0$ と $E_0$ を結ぶ．$\bm{E} = -(1/c)\partial\bm{A}/\partial t$ より

$$ \bm{E} = -\frac{1}{c}\,A_0\hat{\bm{e}}\,\omega\sin(\bm{q}\cdot\rr-\omega t) \qquad\Longrightarrow\qquad \abs{\bm{E}_0} = \frac{\omega A_0}{c} . $$

ステップ6:等置して解く．\eqref{eq:18-pmicro}$=$\eqref{eq:18-pmacro}とおく:

\begin{align} \frac{\omega}{8\pi}\varepsilon_2\cdot\frac{\omega^2A_0^2}{c^2} &= \frac{\pi\omega A_0^2}{c^2\Omega}\sum_{\kk,v,c}\abs{M_{cv}}^2\delta(\cdots) \label{eq:18-eq1}\\[2pt] \frac{\omega^3A_0^2}{8\pi c^2}\,\varepsilon_2 &= \frac{\pi\omega A_0^2}{c^2\Omega}\sum_{\kk,v,c}\abs{M_{cv}}^2\delta(\cdots) \label{eq:18-eq2}\\[2pt] \varepsilon_2(\omega) &= \frac{8\pi^2}{\omega^2\,\Omega}\sum_{\kk,v,c}\abs{M_{cv}}^2\delta(\varepsilon_{c\kk}-\varepsilon_{v\kk}-\omega) . \label{eq:18-eq3} \end{align}

\eqref{eq:18-eq2}から\eqref{eq:18-eq3}へは，両辺を $\omega^3A_0^2/(8\pi c^2)$ で割った($A_0$, $c$ が完全に消えることに注意——物理量が入射強度によらないという当然の帰結である)．ここで $M_{cv} = \braket{c\kk|\hat{\bm{e}}\cdot\hat{\bm{p}}|v\kk}$ と略記した．∎

原子単位系($e = m_e = \hbar = 1$)を解いて，SI・Gauss単位で使われる形に戻すと次のようになる．

$$ \begin{equation} \varepsilon_2(\omega) = \frac{4\pi^2e^2}{m_e^2\omega^2}\cdot\frac{2}{\Omega} \sum_{\kk}\sum_{v}\sum_{c} \abs{\hat{\bm{e}}\cdot\braket{c\kk|\hat{\bm{p}}|v\kk}}^2\, \delta\bigl(\varepsilon_{c\kk}-\varepsilon_{v\kk}-\hbar\omega\bigr) \label{eq:18-eps2-key} \end{equation} $$

因子 $2$ はスピン縮重度である．この式を読み解いておこう．

$\delta$ 関数:エネルギー保存．光子1個で $\varepsilon_c-\varepsilon_v = \hbar\omega$ を満たす遷移だけが起こる．
同じ $\kk$ どうし:双極子近似(光子の運動量が無視できる)の帰結．バンド図の上では「縦向きの矢印」だけが許される．
$1/\omega^2$:運動量形式から来る因子．長さ形式に直せば $\abs{M_p}^2=\omega^2\abs{M_r}^2$ なので $\varepsilon_2\propto\abs{M_r}^2$ となり，$\omega$ 依存性は行列要素側に移る(18.6.4節)．
$1/\Omega$ とセル体積:$\sum_{\kk}$ が体積に比例して増えるので，この因子で示強変数になる．連続極限では $(2/\Omega)\sum_{\kk}\to 2\int_{\mathrm{BZ}}\dd^3k/(2\pi)^3$ とすればよく，これを使うと

$$ \begin{equation} \varepsilon_2(\omega) = \frac{e^2}{\pi m_e^2\omega^2}\sum_{v,c}\int_{\mathrm{BZ}}\dd^3k\; \abs{\hat{\bm{e}}\cdot\bm{p}_{cv}(\kk)}^2\,\delta\bigl(\varepsilon_{c\kk}-\varepsilon_{v\kk}-\hbar\omega\bigr) \label{eq:18-eps2-cont} \end{equation} $$

となる(係数は $4\pi^2\cdot2/(2\pi)^3 = 1/\pi$ から)．部分占有(金属や有限温度)を扱う場合は，$\delta$ 関数の前に占有因子の差 $\bigl(f_{v\kk}-f_{c\kk}\bigr)$ を掛ける．これは「吸収 $-$ 誘導放出」の正味を表しており，$f_v=1$, $f_c=0$ なら $1$ になって上式に帰着する．

物理的意味:結合状態密度とvan Hove特異点

行列要素 $\abs{\bm{p}_{cv}}^2$ が $\kk$ に鈍感だと仮定すると，\eqref{eq:18-eps2-cont}は

$$ \varepsilon_2(\omega)\;\propto\;\frac{\abs{\bm{p}_{cv}}^2}{\omega^2}\,J(\omega), \qquad J(\omega) = \frac{2}{(2\pi)^3}\sum_{v,c}\int\dd^3k\,\delta\bigl(\varepsilon_{c\kk}-\varepsilon_{v\kk}-\omega\bigr) $$

と書ける．$J(\omega)$ を結合状態密度(joint density of states, JDOS)と呼ぶ．第3章・第13章で状態密度を扱ったときと同じ論法により，$J(\omega)$ は

$$ \nabla_{\kk}\bigl(\varepsilon_{c\kk}-\varepsilon_{v\kk}\bigr) = \bm{0} $$

を満たす $\kk$(臨界点)で特異性を持つ．この点では価電子帯と伝導帯が「平行」になり，多数の $\kk$ が同じ遷移エネルギーに集まる．シリコンの $\varepsilon_2$ に現れる $3.4$ eV の $E_1$ ピークと $4.3$ eV の $E_2$ ピークは，いずれもこのvan Hove特異点である．光学スペクトルのピーク位置は，特定の1点のバンドギャップではなく，バンドが平行になる領域の遷移エネルギーを表している——この点は誤解されやすいので注意したい．

実装上の注意:$\delta$ 関数のスメアリングと $\kk$ 点収束

数値計算では $\delta$ 関数を有限幅の関数に置き換える．Gauss型 $\delta(x)\to(\eta\sqrt{2\pi})^{-1}e^{-x^2/2\eta^2}$ かLorentz型 $\delta(x)\to\eta/[\pi(x^2+\eta^2)]$ が使われる．ここで決定的に重要なのは，スメアリング幅 $\eta$ と $\kk$ 点メッシュの細かさを独立に決めてはいけないことである．

$\eta$ が細かすぎると，有限個の $\kk$ 点に対応する離散的なスパイクがそのまま残り，スペクトルがギザギザになる．
$\eta$ が粗すぎると，本来の構造(臨界点のピーク)がなまされて消える．

目安は「隣り合う $\kk$ 点の間での遷移エネルギーの変化 $\abs{\nabla_{\kk}(\varepsilon_c-\varepsilon_v)}\Delta k$ と同程度に $\eta$ を取る」ことである．したがって $\kk$ 点を倍に増やしたら $\eta$ も半分にして，両者を同時に収束させなければならない．バンド構造の計算より1桁多い $\kk$ 点が必要になるのが普通である(シリコンで $20\times20\times20$ 以上)．四面体法(本書では扱わない)を使えばスメアリングなしで積分できるが，実装は複雑になる．

18.8.3 Kubo–Greenwood形式との対応

誘電関数と等価な情報を，電流応答の言葉で書いたものが光学伝導度 $\sigma(\omega)$ である．Maxwell方程式で「変位電流 $\partial\bm{D}/\partial t$」と「伝導電流 $\bm{j} = \sigma\bm{E}$」をどちらに数えるかの違いにすぎず，両者は

$$ \begin{equation} \varepsilon(\omega) = 1 + \frac{4\pi i}{\omega}\sigma(\omega) \qquad\Longleftrightarrow\qquad \sigma_1(\omega) = \frac{\omega}{4\pi}\varepsilon_2(\omega), \quad \sigma_2(\omega) = -\frac{\omega}{4\pi}\bigl[\varepsilon_1(\omega)-1\bigr] \label{eq:18-sigma-eps} \end{equation} $$

で結ばれる(Gauss単位)．したがって\eqref{eq:18-eps2-key}から直ちに

$$ \begin{equation} \sigma_1(\omega) = \frac{\pi e^2}{m_e^2\,\omega}\cdot\frac{2}{\Omega} \sum_{\kk,v,c}\abs{\hat{\bm{e}}\cdot\braket{c\kk|\hat{\bm{p}}|v\kk}}^2\, \delta\bigl(\varepsilon_{c\kk}-\varepsilon_{v\kk}-\hbar\omega\bigr) \label{eq:18-sigma1} \end{equation} $$

が得られる．これがKubo–Greenwoodの式の実部である．より一般には，線形応答理論(Kuboの公式)から，テンソル形の複素伝導度が

$$ \begin{equation} \sigma_{\alpha\beta}(\omega) = \frac{ie^2\hbar}{m_e^2\,\Omega} \sum_{\kk}\sum_{J,J'} \frac{\bigl(f_{J'\kk}-f_{J\kk}\bigr)\, \braket{J\kk|\hat{p}_\alpha|J'\kk}\braket{J'\kk|\hat{p}_\beta|J\kk}} {\bigl(\varepsilon_{J\kk}-\varepsilon_{J'\kk}\bigr) \bigl(\varepsilon_{J\kk}-\varepsilon_{J'\kk}+\hbar\omega+i\eta\bigr)} \label{eq:18-kubo-greenwood} \end{equation} $$

という形で書ける($J, J'$ はバンド指標，$\eta\to0^+$)．この式の実部が\eqref{eq:18-sigma1}に帰着することは，Sokhotski–Plemeljの公式

$$ \lim_{\eta\to0^+}\frac{1}{x+i\eta} = \mathcal{P}\frac{1}{x} - i\pi\delta(x) $$

を分母の第2因子に適用すればよい．$-i\pi\delta$ の項が前置因子 $i$ と掛かって実数になり，$\delta$ 関数によって $\varepsilon_{J}-\varepsilon_{J'} = -\hbar\omega$ が課され，残りの分母 $(\varepsilon_J-\varepsilon_{J'}) = -\hbar\omega$ と占有因子の符号 $(f_{J'}-f_J) = -1$($J$ 占有・$J'$ 非占有のとき)が打ち消しあって\eqref{eq:18-sigma1}が出る．

物理的意味:Kubo公式のご利益

\eqref{eq:18-kubo-greenwood}のような形の式を導くのがKuboの線形応答理論である．その核心は，「摂動に対する応答」を「平衡状態における揺らぎの相関関数」で表すことにある(揺動散逸定理)．式\eqref{eq:18-kubo-greenwood}では，電流演算子の行列要素と，平衡分布 $f$ だけで応答が書けている．非平衡の計算をしなくてよいという点が実用上の最大の利点である．第17章のLandauer公式も，線形応答領域ではKubo公式から導ける．$\varepsilon_2$ を「摂動 $\to$ 吸収パワー」の道筋で導いた本節の計算は，Kubo公式の特別な場合を手で追ったことになる．

18.8.4 Kramers–Kronig関係:因果律が実部と虚部を結ぶ

ここまでで得られたのは $\varepsilon_2$(散逸)だけである．反射率や屈折率を計算するには実部 $\varepsilon_1$ も必要だが，これは独立に計算しなくてもよい．因果律だけから，$\varepsilon_1$ と $\varepsilon_2$ は互いに決定しあう．これがKramers–Kronig(KK)関係である．

数学ノート:因果律 $\to$ 上半平面での解析性 $\to$ Cauchyの積分公式

ステップ1:因果律を式にする．線形応答を時間領域で書くと

$$ P(t) = \int_{-\infty}^{\infty}\chi(\tau)\,E(t-\tau)\,\dd\tau . $$

「原因は結果に先立つ」——時刻 $t$ の分極は，それより前の電場にしか依存しない．すなわち

$$ \chi(\tau) = 0 \qquad (\tau \lt 0) . $$

ステップ2:上半平面での解析性を示す．Fourier変換

$$ \tilde\chi(\omega) = \int_{-\infty}^{\infty}\chi(\tau)e^{i\omega\tau}\dd\tau = \int_{0}^{\infty}\chi(\tau)e^{i\omega\tau}\dd\tau $$

において，$\omega$ を複素数 $\omega = \omega' + i\omega''$ に拡張する．すると

$$ \abs{e^{i\omega\tau}} = \abs{e^{i\omega'\tau}}\,e^{-\omega''\tau} = e^{-\omega''\tau} $$

である．積分範囲が $\tau\ge0$ に限られているので，$\omega''\gt0$(上半平面)では被積分関数が指数的に減衰し，積分は絶対収束する．よって $\tilde\chi(\omega)$ は上半平面全体で解析的である．もし因果律がなく $\tau\lt0$ からも寄与があれば，そこで $e^{-\omega''\tau}$ は発散し，この結論は得られない．因果律と解析性は同じことの2つの言い方である．

ステップ3:Cauchyの積分公式を使う．実軸上の点 $\omega$ を選び，複素平面上で関数 $f(z) = \tilde\chi(z)/(z-\omega)$ を考える．積分路 $C$ を次のように取る:実軸を $-R$ から $\omega-\delta$ まで，次に $z=\omega$ を上側から小さく回り込む半径 $\delta$ の半円(時計回り)，次に実軸を $\omega+\delta$ から $R$ まで，最後に上半平面の半径 $R$ の大半円(反時計回り)．この閉曲線の内部に $f$ の特異点はない($\tilde\chi$ は上半平面で解析的，極 $z=\omega$ は小半円で除外した)．よってCauchyの定理より $\oint_C f(z)\dd z = 0$ である．

各部分の寄与は次のとおり．

実軸部分:$\delta\to0$，$R\to\infty$ で Cauchy の主値積分 $\displaystyle\mathcal{P}\!\int_{-\infty}^{\infty}\frac{\tilde\chi(\omega')}{\omega'-\omega}\dd\omega'$ になる．
小半円:$z = \omega+\delta e^{i\theta}$，$\theta$ が $\pi$ から $0$ まで(時計回り)．$\dd z = i\delta e^{i\theta}\dd\theta$ だから $$\int_{\pi}^{0}\frac{\tilde\chi(\omega)}{\delta e^{i\theta}}i\delta e^{i\theta}\dd\theta = i\tilde\chi(\omega)\int_\pi^0\dd\theta = -i\pi\tilde\chi(\omega) .$$ 半円は「留数の半分」を，時計回りなので符号を変えて拾う．
大半円:物理的な応答関数は $\abs{z}\to\infty$ で $\tilde\chi\to0$ になる(無限に速い振動には物質は応答できない)ので，Jordanの補題によりゼロになる．

合計をゼロと置いて

$$ \mathcal{P}\!\int_{-\infty}^{\infty}\frac{\tilde\chi(\omega')}{\omega'-\omega}\dd\omega' - i\pi\tilde\chi(\omega) = 0 \qquad\Longrightarrow\qquad \tilde\chi(\omega) = \frac{1}{i\pi}\mathcal{P}\!\int_{-\infty}^{\infty}\frac{\tilde\chi(\omega')}{\omega'-\omega}\dd\omega' . $$

ステップ4:実部と虚部に分ける．$1/i = -i$ を使い，$\tilde\chi = \chi_1+i\chi_2$ を代入すると

$$ \chi_1 + i\chi_2 = -\frac{i}{\pi}\mathcal{P}\!\int\frac{\chi_1(\omega')}{\omega'-\omega}\dd\omega' + \frac{1}{\pi}\mathcal{P}\!\int\frac{\chi_2(\omega')}{\omega'-\omega}\dd\omega' . $$

両辺の実部・虚部を比較して

$$ \chi_1(\omega) = \frac{1}{\pi}\mathcal{P}\!\int_{-\infty}^{\infty}\frac{\chi_2(\omega')}{\omega'-\omega}\dd\omega', \qquad \chi_2(\omega) = -\frac{1}{\pi}\mathcal{P}\!\int_{-\infty}^{\infty}\frac{\chi_1(\omega')}{\omega'-\omega}\dd\omega' . $$

あとは，負の周波数の積分を正の周波数だけに畳み込めば実用形になる．そのために使うのが応答の実数性である．

導出:誘電関数のKramers–Kronig関係(正の周波数のみの形)

ステップ1:偶奇性を導く．$\chi(\tau)$ は実関数である(実の電場が実の分極を生む)．したがって

$$ \tilde\chi(-\omega) = \int_0^\infty\chi(\tau)e^{-i\omega\tau}\dd\tau = \left[\int_0^\infty\chi(\tau)e^{i\omega\tau}\dd\tau\right]^* = \tilde\chi^*(\omega) \qquad(\omega\ \text{は実}) $$

である($\chi(\tau)$ が実なので複素共役を取ると指数の符号だけが変わる)．実部・虚部で書けば

$$ \varepsilon_1(-\omega) = \varepsilon_1(\omega)\ (\text{偶関数}), \qquad \varepsilon_2(-\omega) = -\varepsilon_2(\omega)\ (\text{奇関数}) . $$

ステップ2:$\varepsilon$ の関係式に書き換える．Gauss単位で $\varepsilon = 1+4\pi\tilde\chi$ だから，$\varepsilon_1-1 = 4\pi\chi_1$，$\varepsilon_2 = 4\pi\chi_2$ を代入して

$$ \varepsilon_1(\omega) - 1 = \frac{1}{\pi}\mathcal{P}\!\int_{-\infty}^{\infty}\frac{\varepsilon_2(\omega')}{\omega'-\omega}\dd\omega' . $$

ステップ3:負の領域を折り返す．積分を $\int_{-\infty}^{0}+\int_{0}^{\infty}$ に分け，前者で $\omega' = -u$($\dd\omega' = -\dd u$，積分範囲は $u:\infty\to0$)と置換する:

\begin{align} \int_{-\infty}^{0}\frac{\varepsilon_2(\omega')}{\omega'-\omega}\dd\omega' &= \int_{\infty}^{0}\frac{\varepsilon_2(-u)}{-u-\omega}(-\dd u) \label{eq:18-kk1}\\[2pt] &= \int_{0}^{\infty}\frac{\varepsilon_2(-u)}{-u-\omega}\dd u \label{eq:18-kk2}\\[2pt] &= \int_{0}^{\infty}\frac{-\varepsilon_2(u)}{-(u+\omega)}\dd u = \int_{0}^{\infty}\frac{\varepsilon_2(u)}{u+\omega}\dd u . \label{eq:18-kk3} \end{align}

\eqref{eq:18-kk2}では積分の上下限を入れ替えて符号を戻し，\eqref{eq:18-kk3}では $\varepsilon_2$ の奇関数性 $\varepsilon_2(-u)=-\varepsilon_2(u)$ を使って分子と分母の符号を整理した．

ステップ4:2つを足す．積分変数を $\omega'$ に統一して

\begin{align} \varepsilon_1(\omega)-1 &= \frac{1}{\pi}\mathcal{P}\!\int_0^\infty\varepsilon_2(\omega')\left[\frac{1}{\omega'-\omega}+\frac{1}{\omega'+\omega}\right]\dd\omega' \label{eq:18-kk4}\\[2pt] &= \frac{1}{\pi}\mathcal{P}\!\int_0^\infty\varepsilon_2(\omega')\cdot\frac{(\omega'+\omega)+(\omega'-\omega)}{(\omega'-\omega)(\omega'+\omega)}\dd\omega' \label{eq:18-kk5}\\[2pt] &= \frac{2}{\pi}\mathcal{P}\!\int_0^\infty\frac{\omega'\,\varepsilon_2(\omega')}{\omega'^2-\omega^2}\dd\omega' . \label{eq:18-kk6} \end{align}

\eqref{eq:18-kk5}で通分し，\eqref{eq:18-kk6}で分子 $2\omega'$ と分母 $\omega'^2-\omega^2$ に整理した．同じ手続きを逆向きに行えば

$$ \varepsilon_2(\omega) = -\frac{2\omega}{\pi}\mathcal{P}\!\int_0^\infty\frac{\varepsilon_1(\omega')-1}{\omega'^2-\omega^2}\dd\omega' $$

も得られる．∎

$$ \begin{equation} \varepsilon_1(\omega) = 1 + \frac{2}{\pi}\,\mathcal{P}\!\int_0^\infty \frac{\omega'\,\varepsilon_2(\omega')}{\omega'^2-\omega^2}\,\dd\omega' \label{eq:18-kk-key} \end{equation} $$

実務上の使い方はこうである．まず\eqref{eq:18-eps2-key}で $\varepsilon_2(\omega)$ を十分高い周波数まで計算し，次に\eqref{eq:18-kk-key}の数値積分で $\varepsilon_1(\omega)$ を得る．注意点は積分の打ち切りである．$\varepsilon_2$ を有限の $\omega_{\max}$ までしか計算していないと，$\varepsilon_1$ が系統的に過小評価される．$\omega_{\max}$ は興味のある領域の $3$–$5$ 倍は取り，十分な数の空バンドを含めなければならない．

例:KK関係から出る2つのチェック

(a) 静的誘電率．\eqref{eq:18-kk-key}で $\omega=0$ とおくと

$$ \varepsilon_1(0) = 1 + \frac{2}{\pi}\int_0^\infty\frac{\varepsilon_2(\omega')}{\omega'}\dd\omega' $$

となる．すなわち静的誘電率は吸収スペクトルの重み付き積分で決まる．ギャップの小さい物質ほど低周波側に吸収が集まって $1/\omega'$ の重みが効き，$\varepsilon_1(0)$ が大きくなる．Siで $\varepsilon_1(0)\simeq11.9$，ダイヤモンドで $5.7$ という違いは，この積分の重心の違いである．

(b) f-総和則．同様の議論から

$$ \int_0^\infty\omega\,\varepsilon_2(\omega)\,\dd\omega = \frac{\pi}{2}\,\omega_p^2, \qquad \omega_p^2 = \frac{4\pi n e^2}{m_e} $$

が導かれる($n$ は価電子密度，$\omega_p$ はプラズマ振動数)．計算した $\varepsilon_2$ をこの式に入れて，価電子数が正しく再現されるかを見るのは，基底関数や空バンド数の収束を確かめる強力なテストである．

18.8.5 光学定数への変換

$\varepsilon_1$ と $\varepsilon_2$ が揃えば，測定量のすべてが計算できる．複素屈折率を $N = n+i\kappa$ と定義すると $N^2 = \varepsilon$ だから

$$ \begin{equation} \varepsilon_1 = n^2-\kappa^2,\qquad \varepsilon_2 = 2n\kappa . \label{eq:18-nk} \end{equation} $$

これを $n,\kappa$ について解こう．2式から

$$ (n^2+\kappa^2)^2 = (n^2-\kappa^2)^2 + (2n\kappa)^2 = \varepsilon_1^2+\varepsilon_2^2 \qquad\Longrightarrow\qquad n^2+\kappa^2 = \sqrt{\varepsilon_1^2+\varepsilon_2^2} $$

(平方根は正)．これと $n^2-\kappa^2 = \varepsilon_1$ を足し引きして

$$ \begin{equation} n = \sqrt{\frac{\sqrt{\varepsilon_1^2+\varepsilon_2^2}+\varepsilon_1}{2}}, \qquad \kappa = \sqrt{\frac{\sqrt{\varepsilon_1^2+\varepsilon_2^2}-\varepsilon_1}{2}} \label{eq:18-nk-solve} \end{equation} $$

を得る．以下，主要な光学量をまとめる．

量	表式	由来
吸収係数 $\alpha$	$\alpha = \dfrac{2\omega\kappa}{c}$	媒質中の電場 $\propto e^{i\omega Nz/c}$ の強度 $\abs{E}^2\propto e^{-2\omega\kappa z/c}$
垂直入射反射率 $R$	$R = \dfrac{(n-1)^2+\kappa^2}{(n+1)^2+\kappa^2}$	Fresnel係数 $r = (1-N)/(1+N)$ の $\abs{r}^2$
エネルギー損失関数 $L$	$L(\omega) = \mathrm{Im}\!\left[\dfrac{-1}{\varepsilon}\right] = \dfrac{\varepsilon_2}{\varepsilon_1^2+\varepsilon_2^2}$	高速電子が失うエネルギーの分布(EELS)．$\varepsilon\to0$ でプラズモンピーク
誘電正接 $\tan\delta$	$\tan\delta = \varepsilon_2/\varepsilon_1$	誘電体損失の指標

損失関数が登場したことに注目してほしい．式\eqref{eq:18-imfp}で，XPSの表面敏感性を支配する非弾性平均自由行程が $\mathrm{Im}[-1/\varepsilon]$ に比例すると述べた．本節で計算した誘電関数が，18.2節の話に戻ってきたわけである．本章の4本の矢印は，こうして一つの理論で結ばれている．

18.8.6 近似の限界と，その先

式\eqref{eq:18-eps2-key}は，Kohn–Sham固有値と固有関数だけから作られている．すなわち次の3つを含んでいない．

準粒子補正:KS固有値は準粒子エネルギーではない(第10章)．とくにバンドギャップは系統的に $30$–$50\,\%$ 過小評価される．$\varepsilon_2$ の立ち上がりとすべての構造が，実験より低エネルギー側にずれる．
電子-正孔相互作用(励起子効果):光吸収で生成した電子と正孔はCoulomb引力で束縛される．独立粒子近似ではこの2体効果が完全に欠落している．その結果，吸収端近傍のピーク強度が過小評価され，ギャップ内の励起子ピークが現れない．
局所場効果:実際の物質中では，外場に対して誘起された微視的な分極が不均一で，それが作る局所電場が外場と異なる．$\bm{G}\ne\bm{G}'$ の非対角成分(逆誘電行列)を扱わないと取り込めない．

図18.10 シリコンの複素誘電関数の模式図．$\varepsilon_2$(青)の2つのピークは，それぞれ $E_1$($3.4$ eV 付近)と $E_2$($4.3$ eV 付近)と呼ばれるvan Hove臨界点に対応する．$\varepsilon_1$(赤)は $\varepsilon_2$ のピークの直前で最大となり，その後急降下してゼロを切る($\varepsilon_1=0$ の点はプラズモン的な集団励起に対応する)．LDA/GGA計算はこの形を定性的に再現するが，全体が $0.5$ eV ほど低エネルギー側にずれ，$E_1$ ピークの強度が励起子効果の欠如により過小評価される．

これらに対する処方を，コストの低い順に並べておく．

定義:シザーズ補正

すべての伝導帯を一律に $\Delta$ だけ上へずらす:$\varepsilon_{c\kk}\to\varepsilon_{c\kk}+\Delta$．$\Delta$ は実験のギャップまたはGW計算の結果に合わせる．$\varepsilon_2$ のスペクトルは丸ごと $\Delta$ だけ高エネルギー側へ平行移動する．ただし\eqref{eq:18-eps2-key}の前置因子 $1/\omega^2$ はシフトしないので，強度は $(\omega/(\omega+\Delta))^2$ 程度の補正を受ける(運動量行列要素をそのまま使う場合)．最も安価で，スペクトル形状が定性的に正しいときには実用的だが，価電子帯と伝導帯の分散そのものの誤りは直せない．

GW近似(Hedin, 1965):自己エネルギーを $\Sigma = iGW$($G$ はGreen関数，$W$ は遮蔽Coulomb相互作用)で近似し，準粒子方程式を解く．バンドギャップは実験値の $0.1$ eV 程度まで改善する．第13章のGreen関数の枠組みの自然な発展であり，計算コストはDFTの $10$–$100$ 倍．
Bethe–Salpeter方程式(BSE):GWで得た準粒子の上で，電子-正孔対の2体方程式を解く．励起子束縛エネルギーと吸収端近傍のピーク強度の再分配が正しく再現され，Siの $E_1$ ピークもきちんと立ち上がる．コストはさらに $1$–$2$ 桁上がる．
時間依存DFT(TDDFT):時間依存Kohn–Sham方程式の線形応答から $\varepsilon(\omega)$ を求める．分子系では標準的手法として確立している．ただし固体では，励起子を出すために交換相関カーネル $f_{xc}$ が長距離 $\sim-\alpha/q^2$ の振る舞いを持つ必要があり，標準的なLDA/GGAカーネル(短距離)ではRPAとほとんど変わらない結果しか出ない．長距離補正カーネルの開発が続いている．

物理的意味:「バンドギャップが合わない」の3つの層

初学者が混乱しやすいのは，「DFTのバンドギャップが合わない」という言い方が3つの異なる問題を指しうる点である．整理しておこう．

KS固有値差 $\ne$ 準粒子ギャップ．これは近似汎関数の欠陥ではなく，厳密な汎関数を使っても残る原理的な差である(第10章の汎関数微分の不連続)．GWがこれを直す．
近似汎関数の誤差．LDA/GGAの自己相互作用誤差などにより，KSギャップそのものも厳密KSギャップからずれる．ハイブリッド汎関数やDFT+Uが部分的に改善する(第11章)．
光学ギャップ $\ne$ 準粒子ギャップ．光吸収で測るのは「電子と正孔を同時に作る」エネルギーであり，両者が束縛されるぶん(励起子束縛エネルギー)だけ準粒子ギャップより小さい．GWだけでは直らず，BSEが必要である．

したがって，計算した $\varepsilon_2$ の立ち上がりを実験の吸収端と比べるときには，上の3つがすべて絡んでいることを意識しなければならない．偶然の誤差相殺で「LDAの $\varepsilon_2$ が実験に合ってしまう」ことすらある．

18.9 事例研究:理論と実験の協働で構造を決める

ここまでに4種類の分光量の計算法を整備した．最後に，それらを総動員して一つの物質の構造を確定していく実際のワークフローを追いかける．題材は二硼化ジルコニウム表面に自発形成されるシリセンである．この例は，単一の観測量では構造が決まらないこと，そして複数の分光量を同時に説明できるモデルだけが生き残ることを，はっきりと示している．

18.9.1 問題設定:なぜ構造同定が難しいのか

シリセンとは，シリコン原子が蜂の巣格子を組んだ二次元物質である．炭素のグラフェンの類推から，線形分散(Diracコーン)を持つことが期待され，しかもシリコン基盤の半導体技術と親和性が高いことから，大きな注目を集めてきた．

ところが，シリセンはグラフェンとは決定的に違う点を持つ．平坦な構造が安定ではない．理由は結合論から理解できる(第2章)．Si–Si の結合長は C–C より $30\,\%$ 以上長く，そのぶん $p_z$ 軌道どうしの重なりが小さいので $\pi$ 結合が弱い．$\pi$ 結合による安定化が小さければ，$sp^2$ を保つ理由がなくなり，$sp^3$ 的な混成(すなわち原子が面から上下にずれたバックリング構造)のほうが得になる．実際，自立シリセンの最安定構造は $0.44$ Å 程度バックルしている．

この「平坦と非平坦のエネルギー差が小さい」という事情が，構造同定を難しくする．基板の上に載せると，基板との相互作用がこの微妙なバランスを容易にひっくり返し，複数の構造多型が近いエネルギーで共存しうる．しかも，走査トンネル顕微鏡(STM)で見える蜂の巣模様は，どの多型でも似たように見える．「蜂の巣が見えた」だけでは構造は決まらないのである．

例:発見の経緯

この系は，そもそもシリセンを狙って作られたものではない．窒化ガリウム($\mathrm{GaN}$)の成長基板として二硼化ジルコニウム $\mathrm{ZrB_2}$ が注目されたのは，格子定数がよく一致するからである:

$$ \mathrm{GaN}(0001):\ a = 3.189\ \text{Å}, \qquad \mathrm{ZrB_2}(0001):\ a = 3.187\ \text{Å} . $$

$\mathrm{Si}(111)$ 基板上に $\mathrm{ZrB_2}$ の薄膜($200$ Å)を成長させる過程で，下地のシリコンが $\mathrm{ZrB_2}$ 表面まで拡散し，そこで自発的に周期構造を作ることがSTMで見出された．$\mathrm{ZrB_2}$ の $2\times2$ 表面胞に対して $\sqrt{3}\times\sqrt{3}$ の周期を持つ蜂の巣格子——すなわちシリセンである．狙って作ったのではなく，望まぬ拡散として現れたものが新しい二次元物質だった，という点がこの発見の面白さである．

図18.11 第一原理計算による構造同定のワークフロー．実験から候補構造を列挙し，全エネルギーとフォノンで絞り込み，複数の分光量をシミュレートして実験と突き合わせる．決め手になるのは「1つの量が合うこと」ではなく「すべての量が同時に合うこと」である．1つでも矛盾すれば候補は棄却され，ループの最初に戻る．

18.9.2 段階1:候補構造の列挙と全エネルギー比較

最初にすべきことは，実験から分かる制約(周期・組成)を満たす構造を，可能な限り網羅的に作って全エネルギーを比べることである(第16章の構造最適化)．

まず，Si原子1個を $\mathrm{ZrB_2}(0001)$ 表面に置いたときの吸着サイトを調べる．高対称なサイトは3つある．

吸着サイト	結合エネルギー	表面からの高さ
hollow(3回対称の窪み)	$6.08$ eV	$1.97$ Å
bridge(結合の中点)	$5.79$ eV	$2.07$ Å
on-top(Zr原子の真上)	$5.05$ eV	$2.57$ Å

hollowサイトが最安定で，on-topは $1$ eV 以上不利である．この情報は，蜂の巣格子を表面に載せるときの「どのSi原子がどこに落ち着きたがるか」の指針になる．

次に，$\sqrt{3}\times\sqrt{3}$ シリセン $/$ $2\times2$ $\mathrm{ZrB_2}$ という周期の中でSi原子6個の位置を変えながら構造最適化すると，2つの局所安定構造が見つかる．

モデル1(buckled-like):自立シリセンやダイヤモンド構造の $\mathrm{Si}(111)$ 一層に似た，大きくバックルした構造．
モデル2(planar-like):on-topサイトに乗った1個のSi原子だけが高く，残りはほぼ同じ高さに並んだ，平坦に近い構造．

どちらも蜂の巣格子のトポロジーは保っており，STM像だけでは区別しにくい．全エネルギーを比べると

$$ E(\text{モデル2}) - E(\text{モデル1}) = -1.589\ \mathrm{eV/cell} $$

すなわち planar-like のほうが明確に安定である．ここでいったん「答えが出た」ように見えるが，この段階での結論は暫定的である．交換相関汎関数の誤差，基板のモデル化(スラブ厚)，格子定数の取り方——いずれもエネルギー差に効きうる．他の観測量で独立に検証しなければならない．

18.9.3 段階2:フォノンによる動的安定性と，ストライプ構造の起源

全エネルギーの極小は，あくまで「その配置から少し動かすと戻ってくる」ことを保証するだけである．しかも構造最適化は，対称性を仮定したまま行われることが多い．仮定した対称性を破る方向に不安定でないかを確かめるには，フォノン分散を計算しなければならない(第16章のHessian行列を，周期系のすべての波数について計算することに相当する)．調和近似から動的行列 $D_{\alpha\beta}^{\kappa\kappa'}(\bm{q})$ を導き，その固有値 $\omega^2(\bm{q})$ が負になることが構造の不安定性を意味することを証明したのが第16章16.4.5項である．以下ではそこで確立した判定基準をそのまま用いる．

planar-like 構造のフォノン分散を計算すると，ブリルアン域のM点に虚の振動数(ソフトモード)が現れる．虚の振動数とは，その変位方向にエネルギーが下がる——すなわち構造が不安定である——ことを意味する(16.4.5項(7))．その固有ベクトルを調べると，on-topに乗った高いSi原子の列に垂直な方向の変位であった．

これは実験と符合する．STMでは，シリセン膜が自発的にストライプ状のドメインを形成することが観測されていた．均一な $\sqrt{3}\times\sqrt{3}$ 構造がM点で不安定なら，その不安定性を解消するために周期構造が変調される——というのが自然な解釈である．

導出:ストライプ幅の下限

ストライプ構造ができると，その周期に対応する新しい逆格子ベクトルが生まれ，ブリルアン域が縮む．ソフトモードを解消するには，この新しい逆格子が十分細かくなければならない．定量化しよう．

ステップ1:ソフトモードの波数．六方格子(格子定数 $a$)のM点は，$\Gamma$ から

$$ q_M = \frac{2\pi}{\sqrt{3}\,a} $$

の距離にある(M点は隣接する逆格子点を結ぶ線分の中点であり，六方逆格子の最短逆格子ベクトルの長さ $4\pi/(\sqrt3 a)$ の半分)．

ステップ2:ストライプの逆格子ベクトル．ストライプは，原子列の並ぶ方向に垂直に幅 $N$ 列ごとに繰り返す．列間隔は $\frac{\sqrt3}{2}a$ だからストライプの実空間周期は $N\cdot\frac{\sqrt3}{2}a$ であり，対応する逆格子ベクトルの大きさは

$$ q' \cdot \frac{\sqrt{3}}{2}a\cdot N = 2\pi \qquad\Longrightarrow\qquad q' = \frac{4\pi}{\sqrt{3}\,a\,N} . $$

ステップ3:条件を課す．ソフトモードが縮んだブリルアン域の中で $\Gamma$ 点(振動数ゼロのモード)に折り畳まれてしまわないためには，新しい逆格子の刻みがソフトモードの波数の半分より細かくなければならない:

$$ q' \lt \frac{q_M}{2} . $$

ステップ4:$N$ について解く．

\begin{align} \frac{4\pi}{\sqrt{3}\,a\,N} &\lt \frac{1}{2}\cdot\frac{2\pi}{\sqrt{3}\,a} \label{eq:18-stripe1}\\[2pt] \frac{4}{N} &\lt 1 \label{eq:18-stripe2}\\[2pt] N &\gt 4 . \label{eq:18-stripe3} \end{align}

\eqref{eq:18-stripe2}では両辺に $\sqrt3 a/\pi$ を掛けて共通因子を消した．すなわちストライプの幅は5列以上でなければならない．∎

実際に大規模な第一原理計算で，ストライプ幅 $N$ を変えながら終端Zr原子あたりのエネルギーを求めると，$N=5$ で最小になる．この値は，STMで観測されたストライプ周期と一致する．フォノン計算という「エネルギーの2階微分」から出発して，実験で見えている表面のドメイン構造の周期まで説明できたことになる．

18.9.4 段階3:XPS内殻シフトによる決定的な判別

ここで18.3〜18.5節の道具が効く．高分解能XPSでSi $2p$ を測ると，実験スペクトルは3つの成分 $\alpha, \beta, \gamma$ に分解された(面積比はおよそ $2.5:3.7:1$)．$\sqrt{3}\times\sqrt{3}$ 胞にはSi原子が非等価な3サイト存在するので，この3成分は3サイトに対応するはずである．

そこで，2つの候補構造それぞれについて，3サイトのSi $2p_{1/2}$, $2p_{3/2}$ の絶対束縛エネルギーを式\eqref{eq:18-eb-general}で計算し，実験と同じブロードニング関数をかけてスペクトルを合成する．結果は明快であった．

buckled-like 構造:3ピークのうち $\gamma$ の相対位置を再現できない．
planar-like 構造:3ピークの相対位置も絶対値も，定量的に実験と一致する．

ここで絶対値まで比較できたことが決め手になった．相対シフトだけの比較では，基準の取り方に自由度が残り，「どちらの構造でも合わせられる」余地がある．絶対値を第一原理で出せるからこそ，候補を一つに絞れたのである．実際，Si $2p$ の Kohn–Sham 固有値から作った状態密度をそのまま実験と比べても，絶対値が数十 eV ずれてしまい，判別の役に立たない．これは18.3.3節で見た「Koopmans型評価の誤差は $\varepsilon_c'$ に比例する」という結論の，実戦での現れである．

18.9.5 段階4:ARPESとアンフォールディング

4つ目の検証がバンド構造である．ARPESで測ったバンド分散と，計算したバンドを比べる．ここで18.7節のアンフォールディングが必須になる．というのも，計算に使うセルは $2\times2$ の $\mathrm{ZrB_2}$ スラブ(しかも数層)であり，そのままバンドを描いてもバンドが折り畳まれてARPESと比較できないからである．

$\mathrm{ZrB_2}$ のバルク単位胞を参照セルとしてアンフォールドしたスペクトル関数を計算すると，planar-like 構造について，表面状態 $S_1$, $S_2$(Zr $d$ 軌道由来)とシリセン由来のバンド $X_2$, $X_3$ の分散が，いずれもARPESの観測と一致する．軌道分解した重み $W^{\kk}_{\bm{K}JM}$ を使えば，どのバンドが基板由来でどれがシリセン由来かも同時に示せる．

補足:STM像のシミュレーション(Tersoff–Hamann近似)

5つ目の観測量としてSTM像も比較できる．トンネル電流を第一原理で計算する最も簡便な処方がTersoff–Hamann近似である．探針の先端を $s$ 波軌道で近似すると，Bardeenのトンネル行列要素が探針位置 $\bm{r}_t$ における試料波動関数の値に比例することが示され，バイアス電圧 $V$ でのトンネル電流は

$$ \begin{equation} I(\bm{r}_t, V) \;\propto\; \int_{E_F}^{E_F+eV}\rho(\bm{r}_t, E)\,\dd E, \qquad \rho(\bm{r}_t,E) = \sum_n\abs{\psi_n(\bm{r}_t)}^2\delta(E-\varepsilon_n) \label{eq:18-tersoff} \end{equation} $$

となる．すなわち定電流STM像は，エネルギー積分した局所状態密度の等値面である．この事実には2つの含意がある．第一に，STM像は原子位置そのものではなく電子密度を見ている——「明るい点」が必ずしも原子の位置とは限らない．第二に，バイアス電圧を変えれば見えるエネルギー窓が変わるので，正負のバイアスで像が反転することさえある．したがって，STM像の比較は必ず実験と同じバイアス条件で行わなければならない．

18.9.6 方法論の要点

最終的に，planar-like 構造は次の4つの独立な実験すべてと整合した．

高分解能XPS(Si $2p$ の3成分の絶対束縛エネルギーと相対位置)
ARPES(アンフォールドしたバンド分散)
STMで観測されたストライプドメインの周期(フォノンのソフトモードから予言された $N\gt4$，および $N=5$ でのエネルギー最小)
全反射高速陽電子回折による表面原子配置の解析

物理的意味:なぜ「複数の観測量」でなければならないのか

この事例から取り出すべき教訓は，具体的な構造そのものではなく，方法論である．

(1) 単一の観測量は構造を決めない．STMの蜂の巣模様は両方の候補で再現される．全エネルギー差 $1.589$ eV/cell は大きく見えるが，汎関数の系統誤差や基板モデルの不完全さを考えれば，それだけで断定はできない．XPSの3ピークも，相対位置だけなら「合わせにいく」余地がある．1つの量に対して1つのモデルを合わせるのは，常に可能に近い．

(2) 独立な観測量を重ねると自由度が枯渇する．XPSは内殻の静電環境，ARPESは占有バンドの分散，フォノンは力の定数，回折は原子位置——それぞれ電子構造の異なる側面を測っている．1つのモデルがこれらを同時に説明する確率は，偶然では極めて低い．だから「すべてが合った」ときには，その構造が正しいと強く信じられる．

(3) 絶対値で比較できることの価値．相対量(化学シフト，バンド分散の形)だけの比較では，基準の自由度が残る．本章で内殻束縛エネルギーの絶対値を計算する技術に多くの紙面を割いたのは，それが「合わせられる自由度」を一つ潰すからである．理論の価値は，フィッティングパラメータを減らすところにある．

(4) 計算は「反証」に使うのが最も強い．この事例で決定的だったのは，planar-like が合ったことより，buckled-like が $\gamma$ ピークを再現できなかったことである．候補を積極的に殺せる観測量を探すこと——それが理論と実験の協働における理論側の最大の仕事である．

18.10 まとめ・展望・演習・参考文献

18.10.1 本章のまとめ

4種類の電子遷移．始状態(内殻/価電子帯)と終状態(真空/物質内の非占有状態)の組合せで，XPS・XANES・ARPES・UV-Vis-NIRが整理される(図18.1)．XPSだけが全エネルギー差を要求し，他は一電子的な量(固有値と行列要素)で書ける．内殻に空孔を作ると有効ポテンシャルが数十 eV も変わるからである．
XPSのエネルギー保存則．試料と測定器がオーミック接触して共通の化学ポテンシャル $\mu$ を持つとき，測定量 $E_B = h\nu - K_{\mathrm{spec}} - \phi_{\mathrm{spec}}$ には測定器の仕事関数だけが現れる．基準の任意シフト $\Delta\mu$ は係数 $(N-1)-N+1=0$ で厳密に相殺し， $$E_B = E_f^{(0)}(N-1) - E_i^{(0)}(N) + \mu_0$$ という，金属・半導体・絶縁体・分子に共通な絶対束縛エネルギーの式\eqref{eq:18-eb-general}が得られる．$E_f^{(0)}$, $E_i^{(0)}$, $\mu_0$ を同じ電位基準で計算することが決定的に重要である．
金属の簡略式．金属ではフェルミ準位の固有値が電子数によらないので，Janakの定理を積分して $E_f^{(0)}(N-1)-E_f^{(0)}(N) = -\mu_0$ が示され，$\mu_0$ が相殺して $E_B = E_f^{(0)}(N)-E_i^{(0)}(N)$ となる(式\eqref{eq:18-eb-metal})．両方の計算が中性セルで済むので収束が速い．有限ギャップ系にこの式を流用すると $1$ eV 規模の系統誤差が残る．
Janak積分と遷移状態．$E_B = -\int_0^1\varepsilon_c(n)\dd n$(式\eqref{eq:18-eb-integral})．被積分関数を $n=1/2$ のまわりで展開して項別積分すると，奇数次が消えて中点則 $\int_0^1\varepsilon_c\dd n = \varepsilon_c(1/2)+\frac{1}{24}\varepsilon_c''+\cdots$ が出る．ゆえに $E_B\simeq-\varepsilon_c(1/2)$(Slaterの遷移状態)．誤差はKoopmans型が $\frac12\varepsilon_c'$(1階微分!)，台形則が $\frac{1}{12}\varepsilon_c''$，中点則が $\frac{1}{24}\varepsilon_c''$ である．内殻では $\varepsilon_c'$ が数十 eV なのでKoopmans型は使えない．
始状態効果と終状態効果．$E_B = -\varepsilon_c(1) + R$，$R = \int_0^1[\varepsilon_c(1)-\varepsilon_c(n)]\dd n \ge 0$(式\eqref{eq:18-eb-init-final})．化学シフトはこの2項の差として分解され，さらに電荷ポテンシャル模型 $E_B = E_B^0 + kq_A + \sum_{B}q_B/R_{AB} - R_A$(式\eqref{eq:18-potmodel-key})で「自分の電荷 $+$ Madelung $-$ 緩和」に翻訳できる．観測される傾きが $k\simeq17$ eV/電子より小さいのは，実効電荷移動が形式酸化数より小さいことと，Madelung相殺による．
コアホール計算の2つの技術．(i) ペナルティ汎関数 $\hat{P}=\ket{\chi_c}\Delta\bra{\chi_c}$ を加えると，修正KS方程式 $(\hat{h}_{\mathrm{KS}}+\hat{P})\psi=\varepsilon\psi$ により内殻準位が $\Delta$ だけ押し上げられ，空孔が狙った原子に局在したまま自己無撞着に維持される．自己無撞着解では $E_{\mathrm{pen}}\to0$ なのでエネルギーを歪めない．(ii) 打ち切りCoulomb核 $\tilde{v}_c(\bm{G}) = (4\pi/G^2)(1-\cos GR_c)$ は $\bm{G}\to\bm{0}$ で $2\pi R_c^2$ と有限になるので，荷電セルを補償背景なしに厳密に扱える．条件は $R_c=2R$，$L\gt4R$．補償背景では $\alpha_MQ^2/(2L)$ の誤差($L=20$ Bohr で $1.93$ eV)が $1/L$ でしか消えない．
精度．固体8例で平均絶対誤差 $0.4$ eV(相対 $0.16\,\%$)，気相分子23例で $0.5$ eV(相対 $0.22\,\%$)．Si $2p$ のスピン軌道分裂 $0.73$ eV も再現される．一方，$\mathrm{O_2}$ の $1s$ 交換分裂はGGAで半分しか出ず，汎関数の限界が露呈する場合もある．
励起状態DFTの正当化．$[\hat{H},\hat{O}]=0$ なる観測量で状態をクラス分けすると，各クラスの最低状態には変分原理が成り立つ(Gunnarsson–Lundqvist)．Levyの制約付き探索をクラス $A$ に制限すれば $E^{(A)}_{\min}=\min_\rho\{F_A[\rho]+\int v\rho\}$ と密度汎関数形式が組める．ただし $F_A$ はクラス依存で，基底状態用の汎関数を流用するのは制御されていない近似である．コアホールには $\hat{A}=1-\hat{n}_c$ を取ればよく，これがペナルティ汎関数法の理論的裏付けになる．
XANESの公式と選択則．黄金律と双極子近似から $\sigma(\omega)=4\pi^2\alpha\omega\sum_f\abs{\braket{f|\hat{\bm{e}}\cdot\rr|c}}^2\delta(\varepsilon_f-\varepsilon_c-\omega)$．角度積分のパリティ($l-l_c$ が奇数)と三角不等式($\abs{l-l_c}\le1$)から $\Delta l=\pm1$ が出る．内殻軌道が局在しているため，断面積はコアホール存在下の局所 $l\pm1$ 部分状態密度に比例する(終状態則)．K端は $p$，L$_{2,3}$端は $d$ を映す．偏光依存性から吸着分子の配向が決まる．
ARPESとアンフォールディング．$k_\parallel = \sqrt{2mE_{\mathrm{kin}}}\sin\theta/\hbar$ で面内波数が決まる($0.5123\sqrt{E[\mathrm{eV}]}\sin\theta$ Å$^{-1}$)．$k_\perp$ は保存しないので内部ポテンシャル $V_0$ を仮定する．スーパーセル計算ではバンドが折り畳まれるが，基本並進群の既約表現への射影 $\hat{P}_{\kk}=N^{-1}\sum_{\bm{n}}e^{i\kk\cdot\bm{R}_{\bm{n}}}\hat{T}_{\bm{R}_{\bm{n}}}$ を使えば，スペクトル重み $W_{\bm{K}J}(\kk)=\sum_{\bm{G}_p}\abs{\braket{\kk+\bm{G}_p|\Psi_{\bm{K}J}}}^2$($\sum_{\kk}W=1$)によって基本セルのブリルアン域に戻せる．摂動があると $W$ が $\kk$ 方向に分散し，これがARPESの「ぼやけたバンド」に対応する．
誘電関数．黄金律で計算した吸収パワーと，マクロな散逸 $\braket{\bm{j}\cdot\bm{E}}=(\omega/8\pi)\varepsilon_2\abs{\bm{E}_0}^2$ を等置すると，独立粒子近似の $$\varepsilon_2(\omega)=\frac{4\pi^2e^2}{m_e^2\omega^2}\frac{2}{\Omega}\sum_{\kk,v,c}\abs{\hat{\bm{e}}\cdot\bm{p}_{cv}}^2\delta(\varepsilon_{c\kk}-\varepsilon_{v\kk}-\hbar\omega)$$ が得られる(式\eqref{eq:18-eps2-key})．ピークは結合状態密度のvan Hove臨界点($\nabla_{\kk}(\varepsilon_c-\varepsilon_v)=0$)に対応する．$\sigma_1=\omega\varepsilon_2/4\pi$ がKubo–Greenwood形式との橋渡しである．
Kramers–Kronig．因果律 $\chi(\tau\lt0)=0$ $\Rightarrow$ $\tilde\chi(\omega)$ が上半平面で解析的 $\Rightarrow$ Cauchyの積分公式(極を小半円で避ける)$\Rightarrow$ $\varepsilon_1(\omega)=1+\frac{2}{\pi}\mathcal{P}\int_0^\infty\omega'\varepsilon_2(\omega')/(\omega'^2-\omega^2)\dd\omega'$．静的誘電率 $\varepsilon_1(0)=1+\frac{2}{\pi}\int_0^\infty\varepsilon_2/\omega'\dd\omega'$ とf-総和則が実装のチェックに使える．得られた $\varepsilon$ から $n,\kappa,\alpha,R$，そして損失関数 $\mathrm{Im}[-1/\varepsilon]$(XPSの平均自由行程に戻る量)がすべて導かれる．
限界とその先．独立粒子近似は準粒子補正・励起子効果・局所場効果を含まない．処方はコストの低い順にシザーズ補正，GW，BSE，TDDFTである．「バンドギャップが合わない」は，(i) KS固有値差 $\ne$ 準粒子ギャップ，(ii) 近似汎関数の誤差，(iii) 光学ギャップ $\ne$ 準粒子ギャップ，の3層に分けて考えなければならない．
方法論．単一の観測量では構造は決まらない．独立な複数の分光量(XPS・ARPES・XANES・フォノン・回折・STM)を同時に説明できるモデルだけが生き残る．理論の仕事は，フィッティングの自由度を潰すこと(絶対値を出すこと)と，候補を積極的に反証できる観測量を見つけることである．

18.10.2 本書全体を振り返って

ここまでの道のりを，一つの物語として振り返っておきたい．

本書は「多電子のSchrödinger方程式は解けない」という絶望から始まった(第1章)．$N$ 個の電子の波動関数は $3N$ 次元の関数であり，$N=100$ でもう宇宙のどんな計算機にも収まらない．にもかかわらず，化学結合の起源はビリアル定理という驚くほど一般的な関係で捉えられ(第2章)，金属中の電子の集団的な振る舞いは自由電子ガスという単純な模型でよく記述できた(第3章)．

第II部では，この絶望に対する最初の突破口を追った．Thomas–Fermi模型は「エネルギーを電子密度の汎関数で書く」という発想を初めて形にし(第4章)，Hartree–Fock法は反対称性から交換相互作用という新しい物理を引き出した(第5章)．第二量子化という言語を手にし(第6章)，ジェリウム模型で交換・相関エネルギーを最後まで計算し(第7章)，線形応答と遮蔽で電子系の集団的な応答を理解した(第8章)．ここまでで，密度汎関数という考え方の「材料」がすべて揃った．

第III部が核心である．Hohenberg–Kohnの定理は，基底状態密度が外部ポテンシャルを——したがって系のすべてを——一意に決めることを示した(第9章)．$3N$ 次元の波動関数ではなく，$3$ 次元の密度で十分だという，驚くべき主張である．Kohn–Sham方程式は，この主張を実際に解ける形に翻訳した(第10章)．相互作用しない参照系という発想が，運動エネルギーの大部分を正確に扱う道を開いた．そして交換相関汎関数の構成(第11章)によって，理論は実用的な計算手法になった．

第IV部・第V部では，この理論を現実の物質と計算機に接続した．Blochの定理が周期系の扱いを可能にし(第12章)，Green関数と状態密度が電子構造を局所的に読み解く言語を与え(第13章)，基底関数・SCF・電荷混合という数値の技術(第14章)と擬ポテンシャル(第15章)が，実際に手を動かすための道具を提供した．

第VI部が本書の到達点である．力と応力によって原子が動くようになり(第16章)，非平衡Green関数によって電流が流れるようになり(第17章)，そして本章で，計算結果が実験と直接比較できる数値になった．理論は，閉じた体系から，外界と対話する道具になったのである．

物理的意味:本書を貫く3つのモチーフ

振り返ると，本書には繰り返し現れた3つのモチーフがある．

(1) 変分原理．Hartree–Fock方程式，Kohn–Sham方程式，Hellmann–Feynmanの定理，Pulay力，応力テンソル，そして本章のGunnarsson–Lundqvistの議論とペナルティ汎関数——すべてが「あるエネルギー汎関数を停留させる」という一つの原理から出ている．変分原理は，近似を作る指針でもあり，その近似が持つ誤差の次数を保証する道具でもある．

(2) 応答と相関．誘電関数(第8章・第18章)，Green関数と状態密度(第13章)，Kubo公式，揺動散逸定理——「摂動に対する応答」と「平衡状態の揺らぎ」が同じものであるという構造は，遮蔽・輸送・分光のすべてに現れた．

(3) 微分量と差分量の区別．Kohn–Sham固有値は微分量($\partial E/\partial n$)であって，電子1個分のエネルギーではない．この区別を曖昧にしたまま固有値を実験と比べると，本章18.3.3節で見たように数十 eV も外す．逆に，Janakの定理を積分すれば正しい差分量が得られる．物理量が微分なのか差分なのかを常に問うこと——これは本書を通じて最も実用的な教訓かもしれない．

18.10.3 その先の地図

本書が扱ったのは，基底状態密度汎関数理論とその直接の応用である．ここから先には，いくつもの発展方向が広がっている．読者が次に何を学ぶべきかの地図として，簡単に整理しておく．

方向	解こうとしている問題	鍵となる考え方	本書のどこから接続するか
GW近似	準粒子エネルギー(バンドギャップ，光電子スペクトルの絶対位置)	自己エネルギー $\Sigma=iGW$．遮蔽Coulomb相互作用 $W$ を使う	第13章のGreen関数，第8章の遮蔽，本章18.8.6節
Bethe–Salpeter方程式	光学スペクトル，励起子	電子-正孔対の2体方程式．$\varepsilon_2$ のピーク強度の再分配	本章18.8節の $\varepsilon_2$
時間依存DFT	励起エネルギー，実時間の電子ダイナミクス	時間依存KS方程式(Runge–Grossの定理)，交換相関カーネル $f_{xc}$	第10章のKS方程式，本章18.6節の励起状態
動的平均場理論(DMFT)	強相関電子系(Mott絶縁体，遷移金属酸化物)	局所自己エネルギー近似．不純物問題への写像	第11章のDFT+U，第13章のGreen関数
量子モンテカルロ	相関エネルギーのベンチマーク精度	試行波動関数の確率的サンプリング，拡散モンテカルロ	第7章の相関エネルギー(LDAの基礎データはQMC由来)
機械学習ポテンシャル	第一原理精度で $10^5$–$10^9$ 原子・ナノ秒以上のMD	局所環境の記述子とニューラルネット/回帰．DFTデータで学習	第16章の第一原理MD
電子-フォノン結合	超伝導，輸送の温度依存性，バンドの繰り込み	Eliashberg関数 $\alpha^2F(\omega)$，超伝導DFT	第16章の力とHessian行列，第17章の輸送
Berry位相とトポロジー	電気分極，トポロジカル絶縁体，異常Hall効果	Bloch波動関数の位相構造，Wannier中心	第12章のBloch定理，第14章の局在基底

物理的意味:第一原理計算は何のためにあるのか

最後に一言だけ，この分野の意味について述べておきたい．第一原理計算の価値は「実験を計算で置き換えること」ではない．実験には常に，計算では手の届かない現実の複雑さがある．逆に計算には，実験では決して測れない情報——波動関数そのもの，仮想的な構造のエネルギー，起こらなかった反応経路——がある．

本章18.9節の事例が示したように，両者が本当に力を発揮するのは互いに問いを立て合うときである．実験が「この蜂の巣は何か」と問い，計算が「2つの候補がある，XPSの $\gamma$ ピークを測れば区別できる」と答え，実験がそれを測り，計算が予言と照合する．この往復のなかでしか，物質の姿は確定しない．

そして，その往復を成立させるのは定量性である．数百 eV の量を $0.5$ eV で当てられるからこそ，候補を殺すことができる．本書が数式の変形を一行も飛ばさずに書かれてきたのは，この定量性がどこから来るのかを——どの近似がどの精度を支えているのかを——読者自身が判断できるようになるためである．式の中身を知らずに数値を出すことは誰にでもできる．式の中身を知って初めて，その数値をどこまで信じてよいかが分かる．

18.10.4 演習問題

演習18.1 内殻束縛エネルギーの定式化

(a) 18.3.1節の導出を，終状態が $N-2$ 電子系(たとえばAuger過程の終状態のように，電子が2個抜けた状態)である場合に一般化せよ．すなわち，基準シフト $\Delta\mu$ を含めて計算し，$\Delta\mu$ が相殺する条件と，最終的な式に $\mu_0$ が何倍で現れるかを求めよ．

(b) Al K$\alpha$ 線($h\nu=1486.7$ eV)で炭素($M = 12\times1836\,m_e$)の $1s$($E_B=285$ eV)を測るときの反跳エネルギーを見積もれ．また，同じ条件で水素原子核が反跳する場合と比べよ．$0.1$ eV の精度目標に対してどちらが問題になるか．

(c) $\varepsilon_c(n)$ が $n$ の1次関数 $\varepsilon_c(n)=an+b$ であるとき，式\eqref{eq:18-eb-integral}を直接積分して $E_B$ を求め，Slaterの遷移状態近似 $-\varepsilon_c(1/2)$ が厳密であることを示せ．また，このとき台形則 $-\frac12[\varepsilon_c(0)+\varepsilon_c(1)]$ も厳密になることを示し，なぜ2階微分の誤差項が消えるのかを説明せよ．

(d) 炭素 $1s$ で $\varepsilon_c(1)=-270$ eV，$E_B=285$ eV とする．緩和エネルギー $R$ を式\eqref{eq:18-eb-init-final}から求めよ．また，$\varepsilon_c(n)$ を1次関数と仮定して $\varepsilon_c(0)$ と $\varepsilon_c(1/2)$ を推定せよ．

ヒント: (a) 電子が2個抜けるので $E_f(N-2)=E_f^{(0)}(N-2)+(N-2)\Delta\mu$ とし，測定される運動エネルギーが2個分であることに注意して保存則を書き直す．(b) $E_{\mathrm{rec}}\simeq K\,m_e/M$．(c) 1次関数の積分は台形則でも中点則でも厳密である($\varepsilon_c''=0$)．(d) $R=E_B+\varepsilon_c(1)$，1次関数なら $\varepsilon_c(1/2)=-E_B$．

演習18.2 荷電セルの静電問題と厳密クーロンカットオフ

(a) 単純立方セル($\alpha_M=2.8373$)，$Q=1$ について，Makov–Payne展開\eqref{eq:18-makov-payne}の主要項 $\alpha_MQ^2/(2L)$ を $L=15,20,30,60$ Bohr について eV 単位で計算し，$1/L$ で減衰することを確認せよ．誤差を $0.05$ eV 以下にするのに必要な $L$ は何 Bohr か．またそれは何 Å か．

(b) 打ち切りCoulomb核\eqref{eq:18-vc-G}を $G\to0$ で4次まで展開し，$\tilde{v}_c(\bm{G}) = 2\pi R_c^2 - (\pi/6)R_c^4G^2+\cdots$ となることを示せ．$\cos x$ の展開を $x^4$ の項まで使うこと．

(c) シリコン(ダイヤモンド構造，格子定数 $a=5.43$ Å，単位胞あたり8原子)で $\Delta\rho$ の有効半径が $R=7$ Å のとき，条件\eqref{eq:18-cutoff-cond}を満たすのに必要な立方スーパーセルの一辺は何 Å か．また，そのセルに含まれる原子数はおよそいくつか．

(d) 金属用の式\eqref{eq:18-eb-metal}を有限ギャップ系に適用したときに生じる系統誤差を，物理的に議論せよ．とくに，「抜いた電子をどこに置くか」という観点から，誤差の符号(過大評価か過小評価か)を予想し，図18.5(b)の振る舞いと整合するか確かめよ．

ヒント: (a) $1$ Hartree $=27.211$ eV．(b) $1-\cos x = x^2/2 - x^4/24+\cdots$ を $4\pi/G^2$ に掛ける．(c) $L\gt4R$．原子数は $(L/a)^3\times8$．(d) 有限ギャップ系では，抜いた電子を伝導帯の底に置くことになり，$\mu_0$ との差だけ余分なエネルギーを支払う．

演習18.3 双極子選択則とXANES

(a) K端($l_c=0$)から $d$ 状態($l=2$)への双極子遷移が禁制であることを，18.6.4節のパリティ議論と三角不等式の両方から独立に示せ．実際に $3d$ 遷移金属のK端には弱い前縁ピークが観測されることがあるが，その起源として考えられる機構を2つ挙げよ．

(b) L$_{2,3}$端($l_c=1$)では $l=0$ と $l=2$ の両方が許される．それにもかかわらず $d$ 成分が圧倒的に強い理由を，動径積分 $\int R_l(r)\,r\,R_{2p}(r)r^2\dd r$ の観点から説明せよ．

(c) 表面に平らに寝た芳香環分子の $\pi^*$ 軌道は表面法線方向($z$)を向いている．直線偏光の電場ベクトルが法線から角度 $\theta$ をなすとき，$1s\to\pi^*$ 遷移の強度が $\cos^2\theta$ に比例することを示せ．$\theta=0^\circ, 45^\circ, 90^\circ$ での相対強度を求めよ．

(d) ジルコニウムのK端($18.0$ keV)における双極子近似の妥当性を，18.6.4節と同じ方法で評価せよ．Zr $1s$ 軌道の半径を $a_0/Z = 0.529/40$ Å と見積もってよい．四重極遷移が観測にかかるかどうかを論ぜよ．

ヒント: (a) $l-l_c=2$ は偶数．前縁ピークの起源は，四重極遷移と，配位子との混成による $p$ 成分の混入．(b) $2p$ 軌道と重なりの大きい終状態はどちらか．(c) $\hat{\bm{e}}\cdot\rr$ のうち $p_z$ 終状態と結合するのは $z$ 成分だけである．(d) $\lambda[\text{Å}]=12398/E[\mathrm{eV}]$，$q=2\pi/\lambda$．

演習18.4 アンフォールディングと誘電関数

(a) 一次元鎖(格子定数 $a$)を2倍にしたスーパーセル($N=2$)を考える．射影演算子\eqref{eq:18-projector-k}を具体的に $\hat{P}_{\kk}=\frac12\bigl(\hat{1}+e^{i\kk a}\hat{T}_{a}\bigr)$ と書き下し，スーパーセル固有状態が2つの平面波 $\ket{K}$ と $\ket{K+\pi/a}$ の重ね合わせ $c_1\ket{K}+c_2\ket{K+\pi/a}$ であるとき，$\kk=K$ と $\kk=K+\pi/a$ に対するスペクトル重みをそれぞれ求めよ．和則 $\sum_{\kk}W=1$ を確認せよ．

(b) (a)で摂動がまったくない場合($c_1=1, c_2=0$)と，摂動が強く2つの平面波が等しく混ざる場合($c_1=c_2=1/\sqrt2$)について，アンフォールドしたバンド図がどう見えるかを述べよ．実験のARPESでこの違いはどのように現れるか．

(c) 減衰のないLorentz振動子の吸収 $\varepsilon_2(\omega)=\frac{\pi}{2}\frac{\omega_p^2}{\omega_0}\delta(\omega-\omega_0)$($\omega\gt0$)を考える．Kramers–Kronig関係\eqref{eq:18-kk-key}を使って $\varepsilon_1(\omega)$ を求め，古典的なLorentz模型の結果 $\varepsilon_1(\omega)=1+\omega_p^2/(\omega_0^2-\omega^2)$ が再現されることを示せ．またf-総和則 $\int_0^\infty\omega\varepsilon_2\dd\omega=\frac{\pi}{2}\omega_p^2$ が成り立つことを確認せよ．

(d) (c)の模型で $\omega_0=3.4$ eV，$\varepsilon_1(0)=11.9$(シリコンの値)とする．$\omega_p$ を求めよ．次に，シザーズ補正で $\omega_0$ を $0.6$ eV だけ上へずらしたとき($\omega_p$ は不変)，$\varepsilon_1(0)$ が何%変化するかを計算せよ．この結果は，静的誘電率の計算にバンドギャップ誤差がどう効くかについて何を教えるか．

ヒント: (a) $\hat{T}_a\ket{q}=e^{-iqa}\ket{q}$．$e^{i(K-K)a}=1$，$e^{i(K-K-\pi/a)a}=e^{-i\pi}=-1$．(c) $\delta$ 関数があるので積分は一瞬で終わる．(d) $\varepsilon_1(0)=1+\omega_p^2/\omega_0^2$ から $\omega_p$ を出し，$\omega_0\to4.0$ eV で再計算する．

18.10.5 参考文献

光電子分光の基礎

K. Siegbahn et al., ESCA: Atomic, Molecular and Solid State Structure Studied by Means of Electron Spectroscopy (Almqvist & Wiksells, Uppsala, 1967). ——化学シフトを分析手法にまで高めた古典．
S. Hüfner, Photoelectron Spectroscopy: Principles and Applications, 3rd ed. (Springer, 2003). ——XPS/UPSの標準教科書．
M. P. Seah and W. A. Dench, "Quantitative electron spectroscopy of surfaces," Surf. Interface Anal. 1, 2 (1979). ——非弾性平均自由行程のユニバーサル曲線．
C. S. Fadley, "X-ray photoelectron spectroscopy: Progress and perspectives," J. Electron Spectrosc. Relat. Phenom. 178–179, 2 (2010).
U. Gelius, "Binding energies and chemical shifts in ESCA," Phys. Scr. 9, 133 (1974). ——電荷ポテンシャル模型．

内殻束縛エネルギーの第一原理計算

J. F. Janak, "Proof that $\partial E/\partial n_i=\varepsilon_i$ in density-functional theory," Phys. Rev. B 18, 7165 (1978).
J. C. Slater, "Statistical exchange-correlation in the self-consistent field," Adv. Quantum Chem. 6, 1 (1972). ——遷移状態近似．
T. Ozaki and C.-C. Lee, "Absolute binding energies of core levels in solids from first principles," Phys. Rev. Lett. 118, 026401 (2017). ——本章18.3–18.5節の定式化・ペナルティ汎関数法・ベンチマークの原論文．
G. Makov and M. C. Payne, "Periodic boundary conditions in ab initio calculations," Phys. Rev. B 51, 4014 (1995). ——荷電セルの有限サイズ補正．
M. R. Jarvis, I. D. White, R. W. Godby, and M. C. Payne, "Supercell technique for total-energy calculations of finite charged and polar systems," Phys. Rev. B 56, 14972 (1997). ——厳密クーロンカットオフ法．
L. Köhler and G. Kresse, "Density functional study of CO on Rh(111)," Phys. Rev. B 70, 165405 (2004). ——コアホール擬ポテンシャルによる内殻シフト計算の実例．

励起状態のDFTとX線吸収

O. Gunnarsson and B. I. Lundqvist, "Exchange and correlation in atoms, molecules, and solids by the spin-density-functional formalism," Phys. Rev. B 13, 4274 (1976). ——対称性クラスごとの変分原理．18.6.1節の出発点．
M. Levy, "Universal variational functionals of electron densities…," Proc. Natl. Acad. Sci. USA 76, 6062 (1979). ——制約付き探索(第9章も参照)．
U. von Barth and G. Grossmann, "Dynamical effects in x-ray spectra and the final-state rule," Phys. Rev. B 25, 5150 (1982).
J. Stöhr, NEXAFS Spectroscopy (Springer, 1992). ——選択則・偏光依存性・分子配向解析の標準的参考書．
J. J. Rehr and R. C. Albers, "Theoretical approaches to x-ray absorption fine structure," Rev. Mod. Phys. 72, 621 (2000). ——多重散乱理論を含む総説．
M. Taillefumier, D. Cabaret, A.-M. Flank, and F. Mauri, "X-ray absorption near-edge structure calculations with the pseudopotentials," Phys. Rev. B 66, 195107 (2002).
N. A. Besley, A. T. B. Gilbert, and P. M. W. Gill, "Self-consistent-field calculations of core excited states," J. Chem. Phys. 130, 124308 (2009). ——分子の内殻励起に対する $\Delta$SCF 的手法．

ARPESとバンドアンフォールディング

A. Damascelli, Z. Hussain, and Z.-X. Shen, "Angle-resolved photoemission studies of the cuprate superconductors," Rev. Mod. Phys. 75, 473 (2003). ——ARPESの運動学・行列要素効果・スペクトル関数の標準的総説．
W. Ku, T. Berlijn, and C.-C. Lee, "Unfolding first-principles band structures," Phys. Rev. Lett. 104, 216401 (2010). ——アンフォールディングの原論文．
V. Popescu and A. Zunger, "Extracting E versus k effective band structure from supercell calculations on alloys and impurities," Phys. Rev. B 85, 085201 (2012).
C.-C. Lee, Y. Yamada-Takamura, and T. Ozaki, "Unfolding method for first-principles LCAO electronic structure calculations," J. Phys.: Condens. Matter 25, 345501 (2013). ——非直交局在基底での定式化(式\eqref{eq:18-weight-lcao})．

誘電関数・光学応答

R. Kubo, "Statistical-mechanical theory of irreversible processes I," J. Phys. Soc. Jpn. 12, 570 (1957). ——線形応答理論．
D. A. Greenwood, "The Boltzmann equation in the theory of electrical conduction in metals," Proc. Phys. Soc. 71, 585 (1958).
H. Ehrenreich and M. H. Cohen, "Self-consistent field approach to the many-electron problem," Phys. Rev. 115, 786 (1959). ——RPA誘電関数(式\eqref{eq:18-eps2-key})．
F. Wooten, Optical Properties of Solids (Academic Press, 1972). ——Kramers–Kronig関係と光学定数の導出が丁寧．
D. E. Aspnes and A. A. Studna, "Dielectric functions and optical parameters of Si, Ge, GaP, GaAs, GaSb, InP, InAs, and InSb," Phys. Rev. B 27, 985 (1983). ——図18.10の実験の基準データ．
L. Hedin, "New method for calculating the one-particle Green's function with application to the electron-gas problem," Phys. Rev. 139, A796 (1965). ——GW近似．
G. Onida, L. Reining, and A. Rubio, "Electronic excitations: density-functional versus many-body Green's-function approaches," Rev. Mod. Phys. 74, 601 (2002). ——GW/BSE/TDDFTを統一的に比較する総説．
E. Runge and E. K. U. Gross, "Density-functional theory for time-dependent systems," Phys. Rev. Lett. 52, 997 (1984). ——TDDFTの基礎定理．

STMと事例研究

J. Tersoff and D. R. Hamann, "Theory of the scanning tunneling microscope," Phys. Rev. B 31, 805 (1985). ——式\eqref{eq:18-tersoff}．
A. Fleurence, R. Friedlein, T. Ozaki, H. Kawai, Y. Wang, and Y. Yamada-Takamura, "Experimental evidence for epitaxial silicene on diboride thin films," Phys. Rev. Lett. 108, 245501 (2012).
C.-C. Lee, A. Fleurence, R. Friedlein, Y. Yamada-Takamura, and T. Ozaki, "First-principles study on competing phases of silicene," Phys. Rev. B 88, 165404 (2013);同 Phys. Rev. B 90, 075422 (2014). ——構造候補の全エネルギー比較とARPESとの比較．
C.-C. Lee, A. Fleurence, Y. Yamada-Takamura, and T. Ozaki, "Hidden mechanism for embedding the flat bands of Landau levels…," Phys. Rev. B 90, 241402(R) (2014). ——フォノンのソフトモードとストライプドメイン．
C.-C. Lee et al., "Core-level binding energy shifts as a probe of the structure of epitaxial silicene," Phys. Rev. B 95, 115437 (2017). ——18.9.4節のXPSによる構造判別．

本書の先へ

R. M. Martin, L. Reining, and D. M. Ceperley, Interacting Electrons: Theory and Computational Approaches (Cambridge University Press, 2016). ——GW・BSE・DMFT・QMCを一冊で見渡せる．本書の直接の続編にあたる．
A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, "Dynamical mean-field theory of strongly correlated fermion systems," Rev. Mod. Phys. 68, 13 (1996).
W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Rajagopal, "Quantum Monte Carlo simulations of solids," Rev. Mod. Phys. 73, 33 (2001).
J. Behler and M. Parrinello, "Generalized neural-network representation of high-dimensional potential-energy surfaces," Phys. Rev. Lett. 98, 146401 (2007). ——機械学習ポテンシャルの出発点．
F. Giustino, "Electron-phonon interactions from first principles," Rev. Mod. Phys. 89, 015003 (2017).
R. Resta and D. Vanderbilt, "Theory of polarization: A modern approach," in Physics of Ferroelectrics (Springer, 2007). ——Berry位相と電気分極．

本書はここで終わる．しかし読者にとっては，ここからが始まりであってほしい．第1章で「$3N$ 次元の波動関数は扱えない」と書いたところから，本章で「数百 eV の束縛エネルギーを $0.5$ eV で当てる」ところまで来た．その間にあったのは，いくつもの巧妙な近似と，その近似がなぜ許されるかという丁寧な吟味である．次に読者が新しい論文を読み，新しい計算を走らせるとき，そこにある式の一行一行が「どこから来て，何を仮定していて，どこで破れるのか」を問えるようになっていれば，本書の目的は果たされたことになる．