密度汎関数理論入門 — 目次 第III部 密度汎関数理論の核心 / 第9章

第9章ホーヘンベルグ・コーンの定理

第II部では,自由電子ガス(第3章)から出発して,トーマス・フェルミ模型(第4章),ハートリー・フォック法(第5章),ジェリウムモデル(第7章),線形応答と遮蔽(第8章)へと進み,「多電子系のエネルギーを電子密度 $\rho(\bm{r})$ で書きたい」という動機がどこから生まれ,どこまで成功し,どこで挫折したかを見てきた.トーマス・フェルミ模型は密度だけでエネルギーを書き下すが,殻構造も化学結合も再現できない.ハートリー・フォック法は精度は高いが,密度ではなく $N$ 電子波動関数(スレーター行列式)を担ぐ.それでは,密度だけで多電子系を厳密に記述することは,そもそも原理的に可能なのか.この問いに肯定的な答えを与えたのが,1964年のホーヘンベルグ(P. Hohenberg)とコーン(W. Kohn)による2つの定理である.本章ではこの定理を,証明の一行一行まで省略せずにたどる.さらに,証明が暗黙に仮定している「$v$ 表示可能性」の問題を洗い出し,レヴィ(M. Levy)の制約付き探索とハリマン(J. E. Harriman)の軌道構成法によって理論の土台が完全に固められる過程を追う.次章(第10章)では,この土台の上にコーン・シャム方程式という実用的な計算体系が建てられる.

この章で学ぶこと
  • 前史:トーマス・フェルミ・ディラック模型の失敗の原因と,スレーターの X$\alpha$ 法における「交換ポテンシャルの平均化」の完全な導出
  • ホーヘンベルグ・コーンの第1定理(密度 $\rho$ が外部ポテンシャル $v$ を一意に決める)の背理法による証明
  • 第2定理(エネルギー汎関数の変分原理)の証明と,普遍汎関数 $F[\rho]$ の意味
  • $v$ 表示可能性・$N$ 表示可能性の区別と,ハリマンの軌道構成法(任意の妥当な密度を再現する直交軌道の明示的構成)
  • レヴィの制約付き探索による再定式化と,それが $v$ 表示可能性問題を解消する仕組み
  • 定理が保証すること(存在)と保証しないこと($F[\rho]$ の具体形),そして第10章への戦略

本章でも第1章で導入したハートリー原子単位系($\hbar = m_e = e = 4\pi\varepsilon_0 = 1$)を用いる.エネルギーの単位は Hartree(Ha)で,$1\,\mathrm{Ha} = 2\,\mathrm{Ry} = 27.211\,\mathrm{eV}$,長さの単位はボーア半径 $a_0 = 0.5292\,\text{Å}$ である.歴史的な文脈で Rydberg 単位系(Ry)の式形に触れる箇所では,その都度明示する.

9.1 前史:TFD模型とスレーターのX$\alpha$法

ホーヘンベルグ・コーンの定理の価値を理解するには,定理以前に「密度からエネルギーやポテンシャルを作る」試みがどこまで進み,何に悩んでいたかを知るのが早道である.この節では,第4章で学んだトーマス・フェルミ・ディラック(TFD)模型の限界を復習し,次にスレーター(J. C. Slater)が1951年に提案した X$\alpha$ 法を導出する.X$\alpha$ 法は「ハートリー・フォックの非局所交換ポテンシャルを,密度だけで書ける局所ポテンシャルに平均化する」という発想であり,後のコーン・シャム法の直接の先駆けである.

9.1.1 トーマス・フェルミ・ディラック模型の復習と限界

第4章で導出したように,TFD模型は運動エネルギーと交換エネルギーの両方に一様電子ガスの結果を局所的に流用する.全エネルギー汎関数は(原子核間反発を除いて)

\begin{equation} E_{\mathrm{TFD}}[\rho] = \int \dd^3 r\, \rho(\bm{r})\, t(\rho(\bm{r})) + \int \dd^3 r\, \rho(\bm{r})\, \varepsilon_x(\rho(\bm{r})) + \int \dd^3 r\, \rho(\bm{r})\, v_{\mathrm{ext}}(\bm{r}) + \frac{1}{2}\iint \dd^3 r\, \dd^3 r'\, \frac{\rho(\bm{r})\rho(\bm{r}')}{|\bm{r}-\bm{r}'|} \label{eq:9-tfd} \end{equation}

であり,電子1個あたりの運動エネルギー $t(\rho)$ と交換エネルギー $\varepsilon_x(\rho)$ は,それぞれ第3章のフェルミ球の計算と第7章のジェリウムの交換エネルギーの計算から

\begin{equation} t(\rho) = \frac{3}{10}(3\pi^2)^{2/3}\rho^{2/3}, \qquad \varepsilon_x(\rho) = -\frac{3}{4}\left(\frac{3}{\pi}\right)^{1/3}\rho^{1/3} \label{eq:9-tf-t} \end{equation}

である.電子数の保存 $\int \dd^3 r\,\rho = N$ をラグランジュ未定乗数 $\mu$ で課して $\rho$ について変分すると(汎関数微分の規則は第4章で導入した),各項の汎関数微分は次のように計算できる.

\begin{align} \frac{\delta}{\delta\rho(\bm{r})}\int \dd^3 r'\, \rho\, t(\rho) &= \frac{\dd}{\dd\rho}\left[\frac{3}{10}(3\pi^2)^{2/3}\rho^{5/3}\right] = \frac{5}{3}\cdot\frac{3}{10}(3\pi^2)^{2/3}\rho^{2/3} = \frac{1}{2}(3\pi^2)^{2/3}\rho^{2/3}(\bm{r}) \label{eq:9-dtf}\\ \frac{\delta}{\delta\rho(\bm{r})}\int \dd^3 r'\, \rho\, \varepsilon_x(\rho) &= \frac{\dd}{\dd\rho}\left[-\frac{3}{4}\left(\frac{3}{\pi}\right)^{1/3}\rho^{4/3}\right] = -\frac{4}{3}\cdot\frac{3}{4}\left(\frac{3}{\pi}\right)^{1/3}\rho^{1/3} = -\left(\frac{3}{\pi}\right)^{1/3}\rho^{1/3}(\bm{r}) \label{eq:9-dex} \end{align}

1行目では被積分関数が $\rho(\bm{r}')$ の局所関数なので汎関数微分が通常の微分に帰着すること(第4章の局所汎関数の公式)を使い,$\rho^{5/3}$ を微分した.2行目も同様に $\rho^{4/3}$ を微分した.ハートリー項の汎関数微分は,第4章で示したように二重積分の対称性から因子 $\tfrac{1}{2}$ が 2 と相殺して $\int \dd^3 r' \rho(\bm{r}')/|\bm{r}-\bm{r}'|$ となる.以上を集めると TFD のオイラー方程式

\begin{equation} \frac{1}{2}(3\pi^2)^{2/3}\rho^{2/3}(\bm{r}) - \left(\frac{3}{\pi}\right)^{1/3}\rho^{1/3}(\bm{r}) + v_{\mathrm{ext}}(\bm{r}) + \int \dd^3 r'\, \frac{\rho(\bm{r}')}{|\bm{r}-\bm{r}'|} = \mu \label{eq:9-tfd-euler} \end{equation}

が得られる.$\mu$ は化学ポテンシャル(電子を1個加えるときのエネルギー変化)である.ここまでが第4章の復習である.

問題は,この模型の実際の性能である.TFD模型には次の3つのよく知られた失敗がある.

r D(r) = 4πr²ρ(r) K殻 L殻 M殻 ハートリー・フォック トーマス・フェルミ系
図9.1 アルゴン原子の動径密度分布 $D(r)=4\pi r^2\rho(r)$ の模式図.ハートリー・フォック(実線)は K・L・M の3つの殻に対応するピークを示すが,トーマス・フェルミ系の密度(破線)は滑らかで殻構造を持たない.

失敗の主因は交換項ではなく運動エネルギー汎関数にある.このことは,アルゴン原子の運動エネルギーを方法別に比べると定量的に見える.

例:アルゴン原子の運動エネルギーの比較

方法運動エネルギー (Ha)誤差
ハートリー・フォック(参照値)526.82—
トーマス・フェルミ汎関数 \eqref{eq:9-tf-t} で評価489.95約 −7%
コーン・シャム法(LDA,第10章)525.95約 −0.2%

ハートリー・フォックの値は,ビリアル定理(第2章)$T = -E$ を使えば全エネルギー計算から直ちに得られる.TF 汎関数の誤差は約 37 Ha,eV に直すと約 1000 eV にも達する.化学結合のエネルギースケールは 0.1 Ha 程度だから,この誤差の中に化学が丸ごと埋もれてしまう.一方,軌道を使って運動エネルギーを評価するコーン・シャム法(第10章)では誤差が 1 Ha 以下に収まる.密度汎関数理論の実用化の鍵は運動エネルギーの扱いにあることを,この表は端的に示している.

9.1.2 スレーターの発想:交換ポテンシャルの平均化

次に,TFDとは別の方向から「密度でポテンシャルを書く」ことに接近したスレーターの X$\alpha$ 法を導出する.出発点は第5章のハートリー・フォック(HF)方程式である.HF方程式の交換項は,スピン $\sigma$ の軌道 $\phi_{i\sigma}$ に対して

\begin{equation} (\hat{V}_x^{\sigma}\phi_{i\sigma})(\bm{r}) = -\sum_j^{\mathrm{occ}} \phi_{j\sigma}(\bm{r}) \int \dd^3 r'\, \frac{\phi_{j\sigma}^*(\bm{r}')\,\phi_{i\sigma}(\bm{r}')}{|\bm{r}-\bm{r}'|} \label{eq:9-hfx} \end{equation}

と働く(和は同じスピン $\sigma$ の占有軌道について取る.異なるスピン間に交換項がないことは第5章で示した).この項の困った点は2つある.第一に,$\phi_{i\sigma}(\bm{r})$ に「掛かる」のではなく積分核として作用する非局所演算子であること.第二に,作用の結果が軌道ごとに異なることである.もし交換項が,全軌道に共通の局所ポテンシャル $v_x(\bm{r})$ の掛け算で置き換えられれば,HF方程式はただの1粒子シュレーディンガー方程式になり,計算は劇的に簡単になる.

スレーターの発想は2段階からなる.まず,式 \eqref{eq:9-hfx} を形式的に $\phi_{i\sigma}(\bm{r})$ で割って,軌道ごとの局所ポテンシャル

\begin{equation} v_x^{i\sigma}(\bm{r}) \equiv \frac{(\hat{V}_x^{\sigma}\phi_{i\sigma})(\bm{r})}{\phi_{i\sigma}(\bm{r})} \label{eq:9-vxi} \end{equation}

を定義する($\phi_{i\sigma}$ の節では発散するが,いまは平均化への中間段階なので気にしない).次に,軌道依存性を消すために,位置 $\bm{r}$ における軌道 $i$ の「存在割合」

\begin{equation} w_{i}^{\sigma}(\bm{r}) = \frac{|\phi_{i\sigma}(\bm{r})|^2}{\rho^{\sigma}(\bm{r})}, \qquad \rho^{\sigma}(\bm{r}) = \sum_i^{\mathrm{occ}} |\phi_{i\sigma}(\bm{r})|^2, \qquad \sum_i^{\mathrm{occ}} w_i^{\sigma}(\bm{r}) = 1 \label{eq:9-weight} \end{equation}

を重みとして $v_x^{i\sigma}$ を平均する.これから,この重み付き平均が第5章で導入した交換ホールと電子のクーロン相互作用というきれいな形にまとまることを示す.

導出:交換ポテンシャルの重み付き平均と交換ホール

重み付き平均を定義どおりに書き下ろし,1行ずつ変形する.

\begin{align} \bar{v}_x^{\sigma}(\bm{r}) &\equiv \sum_i^{\mathrm{occ}} w_i^{\sigma}(\bm{r})\, v_x^{i\sigma}(\bm{r}) \label{eq:9-avg1}\\ &= \sum_i^{\mathrm{occ}} \frac{|\phi_{i\sigma}(\bm{r})|^2}{\rho^{\sigma}(\bm{r})} \cdot \frac{(\hat{V}_x^{\sigma}\phi_{i\sigma})(\bm{r})}{\phi_{i\sigma}(\bm{r})} \label{eq:9-avg2}\\ &= \frac{1}{\rho^{\sigma}(\bm{r})} \sum_i^{\mathrm{occ}} \phi_{i\sigma}^*(\bm{r})\,(\hat{V}_x^{\sigma}\phi_{i\sigma})(\bm{r}) \label{eq:9-avg3}\\ &= -\frac{1}{\rho^{\sigma}(\bm{r})} \sum_i^{\mathrm{occ}}\sum_j^{\mathrm{occ}} \phi_{i\sigma}^*(\bm{r})\,\phi_{j\sigma}(\bm{r}) \int \dd^3 r'\, \frac{\phi_{j\sigma}^*(\bm{r}')\,\phi_{i\sigma}(\bm{r}')}{|\bm{r}-\bm{r}'|} \label{eq:9-avg4}\\ &= \int \dd^3 r'\, \frac{1}{|\bm{r}-\bm{r}'|} \left[-\frac{1}{\rho^{\sigma}(\bm{r})} \sum_{i}^{\mathrm{occ}} \phi_{i\sigma}^*(\bm{r})\phi_{i\sigma}(\bm{r}') \sum_{j}^{\mathrm{occ}} \phi_{j\sigma}(\bm{r})\phi_{j\sigma}^*(\bm{r}')\right] \label{eq:9-avg5}\\ &= \int \dd^3 r'\, \frac{\rho_x^{\sigma}(\bm{r},\bm{r}')}{|\bm{r}-\bm{r}'|} \label{eq:9-avg6} \end{align}

各ステップの内容は次のとおり.

交換ホールの総和則も1行で確かめられる.軌道の規格直交性 $\int \dd^3 r'\, \phi_{j\sigma}^*(\bm{r}')\phi_{i\sigma}(\bm{r}') = \delta_{ij}$ を使うと

\begin{align} \int \dd^3 r'\, \rho_x^{\sigma}(\bm{r},\bm{r}') &= -\frac{1}{\rho^{\sigma}(\bm{r})}\sum_{ij}^{\mathrm{occ}} \phi_{i\sigma}^*(\bm{r})\phi_{j\sigma}(\bm{r}) \int \dd^3 r'\,\phi_{j\sigma}^*(\bm{r}')\phi_{i\sigma}(\bm{r}') \notag\\ &= -\frac{1}{\rho^{\sigma}(\bm{r})}\sum_{ij}^{\mathrm{occ}} \phi_{i\sigma}^*(\bm{r})\phi_{j\sigma}(\bm{r})\,\delta_{ij} = -\frac{1}{\rho^{\sigma}(\bm{r})}\sum_{i}^{\mathrm{occ}} |\phi_{i\sigma}(\bm{r})|^2 = -1 \label{eq:9-sumrule} \end{align}

1行目では \eqref{eq:9-xhole} の絶対値2乗を二重和に展開して $\bm{r}'$ 積分を中に入れ,2行目で規格直交性により $j=i$ の項だけを残し,最後に $\rho^{\sigma}$ の定義 \eqref{eq:9-weight} を使った.すなわち,どの位置 $\bm{r}$ から見ても,交換ホールはちょうど電子1個分の電荷の欠損を表す.

以上で,平均化された交換ポテンシャルが「電子とその交換ホールの静電相互作用」$\bar{v}_x^{\sigma}(\bm{r}) = \int \dd^3 r' \rho_x^{\sigma}(\bm{r},\bm{r}')/|\bm{r}-\bm{r}'|$ という物理的に透明な形に書けることが分かった.しかしこのままでは $\rho_x^{\sigma}$ が軌道の集合に依存するので,まだ「密度の関数」ではない.スレーターの第2の近似は,交換ホールの形をジェリウム(一様電子ガス)のもので置き換えることである.

9.1.3 ジェリウムの交換ホールと X$\alpha$ ポテンシャル

これから,一様電子ガスの1体密度行列を計算し,交換ホールを閉じた形で求め,それを $\bar{v}_x^{\sigma}$ に代入して動径積分を実行する.この3段の計算で X$\alpha$ ポテンシャルの $\rho^{1/3}$ 則が出てくる.

導出:一様電子ガスの1体密度行列と交換ホール

体積 $V$ のジェリウムでは,スピン $\sigma$ の占有軌道は平面波 $\phi_{\kk\sigma}(\bm{r}) = e^{i\kk\cdot\bm{r}}/\sqrt{V}$($|\kk| \le k_F^{\sigma}$)である(第3章).1体密度行列は $\bm{s} = \bm{r}'-\bm{r}$ のみの関数になる:

\begin{align} \gamma^{\sigma}(\bm{r},\bm{r}') &= \sum_{|\kk|\le k_F^{\sigma}} \frac{e^{i\kk\cdot\bm{r}}}{\sqrt{V}}\frac{e^{-i\kk\cdot\bm{r}'}}{\sqrt{V}} = \frac{1}{V}\sum_{|\kk|\le k_F^{\sigma}} e^{-i\kk\cdot\bm{s}} \label{eq:9-dm1}\\ &= \frac{1}{(2\pi)^3}\int_{|\kk|\le k_F^{\sigma}} \dd^3 k\; e^{-i\kk\cdot\bm{s}} \label{eq:9-dm2}\\ &= \frac{1}{(2\pi)^3}\int_0^{k_F^{\sigma}} \dd k\, k^2 \int_0^{\pi} \dd\theta\, \sin\theta \int_0^{2\pi} \dd\varphi\; e^{-iks\cos\theta} \label{eq:9-dm3}\\ &= \frac{1}{(2\pi)^3}\int_0^{k_F^{\sigma}} \dd k\, k^2 \cdot 2\pi \int_{-1}^{1} \dd u\; e^{-iksu} \label{eq:9-dm4}\\ &= \frac{1}{(2\pi)^2}\int_0^{k_F^{\sigma}} \dd k\, k^2\, \frac{e^{-iks}-e^{iks}}{-iks} = \frac{1}{2\pi^2 s}\int_0^{k_F^{\sigma}} \dd k\, k \sin(ks) \label{eq:9-dm5}\\ &= \frac{1}{2\pi^2 s}\left[\left[-\frac{k\cos(ks)}{s}\right]_0^{k_F^{\sigma}} + \frac{1}{s}\int_0^{k_F^{\sigma}}\dd k\,\cos(ks)\right] \label{eq:9-dm6}\\ &= \frac{1}{2\pi^2 s}\left(\frac{\sin(k_F^{\sigma}s)}{s^2} - \frac{k_F^{\sigma}\cos(k_F^{\sigma}s)}{s}\right) = \frac{\sin x - x\cos x}{2\pi^2 s^3}, \qquad x \equiv k_F^{\sigma} s \label{eq:9-dm7} \end{align}

各ステップ:\eqref{eq:9-dm2} では $V\to\infty$ で $\kk$ 和を積分に置き換えた($\kk$ 点密度は $V/(2\pi)^3$,第3章).\eqref{eq:9-dm3} では $\bm{s}$ 方向を極軸に取って球座標に移り,\eqref{eq:9-dm4} で $\varphi$ 積分($2\pi$)を実行し $u=\cos\theta$ と置換した.\eqref{eq:9-dm5} で $u$ 積分を実行し,指数関数の差を正弦にまとめた.\eqref{eq:9-dm6} は $k\sin(ks)$ の $k$ についての部分積分($k$ を微分する側,$\sin(ks)$ を積分する側に取った)で,境界項は $k=0$ で消える.\eqref{eq:9-dm7} で残りの積分 $\int_0^{k_F}\cos(ks)\dd k = \sin(k_F s)/s$ を実行してまとめた.

1次の球ベッセル関数 $j_1(x) = (\sin x - x\cos x)/x^2$ を使い,さらにスピン $\sigma$ の密度がフェルミ球の体積から $\rho^{\sigma} = (k_F^{\sigma})^3/(6\pi^2)$ と書けること(第3章;スピン1成分あたりなので $6\pi^2$)を使うと

\begin{equation} \gamma^{\sigma}(s) = \frac{(k_F^{\sigma})^3}{2\pi^2}\,\frac{j_1(x)}{x} = 3\rho^{\sigma}\,\frac{j_1(x)}{x} \label{eq:9-dm-final} \end{equation}

を得る.$s\to 0$ では $j_1(x)/x \to 1/3$ なので $\gamma^{\sigma}(0) = \rho^{\sigma}$ となり,密度行列の対角成分が密度に一致するという当然の性質を再現している.これを交換ホールの定義 \eqref{eq:9-xhole} に代入すると

\begin{equation} \rho_x^{\sigma}(s) = -\frac{|\gamma^{\sigma}(s)|^2}{\rho^{\sigma}} = -9\rho^{\sigma}\left\{\frac{j_1(k_F^{\sigma}s)}{k_F^{\sigma}s}\right\}^2 \label{eq:9-jellium-hole} \end{equation}

となる.これがジェリウムの交換ホールである.$s=0$ で深さ $-\rho^{\sigma}$(同スピン電子の完全排除:パウリ原理),$s \sim \pi/k_F$ 程度の半径を持ち,総和則 \eqref{eq:9-sumrule} により全電荷は $-1$ である.

導出:スレーター交換ポテンシャルの動径積分

ジェリウムの交換ホール \eqref{eq:9-jellium-hole} を平均化交換ポテンシャル \eqref{eq:9-avg6} に代入する.ホールが等方的なので球座標で書ける($s = |\bm{r}'-\bm{r}|$,立体角積分は $4\pi$,クーロン核は $1/s$):

\begin{align} \bar{v}_x^{\sigma} &= \int_0^{\infty} 4\pi s^2\, \dd s\; \frac{1}{s}\left[-9\rho^{\sigma}\left\{\frac{j_1(k_F^{\sigma}s)}{k_F^{\sigma}s}\right\}^2\right] \label{eq:9-slater1}\\ &= -36\pi \rho^{\sigma} \int_0^{\infty} \dd s\; s\, \frac{\{j_1(k_F^{\sigma}s)\}^2}{(k_F^{\sigma}s)^2} \label{eq:9-slater2}\\ &= -\frac{36\pi \rho^{\sigma}}{(k_F^{\sigma})^2} \int_0^{\infty} \dd x\; \frac{\{j_1(x)\}^2}{x} \label{eq:9-slater3} \end{align}

\eqref{eq:9-slater2} では $s^2 \cdot (1/s) = s$ と定数をまとめ,\eqref{eq:9-slater3} で $x = k_F^{\sigma}s$($\dd s = \dd x/k_F^{\sigma}$,$s = x/k_F^{\sigma}$)と置換した.残る無次元積分

\begin{equation} I \equiv \int_0^{\infty} \dd x\, \frac{\{j_1(x)\}^2}{x} = \int_0^{\infty} \dd x\, \frac{(\sin x - x\cos x)^2}{x^5} \label{eq:9-I-def} \end{equation}

を,部分積分を2回使って厳密に計算する.まず $F(x) \equiv \sin x - x\cos x$ と置くと,微分は

\begin{equation} F'(x) = \cos x - (\cos x - x\sin x) = x\sin x \label{eq:9-Fprime} \end{equation}

という簡単な形になる(積の微分で $-x\cos x$ から $-\cos x + x \sin x$ が出て,$\cos x$ が相殺した).また原点近傍で $\sin x = x - x^3/6 + \cdots$,$x\cos x = x - x^3/2 + \cdots$ より $F(x) = x^3/3 + O(x^5)$ である.第1の部分積分では $1/x^5$ を積分する側($\int x^{-5}\dd x = -x^{-4}/4$),$F^2$ を微分する側に取る:

\begin{align} I &= \left[-\frac{F(x)^2}{4x^4}\right]_0^{\infty} + \frac{1}{4}\int_0^{\infty} \dd x\, \frac{2F(x)F'(x)}{x^4} \label{eq:9-ibp1}\\ &= 0 + \frac{1}{2}\int_0^{\infty} \dd x\, \frac{F(x)\, x\sin x}{x^4} = \frac{1}{2}\int_0^{\infty} \dd x\, \frac{(\sin x - x\cos x)\sin x}{x^3} \label{eq:9-ibp2} \end{align}

境界項は,$x\to 0$ で $F^2/x^4 \sim (x^3/3)^2/x^4 = x^2/9 \to 0$,$x\to\infty$ で $F^2/x^4 \le (1+x)^2/x^4 \to 0$ となるので両端で消える.\eqref{eq:9-ibp2} では $F' = x\sin x$(式 \eqref{eq:9-Fprime})を代入して $x$ をひとつ約分した.次に分子を

\begin{equation} N(x) \equiv (\sin x - x\cos x)\sin x = \sin^2 x - \frac{x}{2}\sin 2x \label{eq:9-Ndef} \end{equation}

と置く($\sin x\cos x = \tfrac{1}{2}\sin 2x$ を使った).微分は

\begin{equation} N'(x) = 2\sin x\cos x - \frac{1}{2}\sin 2x - x\cos 2x = \frac{1}{2}\sin 2x - x\cos 2x \label{eq:9-Nprime} \end{equation}

であり($2\sin x \cos x = \sin 2x$ とまとめた),原点近傍では $\sin^2 x = x^2 - x^4/3 + \cdots$,$\tfrac{x}{2}\sin 2x = x^2 - \tfrac{2}{3}x^4 + \cdots$ より $N(x) = x^4/3 + O(x^6)$ である.第2の部分積分では $1/x^3$ を積分する側($-x^{-2}/2$)に取る:

\begin{align} \int_0^{\infty} \dd x\, \frac{N(x)}{x^3} &= \left[-\frac{N(x)}{2x^2}\right]_0^{\infty} + \frac{1}{2}\int_0^{\infty} \dd x\, \frac{N'(x)}{x^2} \label{eq:9-ibp3}\\ &= 0 + \frac{1}{2}\int_0^{\infty} \dd x\, \frac{1}{x^2}\left(\frac{1}{2}\sin 2x - x\cos 2x\right) \label{eq:9-ibp4} \end{align}

境界項は $x\to 0$ で $N/x^2 \sim x^2/3 \to 0$,$x\to\infty$ で $N/x^2 \to 0$($|N| \le 1 + x/2$)より消える.ここで注意が要る.\eqref{eq:9-ibp4} を2つの積分に分けると,$\int \sin 2x/x^2\,\dd x$ も $\int \cos 2x/x\,\dd x$ も下端 $x=0$ で個別には対数発散する(それぞれ被積分関数が $2/x$,$1/x$ と振る舞う).組み合わせだけが有限なので,下端を $\epsilon > 0$ で切ってから最後に $\epsilon\to 0$ を取る.$\sin 2x/x^2$ の項を部分積分すると($1/x^2$ を積分する側)

\begin{align} \int_{\epsilon}^{\infty} \dd x\, \frac{\sin 2x}{x^2} &= \left[-\frac{\sin 2x}{x}\right]_{\epsilon}^{\infty} + \int_{\epsilon}^{\infty} \dd x\, \frac{2\cos 2x}{x} = \frac{\sin 2\epsilon}{\epsilon} + 2\int_{\epsilon}^{\infty} \dd x\, \frac{\cos 2x}{x} \label{eq:9-ibp5} \end{align}

となる(上端では $\sin 2x/x \to 0$).これを \eqref{eq:9-ibp4} に戻すと,発散する $\int_{\epsilon}^{\infty}\cos 2x/x\,\dd x$ の係数が $\tfrac{1}{2}\cdot\tfrac{1}{2}\cdot 2 - \tfrac{1}{2} = 0$ となって厳密に相殺し,

\begin{equation} \int_0^{\infty} \dd x\, \frac{N(x)}{x^3} = \lim_{\epsilon\to 0}\ \frac{1}{4}\,\frac{\sin 2\epsilon}{\epsilon} = \frac{1}{4}\cdot 2 = \frac{1}{2} \label{eq:9-ibp6} \end{equation}

が残る.したがって \eqref{eq:9-ibp2} より

\begin{equation} I = \frac{1}{2}\cdot\frac{1}{2} = \frac{1}{4} \label{eq:9-I-result} \end{equation}

である.これを \eqref{eq:9-slater3} に代入し,$\rho^{\sigma} = (k_F^{\sigma})^3/(6\pi^2)$ を使って $\rho^{\sigma}$ を消すと

\begin{align} \bar{v}_x^{\sigma} = -\frac{36\pi \rho^{\sigma}}{(k_F^{\sigma})^2}\cdot\frac{1}{4} = -\frac{9\pi \rho^{\sigma}}{(k_F^{\sigma})^2} = -\frac{9\pi}{(k_F^{\sigma})^2}\cdot\frac{(k_F^{\sigma})^3}{6\pi^2} = -\frac{3}{2\pi}k_F^{\sigma} \label{eq:9-slater-kf} \end{align}

という驚くほど簡潔な結果に到達する.最後に $k_F^{\sigma} = (6\pi^2\rho^{\sigma})^{1/3}$ を代入して密度で書き直す:

\begin{align} \bar{v}_x^{\sigma}(\bm{r}) = -\frac{3}{2\pi}(6\pi^2)^{1/3}\left[\rho^{\sigma}(\bm{r})\right]^{1/3} = -3\left(\frac{3\rho^{\sigma}(\bm{r})}{4\pi}\right)^{1/3} \label{eq:9-slater-final} \end{align}

ここで係数の同値性は $\frac{3}{2\pi}(6\pi^2)^{1/3} = \frac{3}{2}\cdot 6^{1/3}\pi^{-1/3} = 3\cdot 3^{1/3}\cdot 2^{-2/3}\pi^{-1/3} = 3(3/4\pi)^{1/3}$ による($6^{1/3} = 2^{1/3}3^{1/3}$ を使い $2$ の冪をまとめた).

これで,局所密度からその場で決まる交換ポテンシャルが得られた.スレーターはこの式に調整パラメータ $\alpha$ を掛けた形

\begin{equation} v_{x\alpha}^{\sigma}(\bm{r}) = -3\alpha\left(\frac{3\rho^{\sigma}(\bm{r})}{4\pi}\right)^{1/3} \label{eq:9-xalpha} \end{equation}

を提案した.これが X$\alpha$ 法である.スピン無偏極の系では $\rho^{\sigma} = \rho/2$ なので,全密度 $\rho$ で書けば

\begin{equation} v_{x\alpha}(\bm{r}) = -3\alpha\left(\frac{3\rho(\bm{r})}{8\pi}\right)^{1/3} \ \text{[Hartree]} \qquad\Longleftrightarrow\qquad v_{x\alpha}(\bm{r}) = -6\alpha\left(\frac{3\rho(\bm{r})}{8\pi}\right)^{1/3} \ \text{[Rydberg]} \label{eq:9-xalpha-unpol} \end{equation}

となる.歴史的な文献の多くは Rydberg 単位(1 Ha = 2 Ry)で書かれているため係数が 6 になっている点に注意する.

パラメータ $\alpha$ の意味を明確にしておこう.上の導出(交換ホールをフェルミ面全体で平均する)は $\alpha = 1$ に対応する.一方,ディラックの交換エネルギー汎関数 $E_x^{\mathrm{LDA}}[\rho] = \int \dd^3 r\, \rho\,\varepsilon_x(\rho)$(式 \eqref{eq:9-tf-t})を変分してポテンシャルを作ると,式 \eqref{eq:9-dex} で計算したとおり

\begin{align} v_x^{\mathrm{var}}(\bm{r}) = \frac{\delta E_x^{\mathrm{LDA}}}{\delta\rho(\bm{r})} = -\left(\frac{3}{\pi}\right)^{1/3}\rho^{1/3} = \frac{2}{3}\times\left[-3\left(\frac{3\rho}{8\pi}\right)^{1/3}\right] \label{eq:9-alpha23} \end{align}

となる.最後の等号は $3(3/8\pi)^{1/3} = 3\cdot 3^{1/3}\cdot 2^{-1}\cdot\pi^{-1/3}\cdot 2^{0}\cdots$ を整理すると確かめられる:$3\left(\tfrac{3}{8\pi}\right)^{1/3} = \tfrac{3}{2}\left(\tfrac{3}{\pi}\right)^{1/3}$ であり,$\left(\tfrac{3}{\pi}\right)^{1/3} = \tfrac{2}{3}\cdot\tfrac{3}{2}\left(\tfrac{3}{\pi}\right)^{1/3}$ だからである.つまり変分から出る係数はスレーターの平均の $2/3$ 倍であり,これはガシュパール(1954)と後のコーン・シャム(1965)が採用した値 $\alpha = 2/3$ に当たる.平均の取り方の違い(スレーターはフェルミ球全体で平均,変分はフェルミ面上の電子に対応)がこの因子の起源である.実際の原子計算では,HF全エネルギーを再現するように $\alpha$ を合わせると軽い原子で $\alpha \approx 0.77$,重い原子で $\alpha \approx 0.70$ 程度となり,$2/3 \le \alpha \le 1$ の間の値が使われた.$\alpha$ が $2/3$ より大きめに選ばれるのは,交換だけでなく相関効果の一部を経験的に繰り込んでいるためと解釈できる.

物理的意味:「密度からポテンシャルを作る」系譜と,その不安

TFD模型(1927–30)は「エネルギーを密度で書く」試み,X$\alpha$ 法(1951)は「交換ポテンシャルを密度で書く」試みである.どちらも $\rho^{1/3}$ という同じ組合せに到達しており,一様電子ガスを局所的に流用するという同じ思想(局所密度近似)の産物である.X$\alpha$ 法は実際の原子・分子・固体でそれなりに機能し,1960年代の物質科学計算の主力だった.しかし,これらはあくまで近似の積み重ねである.HF交換の軌道依存性を重みで平均した段階でも,交換ホールをジェリウムで置き換えた段階でも,制御されない誤差が入る.そして最大の問題はもっと根本的なところにある:「多電子系の基底状態が密度 $\rho(\bm{r})$ だけで厳密に決まる」という保証が,当時どこにもなかったことである.波動関数 $\Psi(\bm{r}_1,\ldots,\bm{r}_N)$ は $3N$ 個の変数を持つのに,$\rho(\bm{r})$ はたった3変数である.これほど激しい情報の圧縮をしても基底状態の物理が失われない,と信じてよい理由はあるのか.この問いに「よい」と答えたのがホーヘンベルグ・コーンの定理である.次節から本題に入る.

9.2 定理の主張

この節では,まず問題の設定(どんなハミルトニアンの族を考えるのか)を固定し,密度の定義と基本的な恒等式 $\braket{\hat{V}_{\mathrm{ext}}} = \int v\rho\,\dd^3r$ を証明してから,ホーヘンベルグ・コーンの2つの定理を正確に述べる.証明は9.3節と9.4節で与える.

9.2.1 問題の設定:何が系ごとに違い,何が共通か

$N$ 個の電子が外部ポテンシャル $v(\bm{r})$(原子核が作る $-\sum_I Z_I/|\bm{r}-\bm{R}_I|$ がその代表例)の中で運動する系を考える.ハミルトニアンは

\begin{equation} \hat{H} = \hat{T} + \hat{V}_{ee} + \hat{V}_{\mathrm{ext}}, \qquad \hat{T} = -\frac{1}{2}\sum_{i=1}^{N}\nabla_i^2, \qquad \hat{V}_{ee} = \sum_{i<j}\frac{1}{|\bm{r}_i-\bm{r}_j|}, \qquad \hat{V}_{\mathrm{ext}} = \sum_{i=1}^{N} v(\bm{r}_i) \label{eq:9-ham} \end{equation}

である.ここで決定的に重要な観察がある.運動エネルギー $\hat{T}$ と電子間相互作用 $\hat{V}_{ee}$ は,電子数 $N$ さえ決めればどんな系でも同一の演算子である.水素分子でもシリコン結晶でもタンパク質でも,違うのは $v(\bm{r})$(と $N$)だけである.したがって,系を指定する情報は組 $(N, v)$ に尽きる.$(N,v)$ を与えればシュレーディンガー方程式 $\hat{H}\Psi = E\Psi$ が決まり,その基底状態 $\Psi_0$ が決まり(縮退がなければ位相を除いて一意),基底状態のあらゆる物性が決まる.

電子密度は,基底状態 $\Psi$ から

\begin{equation} \rho(\bm{r}) = N \sum_{\sigma_1\cdots\sigma_N}\int \dd^3 r_2 \cdots \dd^3 r_N\; \left|\Psi(\bm{r}\sigma_1, \bm{r}_2\sigma_2, \ldots, \bm{r}_N\sigma_N)\right|^2 \label{eq:9-density-def} \end{equation}

で定義される(スピン変数 $\sigma_i$ はすべて和を取る).$\int \dd^3 r\, \rho(\bm{r}) = N$ は $\Psi$ の規格化から従う.この定義で本章を通じて使う基本恒等式を先に証明しておく.

導出:外部ポテンシャルの期待値は $\int v\rho$ で書ける

記法を軽くするため $x_i = (\bm{r}_i, \sigma_i)$ とまとめ,$\int \dd x_i$ で「$\bm{r}_i$ 積分と $\sigma_i$ 和」を表す.示したいのは

\begin{equation} \braket{\Psi|\hat{V}_{\mathrm{ext}}|\Psi} = \int \dd^3 r\, v(\bm{r})\,\rho(\bm{r}) \label{eq:9-vrho} \end{equation}

である.期待値を定義どおり書き下ろす:

\begin{align} \braket{\Psi|\hat{V}_{\mathrm{ext}}|\Psi} &= \sum_{i=1}^{N} \int \dd x_1 \cdots \dd x_N\; \Psi^*(x_1,\ldots,x_N)\, v(\bm{r}_i)\, \Psi(x_1,\ldots,x_N) \label{eq:9-vd1}\\ &= \sum_{i=1}^{N} \int \dd x_1 \cdots \dd x_N\; v(\bm{r}_1)\, \left|\Psi(x_1,\ldots,x_N)\right|^2 \label{eq:9-vd2}\\ &= N \int \dd x_1\; v(\bm{r}_1) \int \dd x_2 \cdots \dd x_N\; \left|\Psi(x_1,\ldots,x_N)\right|^2 \label{eq:9-vd3}\\ &= \int \dd^3 r\; v(\bm{r})\, \rho(\bm{r}) \label{eq:9-vd4} \end{align}

各ステップの根拠:

この恒等式の意味は深い.$\hat{V}_{\mathrm{ext}}$ は $3N$ 次元空間の演算子なのに,その期待値は3次元の密度だけで正確に計算できる.ハミルトニアン \eqref{eq:9-ham} の3つの項のうち,系ごとに違う唯一の項がちょうど密度だけで書ける項である——この幸運な構造が,以下のすべての証明の土台になる.

9.2.2 2つの定理

準備ができたので定理を述べる.以下,「基底状態」は式 \eqref{eq:9-ham} 型のハミルトニアンの規格化された基底状態を指し,当面は縮退がないと仮定する(縮退がある場合は9.3.4節と9.6節で扱う).

定理:ホーヘンベルグ・コーンの第1定理(1964)

非縮退基底状態の電子密度 $\rho(\bm{r})$ は,外部ポテンシャル $v(\bm{r})$ を定数を除いて一意に決定する.すなわち,2つの外部ポテンシャル $v$ と $v'$ が同じ基底状態密度を与えるならば,$v(\bm{r}) - v'(\bm{r}) = \text{const}$ である.

したがって対応の連鎖 $\rho \to v \to \hat{H} \to \Psi_0$ により,基底状態波動関数は密度の汎関数 $\Psi_0 = \Psi[\rho]$ であり,基底状態のあらゆる観測量は $\rho(\bm{r})$ の汎関数である.

定理:ホーヘンベルグ・コーンの第2定理(1964)

普遍汎関数を $$F_{\mathrm{HK}}[\rho] \equiv \braket{\Psi[\rho]|\,\hat{T}+\hat{V}_{ee}\,|\Psi[\rho]}$$ で定義する($\Psi[\rho]$ は第1定理が保証する,$\rho$ を基底状態密度に持つ波動関数).外部ポテンシャル $v$ を持つ系のエネルギー汎関数を $$E_v[\rho] \equiv F_{\mathrm{HK}}[\rho] + \int \dd^3 r\, v(\bm{r})\rho(\bm{r})$$ とすると,$E_v[\rho]$ は $v$ の真の基底状態密度 $\rho_0$ で最小値を取り,その最小値は基底状態エネルギー $E_0$ に等しい: $$E_v[\tilde{\rho}] \ge E_v[\rho_0] = E_0 \qquad (\text{任意の妥当な試行密度 } \tilde{\rho})$$

第2定理の中身を一つの式で書けば,本章の最重要結果の1つである次の変分原理になる.

$$ E_v[\tilde{\rho}] = F_{\mathrm{HK}}[\tilde{\rho}] + \int \dd^3 r\, v(\bm{r})\,\tilde{\rho}(\bm{r}) \;\ge\; E_0, \qquad \text{等号} \iff \tilde{\rho} = \rho_0 $$

物理的意味:普遍汎関数 $F_{\mathrm{HK}}[\rho]$ の「普遍」とは

$F_{\mathrm{HK}}[\rho]$ の定義には $\hat{T} + \hat{V}_{ee}$ しか現れない.前述のとおりこの2つは系に依らない演算子だから,$F_{\mathrm{HK}}$ は外部ポテンシャルに依存しない,すべての $N$ 電子系に共通の汎関数である.水素分子の $F_{\mathrm{HK}}$ とシリコン結晶の $F_{\mathrm{HK}}$ は同じ汎関数である(電子数が違うだけ).もし誰かが $F_{\mathrm{HK}}[\rho]$ の正確な表式を一度でも書き下せたら,あとは系ごとに $\int v\rho$ を足して3次元の密度について最小化するだけで,あらゆる物質の基底状態が解けることになる.$3N$ 次元のシュレーディンガー方程式が3次元の変分問題に化けるのである.もちろん後述のとおり,定理は $F_{\mathrm{HK}}$ の存在を保証するだけで表式は与えない.それでも「そういう汎関数が存在する」こと自体が,TFD や X$\alpha$ のような試みを「厳密な理論の近似」として正当化する,まったく新しい視座だった.

定理の構造を図9.2に整理する.ポテンシャルから密度へ向かう2つの写像(シュレーディンガー方程式を解く写像 C,波動関数を2乗して積分する写像 D)は誰でも認める.第1定理の主張は,この合成写像が逆向きにたどれる(単射である)という点にある.

外部ポテンシャル v(r) 基底状態波動関数 Ψ(r1,…,rN) 基底状態密度 ρ(r) 写像 C シュレーディンガー 方程式を解く 写像 D |Ψ|² を 積分(式 9.24) 段階(i)で単射を証明 段階(ii)で単射を証明 HK第1定理:ρ(r) から v(r) が(定数を除き)一意に決まる ゆえに ρ → v → H → Ψ → すべての基底状態物性
図9.2 $v \to \Psi \to \rho$ の対応図.順方向の写像 C(シュレーディンガー方程式を解く)と D(密度を計算する)は自明に存在する.ホーヘンベルグ・コーンの第1定理は,破線の逆向き写像が(定数シフトを除いて)一意に存在することを主張する.証明は2つの段階に分かれ,段階(i)が C の単射性,段階(ii)が D の単射性(基底状態に限る)に対応する.

9.3 第1定理の証明

これから第1定理を証明する.示すべきことは「異なる $v$(定数差を除く)は異なる基底状態密度を与える」,対偶で言えば「同じ密度を与える2つのポテンシャルは定数しか違わない」である.証明は2段階の背理法で行う.段階(i)で「$v$ が違えば基底状態波動関数 $\Psi$ も違う」(写像 C の単射性)を,段階(ii)で「$\Psi$ が違えば基底状態密度 $\rho$ も違う」(写像 D の基底状態上での単射性)を示す.両方合わせると $v \ne v' + \text{const} \Rightarrow \rho \ne \rho'$ となり,定理が従う.

9.3.1 段階(i):ポテンシャルが違えば波動関数も違う

$v(\bm{r}) - v'(\bm{r}) \ne \text{const}$ である2つのポテンシャルを考え,それぞれのハミルトニアンとシュレーディンガー方程式を

\begin{align} \hat{H} &= \hat{T} + \hat{V}_{ee} + \sum_i v(\bm{r}_i), & \hat{H}\Psi &= E\Psi \label{eq:9-se1}\\ \hat{H}' &= \hat{T} + \hat{V}_{ee} + \sum_i v'(\bm{r}_i), & \hat{H}'\Psi' &= E'\Psi' \label{eq:9-se2} \end{align}

とする($\Psi,\Psi'$ はそれぞれの基底状態).背理法の仮定として $\Psi = \Psi'$ とおく.同じ関数 $\Psi$ が両方の方程式を満たすので,\eqref{eq:9-se1} から \eqref{eq:9-se2} を辺々引くと

\begin{align} (\hat{H} - \hat{H}')\,\Psi &= (E - E')\,\Psi \label{eq:9-sub1}\\ \left[\sum_{i=1}^{N}\left\{v(\bm{r}_i) - v'(\bm{r}_i)\right\}\right]\Psi(x_1,\ldots,x_N) &= (E - E')\,\Psi(x_1,\ldots,x_N) \label{eq:9-sub2} \end{align}

となる.\eqref{eq:9-sub1}→\eqref{eq:9-sub2} では,$\hat{T}$ と $\hat{V}_{ee}$ が両ハミルトニアンで完全に共通なので差し引きで消え,外部ポテンシャルの差だけが残ることを使った.ここが証明の要である:微分演算子 $\hat{T}$ が消えたため,\eqref{eq:9-sub2} は微分方程式ではなく,ただの掛け算の等式(各配置点 $(\bm{r}_1,\ldots,\bm{r}_N)$ ごとに成り立つ数の等式)になっている.

$\Psi(x_1,\ldots,x_N) \ne 0$ となる配置では,両辺を $\Psi$ で割ってよい:

\begin{equation} \sum_{i=1}^{N}\left\{v(\bm{r}_i) - v'(\bm{r}_i)\right\} = E - E' \qquad (\Psi \ne 0 \text{ となるすべての配置で}) \label{eq:9-const} \end{equation}

基底状態波動関数は(まともなポテンシャルに対しては)配置空間の開集合上で恒等的に零になることはない(下の数学ノート).したがって \eqref{eq:9-const} は事実上すべての配置で成り立つ.そこで $\bm{r}_2, \ldots, \bm{r}_N$ を適当な値に固定し,$\bm{r}_1$ だけを動かすと,左辺の第1項以外は定数だから

\begin{equation} v(\bm{r}_1) - v'(\bm{r}_1) = (E - E') - \sum_{i=2}^{N}\left\{v(\bm{r}_i) - v'(\bm{r}_i)\right\} = \text{const} \label{eq:9-vconst} \end{equation}

となる.これは仮定 $v - v' \ne \text{const}$ と矛盾する.ゆえに背理法の仮定が誤りで,$v - v' \ne \text{const}$ ならば $\Psi \ne \Psi'$ が示された.なお,$v' = v + c$(定数差)の場合は $\hat{H}' = \hat{H} + Nc$ なので,固有関数は完全に同じで固有値だけが $Nc$ ずれる.エネルギーの原点をずらしても物理が変わらないのはいつものことで,「定数を除いて」という但し書きはこの自明な自由度を指している.

数学ノート:波動関数はどこで零になれるか(一意接続性)

上の議論では「$\Psi$ が配置空間の開集合(体積を持つ領域)上で恒等的に零にならない」ことを使った.これは一意接続定理(unique continuation theorem)と呼ばれる楕円型偏微分方程式の定理による:局所的に2乗可積分なポテンシャルを持つシュレーディンガー方程式の解が,ある開集合上で恒等的に零ならば,解は全空間で零である.規格化された基底状態は全空間で零ではないので,開集合上で零になることはない.節面(零になる超曲面)は存在してよいが,それは測度零の集合であり,\eqref{eq:9-const} を「ほとんどすべての配置」で成立させるには十分である.両辺とも連続関数なら,ほとんどすべてで等しいことと至るところ等しいことは同値になる.クーロンポテンシャルはこの定理の適用範囲にあることが知られている.ホーヘンベルグとコーンの原論文ではこの点は暗黙に仮定されており,厳密な検討はリープ(E. H. Lieb, 1983)らによって整備された.初読では「基底状態は空間のかたまり全体で消えたりしない」という直観で先へ進んでよい.

9.3.2 段階(ii):波動関数が違えば基底状態密度も違う

次に,段階(i)で確立した $\Psi \ne \Psi'$ を出発点に,背理法の仮定として「$\Psi$ と $\Psi'$ が同じ密度 $\rho(\bm{r})$ を与える」とおく.使う道具は量子力学の変分原理(レイリー・リッツの原理)だけである.念のため復習しておく.

数学ノート:変分原理と,等号が成立する条件

$\hat{H}$ の規格化された固有状態を $\{\Psi_n\}$,固有値を $E_0 \le E_1 \le \cdots$ とする.任意の規格化された状態 $\Phi$ を $\Phi = \sum_n c_n \Psi_n$($\sum_n |c_n|^2 = 1$)と展開すると,期待値は

\begin{align} \braket{\Phi|\hat{H}|\Phi} = \sum_{n} |c_n|^2 E_n = E_0 \sum_n |c_n|^2 + \sum_{n} |c_n|^2 (E_n - E_0) = E_0 + \sum_{n} |c_n|^2 (E_n - E_0) \;\ge\; E_0 \label{eq:9-ritz} \end{align}

となる.1つ目の等号では固有状態の直交性 $\braket{\Psi_m|\Psi_n} = \delta_{mn}$ により交差項が消えることを使い,2つ目では $E_0$ を足し引きし,3つ目で規格化を使った.最後の不等号は $E_n - E_0 \ge 0$ と $|c_n|^2 \ge 0$ による.等号が成り立つのは,$E_n > E_0$ なるすべての $n$ で $c_n = 0$ のとき,すなわち $\Phi$ が基底状態の固有空間に属するときに限る.基底状態が非縮退なら,等号成立は $\Phi = e^{i\theta}\Psi_0$(位相を除いて基底状態そのもの)と同値である.言い換えると:基底状態でない状態で期待値を取ると,不等号は必ず「真の不等号」になる.この「真の」が以下の証明の生命線である.

では計算する.$\Psi'$($\hat{H}'$ の基底状態)を $\hat{H}$ の試行関数として使う.$\Psi' \ne \Psi$ であり $\hat{H}$ の基底状態は非縮退と仮定しているから,変分原理の等号は成立せず

\begin{align} E &= \braket{\Psi|\hat{H}|\Psi} \;\lt\; \braket{\Psi'|\hat{H}|\Psi'} \label{eq:9-ii1}\\ &= \braket{\Psi'|\hat{H}'|\Psi'} + \braket{\Psi'|\,\hat{H}-\hat{H}'\,|\Psi'} \label{eq:9-ii2}\\ &= E' + \left\langle \Psi'\,\middle|\,\sum_{i=1}^N \left\{v(\bm{r}_i) - v'(\bm{r}_i)\right\}\,\middle|\,\Psi'\right\rangle \label{eq:9-ii3}\\ &= E' + \int \dd^3 r\; \rho(\bm{r})\left\{v(\bm{r}) - v'(\bm{r})\right\} \label{eq:9-ii4} \end{align}

各ステップの根拠は次のとおり.

まったく同じ議論をプライムを付け替えて繰り返す.今度は $\Psi$ を $\hat{H}'$ の試行関数にする($\Psi \ne \Psi'$,$\hat{H}'$ の基底状態も非縮退):

\begin{align} E' = \braket{\Psi'|\hat{H}'|\Psi'} \;\lt\; \braket{\Psi|\hat{H}'|\Psi} = E + \int \dd^3 r\; \rho(\bm{r})\left\{v'(\bm{r}) - v(\bm{r})\right\} \label{eq:9-ii5} \end{align}

(途中の変形は \eqref{eq:9-ii2}–\eqref{eq:9-ii4} と同一なので,同じ4段を $v \leftrightarrow v'$,$\Psi \leftrightarrow \Psi'$,$E \leftrightarrow E'$ の置き換えで読めばよい.密度が共通 $\rho$ であることを再び使った.)さて,\eqref{eq:9-ii4} と \eqref{eq:9-ii5} を辺々加えると,積分項は $\int\rho(v-v') + \int\rho(v'-v) = 0$ と正確に打ち消し合い,

\begin{equation} E + E' \;\lt\; E' + E \label{eq:9-contradiction} \end{equation}

という明白な矛盾に到達する.ゆえに背理法の仮定が誤りであり,$\Psi \ne \Psi'$ ならば両者の密度は異なる.段階(i)と合わせて,$v - v' \ne \text{const} \Rightarrow \rho \ne \rho'$,すなわち第1定理が証明された.

背理法の仮定:v − v′ ≠ const なのに同じ密度 ρ を与える (v, Ψ, E) と (v′, Ψ′, E′) の2つの系 段階(i):Ψ ≠ Ψ′ を示す もし Ψ = Ψ′ なら,2つの方程式を引き算 → 運動項と e-e 相互作用が消える → Σi{v(ri)−v′(ri)} = E−E′ → v−v′ = const で仮定と矛盾 段階(ii):変分原理を2回使う(Ψ ≠ Ψ′ かつ非縮退 ⇒ 真の不等号) E < E′ + ∫ρ(v−v′)d³r (Ψ′ を H の試行関数に) E′ < E + ∫ρ(v′−v)d³r (Ψ を H′ の試行関数に) 辺々加えると積分項が相殺:E + E′ < E + E′ —— 矛盾! 結論:同じ基底状態密度を与える2つの v は定数差しかない ρ(r) → v(r) は(定数を除き)一意 = HK第1定理
図9.3 第1定理の証明の論理構造.段階(i)では運動エネルギーと電子間相互作用が両系で共通であることを使って波動関数の相違を示し,段階(ii)では変分原理を対称に2回適用して矛盾を導く.

9.3.3 何が証明されたのか:定理の射程

得られた結果を正確に言い直しておく.$\rho(\bm{r})$ が分かれば $v(\bm{r})$ が(定数を除き)分かり,$N = \int\rho\,\dd^3r$ も分かるから,ハミルトニアン \eqref{eq:9-ham} が丸ごと復元できる.ハミルトニアンが分かれば,原理的にはシュレーディンガー方程式を解いて基底状態だけでなく励起状態もすべて手に入る.つまり,基底状態密度という3変数の関数は,$3N$ 変数の基底状態波動関数と同じ情報を(この意味で)担っている.ただし注意すべきは,第2定理の変分原理が使えるのは基底状態のエネルギーと密度に対してだけという点である.励起状態は「原理的には $\rho$ の汎関数」だが,それを取り出す実用的な処方を定理は与えない(励起状態を扱う理論は時間依存DFTなどに委ねられる.第18章で触れる).

例:密度は原子核の位置と電荷まで知っている

クーロン系 $v(\bm{r}) = -\sum_I Z_I/|\bm{r}-\bm{R}_I|$ では,第1定理を持ち出すまでもなく,密度から直接 $v$ を読み取れることが知られている(Kato のカスプ定理).基底状態密度は原子核の位置 $\bm{R}_I$ で尖り(カスプ),その尖り方が $$\left.\frac{\partial \bar{\rho}}{\partial r}\right|_{r\to R_I} = -2 Z_I\, \bar{\rho}(\bm{R}_I)$$ を満たす($\bar\rho$ は核のまわりの球平均密度,$r$ は核からの距離).つまり,カスプの位置が $\bm{R}_I$ を,カスプの鋭さが $Z_I$ を教えてくれる.核の位置と電荷が分かれば $v$ が完全に決まる.HK第1定理のクーロン系版を「実物」で確かめた形であり,密度という3次元の関数に系の同定情報が実際に埋め込まれていることを実感できる例である.もちろん HK 定理の偉さは,クーロン型に限らない任意の外部ポテンシャルで一意性を保証する点にある.

9.3.4 縮退がある場合・注意すべきこと

証明を振り返ると,非縮退性は段階(ii)の「真の不等号」にだけ使われている.縮退がある場合はどうなるかを,場合分けで丁寧に見よう.$\hat{H}$ と $\hat{H}'$($v - v' \ne \text{const}$)の基底状態(いまや複数あってよい)から,同じ密度 $\rho$ を持つ $\Psi$($\hat{H}$ の基底状態)と $\Psi'$($\hat{H}'$ の基底状態)が取れたと仮定する.

したがって,「基底状態密度がポテンシャルを決める」という結論自体は,縮退があっても生き残る.縮退の本当の問題は別のところにある:縮退した基底状態 $\Psi_1, \Psi_2, \ldots$ は一般に異なる密度 $\rho_1, \rho_2, \ldots$ を持つので,「密度 $\rho$ に対応する基底状態 $\Psi[\rho]$」という対応が一価でなくなる可能性があり,また逆に,どの単一の基底状態の密度でもないような密度(例えば $\rho_1$ と $\rho_2$ の混合)には $\Psi[\rho]$ が存在しない.すると第2定理で使う汎関数 $F_{\mathrm{HK}}[\rho] = \braket{\Psi[\rho]|\hat{T}+\hat{V}_{ee}|\Psi[\rho]}$ の定義が揺らぐ.この不備は9.6節のレヴィの構成で完全に解消されるので,それまでは非縮退を仮定して先へ進む.

9.4 第2定理の証明

第1定理は「密度が系を一意に決める」という一意性の主張であって,それ自体は計算の処方ではない.実際に密度を求めるための原理——すなわち「正しい密度は何を最小にするのか」——を与えるのが第2定理である.この節では,まず第1定理の帰結として普遍汎関数 $F_{\mathrm{HK}}[\rho]$ を厳密に定義し,次にレイリー・リッツの変分原理(9.3.2節の数学ノート)だけを道具として変分原理を導く.証明そのものは4行で終わるが,その4行のどこでどの仮定を使っているかを正確に把握しておくことが,9.5節・9.6節で理論の土台を修理するときの設計図になる.

9.4.1 普遍汎関数 $F_{\mathrm{HK}}[\rho]$ の定義

汎関数を定義するには,まず定義域を決めなければならない.$F_{\mathrm{HK}}$ の定義には $\Psi[\rho]$ が必要で,$\Psi[\rho]$ の存在を保証するのは第1定理だから,$F_{\mathrm{HK}}$ が定義できるのは「そもそも何らかの外部ポテンシャルの基底状態密度になっている」ような $\rho$ に限られる.この当たり前の事実に名前を付けておく.

定義:$v$ 表示可能な密度と普遍汎関数

電子数 $N$ を固定する.ある外部ポテンシャル $v(\bm{r})$ が存在して,ハミルトニアン $\hat{H}=\hat{T}+\hat{V}_{ee}+\sum_i v(\bm{r}_i)$ の非縮退基底状態の密度がちょうど $\rho(\bm{r})$ に一致するとき,$\rho$ は $v$ 表示可能($v$-representable)であるという.$v$ 表示可能な $N$ 電子密度の全体を $\mathcal{A}_N$ と書く.

$\rho \in \mathcal{A}_N$ に対しては,第1定理の連鎖 $\rho \to v \to \hat{H} \to \Psi_0$ により,$\rho$ を密度に持つ基底状態波動関数 $\Psi[\rho]$ が(全体の位相を除いて)一意に定まる.そこで

\begin{equation} F_{\mathrm{HK}}[\rho] \equiv \braket{\Psi[\rho]\,|\,\hat{T}+\hat{V}_{ee}\,|\,\Psi[\rho]}, \qquad \rho \in \mathcal{A}_N \label{eq:9-FHK-def} \end{equation}

と定義する.これをホーヘンベルグ・コーンの普遍汎関数と呼ぶ.

強調しておきたいのは,この定義が第1定理なしには成立しないことである.もし同じ密度 $\rho$ を与える基底状態が2つ以上あったら,$\braket{\Psi[\rho]|\hat{T}+\hat{V}_{ee}|\Psi[\rho]}$ という式は「どの $\Psi$ を使うのか」が決まらず,数値が定まらない.第1定理は一見すると「だから何なのか」と思える一意性の主張だが,実は第2定理を述べるための語彙を作るという決定的な役割を担っている.

この $F_{\mathrm{HK}}$ を使って,外部ポテンシャル $v$ を持つ系のエネルギー汎関数を

\begin{equation} E_v[\rho] \equiv F_{\mathrm{HK}}[\rho] + \int \dd^3 r\, v(\bm{r})\,\rho(\bm{r}), \qquad \rho \in \mathcal{A}_N \label{eq:9-Ev-def} \end{equation}

で定義する.右辺の第1項は $v$ をまったく含まない($\hat{T}$ と $\hat{V}_{ee}$ しか使っていない).系ごとの個性は第2項にしか現れず,しかもその依存の仕方は $\rho$ について線形という最も単純な形である.9.2.1節で「系を指定する情報は $(N,v)$ に尽きる」と述べたことと,9.2節の恒等式 \eqref{eq:9-vrho} が,ここで綺麗に噛み合っている.

9.4.2 変分不等式の証明

これから,任意の試行密度 $\tilde{\rho} \in \mathcal{A}_N$ に対して $E_v[\tilde{\rho}] \ge E_0$ が成り立つことを示す.証明の戦略は単純である:$\tilde{\rho}$ は別のポテンシャル $\tilde{v}$ の基底状態密度なので,それに対応する波動関数 $\Psi[\tilde{\rho}]$ が存在する.この波動関数は「われわれの」ハミルトニアン $\hat{H}$ にとっては基底状態でも何でもないが,規格化された反対称 $N$ 電子波動関数であることに変わりはないので,$\hat{H}$ の試行関数としては完全に合法である.あとはレイリー・リッツを適用するだけである.

導出:第2定理の変分不等式

$\tilde{\rho} \in \mathcal{A}_N$ を任意の試行密度とし,対応する基底状態波動関数を $\Psi[\tilde{\rho}]$ とする.$E_v[\tilde{\rho}]$ を定義から出発して変形する.

\begin{align} E_v[\tilde{\rho}] &= F_{\mathrm{HK}}[\tilde{\rho}] + \int \dd^3 r\, v(\bm{r})\,\tilde{\rho}(\bm{r}) \label{eq:9-var1}\\ &= \braket{\Psi[\tilde{\rho}]\,|\,\hat{T}+\hat{V}_{ee}\,|\,\Psi[\tilde{\rho}]} + \int \dd^3 r\, v(\bm{r})\,\tilde{\rho}(\bm{r}) \label{eq:9-var2}\\ &= \braket{\Psi[\tilde{\rho}]\,|\,\hat{T}+\hat{V}_{ee}\,|\,\Psi[\tilde{\rho}]} + \braket{\Psi[\tilde{\rho}]\,|\,\hat{V}_{\mathrm{ext}}\,|\,\Psi[\tilde{\rho}]} \label{eq:9-var3}\\ &= \braket{\Psi[\tilde{\rho}]\,|\,\hat{T}+\hat{V}_{ee}+\hat{V}_{\mathrm{ext}}\,|\,\Psi[\tilde{\rho}]} = \braket{\Psi[\tilde{\rho}]\,|\,\hat{H}\,|\,\Psi[\tilde{\rho}]} \label{eq:9-var4}\\ &\ge E_0 \label{eq:9-var5} \end{align}

各ステップの根拠は次のとおり.

9.4.3 等号成立条件と定理の完成

不等式 $E_v[\tilde{\rho}] \ge E_0$ が示せた.残るのは「等号がいつ成り立つか」である.9.3.2節の数学ノートで注意したとおり,レイリー・リッツの等号は $\Psi[\tilde{\rho}]$ が $\hat{H}$ の基底状態であるときに限って成り立つ.順を追って読み解こう.

以上で第2定理が完全に証明された.9.2.2節に先取りして掲げた変分原理

\begin{equation} E_v[\tilde{\rho}] \ge E_0 = E_v[\rho_0], \qquad \text{等号} \iff \tilde{\rho} = \rho_0 \label{eq:9-hk2} \end{equation}

が,これで証明済みの命題になった.読み方はこうである.電子密度という3変数の関数を試行関数とする変分問題を解けば,多体系の基底状態エネルギーと基底状態密度が同時に得られる.試行関数が $3N$ 変数の波動関数から3変数の密度に縮んだ点が革命的である.$N=100$ の系で波動関数を数値的に表現しようとすると 300 次元の格子が要るが,密度なら3次元の格子で足りる.

9.4.4 オイラー・ラグランジュ方程式と化学ポテンシャル

変分原理を「方程式」の形に翻訳しておこう.最小化には拘束条件 $\int \dd^3 r\,\rho(\bm{r}) = N$ が付いている(電子数は勝手に変えられない).第4章と同じくラグランジュ未定乗数 $\mu$ を導入し,拘束なしの停留問題

\begin{equation} \Omega[\rho] \equiv E_v[\rho] - \mu\left(\int \dd^3 r\, \rho(\bm{r}) - N\right), \qquad \frac{\delta \Omega}{\delta \rho(\bm{r})} = 0 \label{eq:9-omega} \end{equation}

に置き換える.汎関数微分を項ごとに計算する.

\begin{align} 0 = \frac{\delta \Omega}{\delta\rho(\bm{r})} &= \frac{\delta F_{\mathrm{HK}}}{\delta\rho(\bm{r})} + \frac{\delta}{\delta\rho(\bm{r})}\int \dd^3 r'\, v(\bm{r}')\rho(\bm{r}') - \mu\,\frac{\delta}{\delta\rho(\bm{r})}\int \dd^3 r'\, \rho(\bm{r}') \label{eq:9-el1}\\ &= \frac{\delta F_{\mathrm{HK}}}{\delta\rho(\bm{r})} + v(\bm{r}) - \mu \label{eq:9-el2} \end{align}

\eqref{eq:9-el1} では \eqref{eq:9-Ev-def} と \eqref{eq:9-omega} を代入し,汎関数微分の線形性で項ごとに分けた(定数 $\mu N$ は $\rho$ に依らないので微分して消える).\eqref{eq:9-el1}→\eqref{eq:9-el2} では,汎関数微分の定義 $\delta A[\rho] = \int \dd^3 r\, \frac{\delta A}{\delta\rho(\bm{r})}\delta\rho(\bm{r})$ に照らして次の2つを使った.$A = \int v\rho$ のときは $\delta A = \int v\,\delta\rho$ だから $\delta A/\delta\rho(\bm{r}) = v(\bm{r})$.$A = \int \rho$ のときは $\delta A = \int 1\cdot\delta\rho$ だから $\delta A/\delta\rho(\bm{r}) = 1$.どちらも第4章で扱った局所汎関数の公式の最も簡単な場合である.整理して

$$ \frac{\delta F_{\mathrm{HK}}[\rho]}{\delta \rho(\bm{r})} + v(\bm{r}) \;=\; \mu $$

を得る.これが密度汎関数理論の厳密なオイラー・ラグランジュ方程式である.右辺の $\mu$ が位置 $\bm{r}$ に依らない定数であることが,基底状態密度を特徴づける条件そのものになっている.

物理的意味:$\mu$ は化学ポテンシャル,そしてTFD模型の位置づけ

未定乗数 $\mu$ の意味は,拘束条件が「電子数を $N$ に固定する」ことだった点から読める.一般に,ラグランジュ未定乗数は拘束値に対するエネルギーの感度に等しく,ここでは

$$\mu = \frac{\partial E_0}{\partial N}$$

すなわち電子を1個足すときのエネルギー変化=化学ポテンシャルである(第4章でトーマス・フェルミ模型に対して述べたことと同じ内容が,いまや厳密な理論の言葉で再現された).パー(R. G. Parr)らはこの $\mu$ を電気陰性度 $\chi$ と結び付け,$\mu = -\chi$ という関係を提唱した.電子が余っている側($\mu$ が高い)から足りない側($\mu$ が低い)へ電子が流れ,両者の $\mu$ が揃ったところで平衡になる——化学反応の描像がそのまま出てくる.

もう一つ重要なのは,9.1.1節で書いた TFD のオイラー方程式 \eqref{eq:9-tfd-euler} が,上の厳密な方程式の近似としてぴったり収まることである.実際,$F_{\mathrm{HK}}$ を

$$F_{\mathrm{HK}}[\rho] \;\approx\; \underbrace{\int \dd^3 r\,\rho\, t(\rho)}_{\text{TF 運動エネルギー}} + \underbrace{\int \dd^3 r\,\rho\,\varepsilon_x(\rho)}_{\text{ディラック交換}} + \underbrace{\frac{1}{2}\iint \dd^3 r\,\dd^3 r'\,\frac{\rho(\bm{r})\rho(\bm{r}')}{|\bm{r}-\bm{r}'|}}_{\text{ハートリー}}$$

と近似して汎関数微分すれば,\eqref{eq:9-dtf} と \eqref{eq:9-dex} で計算した3つの項がそのまま現れ,\eqref{eq:9-tfd-euler} に一致する.1927年のトーマス・フェルミ模型は,1964年になってようやく「厳密な密度汎関数理論の,$F_{\mathrm{HK}}$ に対する一つの近似」という正しい居場所を与えられたのである.同様に X$\alpha$ 法(9.1.3節)も,$\delta F_{\mathrm{HK}}/\delta\rho$ の交換部分に対する近似と読み直せる.ホーヘンベルグ・コーンの定理は新しい計算法ではなく,既存の近似法に「何を近似しているのか」を教える座標系だった——これが定理の第一の歴史的意義である.

注意:この証明が黙って使っている2つの仮定

ここまでの証明は論理的には完結しているが,次の2つの前提の上に乗っている.どちらも「証明の誤り」ではなく「定義域の狭さ」の問題である.

  1. 試行密度 $\tilde{\rho}$ は $v$ 表示可能でなければならない.証明の第1歩 \eqref{eq:9-var2} で $\Psi[\tilde{\rho}]$ を持ち出した時点で,$\tilde{\rho} \in \mathcal{A}_N$ を仮定している.ところが変分計算を実際に行うときは,$\rho$ を連続的に動かして最小値を探したい.動かした先の密度が $v$ 表示可能である保証はどこにもない.もし $\tilde\rho \notin \mathcal{A}_N$ なら $F_{\mathrm{HK}}[\tilde\rho]$ は定義されていないので,$E_v[\tilde\rho] \ge E_0$ という不等式は意味をなさない.汎関数微分 $\delta F_{\mathrm{HK}}/\delta\rho$ を書くときも同じ問題が起こる:微分するには $\rho$ のまわりの近傍全体で $F_{\mathrm{HK}}$ が定義されている必要がある.
  2. 基底状態は非縮退でなければならない.縮退があると,9.3.4節で見たように $\rho \mapsto \Psi[\rho]$ の対応が一価でなくなり,\eqref{eq:9-FHK-def} の右辺が定まらない.

この2つの穴を塞ぐのが9.5節と9.6節の仕事である.結論を先に言えば,$\mathcal{A}_N$ という「よく分からない集合」を捨てて,完全に特徴づけられた広い集合($N$ 表示可能な密度の全体)の上で汎関数を定義し直す,というのがレヴィの解決策である.

9.5 $v$ 表示可能性と $N$ 表示可能性

前節の最後で指摘した「定義域の穴」を,この節で正面から扱う.問題を一言でいえばこうである.ホーヘンベルグ・コーンの証明は $v \to \rho$ という向きの写像が単射であることを示したが,その写像の像がどんな集合なのかは何も教えてくれない.ところが変分計算をするには,$\rho$ を自由に動かせる広い定義域が要る.そこで2種類の「表示可能性」を区別し,狭くて素性の知れない集合($v$ 表示可能)と,広くて完全に特徴づけられた集合($N$ 表示可能)の関係を明らかにする.最後に,$N$ 表示可能な密度なら必ずそれを再現する直交軌道の組を明示的に構成できるという,ハリマン(J. E. Harriman)の美しい結果を証明する.

9.5.1 2つの「表示可能性」

定義:$v$ 表示可能性と $N$ 表示可能性

3次元の関数 $\rho(\bm{r})$ について,次の2つの性質を区別する.

定義から直ちに $\mathcal{A}_N \subseteq \mathcal{J}_N$ である($v$ 表示可能なら,その基底状態が $N$ 表示の証拠になる).問題は,この包含が真の包含かどうか,そして $\mathcal{J}_N$ をどう特徴づけるかである.

この2つの差は「どれだけ強い条件を課しているか」の差である.$N$ 表示可能性は「その密度を持つ波動関数がとにかく1つでもあればよい」という緩い条件で,$v$ 表示可能性は「しかもその波動関数が,ある1体ポテンシャルの基底状態でなければならない」という遥かに強い条件である.基底状態であることは,シュレーディンガー方程式という強い拘束を満たすことを意味するから,後者がずっと狭い集合になりそうだという直観は正しい.

3次元の関数 ρ(r) 全体(ρ ≥ 0, ∫ρ d³r = N) N 表示可能な密度 JN 条件: ρ ≥ 0, ∫ρ d³r = N, ∫|∇√ρ|² d³r < ∞ v 表示可能な密度 AN FHK[ρ] が定義できる範囲 (特徴づけは未知) 青(内側): FHK[ρ] = ⟨Ψ[ρ]|T^+V^ee|Ψ[ρ]⟩ ── 第1定理が Ψ[ρ] を保証する範囲 緑(中間): FL[ρ] = min(Ψ→ρ) ⟨Ψ|T^+V^ee|Ψ⟩ ── レヴィの制約付き探索の定義域(9.6節) 青と緑の隙間: N 表示可能だが v 表示可能でない密度が実在する(9.5.2節)
図9.4 密度の集合の階層.外側の破線は非負で電子数 $N$ に規格化された関数の全体,緑は $N$ 表示可能な密度の集合 $\mathcal{J}_N$,青は $v$ 表示可能な密度の集合 $\mathcal{A}_N$ である.$F_{\mathrm{HK}}$ は青の上でしか定義できないが,青の形は分かっていない.レヴィの制約付き探索(9.6節)は汎関数の定義域を緑まで広げ,しかも緑は3つの不等式で完全に特徴づけられる.

9.5.2 $v$ 表示可能でない密度は実在する

「$\mathcal{A}_N$ が $\mathcal{J}_N$ より真に狭い」ことを,具体的な仕掛けで示す.鍵になるのは9.3.4節でも顔を出した縮退である.縮退した基底状態が2つあってそれぞれの密度が違うとき,その2つの密度の平均を取ってみる.平均は当然 $\rho \ge 0$ と $\int\rho = N$ を満たすし,なめらかな関数だから $N$ 表示可能である.ところがこの平均密度は,どんな外部ポテンシャルの非縮退基底状態密度にもなり得ない.

定理:縮退から作られる $v$ 表示可能でない密度

ハミルトニアン $\hat{H} = \hat{T}+\hat{V}_{ee}+\sum_i v(\bm{r}_i)$ が縮退した基底状態 $\Psi_1, \Psi_2$(同じエネルギー $E_0$,規格化済み)を持ち,それぞれの密度が $\rho_1 \ne \rho_2$ であるとする.このとき平均密度

$$\bar{\rho}(\bm{r}) \equiv \tfrac{1}{2}\left[\rho_1(\bm{r}) + \rho_2(\bm{r})\right]$$

は $N$ 表示可能であるが,いかなる外部ポテンシャルの非縮退基底状態密度でもない.すなわち $\bar\rho \notin \mathcal{A}_N$ である.

導出:上の定理の証明(背理法)

記号を用意する.$i=1,2$ に対して

\begin{equation} F_i \equiv \braket{\Psi_i|\hat{T}+\hat{V}_{ee}|\Psi_i}, \qquad \bar{F} \equiv \tfrac{1}{2}(F_1+F_2) \label{eq:9-Fbar} \end{equation}

と置く.背理法の仮定:$\bar\rho \in \mathcal{A}_N$,すなわちある外部ポテンシャル $v'$ が存在して,$\hat{H}' = \hat{T}+\hat{V}_{ee}+\sum_i v'(\bm{r}_i)$ の非縮退基底状態 $\Psi'$(エネルギー $E_0'$)の密度が $\bar\rho$ であるとする.$F' \equiv \braket{\Psi'|\hat{T}+\hat{V}_{ee}|\Psi'}$ と書く.

段階1:$E_0$ を $\bar\rho$ で表す.$\Psi_1,\Psi_2$ はどちらも $\hat{H}$ の基底状態だから,恒等式 \eqref{eq:9-vrho} を使って

\begin{equation} E_0 = \braket{\Psi_i|\hat{H}|\Psi_i} = F_i + \int \dd^3 r\, v(\bm{r})\rho_i(\bm{r}), \qquad i = 1,2 \label{eq:9-degE0} \end{equation}

が成り立つ.$i=1$ と $i=2$ の式を辺々足して2で割ると,左辺は $E_0$ のまま,右辺は積分の線形性($\int v\rho_1 + \int v\rho_2 = 2\int v\bar\rho$)により

\begin{equation} E_0 = \bar{F} + \int \dd^3 r\, v(\bm{r})\,\bar\rho(\bm{r}) \label{eq:9-degA} \end{equation}

となる.縮退しているおかげで,平均を取っても左辺のエネルギーが変わらない点がこの証明の仕掛けである.

段階2:$\Psi'$ を $\hat{H}$ の試行関数にする.$\Psi'$ の密度は $\bar\rho$ だから,\eqref{eq:9-vrho} より $\braket{\Psi'|\sum_i v(\bm{r}_i)|\Psi'} = \int v\bar\rho$ である.レイリー・リッツより

\begin{equation} F' + \int \dd^3 r\, v\,\bar\rho = \braket{\Psi'|\hat{H}|\Psi'} \;\ge\; E_0 \stackrel{\eqref{eq:9-degA}}{=} \bar{F} + \int \dd^3 r\, v\,\bar\rho \quad\Longrightarrow\quad F' \ge \bar{F} \label{eq:9-degB} \end{equation}

(両辺から共通の $\int v\bar\rho$ を引いた).

段階3:$\Psi_1,\Psi_2$ を $\hat{H}'$ の試行関数にする.今度は逆向きに,$\hat{H}'$ に対してレイリー・リッツを使う:

\begin{equation} \braket{\Psi_i|\hat{H}'|\Psi_i} = F_i + \int \dd^3 r\, v'\rho_i \;\ge\; E_0', \qquad i=1,2 \label{eq:9-degC1} \end{equation}

$i=1,2$ を足して2で割ると

\begin{equation} \bar{F} + \int \dd^3 r\, v'\,\bar\rho \;\ge\; E_0' \label{eq:9-degC2} \end{equation}

一方 $\Psi'$ は $\hat{H}'$ の基底状態で密度が $\bar\rho$ だから $E_0' = F' + \int v'\bar\rho$.これを \eqref{eq:9-degC2} に代入し,共通の $\int v'\bar\rho$ を消すと

\begin{equation} \bar{F} \ge F' \label{eq:9-degC3} \end{equation}

段階4:矛盾を出す.\eqref{eq:9-degB} と \eqref{eq:9-degC3} を合わせると $F' = \bar{F}$ である.これを \eqref{eq:9-degC2} に戻すと,この不等式は等号で成立している:

\begin{equation} \tfrac{1}{2}\left[\braket{\Psi_1|\hat{H}'|\Psi_1} + \braket{\Psi_2|\hat{H}'|\Psi_2}\right] = \bar F + \int \dd^3 r\, v'\bar\rho = F' + \int \dd^3 r\, v'\bar\rho = E_0' \label{eq:9-degD} \end{equation}

ところが \eqref{eq:9-degC1} より個々の項はどちらも $E_0'$ 以上である.$E_0'$ 以上の2つの数の平均が $E_0'$ に等しいなら,両方とも $E_0'$ に等しい(片方が $E_0'$ より大きければ平均も $E_0'$ より大きくなってしまう).したがって

\begin{equation} \braket{\Psi_1|\hat{H}'|\Psi_1} = \braket{\Psi_2|\hat{H}'|\Psi_2} = E_0' \label{eq:9-degE} \end{equation}

すなわち $\Psi_1$ と $\Psi_2$ はどちらも $\hat{H}'$ の基底状態である.しかし仮定により $\hat{H}'$ の基底状態は非縮退だったから,$\Psi_1$ と $\Psi_2$ は位相を除いて同じ波動関数でなければならず,密度も一致して $\rho_1 = \rho_2$ となる.これは前提 $\rho_1 \ne \rho_2$ に矛盾する.よって背理法の仮定が誤りで,$\bar\rho \notin \mathcal{A}_N$ である.∎

例:開殻原子の平均密度

上の定理は絵空事ではない.開殻原子の基底状態は普通に縮退している.例えばホウ素原子(電子配置 $1s^2 2s^2 2p^1$)の基底状態 $^2P$ は軌道の3重縮退を持ち,$2p$ 電子を $p_x$ 型の軌道に入れた状態 $\Psi_1$ と $p_y$ 型に入れた状態 $\Psi_2$ は同じエネルギーだが,密度は

$$\rho_{1,2}(\bm{r}) = \rho_{\mathrm{core}}(r) + \frac{3}{4\pi}R_{2p}(r)^2 \times \begin{cases}\sin^2\theta\cos^2\varphi & (p_x)\\ \sin^2\theta\sin^2\varphi & (p_y)\end{cases}$$

と明らかに異なる(前者は $x$ 軸方向に,後者は $y$ 軸方向に伸びた形をしている).その平均は $\cos^2\varphi+\sin^2\varphi=1$ より

$$\bar\rho(\bm{r}) = \rho_{\mathrm{core}}(r) + \frac{3}{8\pi}R_{2p}(r)^2\sin^2\theta$$

となり,なめらかで,非負で,$\int\bar\rho = 5$ を満たし,$\int|\nabla\bar\rho^{1/2}|^2$ も有限な,どこから見てもまっとうな密度である.それでも定理により,この $\bar\rho$ はどんな外部ポテンシャルの非縮退基底状態密度でもない.$F_{\mathrm{HK}}[\bar\rho]$ は定義されていないのである.

公平を期して補足すると,この特定の $\bar\rho$ は,同じ縮退多重項の中の別の状態($2p_{+1}=-(p_x+ip_y)/\sqrt{2}$ に電子を入れた状態)の密度としては実現できる.つまり「$v$ 表示可能」の定義を縮退した基底状態まで許すように広げれば,この例は救われる.しかし9.4節の証明が文字どおり走るのは非縮退の $\mathcal{A}_N$ の上だけであり,その $\mathcal{A}_N$ がこの程度に素朴な密度すら含まないことが問題なのである.さらに,定義を縮退や統計混合まで広げても($\rho$ を基底状態のアンサンブルの密度として実現することを許しても),なお $v$ 表示可能でない $N$ 表示可能密度が存在することがレヴィ(1982)とリープ(1983)によって示されている.

もう一つ,変分問題を立てる立場から見て致命的な事実がある.上の定理は,$\rho_1$ と $\rho_2$ という2つの「まともな」密度の中点が定義域から外れてしまう,と言っている.つまり $\mathcal{A}_N$ は凸集合ではない.最小化問題を扱うのに定義域が凸でないというのは,たとえば「$\rho$ を少しずつ動かしながらエネルギーを下げていく」という素朴な最急降下法すら正当化できないことを意味する.しかも $\mathcal{A}_N$ に属するかどうかを判定する条件は,現在に至るまで一般には知られていない.リープは,アンサンブルまで許した $v$ 表示可能密度が $N$ 表示可能密度の中で稠密であること(いくらでも近くに $v$ 表示可能な密度がある)を示したが,「稠密」は「全部」ではない.汎関数の定義域としては,これでは足場が悪すぎる.

9.5.3 $N$ 表示可能性の条件

対照的に,$N$ 表示可能性の方は完全に解決している.必要十分条件が3つの単純な条件で書けるのである.

定理:$N$ 表示可能性の条件(ギルバート 1975,ハリマン 1981,リープ 1983)

$\rho(\bm{r})$ が,運動エネルギー期待値が有限な規格化された反対称 $N$ 電子波動関数の密度として実現できるための必要十分条件は,次の3つである.

\begin{equation} \text{(i)}\ \ \rho(\bm{r}) \ge 0 \quad\text{(至るところ)}, \qquad \text{(ii)}\ \ \int \dd^3 r\, \rho(\bm{r}) = N, \qquad \text{(iii)}\ \ \int \dd^3 r\, \left|\nabla \rho(\bm{r})^{1/2}\right|^2 \lt \infty \label{eq:9-Nrep} \end{equation}

(i) と (ii) は密度の意味からして当然である(確率密度は非負,全電子数は $N$).説明が要るのは (iii) である.これは「密度の平方根の勾配が2乗可積分」という条件で,物理的には運動エネルギーが発散しないことを要求している.なぜそう言えるのかを,次の数学ノートで完全に導出する.

数学ノート:フォン・ヴァイツゼッカー汎関数と条件 (iii) の必然性

次の量をフォン・ヴァイツゼッカー運動エネルギー汎関数と呼ぶ.

\begin{equation} T_W[\rho] \equiv \frac{1}{2}\int \dd^3 r\, \left|\nabla\rho^{1/2}(\bm{r})\right|^2 = \frac{1}{8}\int \dd^3 r\, \frac{|\nabla\rho(\bm{r})|^2}{\rho(\bm{r})} \label{eq:9-TW} \end{equation}

(2つ目の等号は連鎖律 $\nabla\rho^{1/2} = \nabla\rho/(2\rho^{1/2})$ から,$|\nabla\rho^{1/2}|^2 = |\nabla\rho|^2/(4\rho)$ となることによる.)これから,任意の規格化された反対称波動関数 $\Psi$ とその密度 $\rho$ について

\begin{equation} T_W[\rho] \;\le\; T \equiv \braket{\Psi|\hat{T}|\Psi} \label{eq:9-TW-ineq} \end{equation}

が成り立つことを示す.これが言えれば,$T$ が有限な $\Psi$ から作った密度は必ず (iii) を満たすことになり,(iii) の必要性が証明される.

記号を軽くするため,$\int' \equiv \sum_{\sigma_1}\int \dd x_2\cdots \dd x_N$(スピン $\sigma_1$ の和と,電子 $2,\ldots,N$ の座標積分)と書く.密度の定義 \eqref{eq:9-density-def} は $\rho(\bm{r}) = N\int' |\Psi|^2$ である.

第1段:運動エネルギーを勾配の形に書く.部分積分により(無限遠で $\Psi\to 0$ なので境界項は消える)

\begin{align} T = \braket{\Psi\Big|-\tfrac{1}{2}\sum_{i=1}^{N}\nabla_i^2\Big|\Psi} = \frac{1}{2}\sum_{i=1}^{N}\int \dd x_1\cdots \dd x_N\, |\nabla_i\Psi|^2 = \frac{N}{2}\int \dd x_1\cdots \dd x_N\, |\nabla_1\Psi|^2 \label{eq:9-Tgrad} \end{align}

最後の等号では,$|\Psi|^2$ が完全対称であること(9.2.1節の導出で使ったのと同じ理由)から $\int|\nabla_i\Psi|^2$ が $i$ に依らず等しく,$N$ 個の等しい項の和になることを使った.

第2段:密度の勾配を評価する.密度の定義を $\bm{r}$ で微分すると

\begin{align} \nabla\rho(\bm{r}) = N\int' \nabla_1|\Psi|^2 = N\int' \left(\Psi^*\nabla_1\Psi + \Psi\nabla_1\Psi^*\right) = 2N\,\mathrm{Re}\int' \Psi^*\nabla_1\Psi \label{eq:9-gradrho} \end{align}

(積の微分法と,$z + z^* = 2\,\mathrm{Re}\,z$ を使った.)絶対値を取り,$|\mathrm{Re}\,z| \le |z|$ と三角不等式,さらにコーシー・シュワルツの不等式を順に適用する:

\begin{align} |\nabla\rho| \le 2N\int' |\Psi||\nabla_1\Psi| \le 2N\left(\int'|\Psi|^2\right)^{1/2}\left(\int'|\nabla_1\Psi|^2\right)^{1/2} = 2N\left(\frac{\rho}{N}\right)^{1/2}\left(\int'|\nabla_1\Psi|^2\right)^{1/2} \label{eq:9-cs} \end{align}

最後で $\int'|\Psi|^2 = \rho/N$(密度の定義)を代入した.両辺を2乗して $4\rho$ で割ると

\begin{equation} \frac{|\nabla\rho(\bm{r})|^2}{4\rho(\bm{r})} \le \frac{4N^2 (\rho/N)}{4\rho}\int'|\nabla_1\Psi|^2 = N\int'|\nabla_1\Psi|^2 \label{eq:9-cs2} \end{equation}

第3段:$\bm{r}$ で積分する.両辺を $\bm{r}$ について積分すると,右辺は $\int'$ と合わさって全変数の積分になる:

\begin{equation} \int \dd^3 r\, \frac{|\nabla\rho|^2}{4\rho} \le N\int \dd x_1\cdots \dd x_N\, |\nabla_1\Psi|^2 \stackrel{\eqref{eq:9-Tgrad}}{=} 2T \label{eq:9-cs3} \end{equation}

左辺は \eqref{eq:9-TW} により $2T_W[\rho]$ に等しい.したがって $2T_W \le 2T$,すなわち $T_W[\rho] \le T$ が示された.∎

この不等式には物理的な読み方がある.$T_W$ は「密度分布が空間的にどれだけ急峻か」だけから決まる運動エネルギーの下限である.実際,$N=1$ の場合は $\rho = |\phi|^2$,$\phi = \rho^{1/2}$(位相を除く)なので $T = \tfrac{1}{2}\int|\nabla\phi|^2 = T_W$ となり,等号が成立する.多電子系では,パウリ原理によって電子が互いに異なる軌道を占める分だけ運動エネルギーが $T_W$ より大きくなる.

条件 (iii) が必要であることは示せた.残るのは十分性,すなわち「(i)(ii)(iii) を満たす任意の $\rho$ に対して,それを密度に持つ波動関数が実際に作れるか」である.これに答えるのが次のハリマンの構成法である.しかも「存在する」という抽象的な保証ではなく,軌道を紙の上に書き下ろす形で答える.

9.5.4 ハリマンの構成法(1次元)

アイデアは驚くほど単純である.求める密度 $\rho$ を $N$ で割った $\sigma = \rho/N$ を「$N$ 個の軌道が共通に持つべき密度」とみなし,すべての軌道に同じ振幅 $\sigma^{1/2}$ を与える.すると各軌道の $|\phi_k|^2$ は自動的に $\sigma$ になり,$N$ 本足せば $\rho$ になる.あとは位相だけを軌道ごとに変えて,直交性を実現すればよい.位相をどう選べば直交するか——それを教えてくれるのが,密度そのものから作られる「累積分布関数」である.

導出:ハリマンの等密度軌道(1次元)

1次元で,$\rho(x) \ge 0$,$\int_{-\infty}^{\infty}\rho\,\dd x = N$ を満たす密度が与えられたとする.次の2つの量を定義する.

\begin{equation} \sigma(x) \equiv \frac{\rho(x)}{N}, \qquad q(x) \equiv \int_{-\infty}^{x} \dd x'\, \sigma(x') \label{eq:9-harr-q} \end{equation}

$\sigma$ は $\int\sigma\,\dd x = 1$ に規格化された確率密度であり,$q$ はその累積分布関数である.$\sigma \ge 0$ だから $q$ は単調非減少で,$q(-\infty) = 0$,$q(+\infty) = 1$,そして微積分学の基本定理により

\begin{equation} \frac{\dd q}{\dd x} = \sigma(x) \label{eq:9-harr-dq} \end{equation}

である.つまり $q$ は無限に広い $x$ 軸を有限区間 $[0,1]$ に押し込める座標変換である(図9.5).この $q$ を位相に使って,整数 $k$ で番号付けられた軌道の族

\begin{equation} \phi_k(x) \equiv \sqrt{\sigma(x)}\; e^{2\pi i k\, q(x)}, \qquad k = 0, \pm 1, \pm 2, \ldots \label{eq:9-harr-phi} \end{equation}

を定義する.以下,この族が求めるものであることを4段階で確かめる.

(1) どの軌道も同じ密度を持つ.$q(x)$ は実数なので $e^{2\pi i k q}$ は絶対値1の複素数である.したがって

\begin{equation} |\phi_k(x)|^2 = \sigma(x)\left|e^{2\pi i k q(x)}\right|^2 = \sigma(x) \qquad (\text{すべての } k \text{ で同じ}) \label{eq:9-harr-mod} \end{equation}

この性質から「等密度軌道」(equidensity orbital)の名がある.

(2) 規格直交性.これが構成の心臓部である.定義どおりに内積を書き下ろし,変数変換する.

\begin{align} \int_{-\infty}^{\infty}\dd x\; \phi_k^*(x)\,\phi_{k'}(x) &= \int_{-\infty}^{\infty}\dd x\; \sqrt{\sigma(x)}\,e^{-2\pi i k q(x)}\,\sqrt{\sigma(x)}\,e^{2\pi i k' q(x)} \label{eq:9-harr-o1}\\ &= \int_{-\infty}^{\infty}\dd x\; \sigma(x)\, e^{2\pi i (k'-k) q(x)} \label{eq:9-harr-o2}\\ &= \int_{0}^{1}\dd q\; e^{2\pi i m q}, \qquad m \equiv k'-k \in \mathbb{Z} \label{eq:9-harr-o3}\\ &= \delta_{kk'} \label{eq:9-harr-o4} \end{align}

各ステップの根拠は次のとおり.

(3) 密度の再現.相異なる整数 $k_1, k_2, \ldots, k_N$ を任意に $N$ 個選び,対応する軌道を占有させる.\eqref{eq:9-harr-mod} より

\begin{equation} \sum_{n=1}^{N} |\phi_{k_n}(x)|^2 = \sum_{n=1}^{N}\sigma(x) = N\sigma(x) = \rho(x) \label{eq:9-harr-dens} \end{equation}

となり,目標の密度がぴたりと再現される.これら $N$ 個の規格直交軌道からスレーター行列式 $\Phi$ を作れば,$\Phi$ は反対称で規格化された $N$ 電子波動関数であり,その密度は $\rho$ である.すなわち $\rho$ は $N$ 表示可能である.

(4) 運動エネルギーの有限性.波動関数の存在だけでなく,運動エネルギーが有限であることも確かめておく.積の微分法と \eqref{eq:9-harr-dq} を使って

\begin{align} \frac{\dd \phi_k}{\dd x} &= \frac{\dd \sqrt{\sigma}}{\dd x}\, e^{2\pi i k q} + \sqrt{\sigma}\,\left(2\pi i k \frac{\dd q}{\dd x}\right) e^{2\pi i k q} = \left[\frac{\dd\sqrt{\sigma}}{\dd x} + 2\pi i k\, \sigma^{3/2}\right] e^{2\pi i k q} \label{eq:9-harr-dphi} \end{align}

($\sqrt{\sigma}\cdot\sigma = \sigma^{3/2}$ とまとめた).角括弧の中は「実部 $+\,i\times$ 実部」の形なので,絶対値の2乗は2つの2乗和になる:

\begin{equation} \left|\frac{\dd\phi_k}{\dd x}\right|^2 = \left(\frac{\dd\sqrt{\sigma}}{\dd x}\right)^2 + 4\pi^2 k^2 \sigma^3 \label{eq:9-harr-dphi2} \end{equation}

したがって軌道 $k$ の運動エネルギーは

\begin{equation} t_k = \frac{1}{2}\int \dd x \left|\frac{\dd\phi_k}{\dd x}\right|^2 = \frac{1}{2}\int \dd x \left(\frac{\dd\sqrt{\sigma}}{\dd x}\right)^2 + 2\pi^2 k^2 \int \dd x\, \sigma^3 = \frac{T_W[\rho]}{N} + 2\pi^2 k^2 \int \dd x\, \sigma^3 \label{eq:9-harr-kin} \end{equation}

最後の等号では $\sigma = \rho/N$ より $\dd\sqrt{\sigma}/\dd x = N^{-1/2}\,\dd\sqrt{\rho}/\dd x$ となることを使い,フォン・ヴァイツゼッカー汎関数 \eqref{eq:9-TW} の1次元版を認識した.第1項が有限であることは,$N$ 表示可能性の条件 (iii) そのものである.条件 (iii) は天下りに置かれた技術的条件ではなく,「構成した軌道の運動エネルギーが発散しないための条件」という具体的な意味を持っている.

ハリマンの等密度軌道の構成.上段左は Gaussian の規格化密度 σ(x)=ρ(x)/N,上段右はその累積分布 q(x) で 0 から 1 まで単調に増える.下段は包絡線 ±√σ(灰色の破線)と,k=1(藍)・k=2(緑)の Re φ_k(x)=√σ cos(2πk q(x)).
図9.5 ハリマンの等密度軌道.上段左は与えられた規格化密度 $\sigma = \rho/N$,上段右はその累積分布 $q(x)$ で,無限に広い $x$ 軸を有限区間 $[0,1]$ に写す.下段は $\phi_k(x) = \sqrt{\sigma}\,e^{2\pi ikq(x)}$ の実部.振幅(灰色の包絡線 $\pm\sqrt{\sigma}$)はすべての $k$ で共通なので密度は同じだが,位相が $q$ に沿って $k$ 回巻くため,異なる $k$ どうしは直交する.

9.5.5 3次元への拡張

1次元の構成の核心は「密度そのものを積分して座標を作り替える」ことだった.3次元でも同じことができる.ただし1回の積分では3変数を1変数に潰してしまうので,$x, y, z$ を順に処理して3つの座標を作る.

導出:3次元の等密度軌道と単位立方体への写像

まず $\rho$ を部分的に積分した2つの補助関数を定義する.

\begin{equation} P(x) \equiv \int_{-\infty}^{\infty}\dd y \int_{-\infty}^{\infty}\dd z\, \rho(x,y,z), \qquad Q(x,y) \equiv \int_{-\infty}^{\infty}\dd z\, \rho(x,y,z) \label{eq:9-harr-PQ} \end{equation}

$P(x)$ は密度を $x$ 軸に射影した線密度,$Q(x,y)$ は $xy$ 面に射影した面密度である.定義から $\int P\,\dd x = N$,$\int Q\,\dd y = P(x)$ が成り立つ.これらを使って3つの座標

\begin{align} f_1(x) &\equiv \frac{1}{N}\int_{-\infty}^{x}\dd x'\, P(x') \label{eq:9-harr-f1}\\ f_2(x,y) &\equiv \frac{1}{P(x)}\int_{-\infty}^{y}\dd y'\, Q(x,y') \label{eq:9-harr-f2}\\ f_3(x,y,z) &\equiv \frac{1}{Q(x,y)}\int_{-\infty}^{z}\dd z'\, \rho(x,y,z') \label{eq:9-harr-f3} \end{align}

を定義する.どれも「非負関数の部分積分 ÷ その全積分」の形だから,値は必ず $[0,1]$ に入り,上端で 1,下端で 0 になる.つまり写像 $\bm{r} = (x,y,z) \mapsto \bm{f} = (f_1,f_2,f_3)$ は $\mathbb{R}^3$ 全体を単位立方体 $[0,1]^3$ に写す.

ヤコビアンの計算.変数変換の要は $\dd^3 f = |\det J|\,\dd^3 r$ のヤコビ行列 $J_{ab} = \partial f_a/\partial x_b$ である.定義の作り方から,$f_1$ は $x$ だけの関数,$f_2$ は $(x,y)$ だけの関数,$f_3$ は $(x,y,z)$ の関数なので

\begin{equation} \frac{\partial f_1}{\partial y} = \frac{\partial f_1}{\partial z} = 0, \qquad \frac{\partial f_2}{\partial z} = 0 \label{eq:9-harr-tri} \end{equation}

である.つまり $J$ は下三角行列であり,行列式は対角成分の積に等しい.対角成分は微積分学の基本定理で直ちに計算できる:

\begin{equation} \frac{\partial f_1}{\partial x} = \frac{P(x)}{N}, \qquad \frac{\partial f_2}{\partial y} = \frac{Q(x,y)}{P(x)}, \qquad \frac{\partial f_3}{\partial z} = \frac{\rho(x,y,z)}{Q(x,y)} \label{eq:9-harr-diag} \end{equation}

(各式とも,積分の上端で微分するので被積分関数がそのまま出てくる.$f_2$ の $\partial/\partial y$ では前係数 $1/P(x)$ は $y$ に依らない定数として外に出る.$f_3$ でも同様に $1/Q(x,y)$ が定数扱いになる.)したがって

\begin{equation} \det J = \frac{P(x)}{N}\cdot\frac{Q(x,y)}{P(x)}\cdot\frac{\rho(x,y,z)}{Q(x,y)} = \frac{\rho(\bm{r})}{N} \label{eq:9-harr-jac} \end{equation}

$P$ と $Q$ が見事に約分して,ヤコビアンが密度そのものになった.これは1次元の $\dd q = \sigma\,\dd x$ の完全な3次元版である.式で書けば

\begin{equation} \dd f_1\,\dd f_2\,\dd f_3 = \frac{\rho(\bm{r})}{N}\,\dd^3 r \label{eq:9-harr-measure} \end{equation}

軌道の定義と検証.3つの整数の組 $\bm{k} = (k_1,k_2,k_3) \in \mathbb{Z}^3$ で番号付けて

\begin{equation} \phi_{\bm{k}}(\bm{r}) \equiv \sqrt{\frac{\rho(\bm{r})}{N}}\; \exp\!\left[2\pi i\left(k_1 f_1 + k_2 f_2 + k_3 f_3\right)\right] \label{eq:9-harr-3d} \end{equation}

と定義する.指数の肩は実数なので,1次元のときと同じく

\begin{equation} \left|\phi_{\bm{k}}(\bm{r})\right|^2 = \frac{\rho(\bm{r})}{N} \qquad(\text{すべての }\bm{k}\text{ で共通}) \label{eq:9-harr-3dmod} \end{equation}

である.規格直交性は次のように確かめられる.$\bm{m} = \bm{k}'-\bm{k} \in \mathbb{Z}^3$ と置くと

\begin{align} \int \dd^3 r\; \phi_{\bm{k}}^*(\bm{r})\,\phi_{\bm{k}'}(\bm{r}) &= \int \dd^3 r\; \frac{\rho(\bm{r})}{N}\, e^{2\pi i\,\bm{m}\cdot\bm{f}(\bm{r})} \label{eq:9-harr-3do1}\\ &= \int_0^1\!\!\int_0^1\!\!\int_0^1 \dd f_1\,\dd f_2\,\dd f_3\; e^{2\pi i(m_1f_1+m_2f_2+m_3f_3)} \label{eq:9-harr-3do2}\\ &= \prod_{a=1}^{3}\int_0^1 \dd f_a\, e^{2\pi i m_a f_a} = \prod_{a=1}^{3}\delta_{m_a 0} = \delta_{\bm{k}\bm{k}'} \label{eq:9-harr-3do3} \end{align}

\eqref{eq:9-harr-3do1}→\eqref{eq:9-harr-3do2} では \eqref{eq:9-harr-measure} を使って積分変数を $\bm{r}$ から $\bm{f}$ に変えた(積分領域 $\mathbb{R}^3$ が単位立方体に移る).\eqref{eq:9-harr-3do2}→\eqref{eq:9-harr-3do3} では指数関数が3つの因子に分かれるので積分も分離し,1次元で計算済みの $\int_0^1 e^{2\pi i m f}\dd f = \delta_{m0}$ を各方向に適用した.

最後に,相異なる $N$ 個のベクトル $\bm{k}_1,\ldots,\bm{k}_N$ を選んで軌道を占有させれば

\begin{equation} \sum_{n=1}^{N}\left|\phi_{\bm{k}_n}(\bm{r})\right|^2 = N\cdot\frac{\rho(\bm{r})}{N} = \rho(\bm{r}) \label{eq:9-harr-3ddens} \end{equation}

となり,これら $N$ 個の規格直交軌道から作ったスレーター行列式が,目標の密度 $\rho$ を持つ.∎

物理的意味:密度を固定したまま波動関数を動かせる

ハリマンの構成が教えてくれることは2つある.

第一に,$N$ 表示可能性の条件 \eqref{eq:9-Nrep} を満たす密度なら,それを与える波動関数が必ず作れる(十分性).9.5.3節で条件の必要性を示したので,これで \eqref{eq:9-Nrep} が必要十分条件であることの筋道が通った.$v$ 表示可能性が「判定条件すら分からない」のと対照的に,$N$ 表示可能性はこの3行で完全に片が付く.

第二に,そしてこちらが次節への鍵だが,同じ密度を与える波動関数は1つではなく無数にある.$N$ 個の整数(3次元なら整数ベクトル)の選び方は無限にあり,どの選び方でも密度はぴったり $\rho$ になるからである.式 \eqref{eq:9-harr-kin} が示すように,$|k|$ が大きい軌道ほど位相が速く巻いて運動エネルギーが高くなる.つまり「密度を固定する」という拘束は波動関数をほとんど拘束せず,その拘束面の上で $\braket{\Psi|\hat{T}+\hat{V}_{ee}|\Psi}$ はさまざまな値を取る.ならば,その中で最小のものを選べばよいのではないか——これがレヴィの発想であり,次節の主題である.

9.6 レヴィの制約付き探索

9.4節で作った理論には,9.4.4節の注意で挙げた2つの穴があった.試行密度が $v$ 表示可能でなければ $F_{\mathrm{HK}}$ が定義できないこと,そして縮退があると $\Psi[\rho]$ が一意に決まらないことである.1979年にレヴィ(M. Levy)が提案した制約付き探索(constrained search)は,この2つを同時に,しかもきわめて簡単な仕掛けで解消する.仕掛けは一言で言える:「$\rho$ から $\Psi$ を一意に決めようとするのをやめて,$\rho$ を与える $\Psi$ の中から最小のものを選ぶことにする」.

9.6.1 発想:一意性を捨てて最小性を取る

ホーヘンベルグ・コーンの構成では,$\rho \mapsto \Psi[\rho]$ という写像を「基底状態であること」によって一意に定めていた.だからこそ,$\rho$ が基底状態密度でなければ($v$ 表示可能でなければ)対応する $\Psi$ が存在せず,また縮退があれば対応が多価になった.

ところが9.5節で分かったように,$N$ 表示可能な $\rho$ に対しては,その密度を持つ波動関数は無数に存在する(ハリマンが実際に構成してみせた).ならば「一意に決める」ことにこだわる必要はない.無数の候補があるなら,その中で $\braket{\Psi|\hat{T}+\hat{V}_{ee}|\Psi}$ を最小にするものを選べば,$\rho$ から一つの数が確定する.基底状態であることを要求しないから,$\rho$ が $N$ 表示可能でありさえすれば必ず候補が存在し,汎関数が定義できる.縮退があっても「最小値」は一意に決まるから困らない.

定義:レヴィの制約付き探索汎関数(1979)

$N$ 表示可能な密度 $\rho$(すなわち条件 \eqref{eq:9-Nrep} を満たす密度)に対して

\begin{equation} F_L[\rho] \equiv \min_{\Psi \to \rho}\ \braket{\Psi\,|\,\hat{T}+\hat{V}_{ee}\,|\,\Psi} \label{eq:9-FL-def} \end{equation}

と定義する.記号 $\Psi \to \rho$ は「$\Psi$ は規格化された反対称 $N$ 電子波動関数であって,その密度 \eqref{eq:9-density-def} が $\rho$ に等しい」という制約を表す.最小値を与える波動関数を $\Psi_\rho^{\min}$ と書く.

この定義には外部ポテンシャルがまったく現れない.したがって $F_L$ も $F_{\mathrm{HK}}$ と同じく普遍汎関数である.違いは定義域で,$F_{\mathrm{HK}}$ が $\mathcal{A}_N$($v$ 表示可能)の上でしか定義できないのに対し,$F_L$ は $\mathcal{J}_N$($N$ 表示可能)の全体で定義される.

ここで「探索(search)」という語について一言.$\min_{\Psi\to\rho}$ は,密度が $\rho$ に等しいという制約を課したうえで,その制約を満たす波動関数の空間を探索して最小値を取る,という意味である.9.5節のハリマンの構成は,この探索の対象となる集合が空でないことを保証する役割を果たしている.空集合上の最小値は定義できないから,これは形式的な確認ではなく本質的な前提である.

9.6.2 2段階最小化

これから,レヴィの汎関数を使った密度の変分問題が,もとのシュレーディンガーの変分問題と完全に等価であることを証明する.証明の道具は「最小化を2段階に分ける」という初等的だが強力な操作だけである.まずその操作を数学ノートで正確にしておく.

数学ノート:最小化の2段階分解

集合 $X$ が互いに交わらない部分集合の族に分割されているとする.すなわち添字集合 $C$ があって

$$X = \bigcup_{c\in C} X_c, \qquad c \ne c' \Rightarrow X_c \cap X_{c'} = \varnothing$$

が成り立つとする(どの $x \in X$ も,ちょうど1つのクラスに属する).このとき実数値関数 $f$ について

\begin{equation} \min_{x \in X} f(x) = \min_{c \in C}\ \left[\ \min_{x \in X_c} f(x)\ \right] \label{eq:9-minsplit} \end{equation}

が成り立つ.証明は両向きの不等号を示せばよい.

両方合わせて \eqref{eq:9-minsplit} を得る.∎

さらに,$f$ が「クラス内では定数の部分」を持つとき,その部分は内側の最小化の外に出せる.すなわち,$x \in X_c$ のとき $f(x) = g(x) + h(c)$ と書けるなら

\begin{equation} \min_{x\in X_c} f(x) = \min_{x\in X_c}\left[g(x) + h(c)\right] = \left[\min_{x\in X_c} g(x)\right] + h(c) \label{eq:9-minconst} \end{equation}

である($h(c)$ は $x$ に依らないので,どの $x$ を選んでも同じ値だけ足される).次の証明で本質的に効くのはこの性質である.

導出:シュレーディンガー変分問題の2段階分解

出発点は,量子力学の変分原理そのものである.すなわち基底状態エネルギーは,すべての規格化された反対称 $N$ 電子波動関数についてのハミルトニアン期待値の最小値である:

\begin{equation} E_0 = \min_{\Psi}\ \braket{\Psi|\hat{H}|\Psi} = \min_{\Psi}\ \braket{\Psi|\hat{T}+\hat{V}_{ee}+\hat{V}_{\mathrm{ext}}|\Psi} \label{eq:9-lv1} \end{equation}

ここで波動関数の集合を密度によってクラス分けする.どんな $\Psi$ も密度 $\rho_\Psi$ をただ1つ持つから,これは数学ノートの意味での分割である.クラスの添字は $N$ 表示可能な密度 $\rho \in \mathcal{J}_N$ 全体を走り,9.5節のハリマンの構成によりどのクラスも空でない.分解 \eqref{eq:9-minsplit} を適用すると

\begin{align} E_0 &= \min_{\rho \in \mathcal{J}_N}\ \left[\ \min_{\Psi\to\rho}\ \left\{ \braket{\Psi|\hat{T}+\hat{V}_{ee}|\Psi} + \braket{\Psi|\hat{V}_{\mathrm{ext}}|\Psi}\right\}\ \right] \label{eq:9-lv2}\\ &= \min_{\rho \in \mathcal{J}_N}\ \left[\ \min_{\Psi\to\rho}\ \left\{ \braket{\Psi|\hat{T}+\hat{V}_{ee}|\Psi} + \int \dd^3 r\, v(\bm{r})\rho(\bm{r})\right\}\ \right] \label{eq:9-lv3}\\ &= \min_{\rho \in \mathcal{J}_N}\ \left[\ \left\{\min_{\Psi\to\rho} \braket{\Psi|\hat{T}+\hat{V}_{ee}|\Psi}\right\} + \int \dd^3 r\, v(\bm{r})\rho(\bm{r})\ \right] \label{eq:9-lv4}\\ &= \min_{\rho \in \mathcal{J}_N}\ \left[\ F_L[\rho] + \int \dd^3 r\, v(\bm{r})\rho(\bm{r})\ \right] \label{eq:9-lv5} \end{align}

各ステップの根拠は次のとおり.

これで,$3N$ 変数の波動関数についての最小化 \eqref{eq:9-lv1} が,3変数の密度についての最小化 \eqref{eq:9-lv5} に厳密に書き換えられた.近似はどこにも入っていない.∎

$$ E_0 \;=\; \min_{\rho \in \mathcal{J}_N}\left[\,F_L[\rho] + \int \dd^3 r\, v(\bm{r})\rho(\bm{r})\,\right], \qquad F_L[\rho] = \min_{\Psi\to\rho}\braket{\Psi|\hat{T}+\hat{V}_{ee}|\Psi} $$
⟨Ψ|Ĥ|Ψ⟩ 密度 ρ ρ₁ ρ₂ ρ₀ ρ₃ ρ₄ E₀ ① 内側の最小化 min (Ψ→ρ) ⇒ FL[ρ] ② 外側の最小化 min(ρ) ● 灰色の点 = 密度が ρ に等しい波動関数 Ψ ● 緑の点 = そのクラスの最小 = FL[ρ] + ∫vρ 縦の点線1本1本が「同じ密度を持つ波動関数の集まり」(クラス)
図9.6 制約付き探索による2段階最小化.すべての $N$ 電子波動関数を密度で分類し(縦の点線が1つのクラス),まず各クラスの中でエネルギーを最小化して緑の点を得る.この値が $F_L[\rho]+\int v\rho$ である.次にクラスをまたいで最小化すると基底状態エネルギー $E_0$ と基底状態密度 $\rho_0$ が得られる.クラス内では $\int v\rho$ が定数なので,内側の最小化は $\hat{T}+\hat{V}_{ee}$ だけに効く——これが $F_L$ が普遍汎関数になる理由である.

9.6.3 レヴィの2つの定理

2段階分解の結果を,変分原理の形に整理し直しておく.

定理:レヴィの変分原理(1979)

外部ポテンシャル $v$ を持つ $N$ 電子系のエネルギー汎関数を,$N$ 表示可能な密度 $\rho \in \mathcal{J}_N$ に対して

$$E_v[\rho] \equiv F_L[\rho] + \int \dd^3 r\, v(\bm{r})\rho(\bm{r})$$

と定義する.このとき次の2つが成り立つ.

  1. (不等式) 任意の $N$ 表示可能な密度 $\rho$ について $E_v[\rho] \ge E_0$.
  2. (等号) $\rho_0$ を基底状態密度とすると $E_v[\rho_0] = E_0$.

導出:レヴィの2定理の証明

(1) の証明.$\rho \in \mathcal{J}_N$ を任意に取り,その制約付き探索の最小値を与える波動関数を $\Psi_\rho^{\min}$ とする.定義より $\Psi_\rho^{\min}$ の密度は $\rho$ である.

\begin{align} E_v[\rho] &= F_L[\rho] + \int \dd^3 r\, v\rho \label{eq:9-lvt1}\\ &= \braket{\Psi_\rho^{\min}|\hat{T}+\hat{V}_{ee}|\Psi_\rho^{\min}} + \int \dd^3 r\, v\rho \label{eq:9-lvt2}\\ &= \braket{\Psi_\rho^{\min}|\hat{T}+\hat{V}_{ee}|\Psi_\rho^{\min}} + \braket{\Psi_\rho^{\min}|\hat{V}_{\mathrm{ext}}|\Psi_\rho^{\min}} \label{eq:9-lvt3}\\ &= \braket{\Psi_\rho^{\min}|\hat{H}|\Psi_\rho^{\min}} \;\ge\; E_0 \label{eq:9-lvt4} \end{align}

\eqref{eq:9-lvt1}→\eqref{eq:9-lvt2} は $F_L$ の定義(最小値が $\Psi_\rho^{\min}$ で達成される).\eqref{eq:9-lvt2}→\eqref{eq:9-lvt3} は恒等式 \eqref{eq:9-vrho}($\Psi_\rho^{\min}$ の密度が $\rho$ であることを使う).\eqref{eq:9-lvt3}→\eqref{eq:9-lvt4} は期待値の線形性,最後の不等号はレイリー・リッツである.9.4.2節の証明と同じ形だが,$\Psi_\rho^{\min}$ はどこかのポテンシャルの基底状態である必要がまったくない点が決定的に違う.だから $\rho$ が $v$ 表示可能でなくても議論が通る.

(2) の証明.2つの向きの不等式を示す.

まず $\rho_0$ も $N$ 表示可能な密度なので,(1) より $E_v[\rho_0] \ge E_0$ である.

逆向きを示す.$\Psi_0$ を $\hat{H}$ の基底状態(縮退があればそのうちの1つ)とすると,$\Psi_0$ の密度は $\rho_0$ だから,$\Psi_0$ は $\rho_0$ に対する制約付き探索の候補の1つである.最小値は候補の値以下だから

\begin{equation} F_L[\rho_0] \le \braket{\Psi_0|\hat{T}+\hat{V}_{ee}|\Psi_0} = \braket{\Psi_0|\hat{H}|\Psi_0} - \braket{\Psi_0|\hat{V}_{\mathrm{ext}}|\Psi_0} = E_0 - \int \dd^3 r\, v\rho_0 \label{eq:9-lvt5} \end{equation}

(2つ目の等号でハミルトニアンから外部ポテンシャル項を引き,3つ目で $\braket{\Psi_0|\hat{H}|\Psi_0}=E_0$ と恒等式 \eqref{eq:9-vrho} を使った).両辺に $\int v\rho_0$ を足すと $E_v[\rho_0] \le E_0$ を得る.

2つを合わせて $E_v[\rho_0] = E_0$ である.∎

なお,この証明のどこにも「基底状態が非縮退である」という仮定は使われていない.縮退がある場合は,縮退した基底状態のどれを $\Psi_0$ に選んでも同じ議論が通り,対応する密度 $\rho_0$ のそれぞれで $E_v$ が最小値 $E_0$ を取る(最小値を与える $\rho$ が複数あるだけで,何も破綻しない).

9.6.4 2つの普遍汎関数の関係

$F_L$ と $F_{\mathrm{HK}}$ は別々に定義された.両者が矛盾しないこと——$F_{\mathrm{HK}}$ が定義できる範囲では両者が一致すること——を確かめておこう.

導出:$v$ 表示可能な密度の上では $F_L[\rho] = F_{\mathrm{HK}}[\rho]$

$\rho \in \mathcal{A}_N$ とし,対応する外部ポテンシャルを $v$,その非縮退基底状態を $\Psi[\rho]$,基底状態エネルギーを $E_0$ とする.$\Psi \to \rho$ を満たす任意の波動関数 $\Psi$ について,9.6.3節 (1) とまったく同じ計算をたどると

\begin{equation} \braket{\Psi|\hat{T}+\hat{V}_{ee}|\Psi} + \int \dd^3 r\, v\rho = \braket{\Psi|\hat{H}|\Psi} \;\ge\; E_0 \label{eq:9-FLHK1} \end{equation}

である.ここで $E_0 = \braket{\Psi[\rho]|\hat{H}|\Psi[\rho]} = F_{\mathrm{HK}}[\rho] + \int v\rho$(定義 \eqref{eq:9-FHK-def} と \eqref{eq:9-vrho})だから,\eqref{eq:9-FLHK1} の両辺から $\int v\rho$ を引いて

\begin{equation} \braket{\Psi|\hat{T}+\hat{V}_{ee}|\Psi} \;\ge\; F_{\mathrm{HK}}[\rho] \qquad (\Psi \to \rho \text{ を満たす任意の } \Psi) \label{eq:9-FLHK2} \end{equation}

つまり $F_{\mathrm{HK}}[\rho]$ は,制約付き探索の対象となるすべての候補の値の下界である.しかも $\Psi[\rho]$ 自身が候補の1つ($\Psi[\rho] \to \rho$)であり,そこで値がちょうど $F_{\mathrm{HK}}[\rho]$ になる.下界が実際に達成されるのだから,それが最小値である:

\begin{equation} F_L[\rho] = \min_{\Psi\to\rho}\braket{\Psi|\hat{T}+\hat{V}_{ee}|\Psi} = F_{\mathrm{HK}}[\rho], \qquad \rho \in \mathcal{A}_N \label{eq:9-FLHK3} \end{equation}

さらに副産物として,$\rho \in \mathcal{A}_N$ のときは制約付き探索の最小値を与える波動関数 $\Psi_\rho^{\min}$ が基底状態 $\Psi[\rho]$ そのものであることも分かる.∎

したがって $F_L$ は $F_{\mathrm{HK}}$ の拡張である.$\mathcal{A}_N$ の上では両者は同じ値を取り,$\mathcal{A}_N$ の外($N$ 表示可能だが $v$ 表示可能でない密度)では $F_L$ だけが値を持つ.以後,単に $F[\rho]$ と書けば $F_L[\rho]$ を指すことにする.

物理的意味:制約付き探索が解決したこと

レヴィの定式化が理論の土台をどう固めたかを整理する.

注意:制約付き探索は「計算法」ではない

$F_L[\rho]$ の定義式 \eqref{eq:9-FL-def} は,$F$ を計算する手続きを与えているように見える.実際には,これを額面どおり実行することは不可能に近い.密度が $\rho$ に等しいという制約を満たす $N$ 電子波動関数の集合は途方もなく大きく,その中を探索する計算量は電子数について指数関数的に増大する.$3N$ 次元の波動関数を扱う困難は,探索という形に姿を変えて残っているだけである.

制約付き探索の価値は定義としての明晰さにある.「$F[\rho]$ とは何か」に対して,$\rho$ さえ与えられれば原理的には一意に決まる数として答えを与え,その定義域を完全に確定させた.これによって,近似汎関数を作る人は「何を近似しているのか」を正確に述べられるようになった.実際,制約付き探索の形は近似汎関数が満たすべき厳密な条件(スケーリング則など)を導くときに繰り返し使われる道具になっている.

もう一つ技術的な注意として,\eqref{eq:9-FL-def} で最小値($\min$)と書いたが,厳密には下限($\inf$)が実際に達成されるかどうかは自明ではない.この点はリープ(1983)によって,凸解析の道具(ルジャンドル変換)を使って厳密に扱われ,適切な関数空間の上で最小値が存在することが示されている.

9.7 定理が保証すること・しないこと

本章で証明したことと証明していないことを,明確に切り分けておく.この切り分けを曖昧にしたまま先へ進むと,「DFTは厳密な理論なのに,なぜ計算結果に誤差があるのか」という初学者がほぼ必ず抱く疑問に自分で答えられなくなる.

9.7.1 保証されたこと

9.7.2 保証されなかったこと

物理的意味:「原理的に厳密,実用上は近似」という構図

この2つのリストを並べると,密度汎関数理論という営みの全体像が見える.理論の骨格——変分原理,定義域,汎関数の普遍性——は厳密である.近似が入るのは,ただ一点,$F[\rho]$ の表式だけである.これは他の多体理論の近似の入り方とはずいぶん様子が違う.例えば摂動論では展開の打ち切り次数が誤差を決め,次数を上げれば系統的に改善するが計算量が急増する.配置間相互作用(CI)では行列式の数を増やせば厳密解に収束するが,これも指数関数的コストである.DFTでは,計算コストは近似汎関数の質にほとんど依存せず,誤差はすべて $F[\rho]$ の近似という一箇所に集中している.

この構図には長所と短所がある.長所は,汎関数を改良すれば計算コストをほとんど増やさずに精度が上がることであり,実際1990年代以降の汎関数開発(GGA,ハイブリッド,非経験的制約に基づく構成)はこの利点を活かしてきた.短所は,系統的に厳密解に近づける手続きが存在しないことである.摂動論やCIには「次数を上げる」という収束の梯子があるが,DFTの近似汎関数の列には梯子がない.ある系で汎関数Aが汎関数Bより良い結果を出したとしても,それが別の系でも成り立つ保証はない.この非系統性こそがDFTの実務上の最大の弱点であり,逆に言えば「厳密な制約条件をできるだけ多く満たす汎関数を作る」という現代の汎関数開発の方法論が生まれた理由でもある.

9.7.3 第10章への戦略:運動エネルギーをどう扱うか

では,$F[\rho]$ をどう近似すればよいのか.ここで9.1.1節の表に戻ろう.アルゴン原子の運動エネルギーは,トーマス・フェルミ汎関数で評価すると約 37 Ha も過小評価された.化学結合のエネルギースケールが 0.1 Ha 程度であることを思えば,これは絶望的な誤差である.しかも運動エネルギーは $F[\rho]$ の中で最大の項である(ビリアル定理により $T = -E$,すなわち運動エネルギーの絶対値は全エネルギーと同程度).$F$ の大部分を占める項を,最も精度の悪い方法で近似している——これがTFD模型の失敗の構造だった.

コーン(W. Kohn)とシャム(L. J. Sham)の1965年の発想は,この構造を正面から崩す.運動エネルギーを密度の陽な汎関数で書くのをあきらめ,軌道を再導入して運動エネルギーの大部分を厳密に(第5章のハートリー・フォック法と同じ精度で)計算する.具体的には,$\rho$ と同じ密度を持つ相互作用しない電子系を仮想的に考え,その運動エネルギー

\begin{equation} T_s[\rho] \equiv \min_{\Phi \to \rho}\braket{\Phi|\hat{T}|\Phi} \qquad (\Phi \text{ はスレーター行列式}) \label{eq:9-Ts} \end{equation}

を主役に据える(これはレヴィの制約付き探索で $\hat{V}_{ee}$ を落とした形であり,本章の枠組みがそのまま使える点に注意).そして

\begin{equation} F[\rho] = \underbrace{T_s[\rho]}_{\text{軌道で厳密に計算}} + \underbrace{\frac{1}{2}\iint \dd^3 r\,\dd^3 r'\,\frac{\rho(\bm{r})\rho(\bm{r}')}{|\bm{r}-\bm{r}'|}}_{\text{ハートリー項:密度で厳密に書ける}} + \underbrace{E_{xc}[\rho]}_{\text{残り物:ここだけ近似}} \label{eq:9-KSsplit} \end{equation}

と分割する.第3項の交換相関エネルギー $E_{xc}[\rho]$ は,「$F$ から最初の2項を引いた残り」として定義される量であり,真の運動エネルギー $T$ と $T_s$ の差($T - T_s$)も含んでいる.重要なのは,この残り物が $F$ 全体に比べて小さいことである.アルゴン原子でいえば,$T_s$ が 525.95 Ha を担い,誤差は 1 Ha 以下に収まる(9.1.1節の表).大きくて難しい部分を厳密に扱い,小さい部分だけを近似する——単純だが決定的な戦略の転換である.

厳密な密度汎関数理論(HK 1964 / Levy 1979) E₀ = minρ [ F[ρ] + ∫ v ρ d³r ] ── F[ρ] は存在し,普遍で,一意 定義域は N 表示可能な密度の全体(凸集合) しかし定理は F[ρ] の存在を保証するだけで,その表式は与えない 戦略A:F 全体を密度で近似する TFD模型(第4章)・Xα法(9.1節) F ≈ Tth[ρ] + Ex[ρ] + EH[ρ] Ar 原子の運動エネルギー誤差 ≈ 37 Ha 化学結合のスケール(≈ 0.1 Ha)を 誤差が完全に飲み込む → 殻構造も結合も出ない 戦略B:運動エネルギーは軌道で扱う コーン・シャム法(第10章) F = Ts[ρ] + EH[ρ] + Exc[ρ] Ar 原子の運動エネルギー誤差 ≈ 0.9 Ha 近似は小さな残り物 Exc だけに閉じ込められる → 化学結合・殻構造が定量的に再現される 第10章:コーン・シャム方程式
図9.7 厳密な密度汎関数理論から実用的計算法への2つの分岐.定理は $F[\rho]$ の存在を保証するが表式を与えないので,どこかで近似が必要になる.$F$ 全体を密度の陽な汎関数で近似する戦略A(TFD,X$\alpha$)は運動エネルギーの誤差で破綻した.運動エネルギーの大部分を軌道で厳密に評価し,残った小さな交換相関項だけを近似する戦略Bが第10章のコーン・シャム法である.

ここで一つ注意しておく.式 \eqref{eq:9-Ts} で定義される $T_s[\rho]$ が任意の $N$ 表示可能な $\rho$ に対して存在するか,という問題は,$v$ 表示可能性と同じ性質の問い(非相互作用 $v$ 表示可能性)である.コーン・シャム法はこの点で「定理」ではなく「仮説(Ansatz)」の要素を含んでいる.詳しくは第10章で扱う.9.5節・9.6節の議論は,そこでもそのまま効いてくる.

9.8 まとめ

次章では,いま述べた戦略Bを実行に移す.$F[\rho]$ を $T_s[\rho] + E_H[\rho] + E_{xc}[\rho]$ と分割し,軌道について変分を取ると,有効ポテンシャル $v_{\mathrm{eff}}$ の中を動く1電子の方程式——コーン・シャム方程式——が得られる.この方程式は自己無撞着に解く必要があり,その反復構造(SCFループ)が第一原理計算プログラムの心臓部になる.本章で固めた土台の上に,いよいよ実際に動く計算体系が建つ.

9.9 演習問題

演習 9.1:相互作用が異なる系には第1定理が使えるか

第1定理の証明では,$\hat{T}$ と $\hat{V}_{ee}$ が2つの系で共通であることが本質的だった.いま,相互作用の強さだけが異なる2つのハミルトニアン $$\hat{H}_\lambda = \hat{T} + \lambda\hat{V}_{ee} + \sum_i v_\lambda(\bm{r}_i), \qquad \hat{H}_{\lambda'} = \hat{T} + \lambda'\hat{V}_{ee} + \sum_i v_{\lambda'}(\bm{r}_i) \qquad (\lambda \ne \lambda')$$ を考える.9.3節の証明を段階(i)・段階(ii)それぞれについてなぞり直し,どこで論理が破綻するかを具体的に指摘せよ.また,それでも「$\lambda$ と $\lambda'$ の基底状態密度が等しくなるように $v_\lambda, v_{\lambda'}$ を選ぶことができる」という主張(断熱接続の出発点になる仮定)が,第1定理と矛盾しない理由を説明せよ.

ヒント:段階(i)の引き算では $(\lambda-\lambda')\hat{V}_{ee}$ が残るので,\eqref{eq:9-sub2} が「掛け算の等式」にならない.段階(ii)では $\braket{\Psi'|\hat{V}_{ee}|\Psi'}$ が現れるが,この量は密度だけでは書けない(恒等式 \eqref{eq:9-vrho} が使えるのは1体演算子だけ).第1定理は「$\hat{T}+\hat{V}_{ee}$ を固定したうえでの $v \leftrightarrow \rho$ の対応」を述べているにすぎず,$\lambda$ が違えば別の対応表になる.

演習 9.2:ハリマン構成の具体例(ガウス型密度)

1次元で $N=2$,密度 $\displaystyle \rho(x) = \frac{2}{\sqrt{\pi}}e^{-x^2}$ とする.

  1. $\sigma(x) = \rho(x)/2$ と累積分布 $q(x)$ を求めよ(誤差関数 $\mathrm{erf}$ を使ってよい).
  2. 等密度軌道 $\phi_0(x)$ と $\phi_1(x)$ を明示的に書き下ろせ.
  3. $\int \phi_0^*\phi_1\,\dd x = 0$ を,$q$ への変数変換によって確かめよ.
  4. 軌道 $k$ の運動エネルギー $t_k$ を \eqref{eq:9-harr-kin} から計算し,$k=0,1$ を占有させたときの全運動エネルギーを数値で求めよ.

ヒント:(1) $q(x) = \tfrac{1}{2}\left[1+\mathrm{erf}(x)\right]$.(4) $\int_{-\infty}^{\infty}x^2 e^{-x^2}\dd x = \sqrt{\pi}/2$ と $\int_{-\infty}^{\infty}e^{-3x^2}\dd x = \sqrt{\pi/3}$ を使うと,$t_k = \tfrac{1}{4} + \tfrac{2\pi}{\sqrt{3}}k^2$ という閉じた形になる.第1項が $T_W[\rho]/N$ に一致していることも確かめよ.

演習 9.3:2電子閉殻系では $T_s = T_W$

非相互作用運動エネルギー汎関数を \eqref{eq:9-Ts} で定義する.$N=2$ でスピンが対になった閉殻系(2個の電子が同じ空間軌道を上下スピンで占める)を考えると,任意の $N$ 表示可能な密度 $\rho$ に対して $$T_s[\rho] = T_W[\rho] = \frac{1}{2}\int \dd^3 r\,\left|\nabla\rho^{1/2}\right|^2$$ が成り立つことを示せ.

ヒント:下からの評価は9.5.3節の数学ノート($T_W \le T$ は任意の波動関数に対して成り立つので,最小値に対しても成り立つ).上からの評価は,空間軌道を $\phi(\bm{r}) = \sqrt{\rho(\bm{r})/2}$(実数,節なし)と取ったスレーター行列式が確かに密度 $\rho$ を与えることを確認し,その運動エネルギーを計算する.両者が一致するので最小値が定まる.この結果は,なぜ水素分子やヘリウム原子のような2電子系がDFTのテストケースとして特別なのかを説明する.

演習 9.4:エネルギーの原点と化学ポテンシャル

外部ポテンシャルを定数だけシフトして $v(\bm{r}) \to v(\bm{r}) + c$ とする.

  1. ハミルトニアン,基底状態波動関数,基底状態エネルギー,基底状態密度がそれぞれどう変わるかを述べよ.
  2. エネルギー汎関数 $E_v[\rho]$ と,オイラー方程式に現れる化学ポテンシャル $\mu$ がどう変わるかを示せ.
  3. 「化学ポテンシャルの絶対値には意味がなく,意味があるのは差だけである」という主張を,この結果に基づいて論じよ.孤立系のイオン化エネルギーや電子親和力を議論するときに $c$ をどう固定するのが自然か.

ヒント:(1) $\hat{H} \to \hat{H} + Nc$(定数演算子を足すだけなので固有関数は不変).(2) $\int(v+c)\rho = \int v\rho + cN$ を使う.$F[\rho]$ は $v$ を含まないので不変であることに注意.(3) 通常は $v(\bm{r}) \to 0$($|\bm{r}|\to\infty$)と取る.

参考文献

  1. P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964). — 第1定理・第2定理の原論文.本章9.2〜9.4節の内容にあたる.$v$ 表示可能性の問題は付録で簡潔に触れられているのみである.
  2. M. Levy, Proc. Natl. Acad. Sci. USA 76, 6062 (1979). — 制約付き探索の原論文.9.6節の2段階最小化はここで導入された.
  3. M. Levy, Phys. Rev. A 26, 1200 (1982). — $v$ 表示可能でない密度の構成.9.5.2節の議論の原型.
  4. E. H. Lieb, Int. J. Quantum Chem. 24, 243 (1983). — 凸解析(ルジャンドル変換)による厳密な再定式化.制約付き探索の最小値の存在,アンサンブル $v$ 表示可能密度の稠密性など,本章で「厳密には」と留保した点がここで扱われている.
  5. J. E. Harriman, Phys. Rev. A 24, 680 (1981). — 等密度軌道の構成.9.5.4節・9.5.5節の内容.
  6. T. L. Gilbert, Phys. Rev. B 12, 2111 (1975). — 密度と密度行列の $N$ 表示可能性.
  7. J. C. Slater, Phys. Rev. 81, 385 (1951). — X$\alpha$ 法(9.1.2節・9.1.3節).
  8. E. Teller, Rev. Mod. Phys. 34, 627 (1962). — トーマス・フェルミ理論で分子が結合しないことの証明.
  9. T. Kato, Commun. Pure Appl. Math. 10, 151 (1957). — 核カスプ条件(9.3.3節の例).
  10. W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965). — 次章の主題.
  11. R. G. Parr and W. Yang, Density-Functional Theory of Atoms and Molecules, Oxford University Press (1989). — 本章の内容を化学の立場から詳述した標準的教科書.化学ポテンシャルと電気陰性度の関係も扱う.
  12. R. M. Dreizler and E. K. U. Gross, Density Functional Theory: An Approach to the Quantum Many-Body Problem, Springer (1990). — 表示可能性問題と厳密な定式化を丁寧に扱った教科書.
  13. W. Kohn, Rev. Mod. Phys. 71, 1253 (1999). — ノーベル賞受賞講演.定理の意義と限界が本人の言葉で整理されている.