密度汎関数理論入門 — 目次 第II部 密度汎関数理論への道 / 第7章

第7章ジェリウムモデルと交換・相関エネルギー

第3章では電子間相互作用を完全に無視した自由電子ガスを調べ,運動エネルギーが $E/N = 2.21/r_s^2$ Ry と密度だけで書けることを見た.第5章ではハートリー・フォック(HF)法によって交換相互作用の意味を学び,第6章ではその計算を機械化する道具として第二量子化を整備した.本章ではこれらをすべて合流させる.舞台はジェリウムモデル(jellium model)——一様に塗り広げられた正電荷の背景の中を,互いにクーロン反発しながら動く電子ガス——である.このモデルは単純金属の伝導電子の理想化であると同時に,密度汎関数理論(DFT)全体の参照系でもある.本章の目標は,このモデルのエネルギーを高密度展開の第0次(運動エネルギー)と第1次(交換エネルギー)まで完全に手で計算しきることである.その結果,1電子あたりのエネルギーが $\rho^{2/3}$ と $\rho^{1/3}$ という密度だけの関数で書けることが分かる.これこそが「エネルギーは密度の汎関数である」という第9章のホーエンベルグ・コーンの定理,そして第10章・第11章の局所密度近似(LDA)の原型である.さらに,摂動論の第2次以降に潜む相関エネルギーの問題を概観し,RPA・ウィグナー結晶・量子モンテカルロ法によってそれがどう決着したかを見る.

この章で学ぶこと
  • ジェリウムモデルの定義と,ハミルトニアンの3分割 $\hat{H}=\hat{H}_{\rm el}+\hat{H}_{\rm b}+\hat{H}_{\rm el\text{-}b}$
  • 長距離クーロンの発散を制御する湯川収束因子 $e^{-\mu r}$ と,$V\to\infty$ を先・$\mu\to 0$ を後とする極限の順序
  • 平面波基底での第二量子化ハミルトニアンと,$q=0$ 項が背景電荷と厳密に相殺する仕組み
  • 無次元化により相互作用が $r_s$ に比例する摂動になること(高密度極限で摂動論が正当化される構造)
  • 0次:運動エネルギー $E_0/N = 2.21/r_s^2$ Ry の第二量子化による再導出
  • 1次:交換エネルギー $E_1/N = -0.916/r_s$ Ry の完全導出(角度積分・部分積分・対数積分をすべて明示)
  • ディラック交換汎関数 $E_x^{\rm LDA}[\rho] = -C_x\int\rho^{4/3}\dd^3r$ と,そのスピン分極系への拡張
  • 相関エネルギー:2次摂動の対数発散,RPA(Gell-Mann–Brueckner),ウィグナー結晶,量子モンテカルロとVWNパラメータ化

数学ノート:本章の単位系(ガウス単位・ハートリー・リュードベリ)

本書の基本はハートリー原子単位系($\hbar = m_e = e = 4\pi\varepsilon_0 = 1$)である.しかし電子ガスの理論には,歴史的な事情からガウス単位系で $e^2$ を陽に残し,エネルギーをリュードベリで測る慣習が根強く残っている.本章では両方の流儀を行き来できるように,導出の途中では $e^2$,$\hbar$,$m$ を陽に書き,最終結果を両方の単位で併記する.必要な換算は次の3本である.

ガウス単位系では,2つの点電荷 $q_1, q_2$ の間のクーロンポテンシャルは $q_1q_2/r$ であって $4\pi\varepsilon_0$ は現れない.したがって電子間の相互作用は $e^2/\abs{\rr_1-\rr_2}$ と書かれ,そのフーリエ変換は $4\pi e^2/q^2$ となる.原子単位に移るには単に $e^2 \to 1$,$\hbar\to1$,$m\to1$ と置けばよい.リュードベリ単位系は $\hbar^2/2m = 1$,$e^2 = 2$ と取る流儀であり,この単位では運動エネルギーが $\varepsilon = k^2$ と書ける.本章に現れる有名な数値,$2.21/r_s^2$ と $-0.916/r_s$ は,いずれもリュードベリ単位の値である.ハートリーに直すには2で割って $1.105/r_s^2$ と $-0.458/r_s$ となる.混同しやすい点なので,以後は必ず単位を明記する.

7.1 モデルの設定

まず,これから何を計算する対象とするのかを正確に定義する.本節ではジェリウムモデルの構成要素(電子ガス・一様正背景・電気的中性)を述べ,なぜこのモデルがDFTにとって特別なのかを説明し,第3章で導入した周期境界条件・平面波基底・密度パラメータ $r_s$ を必要な形に整えておく.

7.1.1 なぜジェリウムなのか

Na や K のような単純金属(アルカリ金属)の電子構造を考えよう.Na の原子核の周りには $1s^2 2s^2 2p^6$ の内殻電子が固く束縛されており,その外側に $3s$ 電子が1個ある.この $3s$ 電子は固体中で原子から離れて結晶全体に広がり,伝導電子となる.すると金属は「格子点に固定された $\mathrm{Na}^+$ イオン(内殻電子ごと)」と「その間を動き回る伝導電子ガス」からできていると見なせる.

ここで大胆な近似を行う.イオンの正電荷を,格子点に局在させるのをやめて,空間に一様に塗り広げてしまうのである.イオンの周期的な配列に由来する情報(バンド構造,ブリルアンゾーン,結晶構造の安定性)はすべて失われるが,その代わりに「相互作用する電子ガス」という多体問題の本質だけが純粋な形で残る.この理想化されたモデルをジェリウムモデルと呼ぶ.名前は正電荷が「ゼリー(jelly)」のように一様に広がっている様子に由来する.

+++ +++ +++ 実在の単純金属 周期的に並んだイオン + 伝導電子 正電荷を一様に 塗り広げる +++++ +++++ +++++ +++++ +++++ ジェリウムモデル 一様正背景 + 相互作用する電子ガス 失われるもの:バンド構造・結晶構造 / 残るもの:多体電子問題の本質
図7.1 ジェリウムモデルの発想.イオンの正電荷を一様な背景に置き換えることで,周期ポテンシャルに由来する複雑さを捨て,電子間相互作用そのものだけを取り出す.

このモデルには2つの顔がある.

定義:ジェリウムモデル

一辺 $L$,体積 $V=L^3$ の立方体に周期境界条件を課す.この中に

  1. 質量 $m$,電荷 $-e$ の電子を $N$ 個入れる.電子どうしは完全なクーロン相互作用 $e^2/\abs{\rr_i-\rr_j}$ で反発する.
  2. 電荷密度 $+e\,n_{\rm b}$($n_{\rm b} = N/V$ は定数)の剛体的な正電荷背景を一様に分布させる.背景は動かない(自由度を持たない).

系全体の電荷は $-eN + e\,n_{\rm b}V = 0$ であり,電気的に中性である.この中性条件が本章の議論の要である.

背景電荷を「動かない」としたのは,イオンの質量が電子の $10^4$ 倍以上あってほとんど動かないという断熱(ボルン・オッペンハイマー)近似(第1章)の反映である.したがってジェリウムのハミルトニアンには背景の運動エネルギーは現れず,背景は電子にとって単なる外部ポテンシャル(と,定数のエネルギー)として作用する.

7.1.2 周期境界条件と平面波基底(第3章の復習)

第3章と同様に,体積 $V=L^3$ の立方体にボルン・フォン・カルマンの周期境界条件

$$ \begin{equation} \psi(x+L, y, z) = \psi(x,y,z),\qquad \psi(x, y+L, z) = \psi(x,y,z),\qquad \psi(x, y, z+L) = \psi(x,y,z) \label{eq:7-pbc} \end{equation} $$

を課す.この条件を満たす規格化された1電子固有関数(スピンも含む)は平面波である:

$$ \begin{equation} \varphi_{\kk\lambda}(x) = \frac{1}{\sqrt{V}}\,e^{i\kk\cdot\rr}\,\eta_\lambda(\sigma), \qquad \kk = \frac{2\pi}{L}\,(n_1, n_2, n_3),\quad n_i = 0,\pm1,\pm2,\ldots \label{eq:7-planewave} \end{equation} $$

ここで $x=(\rr,\sigma)$ はスピンを含む合成変数,$\eta_\lambda(\sigma)$($\lambda=\uparrow,\downarrow$)はスピン関数で $\sum_\sigma \eta_\lambda^*(\sigma)\eta_{\lambda'}(\sigma)=\delta_{\lambda\lambda'}$ を満たす.波数が $2\pi n_i/L$ という離散値に制限されるのは,$e^{ik_xL}=1$ すなわち $k_xL = 2\pi n_1$ が要求されるからである.

本章で繰り返し使う直交性の関係を確認しておく.周期境界条件のもとでは

$$ \begin{equation} \frac{1}{V}\int_V \dd^3r\; e^{i(\kk'-\kk)\cdot\rr} = \delta_{\kk\kk'} \label{eq:7-pworth} \end{equation} $$

が成り立つ.実際,$\kk=\kk'$ なら被積分関数は1で積分は $V$,よって左辺は1である.$\kk\neq\kk'$ なら,たとえば $x$ 成分について $\int_0^L \dd x\, e^{i(k_x'-k_x)x} = \big[e^{i(k_x'-k_x)x}/i(k_x'-k_x)\big]_0^L = 0$ となる($k_x'-k_x = 2\pi\Delta n/L$ が代入されて $e^{2\pi i \Delta n}=1$ となり,上端と下端が打ち消し合う).よって左辺は0である.

また,和を積分に置き換える規則(第3章)

$$ \begin{equation} \sum_{\kk} \;\longrightarrow\; \frac{V}{(2\pi)^3}\int \dd^3k \label{eq:7-sum2int} \end{equation} $$

も本章の中心的な道具である.これは $\kk$ 空間で1つの許容点が体積 $(2\pi/L)^3 = (2\pi)^3/V$ を占めることから従う.

7.1.3 密度パラメータ $r_s$ の再確認

第3章で導入した密度パラメータを再掲する.電子1個あたりの体積を半径 $r_0$ の球で表し,それをボーア半径で測ったものが $r_s$ である:

$$ \begin{equation} \frac{4}{3}\pi r_0^3 = \frac{V}{N} = \frac{1}{\rho}, \qquad r_s \equiv \frac{r_0}{a_0}, \qquad a_0 = \frac{\hbar^2}{me^2}. \label{eq:7-rs-def} \end{equation} $$

ここで $\rho = N/V$ は電子数密度である($\rho$ と $n_{\rm b}$ は数値として等しいが,前者は電子,後者は背景の密度を指す).第3章で導いた関係を書き下しておく:

$$ \begin{align} \rho &= \frac{3}{4\pi r_0^3}, \label{eq:7-rho-r0}\\ k_F &= (3\pi^2\rho)^{1/3} = \left(\frac{9\pi}{4}\right)^{1/3}\frac{1}{r_0} = \frac{1.91916}{r_s}\ a_0^{-1}. \label{eq:7-kF-rs} \end{align} $$

1行目は \eqref{eq:7-rs-def} の第1式を $\rho$ について解いたもの,2行目は $k_F^3 = 3\pi^2\rho = 3\pi^2\cdot 3/(4\pi r_0^3) = (9\pi/4)/r_0^3$ の3乗根である.実在の単体金属では $2\lesssim r_s\lesssim 6$ であった(第3章表3.1).この「狭い窓」の中で理論を精密にすることが本章と第11章の課題である.

7.2 ハミルトニアンの3分割と長距離発散

この節では,ジェリウムのハミルトニアンを電子部分・背景部分・電子背景相互作用の3つに分ける.そのうえで,クーロン相互作用が長距離であるために各部分が個別には発散してしまうという困難に直面し,それを制御する湯川収束因子を導入する.最後に,背景に関する2つの項を球座標積分で完全に評価し,それらが $1 : -2$ という比で現れることを確かめる.

7.2.1 3つの部分

系のエネルギーに寄与するのは,(i) 電子の運動エネルギーと電子間クーロン反発,(ii) 背景電荷どうしのクーロン反発,(iii) 電子と背景電荷の間のクーロン引力,の3種類である.それぞれを $\hat{H}_{\rm el}$,$\hat{H}_{\rm b}$,$\hat{H}_{\rm el\text{-}b}$ と書く:

$$ \begin{equation} \hat{H} = \hat{H}_{\rm el} + \hat{H}_{\rm b} + \hat{H}_{\rm el\text{-}b}. \label{eq:7-H3} \end{equation} $$

具体形は次の通りである.背景の電荷密度を $\rho_{\rm b}(\rr) = n_{\rm b} = N/V$(一定)と書くと,

$$ \begin{align} \hat{H}_{\rm el} &= \sum_{i=1}^{N}\frac{\hat{\bm{p}}_i^2}{2m} + \frac{e^2}{2}\sum_{i\neq j}\frac{1}{\abs{\rr_i-\rr_j}}, \label{eq:7-Hel-raw}\\[2pt] \hat{H}_{\rm b} &= \frac{e^2}{2}\int \dd^3r\int \dd^3r'\; \frac{\rho_{\rm b}(\rr)\,\rho_{\rm b}(\rr')}{\abs{\rr-\rr'}}, \label{eq:7-Hb-raw}\\[2pt] \hat{H}_{\rm el\text{-}b} &= -\,e^2\sum_{i=1}^{N}\int \dd^3r\; \frac{\rho_{\rm b}(\rr)}{\abs{\rr-\rr_i}}. \label{eq:7-Helb-raw} \end{align} $$

符号と係数の由来を確認しておく.\eqref{eq:7-Hel-raw} の因子 $1/2$ は,$i\neq j$ の二重和で電子対 $(i,j)$ を2回数えてしまう分の補正である.\eqref{eq:7-Hb-raw} の $1/2$ も同じ理由(背景の微小要素の対を2回数える)による.\eqref{eq:7-Helb-raw} には $1/2$ が付かない.これは電子と背景という異なる種類の電荷の間の相互作用であり,二重計算が起きないからである.符号は,電子の電荷が $-e$,背景が $+e\rho_{\rm b}$ であるから,積 $(-e)(+e\rho_{\rm b}) = -e^2\rho_{\rm b}$ となって引力(負)になる.なお $\hat{H}_{\rm b}$ は電子の座標を含まないので演算子ではなく単なる定数である.

7.2.2 困難:各項が個別に発散する

ここで直ちに問題が生じる.$\hat{H}_{\rm b}$ の積分を実行してみよう.$\rho_{\rm b}$ が定数なので

$$ \hat{H}_{\rm b} = \frac{e^2 n_{\rm b}^2}{2}\int \dd^3r\int \dd^3r'\;\frac{1}{\abs{\rr-\rr'}} $$

となるが,$\rr' $ の積分を $\bm{s}=\rr'-\rr$ に変えると $\int \dd^3s\,/s = 4\pi\int_0^{R} s\,\dd s = 2\pi R^2$ となり,系の大きさ $R$ とともに発散する.同じことが $\hat{H}_{\rm el\text{-}b}$ でも,また $\hat{H}_{\rm el}$ の相互作用項の「平均場的な部分」でも起きる.原因はクーロン相互作用 $1/r$ が長距離である(遠方で十分速く減衰しない)ことである.

物理的には,これは驚くことではない.正電荷だけを集めた系,あるいは負電荷だけを集めた系のエネルギーが体積とともに発散するのは当然である.しかし全体が中性であるジェリウムでは,3つの発散が互いに打ち消し合うはずである.問題は「打ち消し合う」という操作を数学的に正当な手続きで行うことにある.$\infty - 2\infty + \infty$ という式に意味を与えるには,まず各項を有限にしておかねばならない.

数学ノート:湯川型の収束因子と極限の順序

発散を制御する標準的な方法は,クーロン相互作用を湯川(Yukawa)型に置き換えることである:

$$ \frac{1}{r}\;\longrightarrow\;\frac{e^{-\mu r}}{r},\qquad \mu > 0 . $$

$\mu$ は逆長さの次元を持つ正の定数で,$1/\mu$ が相互作用の到達距離を与える.指数関数の減衰のおかげで,すべての積分が絶対収束する.計算をすべて済ませたのち $\mu\to 0$ とすれば,もとのクーロン系に戻る.この形は湯川秀樹が中間子交換による核力を記述するために導入したポテンシャルと同じ形であり,素粒子・物性の両分野で「長距離力を一時的に短距離化する」ための標準的な道具になっている.物理的には,遠方に置かれた誘電体が電荷を遮蔽している状況,あるいは第8章で扱うトーマス・フェルミ遮蔽ポテンシャル $e^{-\kappa r}/r$ と同じ形である.

ここで決定的に重要なのが極限の順序である.本章では

$$ \lim_{\mu\to0}\ \lim_{V\to\infty}\quad(\text{熱力学極限 }V\to\infty\text{ が先,}\mu\to0\text{ が後}) $$

の順で極限をとる.すなわち,$\mu$ を固定したまま,まず密度 $N/V$ を一定に保って $V\to\infty$,$N\to\infty$ とし,そのあとで $\mu\to0$ とする.この順序が意味するのは,相互作用の到達距離 $1/\mu$ が系の大きさ $L$ よりずっと短い($\mu L \gg 1$)という状況を常に保つことである.この仮定のもとでは,系の内部の1点から見て相互作用が届く範囲は完全に系の内部に収まるので,境界の影響を無視して積分の上端を $\infty$ にしてよい.

順序を逆にする($\mu\to0$ を先にする)と,有限の箱の中で無限に長い到達距離を持つ相互作用を扱うことになり,境界の効果が支配的になって物理的に無意味な結果を得る.実際,7.4節で見るように,相殺し残った項は $\propto 1/(\mu^2 V)$ という形をしており,$V\to\infty$ を先にすればゼロ,$\mu\to0$ を先にすれば無限大になる.同じ式が極限の順序で異なる値をとる典型例であり,順序の指定は「約束事」ではなく物理の一部である.

そこで,\eqref{eq:7-Hel-raw}–\eqref{eq:7-Helb-raw} のすべての $1/\abs{\rr-\rr'}$ を $e^{-\mu\abs{\rr-\rr'}}/\abs{\rr-\rr'}$ に置き換える:

$$ \begin{align} \hat{H}_{\rm el} &= \sum_{i=1}^{N}\frac{\hat{\bm{p}}_i^2}{2m} + \frac{e^2}{2}\sum_{i\neq j}\frac{e^{-\mu\abs{\rr_i-\rr_j}}}{\abs{\rr_i-\rr_j}}, \label{eq:7-Hel-mu}\\[2pt] \hat{H}_{\rm b} &= \frac{e^2}{2}\,n_{\rm b}^2\int \dd^3r\int \dd^3r'\; \frac{e^{-\mu\abs{\rr-\rr'}}}{\abs{\rr-\rr'}}, \label{eq:7-Hb-mu}\\[2pt] \hat{H}_{\rm el\text{-}b} &= -\,e^2\,n_{\rm b}\sum_{i=1}^{N}\int \dd^3r\; \frac{e^{-\mu\abs{\rr-\rr_i}}}{\abs{\rr-\rr_i}}. \label{eq:7-Helb-mu} \end{align} $$

7.2.3 基本積分 $\int e^{-\mu r}/r\,\dd^3r = 4\pi/\mu^2$

3つの項をすべて評価するために必要な積分は,実は次の1つだけである.これを球座標で完全に計算する.

導出:湯川ポテンシャルの全空間積分

球座標 $(r,\theta,\phi)$ を用いる.被積分関数が角度に依存しないことに注意して,

$$ \begin{align} \int \dd^3r\;\frac{e^{-\mu r}}{r} &= \int_0^{\infty}\dd r\, r^2 \int_0^{\pi}\dd\theta\,\sin\theta \int_0^{2\pi}\dd\phi\; \frac{e^{-\mu r}}{r} \label{eq:7-yuk1}\\ &= \left(\int_0^{2\pi}\dd\phi\right)\left(\int_0^{\pi}\sin\theta\,\dd\theta\right) \int_0^{\infty}\dd r\; r^2\,\frac{e^{-\mu r}}{r} \label{eq:7-yuk2}\\ &= 2\pi\cdot 2\cdot \int_0^{\infty}\dd r\; r\,e^{-\mu r} \label{eq:7-yuk3}\\ &= 4\pi\cdot\frac{1}{\mu^2} \;=\;\frac{4\pi}{\mu^2}. \label{eq:7-yuk4} \end{align} $$

\eqref{eq:7-yuk1} は3次元の体積要素 $\dd^3r = r^2\sin\theta\,\dd r\,\dd\theta\,\dd\phi$ を書き下しただけである.\eqref{eq:7-yuk2} では,被積分関数が $r$ のみの関数なので角度部分の積分を分離した.\eqref{eq:7-yuk3} では $\int_0^{2\pi}\dd\phi = 2\pi$,$\int_0^\pi\sin\theta\,\dd\theta = [-\cos\theta]_0^\pi = 1-(-1) = 2$ を代入し,$r^2/r = r$ とした.

\eqref{eq:7-yuk4} の $r$ 積分は部分積分で計算できる.$u = r$,$\dd v = e^{-\mu r}\dd r$($v = -e^{-\mu r}/\mu$)とすると

$$ \int_0^{\infty} r\,e^{-\mu r}\,\dd r = \left[-\frac{r\,e^{-\mu r}}{\mu}\right]_0^{\infty} + \frac{1}{\mu}\int_0^{\infty} e^{-\mu r}\,\dd r = 0 + \frac{1}{\mu}\left[-\frac{e^{-\mu r}}{\mu}\right]_0^{\infty} = \frac{1}{\mu}\cdot\frac{1}{\mu} = \frac{1}{\mu^2}. $$

境界項が消えるのは,$\mu>0$ のとき $r\to\infty$ で $re^{-\mu r}\to0$(指数関数の減衰が $r$ の増大に勝つ)であり,$r=0$ では $r$ 自身がゼロだからである.ここで $\mu>0$ が本質的に効いている:$\mu=0$ ならこの積分は発散する.

結果を書き下す:

$$ \begin{equation} \int \dd^3r\;\frac{e^{-\mu r}}{r} = \frac{4\pi}{\mu^2}. \label{eq:7-yukawa-int} \end{equation} $$

$\mu\to0$ でこれが発散することが,まさにクーロン相互作用の長距離性の表れである.

7.2.4 背景・背景項の評価

これから $\hat{H}_{\rm b}$ を評価する.手順は「相対座標に変数変換して \eqref{eq:7-yukawa-int} を使い,残った重心座標の積分が体積 $V$ を与える」というものである.

導出:$\hat{H}_{\rm b}$

$$ \begin{align} \hat{H}_{\rm b} &= \frac{e^2}{2}\,n_{\rm b}^2 \int_V \dd^3r \int_V \dd^3r'\;\frac{e^{-\mu\abs{\rr-\rr'}}}{\abs{\rr-\rr'}} \label{eq:7-Hb1}\\ &= \frac{e^2}{2}\,n_{\rm b}^2 \int_V \dd^3r \left[\int \dd^3s\;\frac{e^{-\mu s}}{s}\right] \label{eq:7-Hb2}\\ &= \frac{e^2}{2}\,n_{\rm b}^2 \int_V \dd^3r \;\frac{4\pi}{\mu^2} \label{eq:7-Hb3}\\ &= \frac{e^2}{2}\,n_{\rm b}^2\, V\,\frac{4\pi}{\mu^2} \label{eq:7-Hb4}\\ &= \frac{1}{2}\,e^2\,\frac{N^2}{V}\,\frac{4\pi}{\mu^2}. \label{eq:7-Hb5} \end{align} $$

\eqref{eq:7-Hb2} では,$\rr$ を固定して $\bm{s} = \rr'-\rr$ と変数変換した.ヤコビアンは1である.積分範囲は本来 $\rr$ に依存する領域になるが,7.2.2 の数学ノートで述べたとおり $\mu L\gg1$ を仮定しているので,被積分関数は $s\gtrsim 1/\mu \ll L$ で既に無視できるほど小さい.したがって積分範囲を全空間に広げてよい(このとき $\rr$ 依存性が消える).

\eqref{eq:7-Hb3} で基本積分 \eqref{eq:7-yukawa-int} を代入した.\eqref{eq:7-Hb4} では,残った被積分関数が定数なので $\int_V\dd^3r = V$ とした.\eqref{eq:7-Hb5} では $n_{\rm b} = N/V$ を代入し,$n_{\rm b}^2 V = N^2/V$ とした.

1電子あたりに直すと

$$ \begin{equation} \frac{\hat{H}_{\rm b}}{N} = \frac{1}{2}\,e^2\,\frac{N}{V}\,\frac{4\pi}{\mu^2} \label{eq:7-Hb-perN} \end{equation} $$

である.正電荷だけを集めた系のエネルギーなので,当然ながら正(反発)である.

7.2.5 電子・背景項の評価

次に $\hat{H}_{\rm el\text{-}b}$ を評価する.今度は電子の位置 $\rr_i$ を中心とする積分になるが,やはり基本積分 \eqref{eq:7-yukawa-int} に帰着する.

導出:$\hat{H}_{\rm el\text{-}b}$

$$ \begin{align} \hat{H}_{\rm el\text{-}b} &= -e^2 n_{\rm b}\sum_{i=1}^{N}\int_V \dd^3r\;\frac{e^{-\mu\abs{\rr-\rr_i}}}{\abs{\rr-\rr_i}} \label{eq:7-Helb1}\\ &= -e^2 n_{\rm b}\sum_{i=1}^{N}\int \dd^3s\;\frac{e^{-\mu s}}{s} \label{eq:7-Helb2}\\ &= -e^2 n_{\rm b}\sum_{i=1}^{N}\frac{4\pi}{\mu^2} \label{eq:7-Helb3}\\ &= -e^2\,n_{\rm b}\,N\,\frac{4\pi}{\mu^2} \;=\; -\,e^2\,\frac{N^2}{V}\,\frac{4\pi}{\mu^2}. \label{eq:7-Helb4} \end{align} $$

\eqref{eq:7-Helb2} では $\bm{s} = \rr - \rr_i$ と変数変換した.積分は電子の位置 $\rr_i$ に依存しなくなることに注意しよう.これは背景が一様であるためであり,「一様な背景は電子に力を及ぼさない」(第3章)という主張の定量的な表れである.\eqref{eq:7-Helb3} で基本積分を代入し,\eqref{eq:7-Helb4} では被加算項が $i$ に依存しないので和が単に $N$ 倍になることを使った.

1電子あたりでは

$$ \begin{equation} \frac{\hat{H}_{\rm el\text{-}b}}{N} = -\,e^2\,\frac{N}{V}\,\frac{4\pi}{\mu^2} \label{eq:7-Helb-perN} \end{equation} $$

となる.電子と正背景の間は引力なので負である.

7.2.6 $1 : -2$ という比の意味

\eqref{eq:7-Hb-perN} と \eqref{eq:7-Helb-perN} を見比べると,両者は同じ発散因子 $4\pi/\mu^2$ を共有し,係数の比が

$$ \hat{H}_{\rm b} : \hat{H}_{\rm el\text{-}b} = \frac{1}{2} : (-1) = 1 : (-2) $$

になっている.この $1:-2$ という比は偶然ではない.相互作用する2種類の一様電荷分布(密度 $+e n_{\rm b}$ と $-e n_{\rm b}$)を考えると,静電エネルギーは「正・正」「正・負」「負・負」の3種類に分かれ,同種どうしには二重計算を防ぐ $1/2$ が付き,異種間には付かない.したがって係数は $\frac12 : -1 : \frac12$,すなわち $1:-2:1$ になる.3つを足せば $\frac12 - 1 + \frac12 = 0$ である.

いま我々の手元には,この3つのうち「正・正」と「正・負」の2つが揃った.合計すると

$$ \begin{equation} \hat{H}_{\rm b} + \hat{H}_{\rm el\text{-}b} = \frac{1}{2}e^2\frac{N^2}{V}\frac{4\pi}{\mu^2} - e^2\frac{N^2}{V}\frac{4\pi}{\mu^2} = -\frac{1}{2}\,e^2\,\frac{N^2}{V}\,\frac{4\pi}{\mu^2} \label{eq:7-Hb-plus-Helb} \end{equation} $$

である.残る「負・負」,すなわち電子どうしの反発の平均場的部分は $\hat{H}_{\rm el}$ の中に隠れており,それを取り出すのが7.4節の仕事である.そこでちょうど $+\frac12 e^2 (N^2/V)(4\pi/\mu^2)$ が現れて \eqref{eq:7-Hb-plus-Helb} と相殺する,というのがこれから示すシナリオである.

背景 — 背景 + + 係数 +1/2 反発(正) 電子 — 背景 + − 係数 −1 引力(負) 電子 — 電子 − − 係数 +1/2(の q=0 部分) 反発(正) 同じ発散因子 e² (N²/V)(4π/μ²) を共有し,係数の和は +1/2 − 1 + 1/2 = 0 中性条件のもとで長距離発散は完全に消える(7.4節で厳密に示す)
図7.2 ジェリウムにおける静電エネルギーの3成分と,その係数 $+\frac12 : -1 : +\frac12$.中性であることが相殺の必要十分な条件である.

物理的意味:なぜ背景電荷が必要なのか

もし背景電荷を導入せず,電子だけを箱に詰めたとすると,電子間反発の平均場的部分 $\frac12 e^2 (N^2/V)(4\pi/\mu^2)$ が残り,$\mu\to0$ で発散する.1電子あたりでは $\frac12 e^2 (N/V)(4\pi/\mu^2)$ であり,これは有限密度を保ったまま発散する.つまり電荷を帯びた電子ガスは熱力学極限を持たない.$N$ 個の電子を集めた系の全電荷は $-Ne$ であり,そのクーロン自己エネルギーは $\sim (Ne)^2/R \sim N^{5/3}$ のように $N$ より速く増大するから,1電子あたりのエネルギーが発散するのは当然である.

背景電荷は「余計な仮定」ではなく,熱力学極限が存在するための必要条件なのである.実際の金属では正イオンがその役割を果たしており,金属が安定に存在すること自体が電気的中性の帰結である.ジェリウムモデルの正背景は,その最も単純な理想化にほかならない.

7.3 電子部分の第二量子化

ここからは第6章で整備した道具を使う.平面波基底 \eqref{eq:7-planewave} を1電子基底に選び,$\hat{H}_{\rm el}$ を生成・消滅演算子で書き直す.運動エネルギーは対角形になり,相互作用項には運動量保存を表すクロネッカーのデルタが現れる.この節の最終目標は,ジェリウムの電子ハミルトニアンを $\hat{H}_{\rm el} = \hat{T} + \hat{V}$ の完全な第二量子化形で書き下すことである.

7.3.1 1体項:運動エネルギー

第6章の結果によれば,1体演算子は

$$ \begin{equation} \sum_{i}\hat{V}_1(x_i) = \sum_{\mu\nu}\langle\varphi_\mu|\hat{V}_1|\varphi_\nu\rangle\,a_\mu^\dagger a_\nu \label{eq:7-onebody-general} \end{equation} $$

と書ける.$\hat{V}_1 = \hat{\bm{p}}^2/2m = -(\hbar^2/2m)\nabla^2$ とし,基底を平面波にとって行列要素を計算しよう.

導出:運動エネルギーの行列要素

まず平面波にラプラシアンを作用させる.$e^{i\kk'\cdot\rr}$ の各成分について $\partial^2 e^{ik_x'x}/\partial x^2 = (ik_x')^2 e^{ik_x'x} = -k_x'^2 e^{ik_x'x}$ であるから,

$$ -\frac{\hbar^2}{2m}\nabla^2\,\frac{e^{i\kk'\cdot\rr}}{\sqrt{V}}\eta_{\lambda'}(\sigma) = -\frac{\hbar^2}{2m}\big(-\abs{\kk'}^2\big)\frac{e^{i\kk'\cdot\rr}}{\sqrt{V}}\eta_{\lambda'}(\sigma) = \frac{\hbar^2\abs{\kk'}^2}{2m}\,\varphi_{\kk'\lambda'}(x). $$

すなわち平面波は運動エネルギー演算子の固有関数である.したがって

$$ \begin{align} \langle\varphi_{\kk\lambda}|-\tfrac{\hbar^2}{2m}\nabla^2|\varphi_{\kk'\lambda'}\rangle &= \frac{\hbar^2\abs{\kk'}^2}{2m}\int\dd x\;\varphi_{\kk\lambda}^*(x)\varphi_{\kk'\lambda'}(x) \label{eq:7-kin1}\\ &= \frac{\hbar^2\abs{\kk'}^2}{2m} \left[\frac{1}{V}\int_V\dd^3r\,e^{i(\kk'-\kk)\cdot\rr}\right] \left[\sum_\sigma \eta_\lambda^*(\sigma)\eta_{\lambda'}(\sigma)\right] \label{eq:7-kin2}\\ &= \frac{\hbar^2\abs{\kk'}^2}{2m}\,\delta_{\kk\kk'}\,\delta_{\lambda\lambda'} \;=\;\frac{\hbar^2\abs{\kk}^2}{2m}\,\delta_{\kk\kk'}\,\delta_{\lambda\lambda'}. \label{eq:7-kin3} \end{align} $$

\eqref{eq:7-kin2} では $\int\dd x = \sum_\sigma\int\dd^3r$ に従って空間部分とスピン部分を分離した.\eqref{eq:7-kin3} では平面波の直交性 \eqref{eq:7-pworth} とスピン関数の直交性を使い,最後に $\delta_{\kk\kk'}$ があるので $\kk'$ を $\kk$ に書き換えた.

これを \eqref{eq:7-onebody-general} に代入すると,二重和のうち $\kk=\kk'$,$\lambda=\lambda'$ の項だけが残り,

$$ \begin{equation} \hat{T} = \sum_{\kk\lambda}\frac{\hbar^2\abs{\kk}^2}{2m}\,a_{\kk\lambda}^\dagger a_{\kk\lambda} = \sum_{\kk\lambda}\frac{\hbar^2\abs{\kk}^2}{2m}\,\hat{n}_{\kk\lambda} \label{eq:7-T-2q} \end{equation} $$

を得る.運動エネルギーは数演算子の重み付き和という完全に対角な形になった.平面波を選んだ最大の御利益がこれである.

7.3.2 2体項:湯川ポテンシャルのフーリエ変換

次が本節の山場である.第6章の結果によれば,2体演算子は

$$ \begin{equation} \frac{1}{2}\sum_{i\neq j}\hat{V}_2(x_i,x_j) = \frac{1}{2}\sum_{\mu\nu\gamma\delta} \langle\varphi_\mu\varphi_\nu|\hat{V}_2|\varphi_\gamma\varphi_\delta\rangle\, a_\mu^\dagger a_\nu^\dagger a_\delta a_\gamma \label{eq:7-twobody-general} \end{equation} $$

であり,2電子積分は

$$ \begin{equation} \langle\varphi_\mu\varphi_\nu|\hat{V}_2|\varphi_\gamma\varphi_\delta\rangle = \int\dd x_1\int\dd x_2\; \varphi_\mu^*(x_1)\varphi_\nu^*(x_2)\,\hat{V}_2(x_1,x_2)\,\varphi_\gamma(x_1)\varphi_\delta(x_2) \label{eq:7-tei-def} \end{equation} $$

と定義されていた(ブラの第1・ケットの第1が変数 $x_1$,ブラの第2・ケットの第2が変数 $x_2$ に対応する,という入れ子の規約).これを平面波基底で評価する.まず,そのために必要なフーリエ変換を計算しておく.

導出:湯川ポテンシャルの3次元フーリエ変換

求めたいのは

$$ \tilde{v}_\mu(\bm{q}) \equiv \int \dd^3r\;e^{-i\bm{q}\cdot\rr}\,\frac{e^{-\mu r}}{r} $$

である.被積分関数は $\rr$ の向きに $\bm{q}\cdot\rr$ を通してのみ依存する.そこで $\bm{q}$ を $z$ 軸にとる球座標を使う(積分は座標軸の取り方によらないので,この選択は自由である).すると $\bm{q}\cdot\rr = qr\cos\theta$ であり,

$$ \begin{align} \tilde{v}_\mu(\bm{q}) &= \int_0^{\infty}\dd r\,r^2\int_0^{\pi}\dd\theta\,\sin\theta\int_0^{2\pi}\dd\phi\; e^{-iqr\cos\theta}\,\frac{e^{-\mu r}}{r} \label{eq:7-ft1}\\ &= 2\pi\int_0^{\infty}\dd r\; r\,e^{-\mu r} \int_0^{\pi}\dd\theta\,\sin\theta\; e^{-iqr\cos\theta} \label{eq:7-ft2}\\ &= 2\pi\int_0^{\infty}\dd r\; r\,e^{-\mu r} \int_{-1}^{1}\dd u\; e^{-iqru} \label{eq:7-ft3}\\ &= 2\pi\int_0^{\infty}\dd r\; r\,e^{-\mu r}\cdot \left[\frac{e^{-iqru}}{-iqr}\right]_{u=-1}^{u=1} \label{eq:7-ft4}\\ &= 2\pi\int_0^{\infty}\dd r\; r\,e^{-\mu r}\cdot \frac{e^{-iqr}-e^{iqr}}{-iqr} \label{eq:7-ft5}\\ &= 2\pi\int_0^{\infty}\dd r\; r\,e^{-\mu r}\cdot \frac{-2i\sin(qr)}{-iqr} \label{eq:7-ft6}\\ &= \frac{4\pi}{q}\int_0^{\infty}\dd r\; e^{-\mu r}\sin(qr). \label{eq:7-ft7} \end{align} $$

\eqref{eq:7-ft2} では $\phi$ 積分($2\pi$)を実行し,$r^2/r=r$ とした.\eqref{eq:7-ft3} では $u=\cos\theta$ と置換した($\dd u = -\sin\theta\,\dd\theta$ であり,積分範囲 $\theta:0\to\pi$ が $u:1\to-1$ に移るので符号が戻って $\int_{-1}^{1}\dd u$ となる).\eqref{eq:7-ft4} は $u$ についての指数関数の積分である.\eqref{eq:7-ft6} ではオイラーの公式 $e^{i\alpha}-e^{-i\alpha} = 2i\sin\alpha$ を使い,\eqref{eq:7-ft7} では $-2i/(-i) = 2$ と $r/r=1$ を整理して $2\pi\cdot 2/q = 4\pi/q$ とした.

残った1次元積分は,収束因子 $e^{-\mu r}$ のおかげで絶対収束する.複素指数に直して計算する:

$$ \begin{align} \int_0^{\infty}e^{-\mu r}\sin(qr)\,\dd r &= \mathrm{Im}\int_0^{\infty}e^{-\mu r}e^{iqr}\,\dd r \label{eq:7-sin1}\\ &= \mathrm{Im}\int_0^{\infty}e^{-(\mu-iq)r}\,\dd r \label{eq:7-sin2}\\ &= \mathrm{Im}\left[\frac{-e^{-(\mu-iq)r}}{\mu-iq}\right]_0^{\infty} = \mathrm{Im}\,\frac{1}{\mu-iq} \label{eq:7-sin3}\\ &= \mathrm{Im}\,\frac{\mu+iq}{(\mu-iq)(\mu+iq)} = \mathrm{Im}\,\frac{\mu+iq}{\mu^2+q^2} = \frac{q}{\mu^2+q^2}. \label{eq:7-sin4} \end{align} $$

\eqref{eq:7-sin1} は $\sin(qr) = \mathrm{Im}\,e^{iqr}$ による.\eqref{eq:7-sin3} で上端が消えるのは $\abs{e^{-(\mu-iq)r}} = e^{-\mu r}\to0$($\mu>0$)だからである.ここでも $\mu>0$ が本質的である.\eqref{eq:7-sin4} では分母を実数化するために共役 $\mu+iq$ を掛けた.

\eqref{eq:7-ft7} に代入すると

$$ \tilde{v}_\mu(\bm{q}) = \frac{4\pi}{q}\cdot\frac{q}{\mu^2+q^2} = \frac{4\pi}{\mu^2+q^2} $$

を得る.結果が $\bm{q}$ の向きによらず大きさ $q=\abs{\bm{q}}$ のみの関数になっているのは,もとのポテンシャルが球対称だからである.

$$ \begin{equation} \int \dd^3r\;e^{-i\bm{q}\cdot\rr}\,\frac{e^{-\mu r}}{r} = \frac{4\pi}{\mu^2+q^2} \qquad\xrightarrow[\ \mu\to0\ ]{}\qquad \frac{4\pi}{q^2}. \label{eq:7-yukawa-ft} \end{equation} $$

$\mu\to0$ の極限がクーロンポテンシャルのフーリエ変換 $4\pi/q^2$ である.$q\neq0$ では極限が問題なく存在するが,$q=0$ では $4\pi/\mu^2$ が発散する.長距離発散は $\bm{q}=\bm{0}$ という1点に集中している——これが本章の議論全体の鍵である.$\bm{q}=\bm{0}$ の成分とは,実空間で言えば「空間的に一様な成分」,すなわち平均電荷密度が作る平均場のことである.中性系ではこの平均場が背景電荷によって打ち消される.

7.3.3 2電子積分の評価と運動量保存

準備が整ったので,平面波基底での2電子積分 \eqref{eq:7-tei-def} を計算する.$\mu = \kk_1\lambda_1$,$\nu = \kk_2\lambda_2$,$\gamma = \kk_3\lambda_3$,$\delta = \kk_4\lambda_4$ とおく(収束因子の $\mu$ と添字の $\mu$ が紛らわしいので,以下では添字を波数とスピンで陽に書く).

導出:平面波基底での2電子積分

$$ \begin{align} &\langle \kk_1\lambda_1\,\kk_2\lambda_2 |\hat{V}_2| \kk_3\lambda_3\,\kk_4\lambda_4\rangle \notag\\ &\quad= \sum_{\sigma_1\sigma_2}\int\dd^3r_1\int\dd^3r_2\; \frac{e^{-i\kk_1\cdot\rr_1}}{\sqrt V}\eta_{\lambda_1}^*(\sigma_1) \frac{e^{-i\kk_2\cdot\rr_2}}{\sqrt V}\eta_{\lambda_2}^*(\sigma_2)\; e^2\frac{e^{-\mu r_{12}}}{r_{12}}\; \frac{e^{i\kk_3\cdot\rr_1}}{\sqrt V}\eta_{\lambda_3}(\sigma_1) \frac{e^{i\kk_4\cdot\rr_2}}{\sqrt V}\eta_{\lambda_4}(\sigma_2) \label{eq:7-tei1}\\ &\quad= \delta_{\lambda_1\lambda_3}\delta_{\lambda_2\lambda_4}\; \frac{e^2}{V^2}\int\dd^3r_1\int\dd^3r_2\; e^{i(\kk_3-\kk_1)\cdot\rr_1}\,e^{i(\kk_4-\kk_2)\cdot\rr_2}\, \frac{e^{-\mu r_{12}}}{r_{12}} \label{eq:7-tei2} \end{align} $$

ここで $r_{12} = \abs{\rr_1-\rr_2}$.\eqref{eq:7-tei2} ではスピン和を実行した:$\sum_{\sigma_1}\eta_{\lambda_1}^*(\sigma_1)\eta_{\lambda_3}(\sigma_1) = \delta_{\lambda_1\lambda_3}$,同様に $\sum_{\sigma_2}\to\delta_{\lambda_2\lambda_4}$.クーロン相互作用はスピンに作用しないので,電子1のスピンと電子2のスピンがそれぞれ保存することが読み取れる.

次に相対座標と「重心」座標に変換する:

$$ \rr = \rr_1 - \rr_2,\qquad \bm{R} = \rr_2 \qquad\Longleftrightarrow\qquad \rr_1 = \rr+\bm{R},\quad \rr_2 = \bm{R}. $$

この変換のヤコビアンは1である($\dd^3r_1\dd^3r_2 = \dd^3r\,\dd^3R$).指数部分は

$$ (\kk_3-\kk_1)\cdot\rr_1 + (\kk_4-\kk_2)\cdot\rr_2 = (\kk_3-\kk_1)\cdot(\rr+\bm{R}) + (\kk_4-\kk_2)\cdot\bm{R} = (\kk_3-\kk_1)\cdot\rr + (\kk_3+\kk_4-\kk_1-\kk_2)\cdot\bm{R} $$

と整理される.したがって

$$ \begin{align} &\langle \kk_1\lambda_1\,\kk_2\lambda_2 |\hat{V}_2| \kk_3\lambda_3\,\kk_4\lambda_4\rangle \notag\\ &\quad= \delta_{\lambda_1\lambda_3}\delta_{\lambda_2\lambda_4}\frac{e^2}{V^2} \left[\int_V\dd^3R\;e^{i(\kk_3+\kk_4-\kk_1-\kk_2)\cdot\bm{R}}\right] \left[\int\dd^3r\;e^{i(\kk_3-\kk_1)\cdot\rr}\frac{e^{-\mu r}}{r}\right] \label{eq:7-tei3}\\ &\quad= \delta_{\lambda_1\lambda_3}\delta_{\lambda_2\lambda_4}\frac{e^2}{V^2} \cdot V\,\delta_{\kk_1+\kk_2,\;\kk_3+\kk_4} \cdot \frac{4\pi}{\mu^2+\abs{\kk_1-\kk_3}^2} \label{eq:7-tei4}\\ &\quad= \delta_{\lambda_1\lambda_3}\delta_{\lambda_2\lambda_4}\; \delta_{\kk_1+\kk_2,\;\kk_3+\kk_4}\; \frac{e^2}{V}\,\frac{4\pi}{\mu^2+\abs{\kk_1-\kk_3}^2}. \label{eq:7-tei5} \end{align} $$

\eqref{eq:7-tei4} の第1因子には平面波の直交性 \eqref{eq:7-pworth} を用いた:$\int_V\dd^3R\,e^{i\bm{G}\cdot\bm{R}} = V\delta_{\bm{G},\bm{0}}$.第2因子はフーリエ変換 \eqref{eq:7-yukawa-ft} を $\bm{q} = \kk_1-\kk_3$ として適用したものである(指数の符号が $e^{i(\kk_3-\kk_1)\cdot\rr} = e^{-i(\kk_1-\kk_3)\cdot\rr}$ と一致することを確認せよ).相対座標の積分範囲を全空間に広げてよい理由は7.2.4と同じ($\mu L\gg1$)である.

結果を書き下すと

$$ \begin{equation} \langle \kk_1\lambda_1\,\kk_2\lambda_2 |\hat{V}_2| \kk_3\lambda_3\,\kk_4\lambda_4\rangle = \delta_{\lambda_1\lambda_3}\delta_{\lambda_2\lambda_4}\, \delta_{\kk_1+\kk_2,\,\kk_3+\kk_4}\, \frac{e^2}{V}\,\frac{4\pi}{\mu^2+\abs{\kk_1-\kk_3}^2} \label{eq:7-tei-final} \end{equation} $$

である.3種類のクロネッカーのデルタが現れた.それぞれの意味は次の通りである.

物理的意味:3つのデルタが語ること

相互作用の強さ $4\pi e^2/(\mu^2+q^2)$ は移行運動量 $\bm{q}=\kk_1-\kk_3$ のみで決まる.$q$ が小さいほど($\mu\to0$ で $4\pi e^2/q^2$)相互作用は強い.これは実空間でのクーロン力が長距離であることの $\kk$ 空間での表現である.

7.3.4 移行運動量による書き換え

\eqref{eq:7-tei-final} を \eqref{eq:7-twobody-general} に代入する.$\delta_{\lambda_1\lambda_3}\delta_{\lambda_2\lambda_4}$ によりスピン添字は2種類($\lambda_1,\lambda_2$)に減り,$\delta_{\kk_1+\kk_2,\kk_3+\kk_4}$ により波数の独立な自由度は4個から3個に減る.そこで新しい変数を

$$ \begin{equation} \kk_3 = \kk,\qquad \kk_4 = \bm{p},\qquad \kk_1 = \kk+\bm{q},\qquad \kk_2 = \bm{p}-\bm{q} \label{eq:7-relabel} \end{equation} $$

と取り直す.この置き換えは運動量保存則を自動的に満たしている:$\kk_1+\kk_2 = (\kk+\bm{q})+(\bm{p}-\bm{q}) = \kk+\bm{p} = \kk_3+\kk_4$.また移行運動量は $\kk_1-\kk_3 = \bm{q}$ となる.$(\kk_1,\kk_2,\kk_3,\kk_4)$ のうち保存則を満たす組と $(\kk,\bm{p},\bm{q})$ とは1対1に対応するので,和の取り替えは正当である.

2体項は次のようになる:

$$ \begin{equation} \hat{V} = \frac{e^2}{2V}\sum_{\kk\bm{p}\bm{q}}\sum_{\lambda_1\lambda_2} \frac{4\pi}{\mu^2+q^2}\; a_{(\kk+\bm{q})\lambda_1}^\dagger\,a_{(\bm{p}-\bm{q})\lambda_2}^\dagger\, a_{\bm{p}\lambda_2}\,a_{\kk\lambda_1}. \label{eq:7-V-2q} \end{equation} $$

演算子の並び順に注意しよう.第6章の規約 $a_\mu^\dagger a_\nu^\dagger a_\delta a_\gamma$ において,$\mu\leftrightarrow\gamma$($\kk_1\leftrightarrow\kk_3$)と $\nu\leftrightarrow\delta$($\kk_2\leftrightarrow\kk_4$)が対応する.すなわち消滅演算子は $a_{\kk_4\lambda_4}a_{\kk_3\lambda_3} = a_{\bm{p}\lambda_2}a_{\kk\lambda_1}$ の順に並ぶ.この「入れ子」の順序を間違えると交換項の符号が反転してしまうので,機械的にでも守る必要がある.

k, λ₁ p, λ₂ k + q, λ₁ p − q, λ₂ q 4πe²/q² 3つの保存則 ・運動量:(k+q)+(p−q) = k+p ・各電子のスピン:λ₁, λ₂ は不変 ・電子数:a† が2個,a が2個 相互作用の強さは q だけで決まる (q → 0 で発散 = 長距離性) 2電子が運動量 q をやりとりする過程
図7.3 ジェリウムの相互作用項が表す素過程.波数 $\kk$,$\bm{p}$ の2電子が運動量 $\bm{q}$ をやりとりして $\kk+\bm{q}$,$\bm{p}-\bm{q}$ に移る.$\bm{q}$ は移行運動量と呼ばれ,相互作用の強さ $4\pi e^2/q^2$ を決める唯一の変数である.

以上をまとめると,ジェリウムの電子部分のハミルトニアンは

$$ \begin{equation} \hat{H}_{\rm el} = \sum_{\kk\lambda}\frac{\hbar^2\abs{\kk}^2}{2m}a_{\kk\lambda}^\dagger a_{\kk\lambda} + \frac{e^2}{2V}\sum_{\kk\bm{p}\bm{q}}\sum_{\lambda_1\lambda_2} \frac{4\pi}{\mu^2+q^2}\, a_{(\kk+\bm{q})\lambda_1}^\dagger a_{(\bm{p}-\bm{q})\lambda_2}^\dagger a_{\bm{p}\lambda_2}a_{\kk\lambda_1} \label{eq:7-Hel-2q} \end{equation} $$

と書ける.何が得られたか.第一量子化での多重積分は,平面波基底を選んだことによってすべて解析的に実行され,残ったのは生成・消滅演算子の代数と,$\kk$,$\bm{p}$,$\bm{q}$ についての和だけになった.次節では,この和のうち $\bm{q}=\bm{0}$ の項を取り出し,それが7.2節で計算した背景項と厳密に相殺することを示す.

7.4 $\bm{q}=\bm{0}$ 項と背景電荷の相殺

7.3節で得た相互作用項 \eqref{eq:7-V-2q} の和には $\bm{q}=\bm{0}$ の項が含まれており,そこでは相互作用の強さが $4\pi/\mu^2$ となって $\mu\to0$ で発散する.本節では,この項を和から分離して演算子の代数で整理し,それが7.2節で求めた背景項 \eqref{eq:7-Hb-plus-Helb} と厳密に相殺することを示す.相殺後に残るのは $V\to\infty$ で消える項だけであり,その結果として有限のハミルトニアンが得られる.

7.4.1 $\bm{q}=\bm{0}$ 項の分離

まず和を2つに分ける.以下,$\sum'$ は $\bm{q}=\bm{0}$ を除外した和を表す記号とする:

$$ \begin{equation} \hat{V} = \underbrace{\frac{e^2}{2V}\frac{4\pi}{\mu^2} \sum_{\kk\bm{p}}\sum_{\lambda_1\lambda_2} a_{\kk\lambda_1}^\dagger a_{\bm{p}\lambda_2}^\dagger a_{\bm{p}\lambda_2}a_{\kk\lambda_1}}_{\displaystyle \hat{V}_{\bm{q}=\bm{0}}} \;+\; \underbrace{\frac{e^2}{2V}\sum_{\kk\bm{p}\bm{q}}{}'\sum_{\lambda_1\lambda_2} \frac{4\pi}{\mu^2+q^2}\, a_{(\kk+\bm{q})\lambda_1}^\dagger a_{(\bm{p}-\bm{q})\lambda_2}^\dagger a_{\bm{p}\lambda_2}a_{\kk\lambda_1}}_{\displaystyle \hat{V}'}. \label{eq:7-Vsplit} \end{equation} $$

第1項では $\bm{q}=\bm{0}$ を代入した結果,生成演算子の添字が $\kk+\bm{0}=\kk$,$\bm{p}-\bm{0}=\bm{p}$ となって消滅演算子の添字と一致していることに注意しよう.$\hat{V}'$ の方は $q\neq0$ なので $\mu\to0$ の極限が問題なくとれ,$4\pi/q^2$ となる.以下では $\hat{V}_{\bm{q}=\bm{0}}$ に含まれる演算子の積を整理する.

7.4.2 演算子の並べ替え:$\hat{N}^2-\hat{N}$

導出:$\sum_{\kk\bm{p}\lambda_1\lambda_2} a_{\kk\lambda_1}^\dagger a_{\bm{p}\lambda_2}^\dagger a_{\bm{p}\lambda_2}a_{\kk\lambda_1} = \hat{N}^2-\hat{N}$

使う道具は第6章の反交換関係 $\{a_\mu,a_\nu\}=0$ と $\{a_\mu,a_\nu^\dagger\}=\delta_{\mu\nu}$ の2つだけである.中央の2つの消滅演算子を入れ替えることから始める.

$$ \begin{align} a_{\kk\lambda_1}^\dagger a_{\bm{p}\lambda_2}^\dagger a_{\bm{p}\lambda_2}a_{\kk\lambda_1} &= -\,a_{\kk\lambda_1}^\dagger a_{\bm{p}\lambda_2}^\dagger\, a_{\kk\lambda_1}a_{\bm{p}\lambda_2} \label{eq:7-q0a}\\ &= -\,a_{\kk\lambda_1}^\dagger \Big(\delta_{\kk\bm{p}}\delta_{\lambda_1\lambda_2} - a_{\kk\lambda_1}a_{\bm{p}\lambda_2}^\dagger\Big) a_{\bm{p}\lambda_2} \label{eq:7-q0b}\\ &= -\,\delta_{\kk\bm{p}}\delta_{\lambda_1\lambda_2}\,a_{\kk\lambda_1}^\dagger a_{\bm{p}\lambda_2} \;+\; a_{\kk\lambda_1}^\dagger a_{\kk\lambda_1}\,a_{\bm{p}\lambda_2}^\dagger a_{\bm{p}\lambda_2} \label{eq:7-q0c}\\ &= -\,\delta_{\kk\bm{p}}\delta_{\lambda_1\lambda_2}\,\hat{n}_{\kk\lambda_1} \;+\; \hat{n}_{\kk\lambda_1}\,\hat{n}_{\bm{p}\lambda_2}. \label{eq:7-q0d} \end{align} $$

\eqref{eq:7-q0a}:消滅演算子どうしの反交換関係 $a_{\bm{p}\lambda_2}a_{\kk\lambda_1} = -a_{\kk\lambda_1}a_{\bm{p}\lambda_2}$ を用いた.
\eqref{eq:7-q0b}:$a_{\bm{p}\lambda_2}^\dagger a_{\kk\lambda_1} = \delta_{\kk\bm{p}}\delta_{\lambda_1\lambda_2} - a_{\kk\lambda_1}a_{\bm{p}\lambda_2}^\dagger$(基本の反交換関係を移項したもの)を,真ん中の $a_{\bm{p}\lambda_2}^\dagger a_{\kk\lambda_1}$ に適用した.
\eqref{eq:7-q0c}:括弧を展開し,第2項では $-(-1)=+1$ とした.
\eqref{eq:7-q0d}:数演算子 $\hat{n}_{\kk\lambda}=a_{\kk\lambda}^\dagger a_{\kk\lambda}$ の定義を使い,第1項ではデルタ関数により $\bm{p}\lambda_2\to\kk\lambda_1$ とした.

両辺を $\kk,\bm{p},\lambda_1,\lambda_2$ について和をとる.第2項は $\big(\sum_{\kk\lambda_1}\hat{n}_{\kk\lambda_1}\big)\big(\sum_{\bm{p}\lambda_2}\hat{n}_{\bm{p}\lambda_2}\big) = \hat{N}^2$ と因数分解でき,第1項はデルタ関数のせいで二重和が一重和に潰れて $-\sum_{\kk\lambda_1}\hat{n}_{\kk\lambda_1} = -\hat{N}$ となる.ここで $\hat{N} = \sum_{\kk\lambda}\hat{n}_{\kk\lambda}$ は粒子数演算子である(第6章).よって

$$ \sum_{\kk\bm{p}}\sum_{\lambda_1\lambda_2} a_{\kk\lambda_1}^\dagger a_{\bm{p}\lambda_2}^\dagger a_{\bm{p}\lambda_2}a_{\kk\lambda_1} = \hat{N}^2 - \hat{N} = \hat{N}(\hat{N}-1). $$

$\hat{N}(\hat{N}-1)$ という組合せは自然である.$N$ 個の電子から異なる2個を選ぶ「順序付きの対」の個数は $N(N-1)$ であり,$\frac12 N(N-1)$ が非順序対の個数である.$\hat{V}_{\bm{q}=\bm{0}}$ の前に $\frac12$ が付いていることを思い出せば,$\bm{q}=\bm{0}$ 項は「異なる電子対すべてに対する一様な($\bm{q}=\bm{0}$ の)クーロン反発の和」を正しく数えていることになる.$-\hat{N}$ は自己相互作用を除く役割を果たしている(第6章で見たように,第二量子化では自己相互作用は反交換関係によって自動的に落ちる).

電子数が $N$ に固定された部分空間では $\hat{N}\to N$ と置き換えてよいから,

$$ \begin{equation} \hat{V}_{\bm{q}=\bm{0}} = \frac{e^2}{2V}\,\frac{4\pi}{\mu^2}\,\big(N^2-N\big). \label{eq:7-Vq0} \end{equation} $$

7.4.3 相殺

いよいよ相殺を実行する.7.2節の結果 \eqref{eq:7-Hb-plus-Helb} と \eqref{eq:7-Vq0} を足し合わせる.

$$ \begin{align} \hat{V}_{\bm{q}=\bm{0}} + \hat{H}_{\rm b} + \hat{H}_{\rm el\text{-}b} &= \frac{e^2}{2V}\frac{4\pi}{\mu^2}\big(N^2-N\big) \;-\;\frac{1}{2}e^2\frac{N^2}{V}\frac{4\pi}{\mu^2} \label{eq:7-canc1}\\ &= \frac{e^2}{2V}\frac{4\pi}{\mu^2}\Big[\big(N^2-N\big) - N^2\Big] \label{eq:7-canc2}\\ &= -\,\frac{e^2}{2V}\frac{4\pi}{\mu^2}\,N \;=\; -\,\frac{2\pi e^2 N}{\mu^2 V}. \label{eq:7-canc3} \end{align} $$

\eqref{eq:7-canc2} では共通因子 $\frac{e^2}{2V}\frac{4\pi}{\mu^2}$ をくくり出した.\eqref{eq:7-canc3} で $N^2$ が完全に消えた.これが長距離発散の相殺である.

残った項を1電子あたりに直すと

$$ \begin{equation} \frac{1}{N}\Big(\hat{V}_{\bm{q}=\bm{0}} + \hat{H}_{\rm b} + \hat{H}_{\rm el\text{-}b}\Big) = -\,\frac{2\pi e^2}{\mu^2 V} \label{eq:7-residual} \end{equation} $$

である.ここで7.2.2の数学ノートで宣言した極限の順序が効く.$\mu$ を固定したまま $V\to\infty$($N/V$ 一定)とすれば,\eqref{eq:7-residual} は $1/V\to0$ によってゼロに収束する.そのあとで $\mu\to0$ としても何も起こらない.したがって

$$ \lim_{\mu\to0}\ \lim_{V\to\infty}\ \left(-\frac{2\pi e^2}{\mu^2 V}\right) = 0 . $$

もし順序を逆にして $\mu\to0$ を先にすれば $-2\pi e^2/(\mu^2 V)\to-\infty$ となり,意味のある結果は得られない.同じ式が極限の順序によって $0$ にも $-\infty$ にもなるという事実は,この手続きが単なる形式ではなく物理的な内容を持つことを示している.

物理的意味:残った項は「1電子の自己エネルギー」である

相殺し残った \eqref{eq:7-residual} の起源をたどると,それは \eqref{eq:7-Vq0} の $-N$ の項,すなわち自己相互作用を除くための補正であった.1電子あたりの大きさ $\frac{1}{2}e^2\cdot\frac{1}{V}\cdot\frac{4\pi}{\mu^2}$ は,「電子1個分の電荷を体積 $V$ 全体に一様に塗り広げたときの静電自己エネルギー」に等しい.電荷が広がる体積が大きくなればこの自己エネルギーは薄まってゼロになる——これが $V\to\infty$ で消える理由である.

逆に言えば,有限の箱の計算(実際の第一原理計算はすべて有限の周期セルで行われる)では,この項は $O(1/V)$ の有限サイズ効果として残る.周期系のDFT計算で「マーデルング補正」「ジェリウム背景の補正」と呼ばれる操作は,まさにこの種の項を扱っている(第14章).

7.4.4 ジェリウムの最終ハミルトニアン

以上をまとめる.$\mu\to0$ の極限では $\hat{V}'$ の相互作用の強さは $4\pi/q^2$ になり,$\bm{q}=\bm{0}$ を除く限りこれは有限である.したがってジェリウムのハミルトニアンは次の形に確定する.

$$ \begin{equation} \hat{H} = \sum_{\kk\lambda}\frac{\hbar^2\abs{\kk}^2}{2m}\,a_{\kk\lambda}^\dagger a_{\kk\lambda} \;+\;\frac{e^2}{2V}\sum_{\kk\bm{p}\bm{q}}{}'\;\sum_{\lambda_1\lambda_2} \frac{4\pi}{q^2}\; a_{(\kk+\bm{q})\lambda_1}^\dagger\,a_{(\bm{p}-\bm{q})\lambda_2}^\dagger\, a_{\bm{p}\lambda_2}\,a_{\kk\lambda_1} \label{eq:7-H-final} \end{equation} $$

ここで $\sum'$ は $\bm{q}=\bm{0}$ を除いた和である.背景電荷は,この $\sum'$ のプライム記号ただ一つの中に姿を変えて残っている.式の見た目からは背景が消えてしまったように見えるが,$\bm{q}=\bm{0}$ を除くという操作こそが背景電荷の効果そのものである.

物理的意味:$\bm{q}=\bm{0}$ を除くとはどういうことか

$\bm{q}=\bm{0}$ のフーリエ成分とは,空間的に一様な成分のことである.電子の密度分布 $\rho(\rr)$ をフーリエ分解したとき,$\bm{q}=\bm{0}$ 成分は平均密度 $N/V$ を,$\bm{q}\neq\bm{0}$ 成分は平均からのゆらぎを表す.ジェリウムでは電子の平均密度と背景の密度がぴったり等しいので,両者が作る一様な静電ポテンシャルは互いに打ち消し合い,電子は密度ゆらぎどうしの相互作用だけを感じる.\eqref{eq:7-H-final} の $\sum'$ はまさにそれを表している.

この見方は,ジェリウムを「一様な部分(古典的静電エネルギー,中性ゆえにゼロ)+ ゆらぎの部分(量子多体効果)」に分解したことに相当する.以後の計算で得られる交換エネルギーも相関エネルギーも,すべてゆらぎが生む効果である.DFTで交換相関エネルギーが「ハートリー項を差し引いた残り」として定義されるのも同じ構造であり(第10章),ジェリウムはその最も透明な例になっている.

例:$\bm{q}=\bm{0}$ を除かないとどうなるか

仮に背景電荷を忘れて $\bm{q}=\bm{0}$ を含めたまま計算すると,\eqref{eq:7-Vq0} により全エネルギーに $\frac{e^2}{2V}\frac{4\pi}{\mu^2}N^2$ が加わる.1電子あたりでは $\frac{1}{2}e^2\frac{N}{V}\frac{4\pi}{\mu^2}$ であり,密度 $N/V$ を一定に保ったまま $\mu\to0$ とすると発散する.数値で見ると,$r_s=4$ の電子ガス($N/V = 3/(4\pi\cdot4^3) = 3.73\times10^{-3}\ a_0^{-3}$)で $1/\mu = 100\,a_0$(まだ十分「短距離」)としても,この項は1電子あたり $\frac12\times 3.73\times10^{-3}\times 4\pi\times 10^4 \approx 234$ Ha $\approx 6.4\times10^{3}$ eV となる.運動エネルギー($0.069$ Ha)や交換エネルギー($-0.115$ Ha)とは比較にならない大きさである.中性条件を無視した電子ガスの計算は,物理的に意味のある量が完全に埋もれてしまうことがこの数値から分かる.

7.5 無次元化と高密度極限

最終ハミルトニアン \eqref{eq:7-H-final} は厳密であるが,そのままでは解けない.そこで摂動論に持ち込みたいのだが,「何が小さいのか」がまだ見えていない.本節では,長さの単位を電子1個あたりの半径 $r_0$ に取り替える(無次元化する)ことで,運動エネルギー項と相互作用項が $r_s$ の異なるべきでスケールすること,したがって高密度極限 $r_s\to0$ で相互作用が摂動になることを明らかにする.

7.5.1 無次元変数の導入

長さの単位として $r_0$($=r_s a_0$,電子1個あたりの球の半径 \eqref{eq:7-rs-def})を選び,次の無次元量を導入する:

$$ \begin{equation} \bar{\rr} = \frac{\rr}{r_0},\qquad \bar{\kk} = r_0\,\kk,\qquad \bar{\bm{q}} = r_0\,\bm{q},\qquad \bar{V} = \frac{V}{r_0^3}. \label{eq:7-dimless} \end{equation} $$

波数を $r_0$ で掛けるのは,波数が長さの逆数の次元を持つからである($\kk\cdot\rr = \bar{\kk}\cdot\bar{\rr}$ が無次元になるようにする).$\bar V = V/r_0^3$ は「電子1個あたりの球いくつ分の体積か」を表す量で,定義 \eqref{eq:7-rs-def} より

$$ \bar{V} = \frac{V}{r_0^3} = \frac{N\cdot\frac{4}{3}\pi r_0^3}{r_0^3} = \frac{4\pi}{3}N $$

である.つまり $\bar V$ は $N$ に比例する純粋な数であり,$r_s$ には依存しない.ここが重要な点である:無次元化したあとで $r_s$ に依存するのは係数だけになる.

7.5.2 運動エネルギー項のスケーリング

導出:$\hat{T}$ の無次元化

$$ \begin{align} \hat{T} &= \sum_{\kk\lambda}\frac{\hbar^2\abs{\kk}^2}{2m}\,a_{\kk\lambda}^\dagger a_{\kk\lambda} \label{eq:7-Tscale1}\\ &= \frac{\hbar^2}{2m r_0^2}\sum_{\bar{\kk}\lambda}\abs{\bar{\kk}}^2\,a_{\bar\kk\lambda}^\dagger a_{\bar\kk\lambda} \label{eq:7-Tscale2}\\ &= \frac{e^2 a_0}{2 r_0^2}\sum_{\bar{\kk}\lambda}\abs{\bar{\kk}}^2\,a_{\bar\kk\lambda}^\dagger a_{\bar\kk\lambda} \label{eq:7-Tscale3}\\ &= \frac{e^2}{a_0 r_s^2}\sum_{\bar{\kk}\lambda}\frac{1}{2}\abs{\bar{\kk}}^2\,a_{\bar\kk\lambda}^\dagger a_{\bar\kk\lambda}. \label{eq:7-Tscale4} \end{align} $$

\eqref{eq:7-Tscale2}:$\abs{\kk}^2 = \abs{\bar\kk}^2/r_0^2$ を代入し,$r_0$ に依存する因子を和の外に出した.
\eqref{eq:7-Tscale3}:ボーア半径の定義 $a_0=\hbar^2/(me^2)$ を $\hbar^2/m = e^2a_0$ と読み替えて代入した.
\eqref{eq:7-Tscale4}:$r_0 = r_s a_0$ より $r_0^2 = r_s^2a_0^2$ であるから $\dfrac{e^2a_0}{2r_0^2} = \dfrac{e^2a_0}{2r_s^2a_0^2} = \dfrac{e^2}{a_0r_s^2}\cdot\dfrac{1}{2}$.$\frac12$ を和の中に入れた.

7.5.3 相互作用項のスケーリング

導出:$\hat{V}'$ の無次元化

$$ \begin{align} \hat{V}' &= \frac{e^2}{2V}\sum_{\kk\bm{p}\bm{q}}{}'\sum_{\lambda_1\lambda_2}\frac{4\pi}{q^2}\, a^\dagger a^\dagger a a \label{eq:7-Vscale1}\\ &= \frac{e^2}{2\bar V r_0^3}\sum{}'\frac{4\pi r_0^2}{\bar q^2}\,a^\dagger a^\dagger a a \label{eq:7-Vscale2}\\ &= \frac{e^2}{2\bar V r_0}\sum{}'\frac{4\pi}{\bar q^2}\,a^\dagger a^\dagger a a \label{eq:7-Vscale3}\\ &= \frac{e^2}{a_0 r_s^2}\cdot\frac{r_s}{2\bar V}\sum{}'\frac{4\pi}{\bar q^2}\,a^\dagger a^\dagger a a . \label{eq:7-Vscale4} \end{align} $$

\eqref{eq:7-Vscale2}:$V = \bar V r_0^3$ と $q^2 = \bar q^2/r_0^2$(すなわち $1/q^2 = r_0^2/\bar q^2$)を代入した.
\eqref{eq:7-Vscale3}:$r_0^2/r_0^3 = 1/r_0$ と約分した.
\eqref{eq:7-Vscale4}:$\dfrac{e^2}{r_0} = \dfrac{e^2}{a_0 r_s} = \dfrac{e^2}{a_0r_s^2}\cdot r_s$ と書き換えた(分母分子に $r_s$ を掛けた).

2つを合わせると,ジェリウムのハミルトニアンは次のように書ける.

$$ \begin{equation} \hat{H} = \frac{e^2}{a_0 r_s^2}\left[ \sum_{\bar\kk\lambda}\frac{1}{2}\abs{\bar\kk}^2\,a_{\bar\kk\lambda}^\dagger a_{\bar\kk\lambda} \;+\;\frac{r_s}{2\bar V}\sum_{\bar\kk\bar{\bm{p}}\bar{\bm{q}}}{}'\sum_{\lambda_1\lambda_2} \frac{4\pi}{\bar q^{\,2}}\, a_{(\bar\kk+\bar{\bm{q}})\lambda_1}^\dagger a_{(\bar{\bm{p}}-\bar{\bm{q}})\lambda_2}^\dagger a_{\bar{\bm{p}}\lambda_2}a_{\bar\kk\lambda_1} \right] \label{eq:7-H-dimless} \end{equation} $$

何が得られたか.角括弧の中に $r_s$ が現れるのは,相互作用項の前の因子 $r_s$ ただ一つだけである(和の範囲や $\bar V$ は $r_s$ に依存しない).すなわち$r_s$ が結合定数の役割を果たす.全体に掛かる $e^2/(a_0r_s^2) = \mathrm{Ha}/r_s^2 = 2\,\mathrm{Ry}/r_s^2$ は単なるエネルギーの目盛りである.

7.5.4 高密度極限で摂動論が正当化される理由

物理的意味:密度が高いほど「弱く相互作用する」

\eqref{eq:7-H-dimless} が主張しているのは次のことである.エネルギーの目盛りを $e^2/(a_0r_s^2)$ にとると,運動エネルギー項の大きさは $r_s$ によらず $O(1)$,相互作用項の大きさは $O(r_s)$ である.したがって

$$ \frac{\text{相互作用エネルギー}}{\text{運動エネルギー}} \sim r_s . $$

$r_s\to0$(高密度)で相互作用は運動エネルギーに比べて相対的に弱くなるから,相互作用項を摂動として扱うことが正当化される.これは次元解析だけでも見て取れる:典型的な長さが $r_0$ の系では,量子的な運動エネルギーは不確定性関係より $\sim\hbar^2/(mr_0^2)$,クーロンエネルギーは $\sim e^2/r_0$ である.両者の比は

$$ \frac{e^2/r_0}{\hbar^2/(mr_0^2)} = \frac{m e^2 r_0}{\hbar^2} = \frac{r_0}{a_0} = r_s . $$

運動エネルギーは長さの $-2$ 乗,クーロンは $-1$ 乗でスケールするので,系を圧縮すると運動エネルギーの方が速く増える.これが「高密度=弱結合」の正体である.

この結論は古典的な直観に反する.古典気体では,密度を上げれば分子どうしが近づいて相互作用は強くなる.電子ガスでこれが逆転するのは,パウリの排他律に由来する量子的な運動エネルギー(第3章)が圧縮に対して非常に強く反発するからである.逆に $r_s\to\infty$(低密度)では相互作用が支配的になり,電子は運動エネルギーを犠牲にしても互いを避けようとする.その極限で電子が結晶を組む——これが7.9節で扱うウィグナー結晶である.

摂動論の枠組みは次のようになる.\eqref{eq:7-H-dimless} を $\hat H = \hat H_0 + \hat H_1$ と分け,

$$ \begin{equation} \hat{H}_0 = \sum_{\kk\lambda}\frac{\hbar^2\abs{\kk}^2}{2m}a_{\kk\lambda}^\dagger a_{\kk\lambda} \ \ (\text{$O(r_s^{-2})$}), \qquad \hat{H}_1 = \frac{e^2}{2V}\sum_{\kk\bm{p}\bm{q}}{}'\sum_{\lambda_1\lambda_2} \frac{4\pi}{q^2}a^\dagger a^\dagger a a \ \ (\text{$O(r_s^{-1})$}) \label{eq:7-H0H1} \end{equation} $$

とする.$\hat H_0$ は自由電子ガス(第3章)そのもので厳密に解ける.エネルギーを摂動展開すると

$$ \begin{equation} \frac{E}{N} = \underbrace{\frac{E_0}{N}}_{O(r_s^{-2})} + \underbrace{\frac{E_1}{N}}_{O(r_s^{-1})} + \underbrace{\frac{E_2}{N} + \cdots}_{O(r_s^{0}\ln r_s)\ \text{以下}} \label{eq:7-Eexpand} \end{equation} $$

という $r_s$ のべき(と対数)による展開が得られる.7.6節で $E_0$,7.7節で $E_1$ を厳密に計算し,7.9節で $E_2$ 以降(相関エネルギー)を議論する.

例:実在金属は「摂動論の境界領域」にある

第3章の表3.1で見たように,実在の単体金属は $2\lesssim r_s\lesssim6$ の範囲にある.摂動パラメータが $2$ から $6$ ということは,これは小さくない.厳密に言えば,高密度展開 \eqref{eq:7-Eexpand} が収束する保証はこの領域にはない.

それにもかかわらず,最初の2項($E_0$ と $E_1$)が実在金属のエネルギーのおおよその大きさを与えるのは,次の2つの事情による.第一に,$r_s\approx4$ でも交換エネルギー($-0.229$ Ry)は運動エネルギー($0.138$ Ry)と同程度であって,まだ「桁が違う」わけではない.第二に,3項目以降(相関エネルギー,$-0.064$ Ry)が交換エネルギーの3割程度に留まる.すなわち展開は「ゆっくり収束しているように見える」.しかし相関エネルギーが化学的精度(1 kcal/mol $\approx 0.0016$ Ha)を左右する場面では,この $-0.064$ Ry $=-0.032$ Ha を精密に知ることが決定的になる.7.9節で見るように,この量を高精度で求めるには摂動論を超えた方法(量子モンテカルロ法)が必要であった.

7.6 摂動の0次:運動エネルギー

本節では \eqref{eq:7-Eexpand} の第1項 $E_0$ を計算する.これは第3章で既に得た結果であるが,ここでは第二量子化の言葉で再導出し,第6章の道具が正しく働くことを確認する.同時に,7.7節の交換エネルギー計算で使う記法と積分技術を準備する.

7.6.1 無摂動基底状態:フェルミ球

$\hat H_0$ は数演算子の重み付き和 \eqref{eq:7-T-2q} であるから,その固有状態は占有数状態そのものであり,固有値は $\sum_{\kk\lambda}(\hbar^2k^2/2m)n_{\kk\lambda}$ である.これを最小にするには,エネルギーの低い状態($\abs{\kk}$ の小さい状態)から順に,パウリの排他律($n_{\kk\lambda}\in\{0,1\}$)が許す限り詰めればよい.結果は第3章のフェルミ球である.

定義:フェルミ球状態 $|\Phi_0\rangle$

$$ \begin{equation} |\Phi_0\rangle = \prod_{\abs{\kk}\le k_F}\prod_{\lambda=\uparrow,\downarrow} a_{\kk\lambda}^\dagger\,|0\rangle, \qquad n_{\kk\lambda} = \langle\Phi_0|\hat{n}_{\kk\lambda}|\Phi_0\rangle = \begin{cases} 1 & (\abs{\kk}\le k_F)\\ 0 & (\abs{\kk} > k_F) \end{cases} \label{eq:7-fermisphere} \end{equation} $$

フェルミ波数 $k_F$ は電子数の条件 $\sum_{\kk\lambda}n_{\kk\lambda}=N$ から決まり,第3章の結果 $k_F=(3\pi^2\rho)^{1/3}$ を与える.

7.6.2 $E_0$ の評価

導出:$E_0 = \langle\Phi_0|\hat{H}_0|\Phi_0\rangle$

$$ \begin{align} E_0 &= \langle\Phi_0|\sum_{\kk\lambda}\frac{\hbar^2\abs{\kk}^2}{2m}\hat{n}_{\kk\lambda}|\Phi_0\rangle \label{eq:7-E0a}\\ &= \sum_{\kk\lambda}\frac{\hbar^2\abs{\kk}^2}{2m}\,n_{\kk\lambda} \label{eq:7-E0b}\\ &= 2\sum_{\abs{\kk}\le k_F}\frac{\hbar^2\abs{\kk}^2}{2m} \label{eq:7-E0c}\\ &= 2\cdot\frac{V}{(2\pi)^3}\int_{\abs{\kk}\le k_F}\dd^3k\;\frac{\hbar^2 k^2}{2m} \label{eq:7-E0d}\\ &= \frac{2V}{8\pi^3}\cdot\frac{\hbar^2}{2m}\cdot 4\pi\int_0^{k_F}\dd k\; k^2\cdot k^2 \label{eq:7-E0e}\\ &= \frac{V\hbar^2}{2\pi^2 m}\int_0^{k_F}k^4\,\dd k \;=\;\frac{V\hbar^2}{2\pi^2 m}\cdot\frac{k_F^5}{5} \;=\;\frac{V\hbar^2 k_F^5}{10\pi^2 m}. \label{eq:7-E0f} \end{align} $$

\eqref{eq:7-E0b}:数演算子の期待値は占有数である(第6章).
\eqref{eq:7-E0c}:占有数はフェルミ球の内側で1,外側で0.スピン和は2種類とも同じ寄与を与えるので因子2になる.
\eqref{eq:7-E0d}:和を積分に置き換える規則 \eqref{eq:7-sum2int} を適用した.
\eqref{eq:7-E0e}:球座標に移り,被積分関数が角度によらないので立体角の積分 $\int\dd\Omega = 4\pi$ を実行した.体積要素の $k^2$ と被積分関数の $k^2$ が掛かって $k^4$ になる.
\eqref{eq:7-E0f}:$\dfrac{2\cdot4\pi}{8\pi^3} = \dfrac{1}{\pi^2}$ と約分し,$\int_0^{k_F}k^4\dd k = k_F^5/5$ を代入した.

1電子あたりに直す.$N = Vk_F^3/(3\pi^2)$(第3章)を使うと

$$ \begin{align} \frac{E_0}{N} &= \frac{V\hbar^2k_F^5}{10\pi^2 m}\cdot\frac{3\pi^2}{Vk_F^3} = \frac{3\hbar^2k_F^2}{10m} = \frac{3}{5}\cdot\frac{\hbar^2k_F^2}{2m} = \frac{3}{5}\varepsilon_F . \label{eq:7-E0-perN} \end{align} $$

第3章の結果が再現された.$r_s$ で書き直すには \eqref{eq:7-kF-rs} を用いる.原子単位($\hbar=m=e=1$,長さは $a_0$)では $\varepsilon_F = k_F^2/2$ であるから

$$ \begin{align} \frac{E_0}{N} &= \frac{3}{10}k_F^2 = \frac{3}{10}\left(\frac{9\pi}{4}\right)^{2/3}\frac{1}{r_s^2} \label{eq:7-E0-rs1}\\ &= \frac{3}{10}\times 3.683169\times\frac{1}{r_s^2} = \frac{1.10495}{r_s^2}\ \mathrm{Ha} = \frac{2.20990}{r_s^2}\ \mathrm{Ry}. \label{eq:7-E0-rs2} \end{align} $$

数値は $(9\pi/4)^{1/3}=1.919158$ を2乗した $(9\pi/4)^{2/3}=3.683169$ から得た.Ry に直すには2倍する($1$ Ha $=2$ Ry).

$$ \frac{E_0}{N} = \frac{3}{5}\varepsilon_F = \frac{3}{10}\left(\frac{9\pi}{4}\right)^{2/3}\frac{1}{r_s^2} = \frac{1.105}{r_s^2}\ \mathrm{Ha} = \frac{2.21}{r_s^2}\ \mathrm{Ry} \qquad\Longleftrightarrow\qquad t(\rho) = \frac{3}{10}(3\pi^2)^{2/3}\rho^{2/3}\ \ [\mathrm{Ha}] $$

最後の等式は $k_F=(3\pi^2\rho)^{1/3}$ を代入して密度で表したものである.この $\rho^{2/3}$ という依存性が,第4章のトーマス・フェルミ模型の運動エネルギー汎関数 $T_{\rm TF}[\rho]=\frac{3}{10}(3\pi^2)^{2/3}\int\rho^{5/3}\dd^3r$ の起源であった(1電子あたりのエネルギー密度 $t(\rho)$ に密度 $\rho$ を掛けて積分するので,被積分関数は $\rho^{5/3}$ になる).

何が得られたか.摂動の0次は自由電子ガスの運動エネルギーそのものであり,$r_s^{-2}$ でスケールする.次節では,これに $r_s^{-1}$ でスケールする1次の補正——交換エネルギー——を加える.

7.7 摂動の1次:交換エネルギー

本節が本章の核心である.摂動の1次補正 $E_1 = \langle\Phi_0|\hat{H}_1|\Phi_0\rangle$ を,途中を一切省かずに計算しきる.手順は次の5段階である.(a) 反交換関係でどの縮約が生き残るかを決める,(b) 残った和を整理する,(c) 和を積分に置き換える,(d) 内側の $\bm{p}$ 積分を角度積分と部分積分で実行する,(e) 外側の $\kk$ 積分を実行する.最終的に得られるのは,DFTの歴史全体を通じて最も重要な式の一つ,ディラックの交換エネルギーである.

7.7.1 どの縮約が生き残るか

出発点は

$$ \begin{equation} E_1 = \langle\Phi_0|\hat{H}_1|\Phi_0\rangle = \frac{e^2}{2V}\sum_{\kk\bm{p}\bm{q}}{}'\sum_{\lambda_1\lambda_2}\frac{4\pi}{q^2}\, \langle\Phi_0| a_{(\kk+\bm{q})\lambda_1}^\dagger a_{(\bm{p}-\bm{q})\lambda_2}^\dagger a_{\bm{p}\lambda_2}a_{\kk\lambda_1} |\Phi_0\rangle \label{eq:7-E1-start} \end{equation} $$

である.第6章6.8節で証明したように,単一のスレーター行列式 $|\Phi_0\rangle$ に対する2体演算子の期待値は

$$ \begin{equation} \langle\Phi_0|a_\mu^\dagger a_\nu^\dagger a_\delta a_\gamma|\Phi_0\rangle = n_\gamma n_\delta\big(\delta_{\mu\gamma}\delta_{\nu\delta} - \delta_{\mu\delta}\delta_{\nu\gamma}\big) \label{eq:7-wick} \end{equation} $$

で与えられる.これは「生成演算子と消滅演算子を対にして結ぶ(縮約する)方法が2通りあり,線が交差する方に $-1$ が付く」という規則である.いま

$$ \mu = (\kk+\bm{q},\lambda_1),\qquad \nu = (\bm{p}-\bm{q},\lambda_2),\qquad \gamma = (\kk,\lambda_1),\qquad \delta = (\bm{p},\lambda_2) $$

であるから,2つの縮約はそれぞれ次の条件を要求する.

ここで決定的なことが起きる.直接縮約は $\bm{q}=\bm{0}$ を要求するが,\eqref{eq:7-E1-start} の和は $\sum'$ すなわち $\bm{q}=\bm{0}$ を除外している.したがって直接項は完全に消える.残るのは交換縮約だけである.

直接縮約 a†k+q a†p−q ap ak k+q = k かつ p−q = p ⇒ q = 0 Σ′ が q=0 を除くので寄与ゼロ 交換縮約 a†k+q a†p−q ap ak k+q = p かつ p−q = k ⇒ q = p−k λ₁ = λ₂ が必要・符号 −1 線が交差する縮約だけが生き残る:ジェリウムの1次補正は純粋な交換エネルギーである
図7.4 摂動1次に現れる2つの縮約.直接(ハートリー)項は $\bm{q}=\bm{0}$ を要求するが,それは背景電荷との相殺によって既に和から除かれている.したがって交換項だけが残る.

物理的意味:ハートリー項はどこへ行ったのか

第5章のハートリー・フォック法では,エネルギーはクーロン積分 $J$(ハートリー項)と交換積分 $K$ の差 $J-K$ という形で現れた.ジェリウムでは $J$ に相当する直接項が完全に消えている.これは $J$ が存在しないからではなく,一様な電子密度が作るハートリーポテンシャルが,一様な正背景が作るポテンシャルと厳密に打ち消し合っているからである(7.4節).実際,7.4節で相殺したのはまさに $\bm{q}=\bm{0}$ の直接項であった.

この構造は第10章のコーン・シャム法にそのまま受け継がれる.そこでも全エネルギーは「ハートリー項 + 交換相関項」に分けられ,一様系ではハートリー項が外部ポテンシャル(核の引力)と相殺して,交換相関項だけがエネルギーを決める.ジェリウムはその最も純粋な実現である.

7.7.2 和の整理

交換縮約だけを残すと,\eqref{eq:7-E1-start} は

$$ \begin{align} E_1 &= \frac{e^2}{2V}\sum_{\kk\bm{p}\bm{q}}{}'\sum_{\lambda_1\lambda_2}\frac{4\pi}{q^2}\, \Big(-\delta_{\kk+\bm{q},\,\bm{p}}\,\delta_{\bm{p}-\bm{q},\,\kk}\,\delta_{\lambda_1\lambda_2}\Big)\, n_{\kk\lambda_1}n_{\bm{p}\lambda_2} \label{eq:7-E1a}\\ &= -\frac{e^2}{2V}\sum_{\lambda}\sum_{\kk\bm{p}}{}^{(\kk\neq\bm{p})}\; \frac{4\pi}{\abs{\bm{p}-\kk}^2}\;n_{\kk\lambda}\,n_{\bm{p}\lambda} \label{eq:7-E1b}\\ &= -\frac{e^2}{V}\sum_{\abs{\kk}\le k_F}\ \sum_{\abs{\bm{p}}\le k_F}\; \frac{4\pi}{\abs{\kk-\bm{p}}^2}. \label{eq:7-E1c} \end{align} $$

\eqref{eq:7-E1a}:\eqref{eq:7-wick} の交換項を代入した.2つの波数のデルタは同じ条件 $\bm{q}=\bm{p}-\kk$ を表すので,実質的に1つとして働く.
\eqref{eq:7-E1b}:$\bm{q}$ についての和をデルタ関数で実行し,$\bm{q}=\bm{p}-\kk$ を代入した.$\sum'$ の制約 $\bm{q}\neq\bm{0}$ は $\kk\neq\bm{p}$ に翻訳される.またスピンのデルタ $\delta_{\lambda_1\lambda_2}$ により $\lambda_2$ の和が消え,$\lambda_1=\lambda_2=\lambda$ と書いた.
\eqref{eq:7-E1c}:占有数 $n_{\kk\lambda}$,$n_{\bm{p}\lambda}$ はどちらもフェルミ球の内側で1,外側で0であるから,和の範囲を $\abs{\kk}\le k_F$,$\abs{\bm{p}}\le k_F$ に制限した.スピン和 $\sum_\lambda$ は2つの値で同じ寄与を与えるので因子2となり,$-e^2/(2V)\times2 = -e^2/V$ となった.$\abs{\bm{p}-\kk}=\abs{\kk-\bm{p}}$ も使った.なお $\kk=\bm{p}$ の項は本来除外されるが,次に和を積分に直したとき,この1点は3次元積分の測度ゼロの集合であり,しかも $1/\abs{\kk-\bm{p}}^2$ の特異性は3次元では可積分($\int \dd^3s/s^2 = 4\pi\int \dd s$ が有限)なので,除外してもしなくても結果は変わらない.

物理的意味:同じスピンの電子どうしだけが交換する

\eqref{eq:7-E1b} に現れたスピンのデルタ $\delta_{\lambda_1\lambda_2}$ は決定的に重要である.交換相互作用は同じスピンを持つ電子の対にしか働かない.これは第5章で交換積分 $K_{ll'}$ がスピン積分 $\sum_\sigma\eta_l^*(\sigma)\eta_{l'}(\sigma)$ を含むために異スピン間でゼロになったことと同じ事実である.

その起源はパウリの排他律である.同じスピンの2電子は同じ場所に来られないので,互いに「避け合う」.避け合えばクーロン反発によるエネルギーの損が減る.この利得が交換エネルギーであり,必ず負である.異なるスピンの2電子にはこの制約がないので,ハートリー・フォックの近似の範囲では互いに無相関に運動する.異スピン間の避け合いを取り入れるのが相関エネルギー(7.9節)である.

7.7.3 和を積分に直す

\eqref{eq:7-E1c} の2つの和に,それぞれ規則 \eqref{eq:7-sum2int} を適用する:

$$ \begin{align} E_1 &= -\frac{4\pi e^2}{V} \left[\frac{V}{(2\pi)^3}\right]^2 \int_{\abs{\kk}\le k_F}\dd^3k\int_{\abs{\bm{p}}\le k_F}\dd^3p\;\frac{1}{\abs{\kk-\bm{p}}^2} \label{eq:7-E1int1}\\ &= -\frac{4\pi e^2}{V}\cdot\frac{V^2}{64\pi^6} \int_{\abs{\kk}\le k_F}\dd^3k\;I(\kk) \label{eq:7-E1int2}\\ &= -\frac{e^2 V}{16\pi^5}\int_{\abs{\kk}\le k_F}\dd^3k\;I(\kk), \label{eq:7-E1int3} \end{align} $$

ただし内側の積分を

$$ \begin{equation} I(\kk) \equiv \int_{\abs{\bm{p}}\le k_F}\dd^3 p\;\frac{1}{\abs{\kk-\bm{p}}^2} \label{eq:7-Ik-def} \end{equation} $$

と定義した.\eqref{eq:7-E1int2} では $[(2\pi)^3]^2 = (8\pi^3)^2 = 64\pi^6$ を使い,\eqref{eq:7-E1int3} では $\dfrac{4\pi}{64\pi^6} = \dfrac{1}{16\pi^5}$ と約分した.$V^2/V = V$ である.

kF k p q = p − k θ 両方の波数がフェルミ球の内側にある対だけが寄与する 計算の段取り 1. |k−p|² = k² + p² − 2kp cos θ と書く 2. 角度 θ の積分 → 対数関数が現れる 3. p の積分(部分積分)→ I(k) 4. k の積分(部分積分)→ E₁ 6次元積分だが,球対称性のおかげで 実質2重の1次元積分に還元される
図7.5 交換エネルギーの積分領域.$\kk$ と $\bm{p}$ はともにフェルミ球の内側を走り,被積分関数は移行運動量 $\bm{q}=\bm{p}-\kk$ の大きさだけに依存する.$\theta$ は $\kk$ と $\bm{p}$ のなす角である.

7.7.4 内側の積分 $I(\kk)$:角度積分

これから $I(\kk)$ を球座標で計算する.積分変数 $\bm{p}$ の極軸を $\kk$ の向きにとる.すると $\kk$ と $\bm{p}$ のなす角が $\theta$ となり,余弦定理により

$$ \begin{equation} \abs{\kk-\bm{p}}^2 = k^2 + p^2 - 2kp\cos\theta \label{eq:7-cosine} \end{equation} $$

と書ける($k=\abs{\kk}$,$p=\abs{\bm{p}}$).被積分関数は方位角 $\phi$ に依存しないので,$\phi$ 積分は $2\pi$ を与える.

導出:$I(\kk)$ の角度積分

$$ \begin{align} I(\kk) &= \int_0^{k_F}\dd p\,p^2\int_0^{\pi}\dd\theta\,\sin\theta\int_0^{2\pi}\dd\phi\; \frac{1}{k^2+p^2-2kp\cos\theta} \label{eq:7-Ia}\\ &= 2\pi\int_0^{k_F}\dd p\,p^2\int_0^{\pi}\frac{\sin\theta\,\dd\theta}{k^2+p^2-2kp\cos\theta} \label{eq:7-Ib}\\ &= 2\pi\int_0^{k_F}\dd p\,p^2\int_{-1}^{1}\frac{\dd u}{k^2+p^2-2kp\,u}. \label{eq:7-Ic} \end{align} $$

\eqref{eq:7-Ic} では $u=\cos\theta$ と置換した($\dd u = -\sin\theta\,\dd\theta$,積分範囲 $\theta:0\to\pi$ は $u:1\to-1$ に対応するので,符号の反転と上下端の入れ替えが打ち消し合って $\int_{-1}^{1}\dd u$ となる).

$u$ 積分を実行する.$A=k^2+p^2$,$B=2kp$ と略記すると,$\dfrac{\dd}{\dd u}\ln(A-Bu) = \dfrac{-B}{A-Bu}$ であるから

$$ \begin{align} \int_{-1}^{1}\frac{\dd u}{A-Bu} &= \left[-\frac{1}{B}\ln(A-Bu)\right]_{u=-1}^{u=1} \label{eq:7-Id}\\ &= \frac{1}{B}\Big[\ln(A+B) - \ln(A-B)\Big] \label{eq:7-Ie}\\ &= \frac{1}{2kp}\ln\frac{k^2+p^2+2kp}{k^2+p^2-2kp} \label{eq:7-If}\\ &= \frac{1}{2kp}\ln\frac{(k+p)^2}{(k-p)^2} \;=\;\frac{1}{kp}\ln\left|\frac{k+p}{k-p}\right|. \label{eq:7-Ig} \end{align} $$

\eqref{eq:7-If} では $A\pm B$ を代入した.\eqref{eq:7-Ig} では完全平方の形 $k^2+p^2\pm2kp = (k\pm p)^2$ を使い,$\ln(X^2/Y^2) = 2\ln\abs{X/Y}$ により係数 $\frac{1}{2kp}\times2 = \frac{1}{kp}$ とした.絶対値が必要なのは $k

したがって

$$ \begin{equation} I(\kk) = 2\pi\int_0^{k_F}\dd p\;p^2\cdot\frac{1}{kp}\ln\left|\frac{k+p}{k-p}\right| = \frac{2\pi}{k}\int_0^{k_F}\dd p\;p\,\ln\left|\frac{k+p}{k-p}\right|. \label{eq:7-Ih} \end{equation} $$

角度積分から対数関数が現れた.この対数こそが,以後ジェリウム理論のあらゆる場面(交換エネルギー,リントハルト関数,フリーデル振動)に顔を出す共通の構造である.

7.7.5 内側の積分 $I(\kk)$:動径積分

次に \eqref{eq:7-Ih} の $p$ 積分を部分積分で実行する.部分積分の使い方に工夫がある.

導出:$I(\kk)$ の動径積分

$L(p)\equiv\ln\left|\dfrac{k+p}{k-p}\right| = \ln\abs{k+p}-\ln\abs{k-p}$ とおく.その導関数は

$$ L'(p) = \frac{1}{k+p} - \frac{-1}{k-p} = \frac{1}{k+p}+\frac{1}{k-p} = \frac{(k-p)+(k+p)}{(k+p)(k-p)} = \frac{2k}{k^2-p^2}. $$

(第2項の微分では $\dfrac{\dd}{\dd p}\ln\abs{k-p} = \dfrac{-1}{k-p}$ に注意する.)

部分積分 $\int u\,\dd v = [uv]-\int v\,\dd u$ を $u=L(p)$,$\dd v = p\,\dd p$ として使う.ここで $v$ として単純な $p^2/2$ ではなく,定数だけずらした

$$ w(p) \equiv \frac{p^2-k^2}{2} $$

を選ぶ.$w'(p)=p$ なので $v$ の資格を満たしており,定数のずれは部分積分の公式を壊さない.この選び方の利点は,$p=0$ で $L(0)=\ln\abs{k/k}=\ln1=0$,そして $w$ の分子が $L'$ の分母 $k^2-p^2$ とちょうど相殺することである.

$$ \begin{align} \int_0^{k_F}p\,L(p)\,\dd p &= \Big[w(p)L(p)\Big]_0^{k_F} - \int_0^{k_F} w(p)\,L'(p)\,\dd p \label{eq:7-Ip1}\\ &= \frac{k_F^2-k^2}{2}\,L(k_F) - \underbrace{\frac{0-k^2}{2}\,L(0)}_{=\,0} - \int_0^{k_F}\frac{p^2-k^2}{2}\cdot\frac{2k}{k^2-p^2}\,\dd p \label{eq:7-Ip2}\\ &= \frac{k_F^2-k^2}{2}\ln\left|\frac{k_F+k}{k_F-k}\right| - \int_0^{k_F}(-k)\,\dd p \label{eq:7-Ip3}\\ &= \frac{k_F^2-k^2}{2}\ln\left|\frac{k_F+k}{k_F-k}\right| + k\,k_F . \label{eq:7-Ip4} \end{align} $$

\eqref{eq:7-Ip2}:境界項を代入した.$p=0$ の項は $L(0)=0$ のためにゼロである.
\eqref{eq:7-Ip3}:被積分関数を整理した.$\dfrac{p^2-k^2}{2}\cdot\dfrac{2k}{k^2-p^2} = \dfrac{-(k^2-p^2)}{2}\cdot\dfrac{2k}{k^2-p^2} = -k$ となり,対数の微分から出た分母が完全に約分されて定数になる.これがこの部分積分の狙いである.また $L(k_F)=\ln\abs{(k+k_F)/(k-k_F)} = \ln\abs{(k_F+k)/(k_F-k)}$ と書き換えた(分子分母の符号を同時に反転しても絶対値は変わらない).
\eqref{eq:7-Ip4}:$-\int_0^{k_F}(-k)\dd p = +k\,k_F$.

これを \eqref{eq:7-Ih} に戻す:

$$ \begin{align} I(\kk) &= \frac{2\pi}{k}\left[k\,k_F + \frac{k_F^2-k^2}{2}\ln\left|\frac{k_F+k}{k_F-k}\right|\right] \label{eq:7-Ires1}\\ &= 2\pi k_F + \frac{\pi(k_F^2-k^2)}{k}\ln\left|\frac{k_F+k}{k_F-k}\right| \label{eq:7-Ires2}\\ &= 2\pi k_F\left[1 + \frac{k_F^2-k^2}{2kk_F}\ln\left|\frac{k_F+k}{k_F-k}\right|\right]. \label{eq:7-Ires3} \end{align} $$

\eqref{eq:7-Ires3} では $2\pi k_F$ をくくり出した.

結果をまとめておく:

$$ \begin{equation} I(\kk) = \int_{\abs{\bm{p}}\le k_F}\frac{\dd^3p}{\abs{\kk-\bm{p}}^2} = 2\pi k_F\left[1 + \frac{k_F^2-k^2}{2kk_F}\ln\left|\frac{k_F+k}{k_F-k}\right|\right] \equiv 2\pi k_F\,F\!\left(\frac{k}{k_F}\right). \label{eq:7-Ik-final} \end{equation} $$

ここで無次元の関数

$$ \begin{equation} F(x) = 1 + \frac{1-x^2}{2x}\ln\left|\frac{1+x}{1-x}\right| \label{eq:7-Fx} \end{equation} $$

を導入した($x=k/k_F$ と置いて $k_F$ をくくり出せばこの形になる).

例:$I(\kk)$ の検算

結果 \eqref{eq:7-Ik-final} を特別な場合で確かめる.

(i) $k=0$.直接計算すると $I(\bm{0}) = \int_{p\le k_F}\dd^3p/p^2 = 4\pi\int_0^{k_F}p^2\dd p/p^2 = 4\pi k_F$.一方 \eqref{eq:7-Fx} で $x\to0$ とすると,$\ln\frac{1+x}{1-x} = 2x + \frac{2x^3}{3}+\cdots$(対数の展開)より

$$ \frac{1-x^2}{2x}\cdot\left(2x+O(x^3)\right) = (1-x^2)\left(1+O(x^2)\right)\to 1, $$

よって $F(0)=1+1=2$,$I=2\pi k_F\cdot2 = 4\pi k_F$.一致する.

(ii) $k=k_F$.\eqref{eq:7-Fx} で $x\to1$ とすると,対数は $\ln\frac{2}{\abs{1-x}}\to\infty$ と発散するが,前の係数 $\frac{1-x^2}{2x}=\frac{(1-x)(1+x)}{2x}\to0$ が線形にゼロに近づく.積は $(1-x)\ln\frac{1}{1-x}\to0$ であるから $F(1)=1$,すなわち $I(k_F)=2\pi k_F$.フェルミ面上の電子が感じる「重み」は,球の中心にいる電子のちょうど半分である.

(iii) 次元.$I$ は $\dd^3p/p^2$ の次元,すなわち波数の1乗の次元を持つ.\eqref{eq:7-Ik-final} は確かに $k_F$ の1乗に比例している.

7.7.6 外側の積分

残るは $\int_{\abs{\kk}\le k_F}\dd^3k\,I(\kk)$ である.$I$ は $\kk$ の向きによらず大きさ $k$ のみの関数なので,立体角積分は $4\pi$ を与える:

$$ \begin{equation} \int_{\abs{\kk}\le k_F}\dd^3k\;I(\kk) = 4\pi\int_0^{k_F}\dd k\;k^2\,I(k). \label{eq:7-outer1} \end{equation} $$

\eqref{eq:7-Ires2} を代入して2つの部分に分ける:

$$ \begin{align} \int_0^{k_F}k^2 I(k)\,\dd k &= 2\pi k_F\int_0^{k_F}k^2\,\dd k + \pi\int_0^{k_F}k\,(k_F^2-k^2)\ln\left|\frac{k_F+k}{k_F-k}\right|\dd k \label{eq:7-outer2}\\ &= 2\pi k_F\cdot\frac{k_F^3}{3} + \pi\,K, \qquad K \equiv \int_0^{k_F}k(k_F^2-k^2)\ln\frac{k_F+k}{k_F-k}\,\dd k . \label{eq:7-outer3} \end{align} $$

第1項では $\int_0^{k_F}k^2\dd k = k_F^3/3$,また第2項では $k^2\cdot\frac{k_F^2-k^2}{k} = k(k_F^2-k^2)$ とした.$0\le k\le k_F$ の範囲では $k_F+k>0$,$k_F-k\ge0$ なので絶対値は外してよい.

$K$ を計算するために $k = k_F x$ と置換する($\dd k = k_F\,\dd x$,$x:0\to1$):

$$ \begin{align} K &= \int_0^1 (k_F x)\big(k_F^2 - k_F^2x^2\big)\ln\frac{k_F(1+x)}{k_F(1-x)}\;k_F\,\dd x \label{eq:7-K1}\\ &= k_F^4\int_0^1 x\,(1-x^2)\,\ln\frac{1+x}{1-x}\,\dd x \;=\; k_F^4\,\mathcal{A}. \label{eq:7-K2} \end{align} $$

対数の中の $k_F$ は分子分母で約分される.残った無次元の定積分 $\mathcal{A}$ を数学ノートで計算する.

数学ノート:定積分 $\displaystyle\mathcal{A}=\int_0^1 x(1-x^2)\ln\frac{1+x}{1-x}\dd x = \frac13$

$L(x) = \ln\dfrac{1+x}{1-x}$ とおく.その導関数は

$$ L'(x) = \frac{1}{1+x} + \frac{1}{1-x} = \frac{(1-x)+(1+x)}{(1+x)(1-x)} = \frac{2}{1-x^2}. $$

部分積分を行う.$x(1-x^2) = x - x^3$ の原始関数のうち,$x=1$ でゼロになるものを選ぶ:

$$ W(x) = \frac{x^2}{2}-\frac{x^4}{4}-\frac14 = -\frac{1}{4}\left(x^4-2x^2+1\right) = -\frac{1}{4}\left(1-x^2\right)^2 . $$

(展開して確かめると $-\frac14(1-2x^2+x^4) = -\frac14+\frac{x^2}{2}-\frac{x^4}{4}$ となり,微分すれば $x-x^3$ に戻る.)この選び方の狙いは2つある:$W(1)=0$ により発散する対数を抑えること,そして $W$ に含まれる $(1-x^2)^2$ が $L'$ の分母 $1-x^2$ と約分することである.

$$ \begin{align} \mathcal{A} &= \Big[W(x)L(x)\Big]_0^1 - \int_0^1 W(x)L'(x)\,\dd x \label{eq:7-A1}\\ &= 0 - \int_0^1\left(-\frac{(1-x^2)^2}{4}\right)\cdot\frac{2}{1-x^2}\,\dd x \label{eq:7-A2}\\ &= \frac{1}{2}\int_0^1 (1-x^2)\,\dd x \label{eq:7-A3}\\ &= \frac{1}{2}\left[x-\frac{x^3}{3}\right]_0^1 = \frac{1}{2}\left(1-\frac13\right) = \frac{1}{2}\cdot\frac{2}{3} = \frac{1}{3}. \label{eq:7-A4} \end{align} $$

境界項が消える理由.$x=0$ では $W(0)=-\frac14$ だが $L(0)=\ln1=0$ なので積はゼロ.$x=1$ では $L(x)\to\infty$ だが,$1-x^2 = (1-x)(1+x)$ より $W(x) \simeq -\frac{(1-x)^2(1+x)^2}{4}\simeq -(1-x)^2$ であり,$L(x)\simeq\ln\frac{2}{1-x}$ である.したがって積は $-(1-x)^2\ln\frac{2}{1-x}$ の形で,$\varepsilon=1-x\to0$ のとき $\varepsilon^2\ln(1/\varepsilon)\to0$(べきが対数に勝つ)なのでゼロになる.

参考:関連する不定積分.同じ手口で,より一般に

$$ \int x\ln\frac{1+x}{1-x}\,\dd x = \frac{x^2-1}{2}\ln\frac{1+x}{1-x} + x + C $$

が得られる($W(x)=(x^2-1)/2$ を使い,$\frac{x^2-1}{2}\cdot\frac{2}{1-x^2} = -1$ が定数になることを利用する).この不定積分は7.7.5でも(変数名を変えて)使ったものであり,電子ガスの計算で繰り返し現れる.

$\mathcal{A}=1/3$ を \eqref{eq:7-K2} に代入して $K = k_F^4/3$.これを \eqref{eq:7-outer3} に戻すと

$$ \begin{equation} \int_0^{k_F}k^2I(k)\,\dd k = \frac{2\pi k_F^4}{3} + \frac{\pi k_F^4}{3} = \pi k_F^4 . \label{eq:7-outer4} \end{equation} $$

2つの寄与の比が $2:1$ で,和がきれいに $\pi k_F^4$ にまとまるのは気持ちがよい.したがって \eqref{eq:7-outer1} より

$$ \begin{equation} \int_{\abs{\kk}\le k_F}\dd^3k\;I(\kk) = 4\pi\cdot\pi k_F^4 = 4\pi^2 k_F^4 . \label{eq:7-outer5} \end{equation} $$

7.7.7 交換エネルギーの完成

\eqref{eq:7-outer5} を \eqref{eq:7-E1int3} に代入する:

$$ \begin{align} E_1 &= -\frac{e^2V}{16\pi^5}\cdot 4\pi^2k_F^4 = -\frac{e^2Vk_F^4}{4\pi^3}. \label{eq:7-E1total} \end{align} $$

1電子あたりに直す.$N = Vk_F^3/(3\pi^2)$ を使うと

$$ \begin{align} \frac{E_1}{N} &= -\frac{e^2Vk_F^4}{4\pi^3}\cdot\frac{3\pi^2}{Vk_F^3} = -\frac{3\,e^2k_F}{4\pi}. \label{eq:7-E1perN} \end{align} $$

$V$ と $k_F^3$ が約分され,$\pi^2/\pi^3 = 1/\pi$ となった.$r_s$ で表すには \eqref{eq:7-kF-rs} を代入する.原子単位($e^2=1$,長さは $a_0$)では

$$ \begin{align} \frac{E_1}{N} &= -\frac{3}{4\pi}\left(\frac{9\pi}{4}\right)^{1/3}\frac{1}{r_s} \label{eq:7-E1rs1}\\ &= -\frac{3}{4\times3.141593}\times 1.919158\times\frac{1}{r_s} \label{eq:7-E1rs2}\\ &= -0.238732\times1.919158\times\frac{1}{r_s} = -\frac{0.458165}{r_s}\ \mathrm{Ha} = -\frac{0.916331}{r_s}\ \mathrm{Ry}. \label{eq:7-E1rs3} \end{align} $$
$$ \begin{equation} \frac{E_1}{N} = -\frac{3}{4}\frac{e^2 k_F}{\pi} = -\frac{3}{4\pi}\left(\frac{9\pi}{4}\right)^{1/3}\frac{e^2}{a_0 r_s} = -\frac{0.4582}{r_s}\ \mathrm{Ha} = -\frac{0.9163}{r_s}\ \mathrm{Ry} \label{eq:7-exchange-key} \end{equation} $$

何が得られたか.ジェリウムの1電子あたりの交換エネルギーが,密度パラメータ $r_s$ に反比例する簡潔な閉じた式として求まった.値は必ず負であり,電子が同じスピンの相手を避けることによるエネルギーの利得を表す.$r_s^{-1}$ というスケーリングは7.5節の予言と一致している.

例:運動エネルギーと交換エネルギーの競合

Ry 単位で $E/N \simeq 2.21/r_s^2 - 0.916/r_s$ である.両者が釣り合う($E/N=0$ となる)のは $2.21/r_s = 0.916$,すなわち $r_s = 2.41$ のときである.$r_s$ をこれより大きくすると全エネルギーは負になり,系は「束縛」される.$\dd(E/N)/\dd r_s = -4.42/r_s^3 + 0.916/r_s^2 = 0$ から極小は $r_s = 4.42/0.916 = 4.83$,そこでの値は $2.21/23.3 - 0.916/4.83 = 0.0948-0.1897 = -0.0949$ Ry である.相関エネルギーを加えると極小はもう少し高密度側に動き,$r_s\simeq4.2$,$E/N\simeq-0.155$ Ry となる(7.10節).

この極小の存在は重要である.運動エネルギーだけ(第3章)では電子ガスは束縛せず,いくらでも膨張してエネルギーを下げようとする.交換エネルギーの引力的な寄与が加わって初めて,有限の密度でエネルギーが極小になる.金属が特定の体積で安定に存在することの,最も単純な説明がここにある.

7.7.8 交換エネルギーの1電子的な読み方

ここまでは全系のエネルギーを計算した.同じ結果を「1電子が感じる交換ポテンシャル」という視点から見直しておくと,第5章のHF方程式や第10章のコーン・シャム方程式との対応が明確になる.

\eqref{eq:7-E1b} を,波数 $\kk$ スピン $\lambda$ の1電子が感じる交換自己エネルギー(フォック項)

$$ \begin{equation} \Sigma_x(\kk) \equiv -\frac{e^2}{V}\sum_{\abs{\bm{p}}\le k_F}\frac{4\pi}{\abs{\kk-\bm{p}}^2} \label{eq:7-sigmax-def} \end{equation} $$

を使って書き直すと,$E_1 = \frac12\sum_{\kk\lambda}n_{\kk\lambda}\Sigma_x(\kk)$ となる(和の中身が $\kk$ と $\bm{p}$ について対称なので,\eqref{eq:7-E1b} をこの形に組み替えられる).$\Sigma_x$ を \eqref{eq:7-sum2int} と \eqref{eq:7-Ik-final} で評価すると

$$ \begin{align} \Sigma_x(\kk) &= -\frac{4\pi e^2}{V}\cdot\frac{V}{(2\pi)^3}\,I(\kk) = -\frac{4\pi e^2}{8\pi^3}\,I(\kk) = -\frac{e^2}{2\pi^2}\,I(\kk) \label{eq:7-sigmax1}\\ &= -\frac{e^2}{2\pi^2}\cdot 2\pi k_F\,F\!\left(\frac{k}{k_F}\right) = -\frac{e^2k_F}{\pi}\,F\!\left(\frac{k}{k_F}\right) \label{eq:7-sigmax2} \end{align} $$

を得る.$F$ は \eqref{eq:7-Fx} で定義した関数である.この関数を図7.6に示す.

交換自己エネルギーの波数依存性を表す関数 F(x) のグラフ.横軸は x=k/k_F,縦軸は F(x).F(0)=2 から単調に減少し,x=1(フェルミ面,赤の縦破線)で F=1 となり,ここで傾きが −∞ になる.x=2 で約 0.18.
図7.6 交換自己エネルギーの波数依存性を表す関数 $F(x)$.$x=1$(フェルミ面)で値は $1$ だが,導関数が対数的に発散する.この特異性がハートリー・フォック近似の病理を生む.

物理的意味:HF近似の病理と第8章への伏線

$\Sigma_x(\kk)$ を1電子エネルギーに加えると,HF近似での準粒子エネルギーは

$$ \varepsilon^{\rm HF}(k) = \frac{\hbar^2k^2}{2m} - \frac{e^2k_F}{\pi}F\!\left(\frac{k}{k_F}\right) $$

となる.ここで重大な問題が生じる.$F(x)$ の導関数は $x\to1$ で $\ln\abs{1-x}$ のように発散するので,$\varepsilon^{\rm HF}(k)$ の傾き $\dd\varepsilon/\dd k$ がフェルミ面で無限大になる.第3章で見たように状態密度は $D(\varepsilon)\propto k^2/(\dd\varepsilon/\dd k)$ に比例するから,HF近似は金属のフェルミ準位における状態密度をゼロにしてしまう.これは電子比熱が $T\ln T$ に比例するという実験と合わない結論を導く.さらにHF近似はバンド幅を大きく見積もりすぎる(Na で実測の約2.5倍).

病気の原因は,$q\to0$ で $4\pi e^2/q^2$ が発散する裸のクーロン相互作用を遮蔽なしで使ったことにある.実際の電子ガスでは,1個の電子のまわりに他の電子が押しのけられて正の「穴」ができ,遠方から見ると電荷が打ち消されて相互作用が短距離化する.この遮蔽(screening)を取り入れると $4\pi e^2/(q^2+\kappa^2)$ のような形になり,$q\to0$ の特異性が消えて状態密度は有限に戻る.遮蔽の理論が第8章の主題であり,それを摂動論の言葉で系統的に行うのが7.9節のRPAである.

なお,係数 $\frac12$ にも注意しておこう.$E_1 = \frac12\sum_{\kk\lambda}n_{\kk\lambda}\Sigma_x(\kk)$ の $\frac12$ は,1電子ごとの交換自己エネルギーを単純に足すと電子対を2回数えてしまうための補正である.「全エネルギー $\neq$ 1電子エネルギーの和」というこの事実は,第10章でコーン・シャム固有値の和が全エネルギーでないことの原型である.

7.8 ディラック交換と局所密度近似の原型

前節で得た交換エネルギーは $r_s$ の関数として書かれていた.本節ではこれを密度 $\rho$ の関数に書き直す.すると,7.6節の運動エネルギーと合わせて,1電子あたりのエネルギーが $\rho$ だけで表される2つの項として並ぶ.この事実が密度汎関数理論の出発点であり,本書全体を通じて最も重要なメッセージである.さらに,交換エネルギーの実空間での姿(交換ホール)とスピン分極系への拡張を扱う.

7.8.1 交換エネルギー密度

\eqref{eq:7-E1perN} に $k_F=(3\pi^2\rho)^{1/3}$ を代入する.

$$ \begin{align} \varepsilon_x(\rho) \equiv \frac{E_1}{N} &= -\frac{3e^2}{4\pi}\,k_F \label{eq:7-ex1}\\ &= -\frac{3e^2}{4\pi}\left(3\pi^2\rho\right)^{1/3} \label{eq:7-ex2}\\ &= -\frac{3e^2}{4}\left(\frac{3\pi^2}{\pi^3}\right)^{1/3}\rho^{1/3} \label{eq:7-ex3}\\ &= -\frac{3e^2}{4}\left(\frac{3}{\pi}\right)^{1/3}\rho^{1/3}. \label{eq:7-ex4} \end{align} $$

\eqref{eq:7-ex3} では,$\pi$ を3乗根の中に入れた:$\dfrac{1}{\pi}=\left(\dfrac{1}{\pi^3}\right)^{1/3}$.\eqref{eq:7-ex4} では $3\pi^2/\pi^3 = 3/\pi$ と約分した.数値は

$$ \left(\frac{3}{\pi}\right)^{1/3} = (0.954930)^{1/3} = 0.984745, \qquad C_x \equiv \frac{3}{4}\left(\frac{3}{\pi}\right)^{1/3} = 0.75\times0.984745 = 0.738559 . $$

したがって原子単位($e^2=1$)では

$$ \begin{equation} \varepsilon_x(\rho) = -C_x\,\rho^{1/3} = -0.7386\,\rho^{1/3}\ \ [\mathrm{Ha}]. \label{eq:7-eps-x} \end{equation} $$

7.8.2 局所密度近似としての交換汎関数

ここまでは一様な電子ガスの話であった.実在の系(原子・分子・固体)では密度 $\rho(\rr)$ は場所によって変わる.そこで次の大胆な仮定を置く.

定義:交換エネルギーの局所密度近似(LDA)

非一様な密度分布 $\rho(\rr)$ を持つ系の交換エネルギーを,空間を微小体積 $\dd^3r$ に分割し,各点でその点の密度 $\rho(\rr)$ を持つ一様電子ガスの交換エネルギー密度を借りてきて足し合わせることで近似する:

$$ \begin{equation} E_x^{\rm LDA}[\rho] = \int \dd^3r\;\rho(\rr)\,\varepsilon_x\big(\rho(\rr)\big). \label{eq:7-Ex-lda-def} \end{equation} $$

被積分関数が「1電子あたりのエネルギー $\varepsilon_x$」×「その点の電子数密度 $\rho$」の形になっていることに注意せよ.

\eqref{eq:7-eps-x} を代入すると,$\rho\cdot\rho^{1/3}=\rho^{4/3}$ より

$$ \begin{equation} E_x^{\rm LDA}[\rho] = -\,C_x\int \dd^3r\;\rho(\rr)^{4/3}, \qquad C_x = \frac{3}{4}\left(\frac{3}{\pi}\right)^{1/3} = 0.7386 \label{eq:7-dirac-exchange} \end{equation} $$

を得る.これがディラック(Dirac)の交換汎関数である(Dirac 1930).密度の $4/3$ 乗を積分するだけという,驚くほど単純な式である.

この汎関数の汎関数微分(第4章)を計算しておこう.$\delta E_x/\delta\rho(\rr)$ は $\rho^{4/3}$ を $\rho$ で微分すればよいので

$$ \begin{equation} v_x(\rr) = \frac{\delta E_x^{\rm LDA}}{\delta\rho(\rr)} = -C_x\cdot\frac{4}{3}\rho(\rr)^{1/3} = -\frac{4}{3}\cdot\frac{3}{4}\left(\frac{3}{\pi}\right)^{1/3}\rho^{1/3} = -\left(\frac{3\rho(\rr)}{\pi}\right)^{1/3}. \label{eq:7-vx} \end{equation} $$

これが交換ポテンシャルであり,第10章のコーン・シャム方程式に現れる.係数がちょうど $1$ になるという美しい形をしている.なお,スレーターは1951年に,HF交換演算子を占有状態について平均する別の議論から $v_x^{\rm Slater} = -\frac{3}{2}(3\rho/\pi)^{1/3}$($\eqref{eq:7-vx}$ の $3/2$ 倍)を導いた.両者を補間する $v_{x\alpha} = -\frac{3\alpha}{2}(3\rho/\pi)^{1/3}$ が$X\alpha$ 法であり,$\alpha=2/3$ が \eqref{eq:7-vx} に対応する.$\alpha$ を $0.7$ 前後に調整すると原子の全エネルギーがよく再現されることが経験的に知られ,DFTが確立する以前の実用的手法として広く使われた.

7.8.3 トーマス・フェルミ・ディラック模型

第4章では,運動エネルギーを一様電子ガスから借りてくるトーマス・フェルミ(TF)模型を扱った.そこに \eqref{eq:7-dirac-exchange} を加えたものがトーマス・フェルミ・ディラック(TFD)模型である.原子核(電荷 $Z$)の場合,エネルギー汎関数は

$$ \begin{equation} E_{\rm TFD}[\rho] = \underbrace{C_F\int\rho^{5/3}\dd^3r}_{\text{運動}} \;\underbrace{-\,Z\int\frac{\rho(\rr)}{r}\dd^3r}_{\text{核との引力}} \;+\;\underbrace{\frac12\iint\frac{\rho(\rr)\rho(\rr')}{\abs{\rr-\rr'}}\dd^3r\,\dd^3r'}_{\text{ハートリー}} \;\underbrace{-\,C_x\int\rho^{4/3}\dd^3r}_{\text{交換}} \label{eq:7-TFD} \end{equation} $$

となる.ここで

$$ C_F = \frac{3}{10}\left(3\pi^2\right)^{2/3} = 0.3\times 9.57080 = 2.87123, \qquad C_x = 0.73856 $$

である($C_F$ の数値は $(3\pi^2)^{1/3} = 29.6088^{1/3} = 3.09367$ を2乗して $9.57080$ から得た).第4章と同じく粒子数保存 $\int\rho\,\dd^3r = N$ の拘束のもとでラグランジュ未定乗数法により変分すれば,TFD方程式

$$ \begin{equation} \frac{5}{3}C_F\rho^{2/3}(\rr) - \left(\frac{3\rho(\rr)}{\pi}\right)^{1/3} + v_{\rm eff}(\rr) = \mu \label{eq:7-TFD-eq} \end{equation} $$

が得られる($v_{\rm eff}$ は核とハートリーのポテンシャルの和,$\mu$ は化学ポテンシャル).第4章では $\frac53 C_F\rho^{2/3} = \frac12(3\pi^2)^{2/3}\rho^{2/3}$ であったから,TFD方程式はTF方程式に交換の項 $-(3\rho/\pi)^{1/3}$ を加えたものである.

例:TFD模型の限界

交換項を加えるとTF模型は改善されるが,根本的な欠陥は残る.TF模型では原子は互いに結合しない(テラーの定理,第4章)が,交換項を加えてもこの結論は変わらない.また,TFDでは原子の電子密度が有限の半径で厳密にゼロになる(TFでは無限に裾を引く)という定性的な変化はあるものの,原子の全エネルギーの誤差は依然として数%残る.

問題は運動エネルギー汎関数 $C_F\int\rho^{5/3}$ の精度にある.原子の内殻のように密度が急激に変化する領域では,一様電子ガスから借りてきた $\rho^{5/3}$ は運動エネルギーを大きく過小評価する.この難点を根本的に解決したのが,運動エネルギーだけは軌道を使って厳密に計算するというコーン・シャムの発想(第10章)である.交換相関だけをLDAで扱い,運動エネルギーは軌道で扱う——この分業がDFTを実用の理論に変えた.

7.8.4 DFTの原型:密度だけで書けた2つの項

ここで立ち止まって,本章がここまでに達成したことを確認しよう.相互作用する多電子系(ジェリウム)の1電子あたりのエネルギーは,摂動の1次までで

$$ \begin{equation} \frac{E}{N} = \underbrace{\frac{3}{10}(3\pi^2)^{2/3}\rho^{2/3}}_{\text{運動エネルギー}} \;\underbrace{-\,\frac{3}{4}\left(\frac{3}{\pi}\right)^{1/3}\rho^{1/3}}_{\text{交換エネルギー}} \;+\;\cdots \qquad [\mathrm{Ha}] \label{eq:7-dft-prototype} \end{equation} $$

と書けた.どちらの項も,波動関数を一切含まず,電子密度 $\rho$ だけの関数である.これが「エネルギーは密度の汎関数である」という主張の,歴史上最初の具体的な証拠であった.

物理的意味:なぜこれが革命的なのか

$N$ 電子系の波動関数 $\Psi(\rr_1,\ldots,\rr_N)$ は $3N$ 個の変数の関数である.電子が100個あれば300次元の関数であり,格子点を各次元10点にとるだけで $10^{300}$ 個の数値が必要になる(第1章).一方,電子密度 $\rho(\rr)$ はどんなに電子が多くても3変数の関数である.もしエネルギーが $\rho$ だけで決まるなら,多体問題の次元は $3N$ から $3$ へと劇的に縮小する.

\eqref{eq:7-dft-prototype} は,少なくとも一様電子ガスについてはそれが実際に可能であることを示した.しかしこの時点では「一様系だから当たり前ではないか」という反論が可能である.一様系では密度が唯一のパラメータなのだから,エネルギーが密度の関数になるのは同語反復に近い.任意の外部ポテンシャルの下にある任意の系について,エネルギーが密度の汎関数であることを証明したのがホーエンベルグ・コーンの定理(第9章)であり,その汎関数を \eqref{eq:7-dft-prototype} の形から構成する具体的な処方が局所密度近似(第11章)である.本章の計算は,その両方の出発点になっている.

もう一つ強調しておきたいのは,\eqref{eq:7-dft-prototype} の2つの項がべきの異なる関数だということである.運動項は $\rho^{2/3}$,交換項は $\rho^{1/3}$.密度を上げると運動項の方が速く増えるので,高密度では運動エネルギーが,低密度では交換(と相関)が支配する.7.5節の $r_s$ スケーリングを密度の言葉で言い換えたものが,まさにこのべきの違いである.

7.8.5 スピン分極系への拡張

これまではスピン $\uparrow$ と $\downarrow$ の電子数が等しい(常磁性の)場合を扱ってきた.磁性を論じるには,両者が異なる場合に拡張する必要がある.幸い,交換エネルギーは同じスピンの対にしか働かない(7.7.2節)ので,拡張は容易である.

\eqref{eq:7-E1b} をスピンごとに分けて書くと,スピン $\lambda$ の寄与は,その成分だけからなる完全分極電子ガスの交換エネルギーである.フェルミ波数を $k_{F\lambda}$ とすると,7.7節の計算をそのまま繰り返して

$$ \begin{align} E_1^{(\lambda)} &= -\frac{e^2}{2V}\cdot4\pi\cdot\left[\frac{V}{(2\pi)^3}\right]^2 \int_{\abs{\kk}\le k_{F\lambda}}\dd^3k\int_{\abs{\bm{p}}\le k_{F\lambda}}\dd^3p\;\frac{1}{\abs{\kk-\bm{p}}^2} \label{eq:7-spin1}\\ &= -\frac{e^2V}{32\pi^5}\cdot4\pi^2k_{F\lambda}^4 = -\frac{e^2Vk_{F\lambda}^4}{8\pi^3} \label{eq:7-spin2} \end{align} $$

を得る(\eqref{eq:7-E1c} でスピン和による因子2を掛けなかった分だけ,\eqref{eq:7-E1total} の半分になっている).1つのスピン成分だけを数えるときの密度とフェルミ波数の関係は

$$ \rho_\lambda = \frac{N_\lambda}{V} = \frac{1}{V}\cdot\frac{V}{(2\pi)^3}\cdot\frac{4\pi}{3}k_{F\lambda}^3 = \frac{k_{F\lambda}^3}{6\pi^2} \qquad\Longleftrightarrow\qquad k_{F\lambda} = \left(6\pi^2\rho_\lambda\right)^{1/3} $$

である(スピン縮重度2が付かないので,\eqref{eq:7-kF-rs} の $3\pi^2$ が $6\pi^2$ になる).これを代入すると

$$ \begin{align} \frac{E_1^{(\lambda)}}{V} &= -\frac{e^2}{8\pi^3}\left(6\pi^2\rho_\lambda\right)^{4/3} = -\frac{e^2\,6^{4/3}\pi^{8/3}}{8\pi^3}\rho_\lambda^{4/3} = -\frac{e^2\,6^{4/3}}{8\,\pi^{1/3}}\rho_\lambda^{4/3} \label{eq:7-spin3}\\ &= -\,2^{1/3}C_x\,e^2\,\rho_\lambda^{4/3}, \qquad 2^{1/3}C_x = 0.930526 . \label{eq:7-spin4} \end{align} $$

最後の等式は $\dfrac{6^{4/3}}{8\pi^{1/3}} = 2^{1/3}\cdot\dfrac{3}{4}\left(\dfrac{3}{\pi}\right)^{1/3}$ を確かめれば分かる(両辺を3乗すると $\dfrac{6^4}{512\pi} = 2\cdot\dfrac{27}{64}\cdot\dfrac{3}{\pi}$,すなわち $\dfrac{1296}{512\pi} = \dfrac{162}{64\pi}$ となり,$1296/512 = 2.53125$,$162/64 = 2.53125$ で一致する).したがってスピン分極系の交換汎関数は

$$ \begin{equation} E_x[\rho_\uparrow,\rho_\downarrow] = -\,2^{1/3}C_x\int\dd^3r\left(\rho_\uparrow^{4/3}+\rho_\downarrow^{4/3}\right). \label{eq:7-Ex-spin} \end{equation} $$

これを全密度 $\rho = \rho_\uparrow+\rho_\downarrow$ とスピン分極率 $\zeta = (\rho_\uparrow-\rho_\downarrow)/\rho$ で書き直す.$\rho_\uparrow = \rho(1+\zeta)/2$,$\rho_\downarrow=\rho(1-\zeta)/2$ を代入すると

$$ \begin{align} \frac{E_x}{V} &= -2^{1/3}C_x\,\rho^{4/3} \left[\left(\frac{1+\zeta}{2}\right)^{4/3}+\left(\frac{1-\zeta}{2}\right)^{4/3}\right] \label{eq:7-zeta1}\\ &= -2^{1/3}C_x\,\rho^{4/3}\cdot 2^{-4/3} \left[(1+\zeta)^{4/3}+(1-\zeta)^{4/3}\right] \label{eq:7-zeta2}\\ &= -C_x\,\rho^{4/3}\cdot\frac{(1+\zeta)^{4/3}+(1-\zeta)^{4/3}}{2}, \label{eq:7-zeta3} \end{align} $$

\eqref{eq:7-zeta2} で $2^{-4/3}$ を括り出し,\eqref{eq:7-zeta3} で $2^{1/3}\cdot2^{-4/3} = 2^{-1} = 1/2$ とした.すなわち

$$ \begin{equation} \varepsilon_x(\rho,\zeta) = \varepsilon_x(\rho,0)\,f(\zeta), \qquad f(\zeta) = \frac{(1+\zeta)^{4/3}+(1-\zeta)^{4/3}}{2} \label{eq:7-fzeta} \end{equation} $$

である.$f(0)=1$(常磁性),$f(1) = \dfrac{2^{4/3}+0}{2} = 2^{1/3} = 1.2599$(完全分極).すなわち完全にスピン分極すると交換エネルギーは $2^{1/3}$ 倍(約26%)大きくなる.

例:ブロッホの強磁性条件

スピン分極は交換エネルギーを得する一方,運動エネルギーを損する.完全分極状態では片方のスピンだけに $N$ 個の電子を詰めるので,フェルミ波数は $2^{1/3}$ 倍になり,運動エネルギーは $(2^{1/3})^2 = 2^{2/3} = 1.5874$ 倍になる.Ry 単位で比べると

$$ \frac{E}{N}\bigg|_{\zeta=0} = \frac{2.21}{r_s^2}-\frac{0.916}{r_s}, \qquad \frac{E}{N}\bigg|_{\zeta=1} = \frac{2^{2/3}\times2.21}{r_s^2}-\frac{2^{1/3}\times0.916}{r_s} = \frac{3.508}{r_s^2}-\frac{1.154}{r_s}. $$

差をとると $\Delta = \dfrac{1.298}{r_s^2}-\dfrac{0.238}{r_s}$ であり,$\Delta<0$(分極が得)となるのは $r_s > 1.298/0.238 = 5.45$ のときである.これがブロッホ(1929)による強磁性の条件である.

しかしこの予言は実在の金属には当てはまらない.Cs($r_s=5.62$)は $r_s>5.45$ を満たすが強磁性ではないし,逆に強磁性を示す Fe,Co,Ni は $d$ 電子の局在性が本質で,自由電子的なジェリウムでは記述できない.ブロッホの議論の欠陥は相関エネルギーを無視したことにある.相関は常磁性状態でより効きやすい(異スピン電子どうしも避け合えるため)ので,分極による利得を打ち消す.量子モンテカルロ計算(7.9節)によれば,3次元ジェリウムが部分分極を起こすのは $r_s\approx50$ 以上という極端な低密度領域である.とはいえ,\eqref{eq:7-fzeta} の形は局所スピン密度近似(LSDA)の骨格としてそのまま使われており,実際の磁性体の第一原理計算の基礎になっている(第11章).

7.8.6 交換ホール:実空間から見た交換エネルギー

交換エネルギーが負である理由を,実空間の描像で確認しておこう.第5章で導入した交換ホール(exchange hole)が,ジェリウムでは解析的に求まる.

出発点は,\eqref{eq:7-E1c} をフーリエ変換で実空間に戻すことである.$\mu\to0$ とした \eqref{eq:7-yukawa-ft} を逆向きに使う:$\dfrac{4\pi}{\abs{\kk-\bm{p}}^2} = \displaystyle\int\dd^3s\;\frac{e^{-i(\kk-\bm{p})\cdot\bm{s}}}{s}$.

導出:交換ホールと交換エネルギーの関係

$$ \begin{align} E_1 &= -\frac{e^2}{V}\sum_{\abs{\kk}\le k_F}\sum_{\abs{\bm{p}}\le k_F} \int\dd^3s\;\frac{e^{-i(\kk-\bm{p})\cdot\bm{s}}}{s} \label{eq:7-hole1}\\ &= -\frac{e^2}{V}\int\dd^3s\;\frac{1}{s} \left(\sum_{\abs{\kk}\le k_F}e^{-i\kk\cdot\bm{s}}\right) \left(\sum_{\abs{\bm{p}}\le k_F}e^{i\bm{p}\cdot\bm{s}}\right) \label{eq:7-hole2}\\ &= -e^2V\int\dd^3s\;\frac{\abs{\rho_1(s)}^2}{s}, \qquad \rho_1(\bm{s}) \equiv \frac{1}{V}\sum_{\abs{\kk}\le k_F}e^{i\kk\cdot\bm{s}} . \label{eq:7-hole3} \end{align} $$

\eqref{eq:7-hole2} では,2つの和が独立なので積の形にまとめた.\eqref{eq:7-hole3} では,2つの和が互いに複素共役であることを使い,$(V\rho_1)^*(V\rho_1) = V^2\abs{\rho_1}^2$ とした.$\rho_1$ は1体密度行列(1スピン成分)であり,$\rho_1(0) = N/(2V) = \rho/2$ が1スピンあたりの密度になる.

$\rho_1$ を計算する.和を積分に直し,$\bm{s}$ を極軸にとった球座標を使う:

$$ \begin{align} \rho_1(s) &= \frac{1}{(2\pi)^3}\int_{\abs{\kk}\le k_F}\dd^3k\;e^{i\kk\cdot\bm{s}} = \frac{2\pi}{(2\pi)^3}\int_0^{k_F}\dd k\,k^2\int_{-1}^{1}\dd u\;e^{iksu} \label{eq:7-rho1a}\\ &= \frac{1}{4\pi^2}\int_0^{k_F}\dd k\,k^2\cdot\frac{2\sin(ks)}{ks} = \frac{1}{2\pi^2 s}\int_0^{k_F}\dd k\;k\sin(ks) \label{eq:7-rho1b}\\ &= \frac{1}{2\pi^2 s}\cdot\frac{\sin(k_Fs)-k_Fs\cos(k_Fs)}{s^2} = \frac{k_F^3}{2\pi^2}\cdot\frac{j_1(k_Fs)}{k_Fs}. \label{eq:7-rho1c} \end{align} $$

\eqref{eq:7-rho1a} では方位角の積分($2\pi$)を実行し,$u=\cos\theta$ と置換した.\eqref{eq:7-rho1b} では $\int_{-1}^1 e^{iksu}\dd u = 2\sin(ks)/(ks)$(7.3.2節と同じ計算)を代入した.\eqref{eq:7-rho1c} の $k$ 積分は部分積分による:$u=k$,$\dd v=\sin(ks)\dd k$,$v=-\cos(ks)/s$ として

$$ \int_0^{k_F}k\sin(ks)\dd k = \left[-\frac{k\cos(ks)}{s}\right]_0^{k_F} + \frac{1}{s}\int_0^{k_F}\cos(ks)\,\dd k = -\frac{k_F\cos(k_Fs)}{s} + \frac{\sin(k_Fs)}{s^2}. $$

最後に球ベッセル関数 $j_1(x) = \dfrac{\sin x - x\cos x}{x^2}$ の定義を使って書き直した.$s\to0$ では $j_1(x)/x\to1/3$(なぜなら $\sin x - x\cos x = x^3/3 + O(x^5)$)であるから $\rho_1(0)=k_F^3/(6\pi^2) = \rho/2$ となり,確かに1スピンあたりの密度に一致する.

交換ホールは,1体密度行列を用いて

$$ \begin{equation} n_x(s) \equiv -\frac{\abs{\rho_1(s)}^2}{\rho/2} = -\frac{9}{2}\,\rho\left[\frac{j_1(k_Fs)}{k_Fs}\right]^2 \label{eq:7-exchange-hole} \end{equation} $$

と定義される(第5章の定義をジェリウムに適用したもの).係数は $\dfrac{(k_F^3/2\pi^2)^2}{\rho/2} = \dfrac{k_F^6}{4\pi^4}\cdot\dfrac{2}{\rho}$ に $\rho=k_F^3/3\pi^2$ を代入して $\dfrac{k_F^6}{4\pi^4}\cdot\dfrac{6\pi^2}{k_F^3} = \dfrac{3k_F^3}{2\pi^2} = \dfrac92\rho$ から得られる.この量には次の2つの重要な性質がある.

総和則は次のように示される.パーセバルの定理(あるいは直接 \eqref{eq:7-hole3} の定義から)

$$ \int\dd^3s\;\abs{\rho_1(s)}^2 = \frac{1}{(2\pi)^3}\int_{\abs{\kk}\le k_F}\dd^3k = \frac{1}{(2\pi)^3}\cdot\frac{4\pi}{3}k_F^3 = \frac{k_F^3}{6\pi^2} = \frac{\rho}{2}, $$

であるから $\int n_x\,\dd^3s = -(\rho/2)/(\rho/2) = -1$ となる.

さらに \eqref{eq:7-hole3} と \eqref{eq:7-exchange-hole} を組み合わせると,$\abs{\rho_1}^2 = -(\rho/2)n_x$ より

$$ \begin{equation} \frac{E_1}{N} = -\frac{e^2V}{N}\int\dd^3s\,\frac{\abs{\rho_1}^2}{s} = \frac{e^2V}{N}\cdot\frac{\rho}{2}\int\dd^3s\,\frac{n_x(s)}{s} = \frac{e^2}{2}\int\dd^3s\;\frac{n_x(s)}{s} \label{eq:7-Ex-hole} \end{equation} $$

という美しい関係が得られる($V\rho/N = 1$ を使った).交換エネルギーとは,電子が自分のまわりに掘った「電子1個分の穴」との静電引力である.穴の実効的な半径を $R$ とすれば $\varepsilon_x\sim -e^2/(2R)$ 程度になるはずで,実際 \eqref{eq:7-exchange-hole} から $R\sim 1/k_F\sim r_0$ であるから $\varepsilon_x\sim -e^2/(2r_0)\propto -1/r_s$ となり,\eqref{eq:7-exchange-key} のスケーリングが直観的に納得できる.

ハートリー・フォック近似におけるジェリウムの対分布関数.横軸は x=k_F r12.異スピンは全域で g=1(緑),同スピンは x=0 で 0 から立ち上がり x≈3 で 1 に近づく(赤),スピン平均は x=0 で 1/2 から 1 に近づく(青の破線).
図7.7 ハートリー・フォック近似におけるジェリウムの対分布関数.同スピン電子は $r_{12}\to0$ で完全に排除される(交換ホール).異スピン電子はこの近似では完全に無相関であり,$g_{\uparrow\downarrow}=1$ のままである.この欠落を埋めるのが相関エネルギー(7.9節)である.

例:交換ホールから $\varepsilon_x$ を再現する

\eqref{eq:7-Ex-hole} に \eqref{eq:7-exchange-hole} を代入すると,$x=k_Fs$ と置換して

$$ \frac{E_1}{N} = \frac{e^2}{2}\cdot4\pi\int_0^{\infty}\dd s\,s\,\left(-\frac92\rho\right)\left[\frac{j_1(k_Fs)}{k_Fs}\right]^2 = -\frac{9\pi e^2\rho}{k_F^2}\int_0^{\infty}\frac{j_1(x)^2}{x}\,\dd x . $$

($\dd^3s = 4\pi s^2\dd s$ を $1/s$ に掛けて $4\pi s\,\dd s$ とし,$s\,\dd s = x\,\dd x/k_F^2$ を使った.)一方 \eqref{eq:7-exchange-key} は $E_1/N = -3e^2k_F/(4\pi)$ を与えるから,両者が等しいことを要求すると

$$ \int_0^{\infty}\frac{j_1(x)^2}{x}\,\dd x = \frac{3k_F^3}{4\pi\cdot9\pi\rho} = \frac{3k_F^3}{36\pi^2}\cdot\frac{3\pi^2}{k_F^3} = \frac{1}{4} $$

となる($\rho=k_F^3/3\pi^2$ を使った).これは球ベッセル関数の公式 $\int_0^\infty j_n(x)^2\dd x/x = \dfrac{1}{2n(n+1)}$ の $n=1$ の場合($=1/4$)と一致する.7.7節の $\kk$ 空間の計算と本節の実空間の計算が同じ答えを与えることの確認になっている.

7.9 相関エネルギー

摂動の1次までは手で計算できた.2次以降はどうか.本節では,素朴に2次摂動を計算すると対数発散が生じることをまず示し,その病理が遮蔽によってどう治癒されるか(RPA)を説明する.続いて,高密度極限のGell-Mann–Brueckner展開,低密度極限のウィグナー結晶,そして実在金属が属する中間密度領域を決着させた量子モンテカルロ法(Ceperley–Alder)とそのパラメータ化(VWN)を紹介する.これらはすべて第11章のLDA汎関数の材料になる.

7.9.1 相関エネルギーの定義

定義:相関エネルギー

与えられた系の厳密な基底状態エネルギー $E_{\rm exact}$ と,ハートリー・フォック近似のエネルギー $E_{\rm HF}$ の差を相関エネルギーという:

$$ \begin{equation} E_c \equiv E_{\rm exact} - E_{\rm HF}. \label{eq:7-Ec-def} \end{equation} $$

$E_{\rm HF}$ は変分原理により厳密解より高いので,$E_c$ は常に負である.ジェリウムでは $E_{\rm HF}/N = E_0/N + E_1/N$ であるから,$\varepsilon_c = E/N - (2.21/r_s^2 - 0.916/r_s)$ Ry が相関エネルギー密度になる.

物理的には,相関エネルギーは「HF近似が取りこぼした電子どうしの避け合い」を表す.図7.7で見たように,HFでは異なるスピンの電子は完全に無相関($g_{\uparrow\downarrow}=1$)であった.実際にはクーロン反発のために異スピン電子も互いを避けるはずで,そのぶんエネルギーが下がる.この効果が相関エネルギーの主要部分である.

7.9.2 2次摂動の対数発散

摂動の2次は一般に

$$ \begin{equation} E_2 = \sum_{n\neq0}\frac{\abs{\langle n|\hat{H}_1|\Phi_0\rangle}^2}{E_0^{(0)}-E_n^{(0)}} \label{eq:7-E2-def} \end{equation} $$

で与えられる.$\hat{H}_1$ は2個の電子をフェルミ球の内から外へ移す演算子であるから,中間状態 $|n\rangle$ は2個の粒子・正孔対を持つ状態である:電子 $(\kk,\lambda_1)$ が $(\kk+\bm{q},\lambda_1)$ へ,電子 $(\bm{p},\lambda_2)$ が $(\bm{p}-\bm{q},\lambda_2)$ へ励起される.行列要素には直接項と交換項が現れるが,小さい $\bm{q}$ で問題を起こすのは直接項の2乗(いわゆるリング項)である:

$$ \begin{equation} E_2^{\rm ring} = -\frac{1}{2}\sum_{\bm{q}\neq\bm{0}}\left(\frac{4\pi e^2}{Vq^2}\right)^2 \sum_{\kk\bm{p}}\sum_{\lambda_1\lambda_2} \frac{n_\kk(1-n_{\kk+\bm{q}})\,n_{\bm{p}}(1-n_{\bm{p}-\bm{q}})} {\varepsilon_{\kk+\bm{q}}+\varepsilon_{\bm{p}-\bm{q}}-\varepsilon_\kk-\varepsilon_{\bm{p}}} . \label{eq:7-E2-ring} \end{equation} $$

分母は励起エネルギー(正)であり,全体の符号は負である.占有数因子 $n_\kk(1-n_{\kk+\bm{q}})$ は「$\kk$ は占有されており $\kk+\bm{q}$ は空である」という条件を表す.この式が $\bm{q}\to0$ でどう振る舞うかを,次数勘定(power counting)で調べよう.

導出:$\bm{q}\to0$ での次数勘定

$\bm{q}$ が小さいとき,\eqref{eq:7-E2-ring} の各因子は次のように振る舞う.

  1. 相互作用の2乗:$\left(4\pi e^2/q^2\right)^2 \propto q^{-4}$.
  2. $\kk$ 和の位相空間:$n_\kk(1-n_{\kk+\bm{q}})=1$ となるのは,$\kk$ がフェルミ球の内側かつ $\kk+\bm{q}$ が外側,すなわちフェルミ面近傍の厚さ $\sim q$ の三日月形の殻に $\kk$ があるときだけである.その体積は $\propto k_F^2\,q \propto q$.$\bm{p}$ についても同様に $\propto q$.合わせて $q^2$.
  3. エネルギー分母:$\varepsilon_{\kk+\bm{q}}-\varepsilon_\kk = \dfrac{\hbar^2}{2m}\left(2\kk\cdot\bm{q}+q^2\right)\simeq \dfrac{\hbar^2}{m}\kk\cdot\bm{q}\propto q$.2つの励起の和も $\propto q$ であるから,分母は $\propto q$,その逆数は $\propto q^{-1}$.
  4. $\bm{q}$ 和:$\sum_{\bm{q}}\to \dfrac{V}{(2\pi)^3}\int\dd^3q = \dfrac{V}{(2\pi)^3}\cdot4\pi\int q^2\dd q \propto \int q^2\,\dd q$.

これらを掛け合わせると,小さい $q$ での被積分関数は

$$ q^2\,\dd q\;\times\;q^{-4}\;\times\;q^{2}\;\times\;q^{-1} = \frac{\dd q}{q} $$

となる.したがって

$$ E_2^{\rm ring}\;\sim\; -\,\mathrm{const}\times\int_0 \frac{\dd q}{q}\;=\;-\infty $$

で,対数発散する.

物理的意味:発散は「遮蔽を忘れた」ことの罰である

発散の原因は,$\bm{q}\to0$ すなわち非常に長い波長のゆらぎに対して,裸のクーロン相互作用 $4\pi e^2/q^2$ を使い続けたことにある.しかし実際の電子ガスでは,長波長の電荷ゆらぎは他の電子によって速やかに埋め立てられる.1個の電子のまわりには他の電子が押しのけられて正に帯電した領域ができ,遠くから見ると電荷が打ち消されて見えなくなる.これが遮蔽である.遮蔽を取り入れると相互作用は

$$ \frac{4\pi e^2}{q^2}\;\longrightarrow\;\frac{4\pi e^2}{q^2+\kappa^2} $$

のようになり(第8章),$q\to0$ で発散しない.つまり $q\lesssim\kappa$ の長波長領域では摂動論の展開そのものが破綻しているのであって,系が本当に不安定なわけではない.

これは物理学でしばしば起こる状況である:各次数の項は発散するが,それらを無限次まで足し上げると有限になる.摂動級数の項別の発散は,正しい「集団的」自由度(ここではプラズマ振動と遮蔽)を最初から取り入れなかったことの現れである.

7.9.3 RPA:リングダイアグラムの部分和

より高次の項を見ると事態はさらに悪い.$n$ 次のリング項(相互作用線が $n$ 本,粒子・正孔ループが $n$ 個)の次数勘定を同様に行うと

$$ q^2\dd q\times \underbrace{q^{-2n}}_{\text{相互作用}}\times \underbrace{q^{n}}_{\text{位相空間}}\times \underbrace{q^{-(n-1)}}_{\text{分母}} = q^{\,3-2n}\,\dd q $$

となる.$n=2$ では $q^{-1}$(対数発散),$n=3$ では $q^{-3}$,$n=4$ では $q^{-5}$ と,次数が上がるほど発散が悪化する.個々の項を見る限り展開はまったく使い物にならない.

Gell-Mann と Brueckner(1957)の卓越した着想は,各次数で最も強く発散する項(リングダイアグラム)だけを選び出し,それらを無限次まで足し上げることであった.リング項は $q$ の各べきについて等比数列を作るので,和は等比級数として閉じた形に書ける.模式的には

$$ \sum_{n=2}^{\infty}(-1)^n\frac{\lambda^n}{n}\;=\;-\ln(1+\lambda)+\lambda, \qquad \lambda \propto \frac{e^2\,\chi_0(q,\omega)}{q^2} $$

という構造である($\chi_0$ は自由電子ガスの応答関数,第8章).左辺の各項は $q\to0$ で発散するが,右辺の対数は $q\to0$ で穏やかにしか増大しないので,$\bm{q}$ 積分が収束する.この処方を乱雑位相近似(random phase approximation, RPA)と呼ぶ.RPAは,裸のクーロン相互作用を遮蔽された相互作用 $4\pi e^2/[q^2\varepsilon(q,\omega)]$ で置き換えることと等価である($\varepsilon$ は誘電関数,第8章).

対数がどこから来るかは,遮蔽波数を使った簡単な見積もりで納得できる.発散していた積分 $\int\dd q/q$ の下限を遮蔽波数 $\kappa_{\rm TF}$ で切ると

$$ \int_{\kappa_{\rm TF}}^{k_F}\frac{\dd q}{q} = \ln\frac{k_F}{\kappa_{\rm TF}} = -\frac{1}{2}\ln\frac{\kappa_{\rm TF}^2}{k_F^2}. $$

第8章で導くように,トーマス・フェルミ遮蔽波数は原子単位で $\kappa_{\rm TF}^2 = 4k_F/\pi$ であるから $\kappa_{\rm TF}^2/k_F^2 = 4/(\pi k_F)\propto r_s$(なぜなら $k_F\propto1/r_s$).したがって

$$ \int_{\kappa_{\rm TF}}^{k_F}\frac{\dd q}{q} = -\frac12\ln r_s + \text{const} $$

となり,$E_2^{\rm ring}$ が負であることと合わせると相関エネルギーに $+\,\mathrm{const}\times\ln r_s$ の形の項が現れる.$\ln r_s$ の出現は遮蔽の直接の帰結なのである.

7.9.4 高密度極限:Gell-Mann–Brueckner 展開

RPAの和を厳密に実行し,さらに2次の交換項を加えると,$r_s\to0$ での漸近展開が得られる.

$$ \begin{equation} \frac{E}{N} = \left[\frac{2.21}{r_s^2} - \frac{0.916}{r_s} + 0.0622\ln r_s - 0.094 + O(r_s\ln r_s)\right]\ \mathrm{Ry} \label{eq:7-GB} \end{equation} $$

すなわち相関エネルギー密度の高密度漸近形は

$$ \begin{equation} \varepsilon_c(r_s) \simeq \big(0.0622\ln r_s - 0.094\big)\ \mathrm{Ry} = \big(0.0311\ln r_s - 0.047\big)\ \mathrm{Ha} \qquad (r_s\to0) \label{eq:7-GB-c} \end{equation} $$

である.$\ln r_s$ の係数は解析的に

$$ \frac{2}{\pi^2}\left(1-\ln 2\right) = \frac{2\times(1-0.693147)}{9.8696} = \frac{0.613706}{9.8696} = 0.062181 $$

と求まっている(Gell-Mann–Brueckner 1957).定数項 $-0.094$ Ry の方は数値積分を要する.

物理的意味:相関は交換よりずっと小さい

\eqref{eq:7-GB} の4項を $r_s\to0$ で比べると,$r_s^{-2} \gg r_s^{-1} \gg \ln r_s \gg 1$ である.相関エネルギーは交換エネルギーより1桁以上ゆっくりとしか増大しないので,高密度では相対的に無視できる.これが7.5節の「高密度極限で摂動論が正当化される」という主張の具体的な確認である.

逆に $r_s$ が大きくなるにつれて相関の相対的な重みは増す.7.9.7の数値で見るように,$r_s=4$ では $\varepsilon_c$ は $\varepsilon_x$ の約 $28\%$ に達する.金属の凝集エネルギーや化学結合エネルギーを議論する精度ではこの $28\%$ は決して無視できない.

7.9.5 低密度極限:ウィグナー結晶

反対の極限 $r_s\to\infty$(低密度)では,7.5節で見たように相互作用が運動エネルギーを圧倒する.電子は運動エネルギーを犠牲にしても互いに最も遠ざかろうとし,その結果電子が周期的に配列した結晶を作る.これをウィグナー結晶(Wigner 1934, 1938)という.

この状態のエネルギーの主要項は,点電荷の格子が一様な負の背景(いまの場合は逆に一様な正背景)の中に置かれたときの静電エネルギー,すなわちマーデルングエネルギーである.ウィグナー・ザイツ球近似でこれを見積もってみよう.

導出:ウィグナー・ザイツ球のマーデルングエネルギー

各電子に半径 $r_0$ の球(ウィグナー・ザイツ球)を割り当てる.球の中には中心に電子(電荷 $-e$)が1個あり,それを中和する一様な正電荷 $+e$ が球全体に広がっているとする.この球1個あたりの静電エネルギーを計算する.

(i) 点電荷と一様電荷の相互作用.一様電荷密度は $\rho_+ = \dfrac{e}{\frac43\pi r_0^3} = \dfrac{3e}{4\pi r_0^3}$.中心の点電荷 $-e$ との相互作用エネルギーは

$$ U_1 = -e\int_0^{r_0}\frac{\rho_+}{r}\,4\pi r^2\dd r = -e\,\rho_+\,4\pi\int_0^{r_0}r\,\dd r = -e\cdot\frac{3e}{4\pi r_0^3}\cdot4\pi\cdot\frac{r_0^2}{2} = -\frac{3e^2}{2r_0}. $$

(ii) 一様電荷の自己エネルギー.半径 $r_0$ に一様に分布した電荷 $+e$ の静電自己エネルギーは,電磁気学の標準的な結果より

$$ U_2 = \frac{3}{5}\frac{e^2}{r_0}. $$

(半径 $r$ まで球殻を積み上げる過程で $\int_0^{r_0}\frac{q(r)\dd q}{r}$,$q(r)=e(r/r_0)^3$ を計算すれば得られる.)

(iii) 合計.

$$ U = U_1+U_2 = \left(-\frac32+\frac35\right)\frac{e^2}{r_0} = -\frac{9}{10}\frac{e^2}{r_0} = -\frac{0.9}{r_s}\ \mathrm{Ha} = -\frac{1.8}{r_s}\ \mathrm{Ry}. $$

球どうしはそれぞれ中性なので,異なる球の間の相互作用は近似的にゼロである(これがウィグナー・ザイツ球近似).

厳密に体心立方(bcc)格子で計算すると係数は $-1.79186$ Ry となり,上の見積もり $-1.8$ Ry とほとんど一致する.さらに電子の零点振動(ウィグナー結晶のフォノンの零点エネルギー)を加えると

$$ \begin{equation} \frac{E}{N} = \left[-\frac{1.792}{r_s} + \frac{2.66}{r_s^{3/2}} + O(r_s^{-2})\right]\ \mathrm{Ry} \qquad (r_s\to\infty) \label{eq:7-wigner} \end{equation} $$

となる(Carr 1961).第2項が $r_s^{-3/2}$ でスケールするのは,調和振動子の零点エネルギー $\frac12\hbar\omega$ が,プラズマ振動数 $\omega_p = \sqrt{4\pi\rho e^2/m}\propto\rho^{1/2}\propto r_s^{-3/2}$ に比例するからである.

例:ウィグナー結晶はどこで実現するか

\eqref{eq:7-GB}(液体)と \eqref{eq:7-wigner}(結晶)を比べると,$r_s$ が十分大きければ結晶の方が低くなる.量子モンテカルロ計算によれば,3次元ジェリウムの結晶化は $r_s\approx100$ 付近で起こる.実在の金属($r_s\lesssim6$)からははるかに遠い.しかし2次元系では結晶化が $r_s\approx30$ 程度で起こり,液体ヘリウム表面上の電子や半導体ヘテロ構造の2次元電子ガスで実際に観測されている.

本書の文脈で重要なのは,ウィグナー結晶が相関エネルギーの $r_s\to\infty$ での厳密な漸近形を与えることである.VWNやPZのようなパラメータ化は,高密度側で \eqref{eq:7-GB-c},低密度側で \eqref{eq:7-wigner} という2つの厳密な極限を満たすように作られている.

7.9.6 中間密度:量子モンテカルロ法

実在金属が属する $2\lesssim r_s\lesssim6$ は,高密度展開も低密度展開も定量的には信用できない「中間領域」である.この領域を決着させたのが Ceperley と Alder による量子モンテカルロ(QMC)計算(1980)であった.

QMCの基本的な考え方は,多体波動関数を解析的に書き下すのをあきらめ,確率過程によって基底状態を数値的にサンプリングすることである.試行波動関数としてスレーター・ヤストロー型

$$ \begin{equation} \Phi_{\rm QMC}(\rr_1,\ldots,\rr_N) = \exp\left(-\sum_{i<j}u(\abs{\rr_i-\rr_j})\right)\; \det\big[\varphi_1\varphi_2\cdots\varphi_N\big] \label{eq:7-slater-jastrow} \end{equation} $$

が用いられる.右側の行列式は反対称性(フェルミ統計)を保証し,左側のヤストロー因子 $\exp(-\sum u)$ は2電子が近づいたとき($u>0$ なら)振幅を小さくして,電子どうしの避け合い=相関を明示的に取り入れる.行列式だけではHF近似に留まるが,ヤストロー因子を掛けることでスレーター行列式で書けない相関が入る.

数学ノート:符号問題と固定節近似

ボース粒子系では,量子モンテカルロ法は原理的に厳密である(虚時間発展 $e^{-\hat H\tau}$ を確率過程として実装できる).しかしフェルミ粒子系では波動関数が正負両方の値をとるため,確率解釈が破綻する.これを符号問題という.異なる符号の寄与が打ち消し合うため,統計誤差が粒子数とともに指数関数的に増大してしまう.

実用的な回避策が固定節近似(fixed-node approximation)である.試行波動関数の節($\Phi_{\rm QMC}=0$ となる $3N$ 次元空間の超曲面)の位置を正しいと仮定して固定し,その節で区切られた各領域の内部では波動関数が定符号であるとして確率過程を走らせる.この近似のもとで得られるエネルギーは,「その節構造を持つ波動関数の中で最良のもの」の変分エネルギーであり,真の基底状態エネルギーの上限を与える.節の位置が正しければ厳密解が得られる.ジェリウムでは自由電子の行列式が与える節が非常によい近似であることが知られており,Ceperley–Alderの結果の精度は $1$ mRy 程度と見積もられている.

7.9.7 VWNパラメータ化

QMCの計算結果は有限個の $r_s$ における数値データである.これを実際のDFT計算で使うには,任意の $r_s$ で微分可能な解析関数として与える必要がある.そこで Vosko,Wilk,Nusair(1980)は,高密度極限 \eqref{eq:7-GB-c} と低密度極限を正しく再現する関数形を選び,QMCデータにフィットした.VWNの解析形は次のとおりである.

定義:VWN相関エネルギー汎関数

$$ \begin{equation} \varepsilon_c(r_s) = A\left\{ \ln\frac{x^2}{X(x)} + \frac{2b}{Q}\arctan\frac{Q}{2x+b} - \frac{b\,x_0}{X(x_0)}\left[ \ln\frac{(x-x_0)^2}{X(x)} + \frac{2(b+2x_0)}{Q}\arctan\frac{Q}{2x+b} \right]\right\} \label{eq:7-vwn} \end{equation} $$

ここで

$$ x = \sqrt{r_s},\qquad X(x) = x^2 + bx + c,\qquad Q = \sqrt{4c-b^2}, $$

であり,常磁性($\zeta=0$)の場合のパラメータは

$$ A = 0.0310907,\quad b = 3.72744,\quad c = 12.9352,\quad x_0 = -0.10498 $$

である.$\varepsilon_c$ の単位はハートリーである.$r_s$ は密度から $r_s = \left(\frac{4\pi}{3}\rho\right)^{-1/3}$ で決まる.

関数形の由来を確認しておくと,$A = (1-\ln2)/\pi^2 = 0.0310907$ は \eqref{eq:7-GB-c} の $\ln r_s$ の係数($0.0311$ Ha)に一致するように選ばれている.実際 $x=\sqrt{r_s}\to0$ では $X(x)\to c$ となり,$\ln(x^2/X)\to \ln r_s - \ln c$ であるから,\eqref{eq:7-vwn} は $A\ln r_s + \text{const}$ の形に漸近する.高密度極限の解析的結果がフィット関数に組み込まれているのである.逆に $r_s\to\infty$ では $\varepsilon_c\to0$ のように $-1/r_s$ で減衰し,ウィグナー極限の構造を持つ.

例:$r_s=4$ での VWN の値を計算する

$r_s=4$ とすると $x=2$.$X(2) = 4 + 3.72744\times2 + 12.9352 = 4 + 7.45488 + 12.9352 = 24.39008$.$Q = \sqrt{4\times12.9352 - 3.72744^2} = \sqrt{51.7408-13.8938} = \sqrt{37.8470} = 6.15199$.$X(x_0) = 0.011021 - 0.391306 + 12.9352 = 12.554915$.$\arctan\dfrac{Q}{2x+b} = \arctan\dfrac{6.15199}{7.72744} = \arctan 0.79612 = 0.67124$.順に

角括弧の中身は $-1.70598+0.76757 = -0.93841$ であり,その前の $-(-0.031166)$ を掛けて $-0.029246$.総和は $-1.80784+0.81338-0.02925 = -1.02371$.よって

$$ \varepsilon_c(4) = 0.0310907\times(-1.02371) = -0.031828\ \mathrm{Ha} = -0.06366\ \mathrm{Ry}. $$

同じ $r_s=4$ での交換エネルギーは $\varepsilon_x = -0.4582/4 = -0.11454$ Ha $= -0.22908$ Ry であるから,比は

$$ \frac{\abs{\varepsilon_c}}{\abs{\varepsilon_x}} = \frac{0.03183}{0.11454} = 0.278 . $$

相関は交換の約28%である.エネルギーの絶対値としては小さいが,化学結合や凝集エネルギーのスケール(1 eV $=0.037$ Ha)と比べると $0.032$ Ha $=0.87$ eV は決して小さくない.

交換エネルギー密度と相関エネルギー密度の r_s 依存性.(a) 同じ目盛で ε_x(赤)と ε_c(青)を比べると ε_c は ε_x よりはるかに小さい.(b) ε_c の拡大図.VWN(青の実線)と高密度漸近形(緑の破線),r_s=4 で −0.0318 Ha.黄色の帯は 2≤r_s≤6.
図7.8 交換エネルギー密度 $\varepsilon_x$ と相関エネルギー密度 $\varepsilon_c$(Hartree単位).(a) 同じ目盛で比べると $\varepsilon_c$ は $\varepsilon_x$ よりはるかに小さい.(b) $\varepsilon_c$ の拡大図.実線はVWNによるQMCのパラメータ化,破線は高密度漸近形 $0.0311\ln r_s - 0.047$ Ha.RPA漸近形は $r_s\to0$ でのみ正しく,$r_s\gtrsim2$ では相関を大きく過小評価する.黄色の帯は実在金属の範囲 $2\le r_s\le6$.

数学ノート:他のパラメータ化(PZ,PW92)

Ceperley–Alderのデータを解析関数にフィットする方法は一意ではなく,複数の広く使われるパラメータ化がある.

これら3者($\mathrm{VWN}$,PZ,PW92)は $2\lesssim r_s\lesssim6$ ではほぼ同じ値を与える(差は $1$ mHa 以下).実用上,どれを選ぶかは結果にほとんど影響しない.重要なのは,いずれも同じQMCデータと同じ2つの解析的極限に基づいているということである.第11章のLDAの精度は,究極的にはCeperley–Alderの数値計算の精度に支えられている.

7.10 ジェリウムで金属を語る

本章で得られた3つの寄与(運動・交換・相関)を合わせると,ジェリウムの全エネルギーが $r_s$ の関数として完全に決まる.本節では,この曲線から平衡密度と体積弾性率を読み取り,実在のアルカリ金属と比較する.うまく説明できることと,できないことをはっきりさせ,後者が第8章・第15章への出発点になることを示す.

7.10.1 全エネルギー曲線

1電子あたりの全エネルギーは

$$ \begin{equation} \frac{E}{N}(r_s) = \frac{2.21}{r_s^2} - \frac{0.9163}{r_s} + 2\,\varepsilon_c^{\rm VWN}(r_s) \qquad [\mathrm{Ry}] \label{eq:7-Etot} \end{equation} $$

である($\varepsilon_c^{\rm VWN}$ はハートリー単位なので2を掛けてRyにした).この3つの寄与と和を図7.9に示す.

ジェリウムの1電子あたりエネルギー.(a) 運動(緑)・交換(赤)・相関(紫)の寄与と合計(青)の r_s 依存性.(b) 合計の拡大図.r_s=4.18 で極小 −0.1548 Ry をとり,Li・Na・K・Cs の r_s(緑の点)がその近くにある.黄色の帯は 2≤r_s≤6.
図7.9 ジェリウムの1電子あたりエネルギー.(a) 運動(緑)・交換(赤)・相関(紫)の3つの寄与と合計(青).黄色の帯は実在金属の範囲 $2\le r_s\le6$.(b) 合計の拡大図.$r_s=4.18$ で極小 $-0.1548$ Ry $=-2.11$ eV をとる.緑の点は代表的なアルカリ金属の $r_s$ を示す.

曲線を数値的に調べると,極小は

$$ r_s^{\rm eq} = 4.18,\qquad \left.\frac{E}{N}\right|_{\rm min} = -0.1548\ \mathrm{Ry} = -0.0774\ \mathrm{Ha} = -2.11\ \mathrm{eV} $$

にある.この値は実在のアルカリ金属の $r_s$(Li 3.25,Na 3.93,K 4.86,Rb 5.20,Cs 5.62)とほぼ同じ範囲にある.原子核も結晶構造も一切含まないモデルが,金属の電子密度を1割程度の精度で言い当てるのである.これはジェリウムモデルの最大の成功と言ってよい.

物理的意味:なぜ平衡密度が決まるのか

図7.9(a) を見ると,極小の位置は本質的に運動エネルギー(正,$r_s^{-2}$)と交換エネルギー(負,$r_s^{-1}$)の綱引きで決まっている.密度を上げれば($r_s$ を小さくすれば)交換で得をするが,パウリの排他律による運動エネルギーの損の方が速く増える.密度を下げれば運動エネルギーは減るが,交換の利得も減る.両者が釣り合う点が平衡密度である.相関エネルギーは極小の位置をわずかに高密度側へずらす($r_s = 4.83\to4.18$)効果を持つ.

これは「金属結合とは何か」に対する最も簡潔な答えでもある.金属では電子が原子から解放されて結晶全体に広がる.広がることで運動エネルギーは下がり(第2章のビリアル定理の議論),同時に電子どうしが交換・相関によって巧みに避け合うことで反発の損を減らす.ジェリウムはこの2つの効果だけを抽出したモデルなのである.

7.10.2 体積弾性率

エネルギー曲線の曲率からは体積弾性率(bulk modulus)$B$ が求まる.$B$ は $B = -V(\partial P/\partial V)$,$P=-\partial E/\partial V$ で定義される.1電子あたりの体積 $v = V/N = \frac43\pi r_s^3 a_0^3$ を使って $r_s$ の微分に書き換えると

$$ \begin{equation} B = \frac{1}{12\pi r_s a_0^3} \left[\frac{\dd^2\varepsilon}{\dd r_s^2} - \frac{2}{r_s}\frac{\dd\varepsilon}{\dd r_s}\right], \qquad \varepsilon\equiv\frac{E}{N}. \label{eq:7-bulk} \end{equation} $$

導出は次のとおりである.$\dd v/\dd r_s = 4\pi r_s^2a_0^3$ より $P = -\dd\varepsilon/\dd v = -\varepsilon'/(4\pi r_s^2a_0^3)$.これを $r_s$ で微分して $\dd P/\dd r_s = \left[-\varepsilon''+2\varepsilon'/r_s\right]/(4\pi r_s^2a_0^3)$ を得る.最後に $B = -v\,(\dd P/\dd r_s)/(\dd v/\dd r_s)$ に代入し,$v/(4\pi r_s^2a_0^3)^2 = \frac43\pi r_s^3a_0^3/(16\pi^2r_s^4a_0^6) = 1/(12\pi r_s a_0^3)$ を使えば \eqref{eq:7-bulk} になる.

例:自由電子とジェリウムの体積弾性率

運動エネルギーだけ($\varepsilon = 1.105/r_s^2$ Ha)の場合,$\varepsilon' = -2.21/r_s^3$,$\varepsilon''=6.63/r_s^4$ であるから角括弧は $6.63/r_s^4+4.42/r_s^4 = 11.05/r_s^4$,よって

$$ B_{\rm free} = \frac{11.05}{12\pi r_s^5}\ \mathrm{Ha}/a_0^3 . $$

$1\ \mathrm{Ha}/a_0^3 = 29421$ GPa である.Na($r_s=3.93$)では $r_s^5 = 937.5$ より $B_{\rm free} = 11.05/(37.70\times937.5) = 3.13\times10^{-4}$ Ha/$a_0^3 = 9.2$ GPa.実測値 $6.42$ GPa と同じ桁である(第3章表3.2と同じ値).

交換と相関を含めた \eqref{eq:7-Etot} で同じ計算をすると,$r_s=3.93$ で $B = 2.5$ GPa となる.交換・相関を入れると弾性率はかえって実験から遠ざかる.理由は,ジェリウムに欠けている電子・イオン間の静電相互作用が,実在の金属では圧縮に対する強い抵抗を生むからである(7.10.3).ジェリウムは平衡密度を当てるが,そこでの曲率までは正しく与えない.

7.10.3 凝集エネルギーは説明できるか

次に,金属が原子に分解するのに必要なエネルギー(凝集エネルギー)を評価してみよう.Na を例にとる.

例:ナトリウムの凝集エネルギーの見積もり

Na は価電子1個($Z=1$),$r_s=3.93$,実験の凝集エネルギーは $1.113$ eV/原子である.孤立 Na 原子の $3s$ 電子のエネルギーは,イオン化エネルギー $5.139$ eV から $-5.14$ eV と見積もれる.以下,内殻電子は両側で同じとして無視し,価電子1個あたりのエネルギーを比べる.

  1. ジェリウムだけ.\eqref{eq:7-Etot} に $r_s=3.93$ を代入すると $$ \frac{E}{N} = \frac{2.21}{15.445} - \frac{0.9163}{3.93} + 2\times(-0.03208) = 0.1431 - 0.2332 - 0.0642 = -0.1543\ \mathrm{Ry} = -2.10\ \mathrm{eV}. $$ これは孤立原子の $-5.14$ eV より $3.04$ eV 高い.すなわちジェリウムだけでは Na 金属は原子に分解してしまうという誤った結論になる.
  2. イオンを点電荷に戻す(マーデルング項).ジェリウムでは正電荷を一様に塗り広げた.実際には点電荷(イオン)に集中しており,電子はイオンの近くに寄れるぶんだけ静電エネルギーが下がる.7.9.5で計算したように,この差は $-1.792/r_s$ Ry $= -0.456$ Ry $= -6.20$ eV である.合計は $-2.10-6.20 = -8.30$ eV となり,今度は原子より $3.16$ eV も低くなって過剰に束縛される.
  3. イオン芯の斥力(擬ポテンシャル).価電子は内殻電子と直交しなければならないため,イオン芯の内部には入り込めない.この効果は,核の引力 $-Z/r$ を半径 $R_c$ の内側で切り落とす「空芯擬ポテンシャル」で近似できる(第15章).切り落としによるエネルギーの増加は,電子密度を一様と見なせば $$ \Delta E = Z\rho\int_{r<R_c}\frac{1}{r}\,\dd^3r = Z\rho\cdot 4\pi\cdot\frac{R_c^2}{2} = 2\pi Z\rho R_c^2 = \frac{3Z^2R_c^2}{2r_0^3}\ \mathrm{Ha} = \frac{3Z^2R_c^2}{r_0^3}\ \mathrm{Ry} $$ である($\rho = 3Z/(4\pi r_0^3)$ を代入した).Na では $R_c\simeq1.67\,a_0$ とされるから $\Delta E = 3\times2.789/60.70 = 0.138$ Ry $= +1.88$ eV.
  4. 合計.$-8.30+1.88 = -6.42$ eV.孤立原子の $-5.14$ eV と比べて $1.28$ eV 低い.実験値 $1.113$ eV と同じ桁であり,粗い見積もりとしては十分な一致である.

教訓.ジェリウムは電子ガスそのもののエネルギー($-2.10$ eV)を与えるが,凝集エネルギーは電子・イオン間の項($-6.20$ eV と $+1.88$ eV)との大きな数どうしの差として現れる.凝集エネルギー $1.1$ eV は,それぞれ数 eV の項の差引きで生じるのである.第一原理計算が高い精度を要求されるのは,まさにこの「大きな量の差として小さな量を求める」という構造のためである.

7.10.4 ウィグナー・ザイツ近似との関係

7.9.5でウィグナー・ザイツ球の静電エネルギーを計算したとき,$-0.9/r_s$ Ha $=-1.8/r_s$ Ry という値を得た.これは低密度極限のウィグナー結晶のマーデルング項 $-1.792/r_s$ Ry とほとんど同じであり,7.10.3で使ったイオンの点電荷化による利得とも同じ量である.

これは偶然ではない.3つはいずれも「一様に広がった電荷を1点に集めることで得られる静電エネルギー」を計算しており,系が中性であるために遠方の球からの寄与が打ち消し合って,自分のウィグナー・ザイツ球の中だけで話が閉じるからである.ウィグナー・ザイツ球という単位胞の近似は,電荷中性という条件と組み合わさることで驚くほどよく働く.この考え方は,第15章の擬ポテンシャル法で金属の全エネルギーを構成する際の基本的な道具になる.

7.10.5 ジェリウム表面と仕事関数

ジェリウムモデルを半無限系に拡張すると,金属表面の物理が扱える.正背景を $z<0$ にだけ置き,$z>0$ を真空とする.電子は背景の端よりわずかに外へ染み出し(トンネル効果),その結果,表面には「外側が負・内側が正」の電気二重層ができる.この二重層を横切るときに電子が失う静電エネルギーが表面双極子障壁である.

仕事関数 $W$(電子を金属内部から真空へ取り出すのに必要な最小エネルギー)は,この双極子障壁 $\Delta\phi$ とバルクの化学ポテンシャル $\bar\mu$(フェルミ準位を背景の平均ポテンシャルから測ったもの)の差

$$ \begin{equation} W = \Delta\phi - \bar{\mu} \label{eq:7-workfunction} \end{equation} $$

として書ける.Lang と Kohn(1970-71)は,コーン・シャム法(第10章)をジェリウム表面に適用してこれを数値的に求め,$r_s=2$ で約 $3.9$ eV,$r_s=6$ で約 $2.4$ eV という値を得た.Na($r_s=3.93$)では約 $3.0$ eV となり,実験値 $2.75$ eV と近い.結晶構造を一切含まないモデルが,単純金属の仕事関数を $10\%$ 程度の精度で再現する.これもジェリウムの成功例の一つである.ただし多価金属(Al など)では結晶面依存性が大きく,ジェリウムだけでは不十分になる.

7.10.6 モデルの限界

物理的意味:ジェリウムで説明できないこと

ジェリウムモデルが原理的に扱えないものを整理しておく.

しかし,これらの「できないこと」の多くは,ジェリウムを出発点として摂動的に取り込める.イオンの周期ポテンシャルを弱い摂動と見なす擬ポテンシャル理論(第15章)は,まさにその戦略である.そして何より,ジェリウムで求めた $\varepsilon_x(\rho)$ と $\varepsilon_c(\rho)$ は,局所密度近似を通じてあらゆる非一様系の計算に持ち込まれる.原子でも分子でも表面でも,LDA計算の各グリッド点では本章の結果が呼び出されているのである.

7.11 まとめと演習

7.11.1 まとめ

7.11.2 次章への橋渡し

本章は2つの宿題を残した.

第一に,遮蔽である.7.7.8で見たように,裸のクーロン相互作用を使ったHF近似はフェルミ面で状態密度をゼロにしてしまう.7.9.2では2次摂動が長波長で発散した.どちらの病理も,電子ガスが自ら電荷を遮蔽するという集団的な効果を取り入れれば治る.第8章では,外部電荷に対する電子ガスの線形応答を定式化し,トーマス・フェルミ遮蔽とリントハルト応答関数を導く.そこには本章の $F(x)$ とほとんど同じ形の対数関数が再び現れ,$2k_F$ での特異性がフリーデル振動を生むことになる(第8章8.6節).

第二に,密度汎関数への飛躍である.本章で得た $\varepsilon_x(\rho)$ と $\varepsilon_c(\rho)$ は,一様系についてのみ正当化された量である.これを非一様な系に持ち込んでよい根拠はまだない.その根拠を与えるのが第9章のホーエンベルグ・コーンの定理(基底状態エネルギーは密度の汎関数である)であり,実用的な処方を与えるのが第10章のコーン・シャム方程式である.そして第11章では,本章の $\varepsilon_x+\varepsilon_c$ をそのまま局所的に使う局所密度近似が導入される.第11章のLDA汎関数の中身は,文字通り本章の \eqref{eq:7-dirac-exchange} と \eqref{eq:7-vwn} なのである.本章の計算は,DFTという建物の土台にそのまま埋め込まれている.

7.11.3 演習問題

演習 7.1:湯川因子と極限の順序

  1. 湯川ポテンシャル $\phi(r)=e^{-\mu r}/r$ が方程式 $\left(\nabla^2-\mu^2\right)\phi = -4\pi\delta^3(\rr)$ を満たすことを,$r\neq0$ では直接微分して,$r=0$ の寄与は微小球上でガウスの定理を使って確かめよ.この式の両辺を全空間で積分することで,$\int\phi\,\dd^3r = 4\pi/\mu^2$ を再導出せよ.
  2. 相殺し残った項 \eqref{eq:7-residual} を,$r_s=4$,一辺 $L=40\,a_0$ の立方体セル,$1/\mu=8\,a_0$ という具体的な値で評価し,1電子あたりの交換エネルギー $-0.115$ Ha と比べよ.
  3. もし $\mu\to0$ を先に取ったら 2. の値はどうなるか.極限の順序が結果を変えることを,この例で確認せよ.

ヒント:1. では $\nabla^2 f(r) = \frac{1}{r}\frac{\dd^2}{\dd r^2}(rf)$ という球対称関数に対する公式を使うとよい.2. では \eqref{eq:7-residual} が $-2\pi e^2/(\mu^2V)$ であり,$V=L^3=64000\,a_0^3$,$\mu^2 = 1/64\ a_0^{-2}$ である.

演習 7.2:交換積分を「2つのフェルミ球の重なり」で計算する

7.7節では,内側の $\bm{p}$ 積分を先に実行した.ここでは順序を変えて,移行運動量 $\bm{q}=\kk-\bm{p}$ を積分変数にとる.

  1. $\displaystyle J \equiv \int_{\abs{\kk}\le k_F}\dd^3k\int_{\abs{\bm{p}}\le k_F}\dd^3p\;\frac{1}{\abs{\kk-\bm{p}}^2}$ において $(\kk,\bm{p})\to(\bm{q},\bm{p})$ と変数変換し,$J = \int\dd^3q\,\dfrac{\Omega(q)}{q^2}$ と書けることを示せ.ここで $\Omega(q)$ は,半径 $k_F$ の球を2つ,中心を $\bm{q}$ だけずらして重ねたときの共通部分の体積である.
  2. 球の重なり体積が $q\le 2k_F$ で $$ \Omega(q) = \frac{4\pi}{3}k_F^3\left[1-\frac32\left(\frac{q}{2k_F}\right)+\frac12\left(\frac{q}{2k_F}\right)^3\right] $$ であることを,球冠の体積の公式から導け($q>2k_F$ では $\Omega=0$).
  3. 1. と 2. から $J = 4\pi^2k_F^4$ を示し,7.7.6の結果 \eqref{eq:7-outer5} を再現せよ.

ヒント:3. では $\int\dd^3q = 4\pi\int q^2\dd q$ が $1/q^2$ と相殺して $4\pi\int_0^{2k_F}\Omega(q)\dd q$ になる.$x=q/(2k_F)$ と置換すると $\int_0^1(1-\frac32x+\frac12x^3)\dd x = 1-\frac34+\frac18=\frac38$ である.1/q² の発散が球面座標の $q^2$ で完全に打ち消されることが,この計算のいちばんの見どころである.

演習 7.3:ハートリー・フォックのバンド幅と状態密度の消失

ジェリウムのHF準粒子エネルギーは,原子単位で $\varepsilon^{\rm HF}(k) = \dfrac{k^2}{2} - \dfrac{k_F}{\pi}F\!\left(\dfrac{k}{k_F}\right)$, $F(x) = 1+\dfrac{1-x^2}{2x}\ln\left|\dfrac{1+x}{1-x}\right|$ である.

  1. 占有帯の幅 $\Delta = \varepsilon^{\rm HF}(k_F)-\varepsilon^{\rm HF}(0)$ を求め,自由電子の値 $\varepsilon_F = k_F^2/2$ との比が $1+\dfrac{2}{\pi k_F}$ になることを示せ.
  2. Na($r_s=3.93$,$k_F=0.4883\,a_0^{-1}$)についてこの比を数値評価し,HFがバンド幅を何倍に見積もるかを述べよ.
  3. $F'(x) = -\left(\dfrac{1}{2x^2}+\dfrac12\right)\ln\left|\dfrac{1+x}{1-x}\right| + \dfrac1x$ を導き,$x\to1$ で $F'\to-\infty$(対数発散)となることを示せ.
  4. 状態密度が $D(\varepsilon)\propto k^2/\abs{\dd\varepsilon/\dd k}$ で与えられることを使い,$\varepsilon\to\varepsilon_F$ で $D\to0$ となることを示せ.これはHF近似のどのような欠陥を意味するか.

ヒント:1. では $F(1)=1$,$F(0)=2$(7.7.5の検算)を使う.3. では $\dfrac{\dd}{\dd x}\dfrac{1-x^2}{2x} = -\dfrac{1}{2x^2}-\dfrac12$ と,対数の微分 $\dfrac{2}{1-x^2}$ が $\dfrac{1-x^2}{2x}$ と約分して $1/x$ になることを使う.4. の結論は「フェルミ準位近傍の電子比熱・輸送が正しく記述できない」ことを意味し,第8章の遮蔽によって修復される.

演習 7.4:水素原子に交換汎関数を適用する(自己相互作用誤差)

水素原子の基底状態の電子密度は原子単位で $\rho(\rr) = \dfrac{e^{-2r}}{\pi}$ である.電子は1個しかないので,電子間相互作用は本来ゼロでなければならない.

  1. ハートリーエネルギー $E_H = \dfrac12\displaystyle\iint\frac{\rho(\rr)\rho(\rr')}{\abs{\rr-\rr'}}\dd^3r\dd^3r' = \dfrac{5}{16}$ Ha を認めたうえで,厳密な交換エネルギーが $E_x^{\rm exact} = -E_H = -\dfrac{5}{16}$ Ha でなければならない理由を述べよ.
  2. ディラック交換汎関数 \eqref{eq:7-dirac-exchange} でこの密度の交換エネルギーを計算せよ.$\int_0^\infty r^2e^{-ar}\dd r = 2/a^3$ を使うとよい.
  3. 2. の結果と $-5/16 = -0.3125$ Ha を比べ,残った差 $E_H + E_x^{\rm LDA}$ を eV 単位で求めよ.

ヒント:2. では $E_x^{\rm LDA} = -C_x\cdot4\pi\int_0^\infty r^2\left(\dfrac{e^{-2r}}{\pi}\right)^{4/3}\dd r = -C_x\cdot4\pi^{-1/3}\cdot\dfrac{2}{(8/3)^3}$ となる.数値は $-0.213$ Ha 程度になるはずである.3. の差 $\approx +0.10$ Ha $\approx 2.7$ eV が,LDAが1電子系に対して残してしまう自己相互作用誤差(self-interaction error)である.この誤差は,原子のイオン化ポテンシャルの過小評価,半導体のバンドギャップの過小評価,電荷移動状態の記述の失敗などの形で現れる.第11章で自己相互作用補正(Perdew–Zunger)やハイブリッド汎関数を扱う際の出発点になる.

参考文献