マテリアル計算科学入門 — 目次 第II部 近似解を求める道具 / 第3章

第3章変分原理

第2章では,水素原子のSchrödinger方程式を極座標で変数分離し,3つの量子数 $n,\ell,m$ と厳密な波動関数を手に入れた.しかし,材料として水素原子ばかりを扱うわけにはいかない.実際に相手にするのは酸素であり,チタンであり,ランタンである.つまり多電子原子の集合体である.ならば多電子原子のSchrödinger方程式はどう解けるのか——そう問うた瞬間に,我々は最初の壁にぶつかる.電子が2個になっただけで,方程式はもう厳密には解けなくなるのである.

本章では,多電子原子のうち最も単純なヘリウム原子を題材に,(1) なぜ厳密に解けないのかを式で確認し,(2) それでも「尤(もっと)もらしい」エネルギーを得るための道具として変分原理を導入する.変分原理は「どんな試行関数を使ってエネルギー期待値を計算しても,その値は真の基底状態エネルギーより低くはならない」という,たった一行の命題である.しかしこの一行があるからこそ,Gaussian・GAMESS・VASP・Quantum ESPRESSOといった世界中の計算ソフトウェアは「エネルギーを下げれば下げるほど良い答えに近づく」という指針のもとで動くことができる.本章の到達点は,ヘリウムの有効核電荷 $Z'=Z-5/16=1.6875$ を自分の手で導き,全エネルギー $-77.49$ eV を実験値 $-79.005$ eV と比べられるようになることである.あわせて,以後の全章で使うブラケット記法を身につける.

この章で学ぶこと
  • ヘリウム原子のハミルトニアン(運動エネルギー2項+核-電子引力2項+電子間反発1項)の書き下し
  • 電子間反発項があると変数分離法が破綻すること,その破綻を式のどこで見抜くか
  • 電子間反発を無視した粗い見積り $E=2Z^2\epsilon_{1s}=-108.8$ eV と実験値 $-79.005$ eV の落差
  • ブラケット記法:$\langle a|b\rangle$,$|b\rangle\langle a|$,規格直交性 $\langle\psi_n|\psi_m\rangle=\delta_{nm}$,完全性条件 $\sum_n|\psi_n\rangle\langle\psi_n|=\mathbb{I}$
  • 有効核電荷 $Z'$ を変分パラメータとする試行関数 $\psi_T=\phi_{Z'}(\bm r_1)\phi_{Z'}(\bm r_2)$ と,期待値の5項分解
  • 平方完成による最小化:$Z'=Z-\dfrac{5}{16}=1.6875$,$\braket{E}=2\epsilon_{1s}\left(Z-\dfrac{5}{16}\right)^2=-77.49$ eV
  • 変分原理 $\braket{\Psi|\Ham|\Psi}\ge E_0$ の完全証明(固有関数の完全系展開を使う)
  • 試行関数=基底関数という発想,局在基底(STO-3G, 6-31G)と第5・6章への橋渡し
  • 寄り道:分散 $(\Delta p)^2=\braket{p^2}-\braket{p}^2$,そして「比熱はエネルギーのゆらぎである」$C_V=\dfrac{\braket{E^2}-\braket{E}^2}{k_{\mathrm B}T^2}$
前提:第2章(水素原子の厳密解,期待値と演算子,$\epsilon_{1s}=-13.606$ eV,Bohr半径 $a_0$).積分公式 $\int_0^\infty r^n e^{-ar}\dd r=n!/a^{n+1}$ と極座標の体積要素は付録A(数学の道具箱)にまとめてある.線形代数の内積・行列の積・固有値の知識も付録Aで復習できる.

本章でもSI単位系を用いる.電気素量を $e_0$,真空の誘電率を $\varepsilon_0$,電子の質量を $m$,Planck定数を $h$($\hbar=h/2\pi$)と書く.第2章で得た2つの基本量,

$$ \begin{equation} a_0=\frac{\varepsilon_0 h^2}{\pi m e_0^2}=0.529\times10^{-10}\ \mathrm{m}, \qquad \epsilon_{1s}=-\frac{m e_0^4}{8\varepsilon_0^2 h^2}=-13.606\ \mathrm{eV} \label{eq:3-constants} \end{equation} $$

を繰り返し使う.この2つは独立ではなく,次の関係で結ばれていることを最初に確認しておこう.これは本章で何度も登場する「換算の要」である.

$$ \begin{align} \frac{e_0^2}{4\pi\varepsilon_0 a_0} &= \frac{e_0^2}{4\pi\varepsilon_0}\cdot\frac{\pi m e_0^2}{\varepsilon_0 h^2} = \frac{m e_0^4}{4\varepsilon_0^2 h^2} = 2\times\frac{m e_0^4}{8\varepsilon_0^2 h^2} = -2\epsilon_{1s} \nonumber\\[4pt] &= 27.211\ \mathrm{eV}. \label{eq:3-hartree} \end{align} $$

数学ノート:本章で使う積分公式

本章の積分はすべて次の一本に帰着する(付録A参照).$a\gt 0$,$n$ を0以上の整数として

$$ \int_0^\infty r^n e^{-ar}\,\dd r=\frac{n!}{a^{n+1}} . $$

実際,$n=0$ なら $\int_0^\infty e^{-ar}\dd r=1/a$ であり,これを $a$ で微分するたびに左辺には $-r$ が,右辺には $-1/a^2,\ +2/a^3,\dots$ が現れることから帰納的に示せる.また球対称関数 $f(r)$ の3次元積分は,極座標の体積要素 $\dd^3r=r^2\sin\theta\,\dd r\,\dd\theta\,\dd\varphi$ を使って

$$ \int f(r)\,\dd^3 r=4\pi\int_0^\infty f(r)\,r^2\,\dd r $$

となる(角度部分の積分が $\int_0^\pi\sin\theta\,\dd\theta\int_0^{2\pi}\dd\varphi=2\cdot2\pi=4\pi$ を与えるため).

3.1 水素原子の次はヘリウム原子

なぜ?:どうして水素原子だけでは足りないのか

第2章で扱った水素原子は,電子が1個しかない特殊な系である.だからこそ厳密に解けた.しかし材料科学の現場で相手にするのは,酸化物なら酸素(8電子)とチタン(22電子),電池材料ならリチウム(3電子)とコバルト(27電子)である.水素イオン伝導体のように水素そのものが主役の材料もあるが,それでも周囲には必ず多電子原子がいる.電子が2個以上になると何が起きるのか——この問いを避けて計算科学は始まらない.そこでまず,多電子原子のうち最も単純な,電子2個のヘリウム原子から始める.

3.1.1 ヘリウム原子の座標系

原点に陽子が2個(電荷 $+2e_0$,一般には $+Ze_0$)ある原子核を置き,その周りに電子が2個ある系を考える.電子1の位置ベクトルを $\bm r_1$,電子2の位置ベクトルを $\bm r_2$ とし,それぞれの原点からの距離を $r_1=\abs{\bm r_1}$,$r_2=\abs{\bm r_2}$ と書く.2つの電子の間の距離は $r_{12}=\abs{\bm r_1-\bm r_2}$ である.

原子核 (+Z e₀, Z = 2) 原点 O 電子1 (−e₀) 電子2 (−e₀) r₁ r₂ r₁₂ = |r₁ − r₂| 引力:核 ↔ 各電子(実線) / 斥力:電子 ↔ 電子(破線)
図3.1 ヘリウム原子の座標系.原子核を原点に固定し,電子1と電子2の位置を $\bm r_1,\bm r_2$ で表す.核と電子の間には引力($-Ze_0^2/4\pi\varepsilon_0 r_i$)が,電子どうしの間には斥力($+e_0^2/4\pi\varepsilon_0 r_{12}$)がはたらく.

注意:原子核は動かないのか

ここでは原子核を原点に固定した.厳密には原子核も運動するが,陽子の質量は電子の約1836倍あるため,核の運動エネルギーは電子に比べて桁違いに小さい.この「核を止めてよい」という近似はBorn–Oppenheimer近似(Born–Oppenheimer approximation)と呼ばれ,第8章で正面から扱う.本章ではこの近似を暗黙に使っている.

3.1.2 ハミルトニアンを書き下す

Schrödinger方程式とは,粒子の運動エネルギーとポテンシャルエネルギーの合計(=ハミルトニアン $\Ham$)に関する固有値方程式であった.したがって,まずヘリウム原子のエネルギーを構成する項をすべて数え上げればよい.項は全部で5つある.

(1) 電子1の運動エネルギー.第2章と同様に,運動量演算子 $\hat{\bm p}_1=-i\hbar\nabla_1$ を使って

$$ \frac{\hat{\bm p}_1^{\,2}}{2m} =\frac{-\hbar^2}{2m}\left(\frac{\partial^2}{\partial x_1^2}+\frac{\partial^2}{\partial y_1^2}+\frac{\partial^2}{\partial z_1^2}\right) =\frac{-h^2}{8\pi^2 m}\nabla_1^2 . $$

ここで $\hbar=h/2\pi$ より $\hbar^2/2m=h^2/8\pi^2m$ を使った.$\nabla_1^2$ は電子1の座標だけに作用するLaplace演算子である.

(2) 電子2の運動エネルギー.まったく同じ形で,添字だけが2になる:$\dfrac{\hat{\bm p}_2^{\,2}}{2m}=\dfrac{-h^2}{8\pi^2 m}\nabla_2^2$.

(3) 電子1が原子核から受けるポテンシャルエネルギー.電荷 $+Ze_0$ と $-e_0$ の間のCoulombポテンシャルであるから,

$$ -\frac{Ze_0^2}{4\pi\varepsilon_0}\frac{1}{r_1} . $$

符号が負なのは引力(近づくほどエネルギーが下がる)だからである.

(4) 電子2が原子核から受けるポテンシャルエネルギー.同様に $-\dfrac{Ze_0^2}{4\pi\varepsilon_0}\dfrac{1}{r_2}$.

(5) 電子1と電子2の反発によるポテンシャルエネルギー.電荷 $-e_0$ どうしであるから積は $(-e_0)(-e_0)=+e_0^2$ となり,符号は正(斥力)である:

$$ +\frac{e_0^2}{4\pi\varepsilon_0}\frac{1}{\abs{\bm r_1-\bm r_2}} . $$

この5つを足し合わせたものがヘリウム原子のハミルトニアンである.

$$ \begin{equation} \Ham=\frac{\hat{\bm p}_1^{\,2}}{2m}-\frac{Ze_0^2}{4\pi\varepsilon_0}\frac{1}{r_1} +\frac{\hat{\bm p}_2^{\,2}}{2m}-\frac{Ze_0^2}{4\pi\varepsilon_0}\frac{1}{r_2} +\frac{e_0^2}{4\pi\varepsilon_0}\frac{1}{\abs{\bm r_1-\bm r_2}} . \label{eq:3-he-hamiltonian} \end{equation} $$

したがってヘリウム原子のSchrödinger方程式は

$$ \begin{equation} \left[\frac{\hat{\bm p}_1^{\,2}}{2m}-\frac{Ze_0^2}{4\pi\varepsilon_0 r_1} +\frac{\hat{\bm p}_2^{\,2}}{2m}-\frac{Ze_0^2}{4\pi\varepsilon_0 r_2} +\frac{e_0^2}{4\pi\varepsilon_0\abs{\bm r_1-\bm r_2}}\right]\psi(\bm r_1,\bm r_2)=E\,\psi(\bm r_1,\bm r_2) \label{eq:3-he-schrodinger} \end{equation} $$

と書ける.多電子波動関数は $\psi$,1電子波動関数は $\phi$ と書き分けるのが本書の約束である(仕様どおり,第7章以降でスピンを含む軌道 $\chi$ も登場する).式 \eqref{eq:3-he-schrodinger} の波動関数 $\psi(\bm r_1,\bm r_2)$ は,極座標で書けば $\psi(r_1,\theta_1,\varphi_1,r_2,\theta_2,\varphi_2)$ という6変数の関数である.水素原子の3変数から一気に倍になった.

物理的意味:最後の1項が問題児である

式 \eqref{eq:3-he-hamiltonian} の5項のうち,はじめの4項は「電子1だけに関わる量」と「電子2だけに関わる量」にきれいに分かれている.ところが第5項 $e_0^2/4\pi\varepsilon_0\abs{\bm r_1-\bm r_2}$ だけは,$\bm r_1$ と $\bm r_2$ の両方が同時に入っている.この「両方が同時に入っている」ことこそ,以下で見るように,厳密解を不可能にする張本人である.物理としては当たり前のことを言っているにすぎない——電子1がどこにいるかによって電子2の感じる力が変わる,つまり電子相関(electron correlation)があるのだ.

3.2 なぜ厳密に解けないのか

3.2.1 変数分離法を試してみる

偏微分方程式を解くときの最も標準的な手段は変数分離法(separation of variables)である.第2章で水素原子を解いたときも,$\psi(r,\theta,\varphi)=R(r)\Theta(\theta)\Phi(\varphi)$ と置いて3本の常微分方程式に分解した.同じ手を打ってみよう.すなわち,ヘリウムの6変数の波動関数を

$$ \begin{equation} \psi(r_1,\theta_1,\varphi_1,r_2,\theta_2,\varphi_2) =\phi_1(r_1,\theta_1,\varphi_1)\,\phi_2(r_2,\theta_2,\varphi_2) \equiv \phi_1\phi_2 \label{eq:3-separation-ansatz} \end{equation} $$

と,電子1だけの関数 $\phi_1$ と電子2だけの関数 $\phi_2$ の積に書けると仮定する.これを式 \eqref{eq:3-he-schrodinger} に代入すると

$$ \left[\frac{\hat{\bm p}_1^{\,2}}{2m}-\frac{Ze_0^2}{4\pi\varepsilon_0 r_1} +\frac{\hat{\bm p}_2^{\,2}}{2m}-\frac{Ze_0^2}{4\pi\varepsilon_0 r_2} +\frac{e_0^2}{4\pi\varepsilon_0\abs{\bm r_1-\bm r_2}}\right]\phi_1\phi_2=E\,\phi_1\phi_2 . $$

ここで重要な観察をする.$\hat{\bm p}_1^{\,2}$ も $1/r_1$ も,電子1の座標にしか作用しない.したがって $\phi_2$ はこれらの演算子から見れば「ただの定数」であり,前に出すことができる.同様に $\phi_1$ は電子2に関する演算子から見れば定数である.よって上式は

$$ \begin{equation} \left(\frac{\hat{\bm p}_1^{\,2}\phi_1}{2m}-\frac{Ze_0^2}{4\pi\varepsilon_0 r_1}\phi_1\right)\phi_2 +\left(\frac{\hat{\bm p}_2^{\,2}\phi_2}{2m}-\frac{Ze_0^2}{4\pi\varepsilon_0 r_2}\phi_2\right)\phi_1 +\frac{e_0^2}{4\pi\varepsilon_0\abs{\bm r_1-\bm r_2}}\phi_1\phi_2=E\,\phi_1\phi_2 \label{eq:3-separation-substituted} \end{equation} $$

と整理できる.次に,変数分離法の常套手段として両辺を $\phi_1\phi_2$ で割る.第1項では $\phi_2$ が約分され,第2項では $\phi_1$ が約分され,第3項は $\phi_1\phi_2$ がまるごと約分される.結果は

$$ \begin{equation} \underbrace{\left(\frac{\hat{\bm p}_1^{\,2}\phi_1}{2m}-\frac{Ze_0^2}{4\pi\varepsilon_0 r_1}\phi_1\right)\Big/\phi_1}_{\text{電子1の座標のみ}} +\underbrace{\left(\frac{\hat{\bm p}_2^{\,2}\phi_2}{2m}-\frac{Ze_0^2}{4\pi\varepsilon_0 r_2}\phi_2\right)\Big/\phi_2}_{\text{電子2の座標のみ}} +\underbrace{\frac{e_0^2}{4\pi\varepsilon_0}\frac{1}{\abs{\bm r_1-\bm r_2}}}_{\text{両方の座標}} =\underbrace{E}_{\text{定数}} \label{eq:3-separation-divided} \end{equation} $$

である.

3.2.2 破綻はどこで起きるか

変数分離法の論理を思い出そう.「$x$ だけの関数」+「$y$ だけの関数」=「定数」という形の等式が成り立つとき,$x$ を動かしても等式が保たれるためには第1項は $x$ に依らない定数でなければならず,同様に第2項も定数でなければならない.これが変数分離が成立する仕組みである.

式 \eqref{eq:3-separation-divided} を見ると,右辺は定数 $E$,第1項は電子1の座標だけの関数,第2項は電子2の座標だけの関数——ここまでは順調である.ところが第3項が電子1と電子2の座標の両方を含んでいる.この項は $\bm r_1$ を動かしても $\bm r_2$ を動かしても値が変わってしまうから,定数として括り出すことができない.したがって「第1項=定数,第2項=定数」という結論を導けず,変数分離法は成立しない.

第1項 r₁ だけの関数 (θ₁, φ₁ も含む) + 第2項 r₂ だけの関数 (θ₂, φ₂ も含む) + 第3項(電子間反発) e₀² / 4πε₀|r₁ − r₂| r₁ と r₂ の両方に依存 定数にできる ✓ 定数にできる ✓ 定数にできない ✗ 「定数 + 定数 + 非定数 = 定数」は矛盾 ⇒ ψ(r₁, r₂) = φ₁(r₁) φ₂(r₂) という積の形は Schrödinger 方程式を満たせない ⇒ 変数分離法による厳密解の構成は不可能
図3.2 変数分離の破綻.式 \eqref{eq:3-separation-divided} の3つの項のうち,電子間反発項だけが2つの電子の座標を同時に含むため定数として括り出せない.これがヘリウム原子が厳密に解けない理由である.

注意:「解けない」の意味

「厳密に解けない」というのは,$\psi=\phi_1\phi_2$ という積の形の解が存在しないという意味であり,より一般には初等関数の有限個の組み合わせで閉じた形の厳密解が書けないという意味である.数値的にいくらでも精度よく解くことは可能で,実際Pekerisらは変分法を極限まで押し進めてヘリウムの基底状態エネルギーを小数点以下10桁以上求めている(文献[6]).また,変数分離法で解けない偏微分方程式そのものを研究している数学者もいる.「解けない」=「何もわからない」ではないことに注意してほしい.

なぜ?:3体問題は古典力学でも解けない

実は「3つ以上の粒子が互いに力を及ぼし合う系は解析的に解けない」という事実は,量子力学に固有のものではない.古典力学でも,太陽・地球・月のような三体問題(three-body problem)には一般的な閉じた解が存在しないことが知られている(Poincaré, 1890年代).ヘリウム原子は「核+電子+電子」というまさに三体問題であり,解けないのはむしろ自然なことである.したがって我々が目指すべきは厳密解ではなく,制御された近似である.

3.3 電子間反発を無視するとどうなるか

3.3.1 まず一番粗い近似を試す

厳密に解けないとわかったら,次に考えるのは「大雑把でもよいから近似値を知りたい」ということである.最も乱暴な近似は,じゃまをしている第3項をゼロと置いてしまうことだ.すなわち式 \eqref{eq:3-separation-divided} で電子間反発項を落とすと

$$ \left(\frac{\hat{\bm p}_1^{\,2}\phi_1}{2m}-\frac{Ze_0^2}{4\pi\varepsilon_0 r_1}\phi_1\right)\Big/\phi_1 +\left(\frac{\hat{\bm p}_2^{\,2}\phi_2}{2m}-\frac{Ze_0^2}{4\pi\varepsilon_0 r_2}\phi_2\right)\Big/\phi_2=E $$

となり,これは「定数+定数=定数」の形なので変数分離が成立する.しかも各項は,核電荷 $Ze_0$ の原子核に電子が1個だけ束縛された系——つまり水素様原子(hydrogen-like atom)のSchrödinger方程式そのものである.第2章で得たとおり,水素様原子のエネルギー準位は

$$ \begin{equation} E_n^{(Z)}=Z^2\,\frac{\epsilon_{1s}}{n^2} =Z^2\left(-\frac{m e_0^4}{8\varepsilon_0^2 h^2}\right)\frac{1}{n^2} =Z^2\left(-\frac{e_0^2}{8\pi\varepsilon_0 a_0}\right)\frac{1}{n^2} \label{eq:3-hydrogenic} \end{equation} $$

である.$Z$ の2乗で効いてくることに注意してほしい.核電荷が大きいほど電子は強く引きつけられ,準位は深く沈む.

ヘリウムでは $Z=2$,1s軌道なので $n=1$ である.式 \eqref{eq:3-hydrogenic} より1電子あたり $Z^2\epsilon_{1s}=4\epsilon_{1s}$,電子が2個あるから

$$ \begin{equation} E=2Z^2\epsilon_{1s}=2\times4\times(-13.606\ \mathrm{eV})=8\times(-13.606\ \mathrm{eV})=-108.8\ \mathrm{eV} \label{eq:3-nonint-energy} \end{equation} $$

を得る.これがヘリウムの全エネルギーの,最も粗い見積りである.

3.3.2 実験値と比べる

ヘリウムの全エネルギーの実験値は $-79.005$ eV である.これは「ヘリウム原子から電子を2個とも無限遠へ引き離すのに必要なエネルギー」であり,第1イオン化エネルギー $24.587$ eV と第2イオン化エネルギー $54.418$ eV の和として測定される($24.587+54.418=79.005$ eV).

計算値 $-108.8$ eV と実験値 $-79.005$ eV の差は

$$ \abs{-108.8-(-79.005)}=29.84\ \mathrm{eV} $$

と,実に約 30 eV もある.化学結合1本のエネルギーが数 eV であることを思えば,これは桁違いに大きい誤差である.

注意:「27%の誤差」と「38%の誤差」

誤差の大きさを百分率で表すとき,何を分母に取るかで数字が変わる.計算値 $-108.8$ eV を基準にすれば $29.84/108.8=27.4\%$,実験値 $-79.005$ eV を基準にすれば $29.84/79.005=37.8\%$ である.講義資料では前者の「約27%」が使われている.いずれにせよ「数十パーセント」の誤差であり,材料の性質を議論するにはまったく使い物にならない.本書では以後,誤差は実験値を基準に取ることを原則とし,必要なら絶対値(eV)を併記する.

さらに悪いことに,計算値は実験値より低い(より安定側).3.7節で見るように,正しく計算された期待値は真の基底状態エネルギーより低くなり得ない.$-108.8$ eV が $-79.005$ eV を下回ったのは,我々がハミルトニアンの一部(電子間反発)を勝手に捨てた——つまり別の系を計算してしまった——からである.

物理的意味:電子間反発は決して無視できない

電子間反発を無視したときの誤差 30 eV は,ヘリウムの全エネルギーの約4割にあたる.つまり電子間反発(electron–electron repulsion)は「小さな補正」ではなく,系のエネルギーを決める主役級の項である.多電子系を扱う理論——変分法(本章),摂動論(第4章),Hartree–Fock法(第9章),密度汎関数理論——は,すべてこの項をどう扱うかをめぐる工夫だと言ってよい.

例題3.1 水素様イオンのイオン化エネルギー

(a) He$^+$(電子1個,$Z=2$)の第1イオン化エネルギーを式 \eqref{eq:3-hydrogenic} から求めよ.(b) その結果と式 \eqref{eq:3-nonint-energy} を使って,電子間反発を無視した近似のもとでヘリウム原子の第1イオン化エネルギーを予測し,実験値 $24.587$ eV と比べよ.

解答 (a) He$^+$ の1s準位は $E=Z^2\epsilon_{1s}=4\times(-13.606)=-54.42$ eV.イオン化エネルギーはこれを0にするのに必要なエネルギーだから $+54.42$ eV である.実験値 $54.418$ eV と3桁一致する.電子が1個の系では厳密解が使えるので当然である.

(b) 電子間反発を無視した近似では He の全エネルギーが $-108.8$ eV,He$^+$ の全エネルギーが $-54.42$ eV であるから,第1イオン化エネルギーは

$$-54.42-(-108.8)=54.42\ \mathrm{eV}$$

と予測される.実験値 $24.587$ eV の2倍以上であり,まったく合わない.原因は明らかで,He には電子が2個いるため電子間反発の分だけエネルギーが持ち上がっている(=結合が弱まっている)のに,それを無視したからである.電子が1個の系は厳密に解けるが,2個になった瞬間に破綻する——この落差が本章の出発点である.

なぜ?:それでも「解いた」ことにしたい

厳密解が書けないなら,どうすればよいか.ここで発想を切り替える.固有値と固有関数を厳密に求めるのをあきらめ,「尤もらしい」エネルギーの値だけを求めるのである.第2章で学んだように,ある波動関数 $\psi$ のもとでの物理量 $\hat A$ の期待値は $\braket{A}=\int\psi^*\hat A\psi\,\dd x$ で与えられた.$\hat A$ としてハミルトニアン $\Ham$ を取れば,それはエネルギーの期待値である.厳密な固有関数でなくても,期待値なら計算できる.あとは「どの $\psi$ を使うか」だけが問題になる.この方針が本章の主題であり,その正当性を保証するのが変分原理である.

3.4 ブラケット記法

3.4.1 積分記号だらけになる問題

第2章で学んだとおり,1次元の系におけるエネルギー期待値は

$$ \begin{equation} \braket{E}=\int_{-\infty}^{\infty}\psi^*\,\Ham\,\psi\,\dd x \label{eq:3-expectation-1d} \end{equation} $$

である.3次元関数なら積分記号が3つ必要で,

$$ \begin{equation} \braket{E}=\int_{-\infty}^{\infty}\!\int_{-\infty}^{\infty}\!\int_{-\infty}^{\infty}\psi^*\,\Ham\,\psi\,\dd x\,\dd y\,\dd z \label{eq:3-expectation-3d} \end{equation} $$

となる.ところがヘリウムの波動関数は6変数なので,積分記号がさらに3つ増えて6重積分になる.多電子系ではこれが $3N$ 重積分に膨れ上がる.毎回これを書いていては,式が積分記号だらけになって本質が見えなくなってしまう.

そこで物理学者が使う略記法がブラケット記法(bra–ket notation,Diracの記法)である.積分をまるごと1つの記号にたたみ込んでしまう.

定義:ブラケット記法

関数 $f(x)$,$g(x)$ に対して,ブラ(bra)$\bra{f}$ とケット(ket)$\ket{g}$ を導入し,その積を次の「関数の内積」と定める:

$$ \begin{equation} \langle f|g\rangle \equiv \int_{-\infty}^{\infty} f^*(x)\,g(x)\,\dd x . \label{eq:3-braket-def} \end{equation} $$

ブラの側には必ず複素共役 $*$ が付くことに注意する.演算子 $\hat A$ を挟んだ形も同様に

$$ \langle f|\hat A|g\rangle \equiv \int_{-\infty}^{\infty} f^*(x)\,\hat A\,g(x)\,\dd x $$

と定める.多次元の場合は右辺の積分が多重積分になるだけで,左辺の書き方は変わらない.すなわちエネルギー期待値は,次元によらず

$$ \braket{E}=\bra{\psi}\Ham\ket{\psi} $$

の一行で書ける.

「ブラケット(bracket,括弧)」を「ブラ(bra)」と「ケット(ket)」に割ったのがこの名前の由来である.Diracの命名によるものだが,なかなか洒落ている.

3.4.2 重なり積分が一行で書ける

ブラケット記法のありがたみを実感するために,第5章で主役になる量を先取りしよう.原子 A と原子 B が隣り合って並んでいるとき,それぞれの波動関数を $\psi_{\mathrm A},\psi_{\mathrm B}$ とすると,両者の重なり積分(overlap integral)$S$ は

$$ S=\iiint \psi_{\mathrm A}^*(\bm r)\,\psi_{\mathrm B}(\bm r)\,\dd x\,\dd y\,\dd z $$

と定義される.これがブラケット記法では

$$ \begin{equation} S=\langle\psi_{\mathrm A}|\psi_{\mathrm B}\rangle \label{eq:3-overlap-braket} \end{equation} $$

のたった一行になる.なんとシンプルなことか.$S$ は「2つの軌道がどれだけ空間的に重なっているか」を表す量で,化学結合の強さを決める最重要量の一つである.第5章では,これを1s軌道どうしについて楕円座標系で手計算する.

3.4.3 ベクトルの内積との対応

ブラケット記法は関数だけでなく,ふつうのベクトルにも使える.行ベクトルをブラ,列ベクトルをケットに対応させればよい:

$$ \vec a=[a_x,\ a_y,\ a_z]=\bra{a}, \qquad \vec b=\begin{bmatrix}b_x\\ b_y\\ b_z\end{bmatrix}=\ket{b} . $$

このとき内積は

$$ \vec a\cdot\vec b=\langle a|b\rangle=a_x b_x+a_y b_y+a_z b_z $$

と表される.ブラを左,ケットを右に並べて掛けると内積(スカラー)になる,というわけである.

では,逆にケットを左,ブラを右に並べたらどうなるか.列ベクトル×行ベクトルであるから,答えは行列である:

$$ \ket{b}\bra{a} =\begin{bmatrix}b_x\\ b_y\\ b_z\end{bmatrix}[a_x,\ a_y,\ a_z] =\begin{bmatrix} b_x a_x & b_x a_y & b_x a_z\\ b_y a_x & b_y a_y & b_y a_z\\ b_z a_x & b_z a_y & b_z a_z \end{bmatrix}. $$

「ブラとケットの順序を入れ替えるだけでスカラーが行列になる」——この事実は3.7節の完全性条件で決定的な役割を果たす.

なぜ?:どうして「掛け算」が「積分」になるのか

ベクトルの内積は $\sum_i a_i b_i$ という和である.一方,関数の内積は $\int f^*g\,\dd x$ という積分である.この2つが同じ「内積」と呼ばれるのは奇妙に見えるかもしれない.しかし,内積の定義には $\sum$ があることを思い出そう.関数 $f(x)$ を,$x$ を細かく刻んだ点 $x_1,x_2,\dots$ での値 $f(x_1),f(x_2),\dots$ を並べた「無限次元のベクトル」だと思えば,

$$ \sum_i f^*(x_i)g(x_i)\,\Delta x \ \xrightarrow[\Delta x\to0]{}\ \int f^*(x)g(x)\,\dd x $$

となり,和が積分に化ける.つまり関数とは無限次元のベクトルである.この見方(関数空間,Hilbert空間)に慣れると,量子力学の記法が一気に見通しよくなる.第11章の永年方程式で,波動関数の展開係数を並べたベクトルの固有値問題を解くことになるが,そこで効いてくるのがこの発想である.

3.4.4 Schrödinger方程式と規格直交性

ブラケット記法を使うと,Schrödinger方程式そのものが

$$ \Ham\ket{\psi}=E\ket{\psi} $$

と書ける.水素原子の主量子数 $n$ や,無限井戸型ポテンシャルの量子数 $n$ を添字として付ければ

$$ \begin{equation} \Ham\ket{\psi_n}=E_n\ket{\psi_n}\qquad(n=0,1,2,\dots) \label{eq:3-schrodinger-braket} \end{equation} $$

である.ここで,異なる $n$ の固有関数どうしを掛けて積分すると0になり(直交性),同じ $n$ の固有関数どうしなら1になる(規格化条件).これをまとめて

$$ \begin{equation} \langle\psi_n|\psi_m\rangle=\delta_{nm} =\begin{cases}1 & (n=m)\\ 0 & (n\neq m)\end{cases} \label{eq:3-orthonormal} \end{equation} $$

と書く.$\delta_{nm}$ はKroneckerのデルタ(Kronecker delta)である.この性質をもつ関数の組を正規直交系(orthonormal system)と呼ぶ.ハミルトニアンのようなエルミート演算子の固有関数は,必ず正規直交系に取れることが知られている.

注意:ベクトルと関数で流儀が違う

ブラ $\bra{a}$ とケット $\ket{b}$ がベクトルのときは $\langle a|b\rangle$ は単なる内積(成分の積の和)である.しかし $\bra{a}=f^*(x)$,$\ket{b}=g(x)$ のように関数のときは,$\langle a|b\rangle$ は積分 $\int f^*(x)g(x)\dd x$ を意味する.同じ記号が文脈によって和にも積分にもなるので,いま何を扱っているのかを常に意識してほしい.なお,第7章で登場するスピン関数 $\bm\alpha(\sigma),\bm\beta(\sigma)$ は2成分の「ベクトル」(スピノル)であり,その内積は和である.空間部分は積分,スピン部分は和,と使い分けることになる.

3.4.5 固有関数を「縦ベクトル」だと思う

正規直交系という言葉が出てきたので,より直感的な見方を紹介しておく.Schrödinger方程式の解 $\ket{\psi_0},\ket{\psi_1},\ket{\psi_2},\dots$ を,無限次元空間の互いに直交する単位ベクトルだと思うのである:

$$ \ket{\psi_0}=\begin{bmatrix}\psi_0\\ 0\\ 0\\ \vdots\end{bmatrix},\quad \ket{\psi_1}=\begin{bmatrix}0\\ \psi_1\\ 0\\ \vdots\end{bmatrix},\quad \ket{\psi_2}=\begin{bmatrix}0\\ 0\\ \psi_2\\ \vdots\end{bmatrix},\quad\dots $$

このとき $\ket{\psi_0}\bra{\psi_0}$ は,$(1,1)$ 成分だけが $\int\psi_0^*\psi_0\,\dd x=1$ で,それ以外の成分がすべて0の行列になる:

$$ \ket{\psi_0}\bra{\psi_0}= \begin{bmatrix} 1 & 0 & 0 & \cdots\\ 0 & 0 & 0 & \cdots\\ 0 & 0 & 0 & \cdots\\ \vdots & \vdots & \vdots & \ddots \end{bmatrix},\qquad \ket{\psi_1}\bra{\psi_1}= \begin{bmatrix} 0 & 0 & 0 & \cdots\\ 0 & 1 & 0 & \cdots\\ 0 & 0 & 0 & \cdots\\ \vdots & \vdots & \vdots & \ddots \end{bmatrix},\ \dots $$

したがって,すべての $n$ について足し合わせると単位行列になる.これを完全性条件(completeness relation)という:

$$ \begin{equation} \sum_{n=0}^{\infty}\ket{\psi_n}\bra{\psi_n} =\begin{bmatrix} 1 & 0 & \cdots\\ 0 & 1 & \cdots\\ \vdots & \vdots & \ddots \end{bmatrix} =\mathbb{I} . \label{eq:3-completeness} \end{equation} $$

「単位行列を掛けても何も変わらない」から,この $\mathbb{I}$ は式の好きなところにこっそり挿入してよい.3.7節の変分原理の証明では,まさにこの技を使う.

例題3.2 重ね合わせ状態のエネルギー期待値

正規直交系 $\{\ket{\psi_n}\}$ を用いて $\ket{\Psi}=\dfrac{1}{\sqrt5}\bigl(\ket{\psi_0}+2\ket{\psi_1}\bigr)$ という状態をつくった.(a) $\ket{\Psi}$ が規格化されていることを確かめよ.(b) この状態のエネルギー期待値を $E_0,E_1$ で表せ.(c) $E_0=-13.606$ eV,$E_1=-3.401$ eV(水素原子の $n=1,2$)のとき数値を求め,$E_0$ と比べよ.

解答 (a) ブラは $\bra{\Psi}=\dfrac{1}{\sqrt5}\bigl(\bra{\psi_0}+2\bra{\psi_1}\bigr)$ であるから(係数が実数なので複素共役はそのまま),

$$ \langle\Psi|\Psi\rangle=\frac{1}{5}\Bigl(\langle\psi_0|\psi_0\rangle+2\langle\psi_0|\psi_1\rangle+2\langle\psi_1|\psi_0\rangle+4\langle\psi_1|\psi_1\rangle\Bigr) =\frac{1}{5}(1+0+0+4)=1 . $$

(b) $\Ham\ket{\psi_n}=E_n\ket{\psi_n}$ と規格直交性 \eqref{eq:3-orthonormal} を使うと

$$ \bra{\Psi}\Ham\ket{\Psi} =\frac{1}{5}\Bigl(E_0\langle\psi_0|\psi_0\rangle+2E_1\langle\psi_0|\psi_1\rangle+2E_0\langle\psi_1|\psi_0\rangle+4E_1\langle\psi_1|\psi_1\rangle\Bigr) =\frac{E_0+4E_1}{5} . $$

(c) 数値を入れると $\dfrac{-13.606+4\times(-3.401)}{5}=\dfrac{-13.606-13.604}{5}=\dfrac{-27.210}{5}=-5.442$ eV.基底状態エネルギー $E_0=-13.606$ eV より高い.これは偶然ではなく,3.7節で証明する変分原理の帰結である.基底状態以外の成分(ここでは $\ket{\psi_1}$)が混ざると,期待値は必ず $E_0$ より持ち上がる.

3.5 試行関数と有効核電荷 — ヘリウムのエネルギー期待値

3.5.1 計算の方針

エネルギー期待値 $\braket{E}=\bra{\psi}\Ham\ket{\psi}$ を計算するには,$\psi$ を具体的に決めなければならない.厳密解は手に入らないのだから,こちらで「これらしい」関数を用意するしかない.そこで次の2段構えで進める.

  1. 水素原子の厳密解をまねた,試し(trial)の意味合いの波動関数——試行関数(trial function)——を用意し,エネルギー期待値の式に代入する.試行関数には調整可能なパラメータを1つ以上含めておく.
  2. 得られた期待値をパラメータの関数と見て,最小化するパラメータを求める.

「系のエネルギー期待値を最小にするという意味で尤もらしい波動関数を選べば,それっぽい答えに辿り着けるのではないか」——これが変分法の発想である.なぜこの発想が正当なのかは3.7節で証明する.

3.5.2 遮蔽の描像とパラメータの選び方

なぜ?:パラメータに何を選ぶべきか

試行関数のパラメータは何でもよいわけではない.物理的に意味のある量を選ぶと,少ないパラメータで大きな改善が得られる.ヘリウムの場合,鍵になるのは遮蔽(screening,shielding)である.電子1は原子核($+2e_0$)に引きつけられるが,その途中に電子2の負電荷の雲が広がっている.すると電子1から見た「実効的な核電荷」は $+2e_0$ より小さくなる.つまり電子は互いに相手をかばい,核の正電荷を薄めてしまう.ならば,核電荷 $Z$ そのものをパラメータ $Z'$ に置き換え,$Z'$ をいくつにすればエネルギーが最も低くなるかを調べればよい.この $Z'$ を有効核電荷(effective nuclear charge)と呼ぶ.

遮蔽を考えない場合 +2 電子 強い引力 感じる核電荷 = Z = 2 → E = 2Z²ε₁ₛ = −108.8 eV(低すぎる) もう1つの電子が遮蔽する場合 +2 もう1つの電子の雲 電子 弱い引力 感じる核電荷 = Z′ < Z → Z′ を変分パラメータにする 遮蔽定数 s ≡ Z − Z′ を求めることが,本節の目標
図3.3 遮蔽の描像.ヘリウムでは,一方の電子の負電荷が原子核の正電荷を部分的に打ち消し,もう一方の電子が感じる実効的な核電荷 $Z'$ は $Z=2$ より小さくなる.この $Z'$ を変分パラメータとして扱う.

3.5.3 試行関数を書き下す

電子間反発を考えない場合,ヘリウムの波動関数は水素様原子の1s軌道の積であった.そこで核電荷を $Z$ から $Z'$ に置き換えた1s軌道

$$ \phi_{Z'}(\bm r)=\frac{1}{\sqrt\pi}\left(\frac{Z'}{a_0}\right)^{3/2}\exp\!\left(-\frac{Z' r}{a_0}\right) $$

を用意し,これを2つ掛けたものを試行関数とする.

定義:ヘリウムの1パラメータ試行関数

$$ \begin{equation} \psi_T(\bm r_1,\bm r_2)=\phi_{Z'}(\bm r_1)\,\phi_{Z'}(\bm r_2) =\frac{1}{\pi}\left(\frac{Z'}{a_0}\right)^{3}\exp\!\left[-\frac{Z'}{a_0}(r_1+r_2)\right]. \label{eq:3-trial-function} \end{equation} $$

添字 $T$ は trial(試行)の意味である.$Z'$ は有効核電荷であり,この段階では値が未定の変分パラメータである.$Z'=Z$ と置けば,電子間反発を無視したときの波動関数に戻る.

数学ノート:試行関数の規格化を確かめる

$\phi_{Z'}$ が規格化されていることを確認しておこう.$\zeta\equiv Z'/a_0$ と置くと

$$ \begin{align} \int\abs{\phi_{Z'}(\bm r)}^2\dd^3 r &=\frac{1}{\pi}\zeta^3\cdot 4\pi\int_0^\infty e^{-2\zeta r}r^2\,\dd r \nonumber\\ &=4\zeta^3\cdot\frac{2!}{(2\zeta)^3} =4\zeta^3\cdot\frac{2}{8\zeta^3}=1 . \nonumber \end{align} $$

ここで積分公式 $\int_0^\infty r^ne^{-ar}\dd r=n!/a^{n+1}$($n=2$,$a=2\zeta$)を使った.$\psi_T$ は $\phi_{Z'}$ の積なので,当然 $\langle\psi_T|\psi_T\rangle=1$ である.規格化されていることは変分原理を適用する前提なので,必ず確認する習慣をつけてほしい.

3.5.4 ハミルトニアンを5つの項に分ける

ここが計算の要である.ハミルトニアン \eqref{eq:3-he-hamiltonian} をそのまま試行関数で挟むと,核-電子引力の項が「核電荷 $Z$」なのに試行関数は「核電荷 $Z'$」でできているため,水素様原子の既知の結果が直接使えない.そこで核-電子引力の項を $Z'$ の部分と残りに分割するという工夫をする.$-Z=-Z'+(Z'-Z)$ と書けることを使って

$$ -\frac{Ze_0^2}{4\pi\varepsilon_0}\frac{1}{r_1} =-\frac{Z'e_0^2}{4\pi\varepsilon_0}\frac{1}{r_1} +\frac{(Z'-Z)e_0^2}{4\pi\varepsilon_0}\frac{1}{r_1} $$

とする(電子2についても同様).これを式 \eqref{eq:3-he-hamiltonian} に代入すると,ハミルトニアンは次の5つの部分に整理される.

$$ \begin{align} \Ham=\ &\underbrace{\left[\frac{\hat{\bm p}_1^{\,2}}{2m}-\frac{Z'e_0^2}{4\pi\varepsilon_0}\frac{1}{r_1}\right]}_{\text{(i-a) 水素様}} +\underbrace{\left[\frac{\hat{\bm p}_2^{\,2}}{2m}-\frac{Z'e_0^2}{4\pi\varepsilon_0}\frac{1}{r_2}\right]}_{\text{(i-b) 水素様}} \nonumber\\[4pt] &+\underbrace{\frac{(Z'-Z)e_0^2}{4\pi\varepsilon_0}\frac{1}{r_1}}_{\text{(ii-a) 核電荷のズレ}} +\underbrace{\frac{(Z'-Z)e_0^2}{4\pi\varepsilon_0}\frac{1}{r_2}}_{\text{(ii-b) 核電荷のズレ}} +\underbrace{\frac{e_0^2}{4\pi\varepsilon_0}\frac{1}{\abs{\bm r_1-\bm r_2}}}_{\text{(iii) 電子間反発}} . \label{eq:3-hamiltonian-split} \end{align} $$

したがってエネルギー期待値は5つの期待値の和になる:

$$ \begin{align} \braket{E}=\bra{\psi_T}\Ham\ket{\psi_T} =\ &\bra{\phi_{Z'}(\bm r_1)}\frac{\hat{\bm p}_1^{\,2}}{2m}-\frac{Z'e_0^2}{4\pi\varepsilon_0 r_1}\ket{\phi_{Z'}(\bm r_1)} +\bra{\phi_{Z'}(\bm r_2)}\frac{\hat{\bm p}_2^{\,2}}{2m}-\frac{Z'e_0^2}{4\pi\varepsilon_0 r_2}\ket{\phi_{Z'}(\bm r_2)} \nonumber\\[4pt] &+\bra{\phi_{Z'}(\bm r_1)}\frac{(Z'-Z)e_0^2}{4\pi\varepsilon_0 r_1}\ket{\phi_{Z'}(\bm r_1)} +\bra{\phi_{Z'}(\bm r_2)}\frac{(Z'-Z)e_0^2}{4\pi\varepsilon_0 r_2}\ket{\phi_{Z'}(\bm r_2)} \nonumber\\[4pt] &+\bra{\psi_T(\bm r_1,\bm r_2)}\frac{e_0^2}{4\pi\varepsilon_0\abs{\bm r_1-\bm r_2}}\ket{\psi_T(\bm r_1,\bm r_2)} . \label{eq:3-five-terms} \end{align} $$

1〜4項目で試行関数が $\psi_T$ から $\phi_{Z'}$ 1個に減っているのは,たとえば第1項では電子2に関する積分が $\int\abs{\phi_{Z'}(\bm r_2)}^2\dd^3r_2=1$ となって落ちるからである.以下,5つの項を順に評価していく.

3.5.5 (i) 水素様項

第1項の演算子は,核電荷が $Z'e_0$ の原子核に電子1個が束縛された系のハミルトニアンそのものである.そして試行関数 $\phi_{Z'}$ は,まさにその系の1s固有関数である.したがって固有値方程式

$$ \left[\frac{\hat{\bm p}_1^{\,2}}{2m}-\frac{Z'e_0^2}{4\pi\varepsilon_0 r_1}\right]\phi_{Z'}(\bm r_1)=Z'^2\epsilon_{1s}\,\phi_{Z'}(\bm r_1) $$

が成り立ち,規格化条件 $\langle\phi_{Z'}|\phi_{Z'}\rangle=1$ と合わせて

$$ \begin{equation} \text{(i-a)}=\text{(i-b)}=Z'^2\epsilon_{1s} \qquad\Longrightarrow\qquad \text{(i-a)}+\text{(i-b)}=2Z'^2\epsilon_{1s} \label{eq:3-term12} \end{equation} $$

を得る.

導出:水素様項を運動エネルギーと引力に分けて検算する

式 \eqref{eq:3-term12} を,運動エネルギーと引力の期待値を別々に計算して確かめてみよう.まず運動エネルギーは,1s軌道(指数部 $e^{-\zeta r}$,$\zeta=Z'/a_0$)に対して

$$ \bra{\phi_{Z'}}\frac{\hat{\bm p}^{\,2}}{2m}\ket{\phi_{Z'}}=\frac{\hbar^2\zeta^2}{2m}=\frac{h^2}{8\pi^2m}\left(\frac{Z'}{a_0}\right)^2 $$

である.ここに $a_0=\varepsilon_0h^2/\pi me_0^2$ を代入すると

$$ \frac{h^2Z'^2}{8\pi^2m}\cdot\frac{\pi^2m^2e_0^4}{\varepsilon_0^2h^4} =Z'^2\frac{me_0^4}{8\varepsilon_0^2h^2}=-Z'^2\epsilon_{1s} $$

となり,正の値(運動エネルギーだから当然)である.次に引力の期待値には,後で使う重要な結果

$$ \begin{equation} \bra{\phi_{Z'}}\frac{1}{r}\ket{\phi_{Z'}}=\frac{Z'}{a_0} \label{eq:3-mean-inv-r} \end{equation} $$

を用いる(次の数学ノートで導出する).すると

$$ \bra{\phi_{Z'}}-\frac{Z'e_0^2}{4\pi\varepsilon_0 r}\ket{\phi_{Z'}} =-Z'\frac{e_0^2}{4\pi\varepsilon_0}\cdot\frac{Z'}{a_0} =-Z'^2\cdot\frac{e_0^2}{4\pi\varepsilon_0a_0} \overset{\eqref{eq:3-hartree}}{=}-Z'^2(-2\epsilon_{1s})=2Z'^2\epsilon_{1s} . $$

両者の和は $-Z'^2\epsilon_{1s}+2Z'^2\epsilon_{1s}=Z'^2\epsilon_{1s}$ となり,確かに式 \eqref{eq:3-term12} と一致する.ついでに $\braket{T}=-\tfrac12\braket{V}$,すなわち $\braket{V}=-2\braket{T}$ という関係(ビリアル定理,virial theorem)が成り立っていることも読み取れる.

数学ノート:$\braket{1/r}=Z'/a_0$ の計算

$\zeta=Z'/a_0$ と置く.$\phi_{Z'}=\pi^{-1/2}\zeta^{3/2}e^{-\zeta r}$ は球対称なので,極座標の体積要素を使って

$$ \begin{align} \bra{\phi_{Z'}}\frac{1}{r}\ket{\phi_{Z'}} &=\int \abs{\phi_{Z'}(\bm r)}^2\frac{1}{r}\dd^3r =\frac{\zeta^3}{\pi}\cdot4\pi\int_0^\infty e^{-2\zeta r}\frac{1}{r}\,r^2\,\dd r \nonumber\\ &=4\zeta^3\int_0^\infty r\,e^{-2\zeta r}\,\dd r =4\zeta^3\cdot\frac{1!}{(2\zeta)^2} =4\zeta^3\cdot\frac{1}{4\zeta^2}=\zeta=\frac{Z'}{a_0} . \nonumber \end{align} $$

$r^2$(体積要素)と $1/r$ が打ち消し合って $r^1$ の積分になるところがポイントである.結果は「$\braket{1/r}$ は軌道の広がり $a_0/Z'$ の逆数」という,きわめて自然な形をしている.なお $\braket{r}$ は $\tfrac32 a_0/Z'$ であって $1/\braket{1/r}$ とは一致しない(第2章参照).期待値の逆数と逆数の期待値は別物である.

3.5.6 (ii) 核電荷のズレの項

第3項・第4項は,式 \eqref{eq:3-mean-inv-r} を使えばただちに計算できる.

$$ \begin{align} \text{(ii-a)}&=\frac{(Z'-Z)e_0^2}{4\pi\varepsilon_0}\bra{\phi_{Z'}}\frac{1}{r_1}\ket{\phi_{Z'}} =\frac{(Z'-Z)e_0^2}{4\pi\varepsilon_0}\cdot\frac{Z'}{a_0} \nonumber\\ &=(Z'-Z)Z'\cdot\frac{e_0^2}{4\pi\varepsilon_0 a_0} \overset{\eqref{eq:3-hartree}}{=}(Z'-Z)Z'\cdot(-2\epsilon_{1s}) =2(Z-Z')Z'\epsilon_{1s} . \label{eq:3-term34} \end{align} $$

電子2についても同じなので,2つ合わせて

$$ \text{(ii-a)}+\text{(ii-b)}=4(Z-Z')Z'\epsilon_{1s}=-4(Z'-Z)Z'\epsilon_{1s} . $$

符号を確認しておこう.遮蔽が効いていれば $Z'\lt Z$ なので $Z-Z'\gt0$,また $\epsilon_{1s}\lt0$ であるから,この項は負である.物理的には「試行関数が本来の核電荷 $Z$ より弱い引力しか感じていない分を,あとから足し戻している」ことに対応する.

3.5.7 (iii) 電子間反発の項

残るは第5項,すなわち電子間反発の期待値

$$ \text{(iii)}=\frac{e_0^2}{4\pi\varepsilon_0}\iint\frac{\abs{\phi_{Z'}(\bm r_1)}^2\abs{\phi_{Z'}(\bm r_2)}^2}{\abs{\bm r_1-\bm r_2}}\,\dd^3r_1\,\dd^3r_2 $$

である.これは6重積分であり,$\abs{\bm r_1-\bm r_2}$ という2つの座標が絡んだ分母を扱うために,Legendre多項式による展開(あるいは古典電磁気学のGaussの法則を使う方法)が必要になる.この積分の完全な計算は第4章4.7節で行うので,ここでは結果だけを引用する:

$$ \begin{equation} \text{(iii)}=\frac{5}{8}\,\frac{Z'e_0^2}{4\pi\varepsilon_0 a_0} \overset{\eqref{eq:3-hartree}}{=}\frac{5}{8}Z'(-2\epsilon_{1s}) =-\frac{5}{4}Z'\epsilon_{1s} . \label{eq:3-term5} \end{equation} $$

$\epsilon_{1s}\lt0$ なので $-\tfrac54 Z'\epsilon_{1s}\gt0$,すなわちこの項は正である.反発なのだからエネルギーを持ち上げる方向に働くはずで,符号は正しい.

注意:$5/8$ という係数はどこから来るのか

もし電子間の反発を「2つの点電荷が距離 $\braket{r}$ だけ離れている」と素朴に見積もると,係数は $5/8=0.625$ にはならない.実際には電子は雲のように広がっており,$1/\abs{\bm r_1-\bm r_2}$ を電子密度で二重に平均する必要がある.その結果として出てくるのが $5/8$ である.この係数は水素様1s軌道どうしのCoulomb積分(Coulomb integral)の値であり,第4章4.7節で導出し,第9章では交換積分と対にして再登場する.いま覚えるべきは,電子間反発の期待値が $Z'$ の1次に比例するという点である.この「1次」が,次節の平方完成で決定的な役割を果たす.

3.5.8 5つの項をまとめる

以上をすべて足し合わせると,エネルギー期待値が有効核電荷 $Z'$ の関数として得られる.

$$ \begin{equation} \braket{E}(Z')=2Z'^2\epsilon_{1s}-4(Z'-Z)Z'\epsilon_{1s}-\frac{5}{4}Z'\epsilon_{1s} . \label{eq:3-energy-sum} \end{equation} $$

$\epsilon_{1s}$ でくくって整理すると

$$ \begin{align} \braket{E}(Z') &=\left[2Z'^2-4(Z'-Z)Z'-\frac{5}{4}Z'\right]\epsilon_{1s} \nonumber\\ &=\left[2Z'^2-4Z'^2+4ZZ'-\frac{5}{4}Z'\right]\epsilon_{1s} \nonumber\\ &=\left(-2Z'^2+4ZZ'-\frac{5}{4}Z'\right)\epsilon_{1s} \label{eq:3-energy-quadratic} \end{align} $$

となる.$Z'$ についての2次関数である.$\epsilon_{1s}\lt0$ で $Z'^2$ の係数は $-2\epsilon_{1s}\gt0$ だから,この放物線は下に凸である.すなわち最小値をもつ.

例題3.3 各項の数値を確かめる

$Z=2$,$\epsilon_{1s}=-13.606$ eV とする.(a) $Z'=2$(遮蔽なし)のときの3種類の項の数値と合計を求めよ.(b) $Z'=1.6875$ のときも同様に求めよ.

解答 (a) $Z'=Z=2$ のとき:

合計は $-108.85+0+34.02=-74.83$ eV.これは3.3節の $-108.8$ eV に電子間反発を足し戻しただけの値であり,第4章で学ぶ1次の摂動論の結果 $\left(2Z^2-\tfrac54Z\right)\epsilon_{1s}=(8-2.5)\times(-13.606)=-74.83$ eV と完全に一致する.

(b) $Z'=1.6875$ のとき:$Z'^2=2.84766$,$Z'-Z=-0.3125$ に注意して

合計は $-77.49-28.70+28.70=-77.49$ eV.(a) より 2.66 eV 低く,実験値 $-79.005$ eV に近づいた.興味深いことに,最適な $Z'$ では (ii) と (iii) が正確に打ち消し合っている.これは偶然ではなく,次節で見る $Z'=Z-5/16$ という条件そのものである.実際 $Z'-Z=-5/16$ を代入すれば $\text{(ii)}=-4\times(-\tfrac{5}{16})Z'\epsilon_{1s}=\tfrac54Z'\epsilon_{1s}$ となり,$\text{(iii)}=-\tfrac54Z'\epsilon_{1s}$ と厳密に相殺する.

3.6 平方完成による最小化

3.6.1 平方完成する

式 \eqref{eq:3-energy-quadratic} の最小値を求める.高校数学以来おなじみの平方完成を使おう.まず $Z'^2$ の係数が $+1$ になるように $-2\epsilon_{1s}$ でくくる:

$$ \begin{align} \braket{E}(Z') &=\left(-2Z'^2+4ZZ'-\frac{5}{4}Z'\right)\epsilon_{1s} \nonumber\\ &=-2\epsilon_{1s}\left(Z'^2-2ZZ'+\frac{5}{8}Z'\right) . \label{eq:3-factored} \end{align} $$

($4ZZ'\epsilon_{1s}=-2\epsilon_{1s}\times(-2ZZ')$,$-\tfrac54Z'\epsilon_{1s}=-2\epsilon_{1s}\times\tfrac58Z'$ である.分数を間違えないよう手を動かして確かめてほしい.)次に $Z'$ の1次の項をまとめる.$-2ZZ'+\tfrac58Z'=-2\left(Z-\tfrac{5}{16}\right)Z'$ であるから

$$ \braket{E}(Z')=-2\epsilon_{1s}\left[Z'^2-2\left(Z-\frac{5}{16}\right)Z'\right] . $$

ここで $\tfrac58$ が $2\times\tfrac{5}{16}$ に化けたのが,$5/16$ という有名な数の出どころである.あとは $\left(Z-\tfrac{5}{16}\right)^2$ を足して引けば平方完成が完了する:

$$ \begin{equation} \braket{E}(Z')=-2\epsilon_{1s}\left[Z'-\left(Z-\frac{5}{16}\right)\right]^2 +2\epsilon_{1s}\left(Z-\frac{5}{16}\right)^2 . \label{eq:3-complete-square} \end{equation} $$

3.6.2 最小値を読み取る

注意:$\epsilon_{1s}\lt0$ だから $-2\epsilon_{1s}\gt0$

平方完成した式から最小値を読むとき,2乗の項の係数の符号を必ず確認しなければならない.ここでは係数は $-2\epsilon_{1s}$ であり,$\epsilon_{1s}=-13.606$ eV は負だから $-2\epsilon_{1s}=+27.212$ eV は正である.したがって式 \eqref{eq:3-complete-square} の第1項は常に0以上であり,$\braket{E}$ は第1項がゼロになるときに最小となる.もし符号を取り違えて $-2\epsilon_{1s}\lt0$ だと思ってしまうと,同じ式から「最大値」を読み取ってしまう.学部生が最もよく間違えるポイントである.

第1項がゼロになるのは $Z'=Z-\dfrac{5}{16}$ のときである.したがって

$$ \begin{equation} Z'_{\min}=Z-\frac{5}{16}, \qquad \braket{E}_{\min}=2\epsilon_{1s}\left(Z-\frac{5}{16}\right)^2 . \label{eq:3-zprime-opt} \end{equation} $$

ヘリウムの陽子数 $Z=2$ を代入すると,$5/16=0.3125$ より

$$ Z'=2-\frac{5}{16}=2-0.3125=1.6875 $$

を得る.そして「尤もらしい」エネルギー期待値は

$$ \begin{align} \braket{E}_{\min} &=2\times(-13.606\ \mathrm{eV})\times(1.6875)^2 \nonumber\\ &=-27.212\times2.84766\ \mathrm{eV} \nonumber\\ &=-77.49\ \mathrm{eV} \label{eq:3-emin} \end{align} $$

である.

導出:微分でも同じ答えが出ることの確認

平方完成に自信がないときは,微分して確かめればよい.式 \eqref{eq:3-energy-quadratic} を $Z'$ で微分すると

$$ \frac{\dd\braket{E}}{\dd Z'}=\left(-4Z'+4Z-\frac{5}{4}\right)\epsilon_{1s} . $$

これがゼロになるのは($\epsilon_{1s}\neq0$ だから)

$$ -4Z'+4Z-\frac{5}{4}=0 \quad\Longrightarrow\quad Z'=Z-\frac{5}{16} $$

のときで,平方完成の結果 \eqref{eq:3-zprime-opt} と一致する.さらに2階微分は $\dfrac{\dd^2\braket{E}}{\dd Z'^2}=-4\epsilon_{1s}=+54.42\ \mathrm{eV}\gt0$ なので,この停留点は確かに極小である.$Z=2$ を数値で入れると $Z'=2-1.25/4=2-0.3125=1.6875$.

3.6.3 放物線を描いてみる

$Z=2$,$\epsilon_{1s}=-13.606$ eV を式 \eqref{eq:3-energy-quadratic} に代入すると

$$ \begin{equation} \braket{E}(Z')=-13.606\left(-2Z'^2+6.75\,Z'\right)\ \mathrm{eV} =\left(27.212\,Z'^2-91.841\,Z'\right)\ \mathrm{eV} \label{eq:3-parabola-numeric} \end{equation} $$

となる($4Z-\tfrac54=8-1.25=6.75$ を使った).この式に具体的な $Z'$ を入れて表にしてみよう.

表3.1 有効核電荷 $Z'$ を変えたときのヘリウムのエネルギー期待値(式 \eqref{eq:3-parabola-numeric}).$Z'=1.6875$ で最小値をとる.
$Z'$1.001.201.401.601.68751.802.002.202.50
$\braket{E}$ (eV)$-64.63$$-71.02$$-75.24$$-77.28$$\mathbf{-77.49}$$-77.15$$-74.83$$-70.34$$-59.53$
有効核電荷 Z′ に対するヘリウムのエネルギー期待値 ⟨E⟩ のグラフ.下に凸の放物線で,Z′ = 1.6875 で最小値 −77.49 eV をとり,曲線は実験値 −79.005 eV の破線より常に上にある.Z′ = 2(遮蔽なし)の点は −74.83 eV.
図3.4 式 \eqref{eq:3-parabola-numeric} が与える放物線.$Z'$ を 1.0 から 2.5 まで動かすと $\braket{E}$ は下に凸の放物線を描き,$Z'=1.6875$ で最小値 $-77.49$ eV をとる.破線は実験値 $-79.005$ eV.変分原理により,この曲線はどこでも実験値(真の基底状態エネルギー)より上にある.

3.6.4 結果を吟味する

得られた結果を整理しよう.

表3.2 ヘリウム原子の全エネルギーの各種近似と実験値の比較.誤差は実験値 $-79.005$ eV を基準とした.
方法波動関数$E$ (eV)誤差 (eV)誤差 (%)
電子間反発を無視(3.3節)$Z'=Z=2$,反発項を捨てる$-108.85$$-29.84$37.8
1次摂動論(第4章)$Z'=Z=2$,反発項は期待値で評価$-74.83$$+4.17$5.3
変分法(本章)$Z'=1.6875$$-77.49$$+1.52$1.9
実験値—$-79.005$——
ヘリウム原子の全エネルギーの比較.1次摂動論 −74.83 eV,変分法 −77.49 eV,実験値 −79.005 eV(破線),電子間反発を無視した −108.85 eV の準位を並べ,変分法で到達できる領域(実験値より上)を緑の帯で示した図.
図3.5 ヘリウム原子の全エネルギー.変分法($-77.49$ eV)は1次摂動論($-74.83$ eV)より真の値に近い.緑の帯が変分法で到達可能な領域であり,実験値の線を下回ることはできない.

物理的意味:遮蔽定数 $5/16$ の読み方

$Z'=Z-5/16$ は「相手の電子が核電荷を $5/16=0.3125$ 個分だけ遮蔽している」と読める.この $s\equiv Z-Z'$ を遮蔽定数(shielding constant)という.興味深いのは,この値が $Z$ に依らないことである.He でも Li$^+$ でも Be$^{2+}$ でも,1s電子どうしの遮蔽定数は $5/16$ である.無機化学で使う経験則Slater則(Slater's rules)では,同じ1s殻の相手電子の遮蔽定数を $0.30$ と取るが,これは変分法から出た $0.3125$ を丸めた値にほかならない.教科書の経験則の裏に,こうした具体的な計算があることを知っておいてほしい.

なぜ?:たった1つのパラメータでここまで来られるのか

誤差は $37.8\%\to1.9\%$ に減った.パラメータはたった1つ,しかも「核電荷を少し小さくする」というだけの単純な工夫である.なぜこれほど効くのか.理由は,遮蔽が電子相関の中で最も平均的で系統的な成分だからである.「電子1がどこにいるかによって電子2が逃げる」という細かい相関(左右相関,角度相関)はまだ全く入っていないが,「相手がいるおかげで核の引力が全体として弱まる」という平均的効果は $Z'$ 一つで表現できる.この考え方——電子どうしの相互作用を平均場に押し込める——を徹底的に推し進めたものが,第9章のHartree–Fock法であり,密度汎関数理論である.逆に,残った 1.5 eV を削るには,$\psi_T$ に $r_{12}$ を陽に含めるなど,質的に違う工夫が必要になる(文献[5][6]).

実物を動かす:変分原理シミュレーター

ここまでの計算は,そのまま画面の上で動かせる.下のシミュレーターのタブ①で $Z'$ のスライダーを動かすと,放物線 $\braket{E}(Z')$ の上を点が走り,最小値 $Z'=Z-5/16$ と,3項(水素様・補正・反発)の内訳が同時に表示される.H$^-$ から B$^{3+}$ までの等電子系列も切り替えられるので,表3.2の数値を自分の手で再現してみてほしい.曲線上には $Z'=Z$ の点にも印が付いている——この値が次章で導く1次摂動のエネルギーとぴたり一致することは,第4章を読んでから確かめると面白い.タブ②は試行関数の「型」を変える実験(Gauss型 vs 指数型)で,3.8節と第6章の主題の予告編である.

別ウィンドウで大きく開く →  シミュレーターの解説ページ →

3.7 変分原理とその証明

3.7.1 主張を述べる

前節では「エネルギー期待値を最小にするパラメータを選ぶ」という操作を,いわば当然のように行った.しかし考えてみれば,これは自明ではない.エネルギーを下げれば下げるほど良い近似になる,という保証はどこにあるのか.その保証を与えるのが変分原理である.

定理:変分原理(variational principle)

ハミルトニアン $\Ham$ の固有値を小さい順に $E_0\le E_1\le E_2\le\cdots$ とし,$E_0$ を基底状態エネルギーとする.このとき,境界条件を満たす任意の規格化された波動関数 $\ket{\Psi}$($\langle\Psi|\Psi\rangle=1$)について

$$ \begin{equation} \bra{\Psi}\Ham\ket{\Psi}\ \ge\ E_0 \label{eq:3-variational} \end{equation} $$

が成り立つ.等号は $\ket{\Psi}$ が基底状態そのもの(縮退があればその線形結合)のときに限る.規格化されていない場合は,左辺を $\dfrac{\bra{\Psi}\Ham\ket{\Psi}}{\langle\Psi|\Psi\rangle}$ に置き換えればよい.

言葉で言えば,「どんな波動関数を持ってきてエネルギー期待値を計算しても,その値は真の基底状態のエネルギーより低くならない」ということである.したがって,試行関数のパラメータを動かしてエネルギーを下げる操作は,必ず真の答えに近づく方向の操作である.前節でヘリウムの $Z'$ を調整してよかった理屈は,これである.

波動関数の空間 境界条件を満たす関数すべて 試行関数の族 ψT(Z′) Z′=1.0 Z′=2.5 Z′=1.6875 (族の中の最良点) 真の基底状態 ψ₀ 届かない差 エネルギー E₀(真の基底状態) この領域には決して入れない ⟨E⟩(Z′=1.0) ⟨E⟩(Z′=1.4) ⟨E⟩(Z′=2.0) ⟨E⟩(Z′=1.6875) 最小化
図3.6 変分原理の概念図.試行関数の族(左図の曲線)は波動関数の空間のごく一部しか覆っていないため,真の基底状態 $\psi_0$ にぴったり届くとは限らない.しかしどの試行関数を使っても,そのエネルギー期待値は必ず $E_0$ 以上(右図)である.したがってエネルギーを下げる操作は必ず改善である.

3.7.2 証明の準備:エルミート性と完全性

証明に必要な材料は2つだけである.

材料1:エルミート性.ハミルトニアンはエルミート演算子(Hermitian operator)である.エルミート演算子とは,その固有値が必ず実数になる演算子のことである.エネルギーは観測できる量(実数)なのだから,$\Ham$ がエルミートであることは物理的要請である.これと式 \eqref{eq:3-schrodinger-braket},\eqref{eq:3-orthonormal} を組み合わせると

$$ \begin{equation} \bra{\psi_n}\Ham\ket{\psi_m}=E_m\langle\psi_n|\psi_m\rangle=E_m\delta_{nm} \label{eq:3-hamiltonian-matrix} \end{equation} $$

が成り立つ.すなわち,固有関数を基底に取ればハミルトニアンは対角行列になり,対角成分が固有値 $E_n$ である.

材料2:完全性.$\Ham$ の固有関数の組 $\{\ket{\psi_n}\}$ は完全系(complete set)をなす.これは「同じ境界条件を満たす任意の関数を,固有関数の線形結合で表せる」という意味である.したがって,境界条件を満たす任意の波動関数 $\ket{\Psi}$ は

$$ \begin{equation} \ket{\Psi}=\sum_{n=0}^{\infty}\ket{\psi_n}c_n \label{eq:3-expansion} \end{equation} $$

と展開できる.これは波の重ね合わせの原理と同じ形である.展開係数 $c_n$ は,$\ket{\Psi}$ に $\bra{\psi_n}$ を作用させて射影した結果として求まる.実際,式 \eqref{eq:3-expansion} の両辺に左から $\bra{\psi_m}$ を掛けると

$$ \langle\psi_m|\Psi\rangle=\sum_{n=0}^{\infty}\langle\psi_m|\psi_n\rangle c_n =\sum_{n=0}^{\infty}\delta_{mn}c_n=c_m $$

となる.つまり

$$ \begin{equation} c_n=\langle\psi_n|\Psi\rangle \qquad\Longrightarrow\qquad \ket{\Psi}=\sum_{n=0}^{\infty}\ket{\psi_n}\langle\psi_n|\Psi\rangle =\left(\sum_{n=0}^{\infty}\ket{\psi_n}\bra{\psi_n}\right)\ket{\Psi} . \label{eq:3-cn} \end{equation} $$

最後の等式は,まさに完全性条件 \eqref{eq:3-completeness}($\sum_n\ket{\psi_n}\bra{\psi_n}=\mathbb{I}$)が「単位行列を挿入してよい」と言っていることの現れである.

なぜ?:完全系という言葉のイメージ

3次元空間の任意のベクトルが $\vec v=v_x\hat e_x+v_y\hat e_y+v_z\hat e_z$ と3本の単位ベクトルで表せるのと同じことを,無限次元でやっているだけである.$\hat e_x,\hat e_y,\hat e_z$ が「空間を張る基底」であるように,$\ket{\psi_0},\ket{\psi_1},\dots$ は「関数空間を張る基底関数(basis function)」である.成分 $v_x=\hat e_x\cdot\vec v$ を取る操作が,係数 $c_n=\langle\psi_n|\Psi\rangle$ を取る操作に対応している.3.8節で登場する STO-3G や 6-31G といった基底関数も,まったく同じ発想の産物である.

3.7.3 証明

定理の証明:変分原理

ステップ1(ブラの展開).式 \eqref{eq:3-cn} の複素共役を取れば,$\ket{\Psi}$ の共役であるブラは

$$ \bra{\Psi}=\sum_{n=0}^{\infty}c_n^*\bra{\psi_n}=\sum_{n=0}^{\infty}\langle\Psi|\psi_n\rangle\bra{\psi_n} $$

となる.ここで $c_n^*=\langle\psi_n|\Psi\rangle^*=\langle\Psi|\psi_n\rangle$ を使った(ブラケットの複素共役は左右を入れ替えたものに等しい).

ステップ2(エネルギー期待値の展開).ブラとケットをそれぞれ展開して代入する.ケット側の添字を $m$,ブラ側の添字を $n$ とすると

$$ \begin{align} \bra{\Psi}\Ham\ket{\Psi} &=\sum_{m=0}^{\infty}\sum_{n=0}^{\infty}\langle\Psi|\psi_n\rangle\,\bra{\psi_n}\Ham\ket{\psi_m}\,\langle\psi_m|\Psi\rangle \nonumber\\ &\overset{\eqref{eq:3-hamiltonian-matrix}}{=}\sum_{m=0}^{\infty}\sum_{n=0}^{\infty}\langle\Psi|\psi_n\rangle\,E_m\delta_{nm}\,\langle\psi_m|\Psi\rangle \nonumber \end{align} $$

となる.$\delta_{nm}$ は $n=m$ のときだけ1なので二重和が一重和につぶれ,

$$ \begin{equation} \bra{\Psi}\Ham\ket{\Psi}=\sum_{m=0}^{\infty}E_m\,\langle\Psi|\psi_m\rangle\langle\psi_m|\Psi\rangle =\sum_{m=0}^{\infty}E_m\,\abs{\langle\psi_m|\Psi\rangle}^2 \label{eq:3-expect-sum} \end{equation} $$

を得る.ここで $\langle\Psi|\psi_m\rangle\langle\psi_m|\Psi\rangle=\langle\psi_m|\Psi\rangle^*\langle\psi_m|\Psi\rangle=\abs{\langle\psi_m|\Psi\rangle}^2$ を使った.

ステップ3(規格化条件).まったく同様に,$\Ham$ の代わりに恒等演算子を挟むと

$$ \begin{align} \langle\Psi|\Psi\rangle &=\sum_{m=0}^{\infty}\sum_{n=0}^{\infty}\langle\Psi|\psi_n\rangle\langle\psi_n|\psi_m\rangle\langle\psi_m|\Psi\rangle =\sum_{m=0}^{\infty}\sum_{n=0}^{\infty}\langle\Psi|\psi_n\rangle\,\delta_{nm}\,\langle\psi_m|\Psi\rangle \nonumber\\ &=\sum_{m=0}^{\infty}\abs{\langle\psi_m|\Psi\rangle}^2=1 . \label{eq:3-norm-sum} \end{align} $$

つまり $\abs{\langle\psi_m|\Psi\rangle}^2$ たちは,足すと1になる非負の数の組——確率分布である.

ステップ4(不等式).ここで $E_0\le E_1\le E_2\le\cdots$ を思い出す.すべての $m$ について $E_m\ge E_0$ であり,係数 $\abs{\langle\psi_m|\Psi\rangle}^2$ は非負だから,項ごとに $E_m\abs{\langle\psi_m|\Psi\rangle}^2\ge E_0\abs{\langle\psi_m|\Psi\rangle}^2$ が成り立つ.和を取れば

$$ \begin{equation} \bra{\Psi}\Ham\ket{\Psi}=\sum_{m=0}^{\infty}E_m\abs{\langle\psi_m|\Psi\rangle}^2 \ \ge\ \sum_{m=0}^{\infty}E_0\abs{\langle\psi_m|\Psi\rangle}^2 =E_0\underbrace{\sum_{m=0}^{\infty}\abs{\langle\psi_m|\Psi\rangle}^2}_{=1\ \text{(ステップ3)}}=E_0 . \label{eq:3-proof-final} \end{equation} $$

すなわち $\bra{\Psi}\Ham\ket{\Psi}\ge E_0$.Q.E.D.

証明を振り返ると,使ったのは「固有関数で展開する」「規格直交性で二重和を潰す」「$E_m\ge E_0$ を項ごとに使う」という3つだけである.こうして示してみれば,当たり前のことのように思える.実際,式 \eqref{eq:3-expect-sum} は「エネルギー期待値は,各固有値 $E_m$ をその出現確率 $\abs{\langle\psi_m|\Psi\rangle}^2$ で重みづけした平均である」と言っているにすぎない.平均が最小値を下回らないのは当然である.

物理的意味:期待値は「サイコロの平均」と同じ

サイコロの出目の期待値は $(1+2+\cdots+6)/6=3.5$ であり,最小の目である1を下回ることは絶対にない.どんなにいかさまサイコロを作っても,出目の平均が1を下回ることはできない.式 \eqref{eq:3-expect-sum} はまったく同じ構造をしている.「観測される可能性のあるエネルギー値」の集合が $\{E_0,E_1,E_2,\dots\}$ で,それぞれの出現確率が $\abs{\langle\psi_m|\Psi\rangle}^2$ である.平均が最小値 $E_0$ を下回れないのは,確率論の当然の帰結なのである.等号が成り立つのは $\abs{\langle\psi_0|\Psi\rangle}^2=1$,すなわち「必ず $E_0$ が出る」=「$\ket{\Psi}$ が基底状態そのもの」のときだけである.

注意:$-108.8$ eV は変分原理の反例ではない

3.3節で得た $-108.8$ eV は実験値 $-79.005$ eV より低い.変分原理に反しているように見えるが,そうではない.あのとき我々は電子間反発を含まない別のハミルトニアンの固有値を求めたのであって,ヘリウムの真のハミルトニアン \eqref{eq:3-he-hamiltonian} の期待値を計算したわけではない.変分原理が保証するのは「同じハミルトニアンの期待値を,いろいろな波動関数で計算した結果」についてである.ハミルトニアンを勝手にいじれば,いくらでも低い値が出てしまう.計算値が実験値を下回ったときは,まず「自分は本当に正しいハミルトニアンを使っているか」を疑うべきである.

例題3.4 規格化されていない試行関数の場合

試行関数 $\ket{\Phi}$ が規格化されていない($\langle\Phi|\Phi\rangle\neq1$)とき,

$$ \frac{\bra{\Phi}\Ham\ket{\Phi}}{\langle\Phi|\Phi\rangle}\ \ge\ E_0 $$

が成り立つことを示せ.

解答 $N\equiv\sqrt{\langle\Phi|\Phi\rangle}\gt0$ と置き,$\ket{\Psi}\equiv\ket{\Phi}/N$ を作る.すると

$$ \langle\Psi|\Psi\rangle=\frac{\langle\Phi|\Phi\rangle}{N^2}=\frac{N^2}{N^2}=1 $$

となり,$\ket{\Psi}$ は規格化されている.したがって定理 \eqref{eq:3-variational} が使えて

$$ E_0\le\bra{\Psi}\Ham\ket{\Psi}=\frac{1}{N^2}\bra{\Phi}\Ham\ket{\Phi}=\frac{\bra{\Phi}\Ham\ket{\Phi}}{\langle\Phi|\Phi\rangle} . $$

これで示された.分母 $\langle\Phi|\Phi\rangle$ をRayleigh商(Rayleigh quotient)の分母と呼ぶ.実際の計算プログラムでは,試行関数を毎回規格化するより,この形(分母付き)で扱うほうが便利なことが多い.第11章の永年方程式は,まさにこのRayleigh商を係数について最小化することから導かれる.

3.8 試行関数の選び方

3.8.1 試行関数は「基底関数」で作る

変分原理が保証するのは「エネルギーを下げれば改善する」ということだけである.どこまで下げられるかは,用意した試行関数の柔軟さで決まる.前節の図3.6でいえば,試行関数の族(曲線)が真の基底状態 $\psi_0$ の近くを通っていなければ,いくら最小化しても $E_0$ には届かない.

そこで実際の計算では,あらかじめ用意した関数の組 $\{\chi_1,\chi_2,\dots,\chi_K\}$ の線形結合

$$ \begin{equation} \phi=\sum_{p=1}^{K}c_p\,\chi_p \label{eq:3-basis-expansion} \end{equation} $$

を試行関数とし,係数 $c_p$ をすべて変分パラメータとして最小化する.この $\chi_p$ を基底関数(basis function),その集合を基底関数系(basis set)と呼ぶ.3.7.2節で見たとおり,無限個の完全系を使えば厳密解に到達できるが,計算機で扱えるのは有限個だけである.したがって「少ない個数で,いかに真の波動関数をよく張るか」が基底関数系の設計思想になる.

なぜ?:どうして試行関数の「型」にこだわるのか

同じ分子を計算しても,どんな試行関数を使うかによって得られるエネルギーも波動関数も変わってくる.だから計算科学の研究者の中には,試行関数の「型」に強いこだわりを持つ人が少なくない.基底関数を増やせばエネルギーは必ず下がる(変分原理)が,計算コストは基底関数の数 $K$ の4乗程度で増大する.精度とコストのトレードオフ——これが実際の計算科学における最大の判断ポイントであり,ソフトウェアのマニュアルに何十種類もの基底関数系が並んでいる理由である.

3.8.2 局在基底:STO-3G と 6-31G

分子の量子化学計算では,水素原子のSchrödinger方程式の厳密解($e^{-\zeta r}$ 型)や,それをまねた関数を基底関数として使うことが多い.関数が各原子の位置に「へばりついている」ので,これを局在基底(localized basis)と呼ぶ.固体を扱う場合でも局在基底を使う手法はある(第14章).

代表的なものを挙げる.

STO-3G の具体的な構成——3本のGaussianの指数と係数をどう決めるか——は第5章で導出する.また,Gauss基底で水素原子を実際に解いて厳密解と比べる作業は第6章で行う.

水素原子 1s 軌道の動径関数の比較.厳密な Slater 型の曲線は核の位置で尖り(カスプ),STO-3G は中間領域でよく一致し,STO-1G は核の位置で尖らず遠方で減衰が速すぎる.
図3.7 基底関数の違いによる波動関数の差.水素原子1s軌道の厳密解 $\phi=\pi^{-1/2}e^{-r/a_0}$(実線)に対し,Gauss関数1本(STO-1G,点線)はカスプも裾も再現できない.3本重ねた STO-3G(破線)は,中間領域でかなりよく一致する.基底関数を増やすほど厳密解に近づき,変分原理によりエネルギーも単調に下がる.

物理的意味:基底関数を変えると何が変わるか

Szabo–Ostlundの教科書には,水素分子 H$_2$ のポテンシャル曲線を STO-3G で計算した場合と 6-31G$^{**}$ で計算した場合を並べた有名な図がある(文献[2]).基底を良くすると,結合距離での曲線の深さ(結合エネルギー)が正確な解(Kolos–Wolniewiczの値)に近づく.しかし,原子間距離を無限に離したときの振る舞い——制限つきHartree–Fock(RHF)が解離極限を大きく誤る一方,非制限Hartree–Fock(UHF)はほぼ正しくなる——は基底関数を良くしても直らない.基底関数系の質と,波動関数の型(RHF/UHF)の選択は,別々の問題である.この区別は第6章と第9章で改めて扱う.

3.8.3 ほとんどすべての計算ソフトウェアが変分原理に立脚している

世界中で使われている量子化学計算ソフトウェア——分子計算用の Gaussian,GAMESS,PySCF,密度汎関数理論に基づく固体向け第一原理計算ソフトウェアの VASP,Quantum ESPRESSO,CASTEP,OpenMX など——は,ほとんどすべてがこの変分原理に従って「尤もらしい」エネルギーと波動関数を求めている.

具体的には,次の流れである.

  1. 計算したい原子配置を決める.
  2. 各原子に基底関数を割り当てる(局在基底なら STO-3G,6-31G,cc-pVDZ など.固体では平面波を使うことも多い).
  3. 展開係数 $c_p$ を変分パラメータとして,全エネルギーが最小になるように反復的に解く(自己無撞着場,self-consistent field, SCF).
  4. 得られたエネルギー・波動関数から,電子密度・状態密度・force などの物理量を評価する.

本章のヘリウム計算は,パラメータがたった1個($Z'$),基底関数がたった1種類という最小構成版だが,論理構造はまったく同じである.ヘリウムを手で解けるようになることは,VASP を使えるようになるための最短経路である.

数学ノート:変分原理が使えない場合

変分原理は基底状態にしか使えないことに注意してほしい.励起状態のエネルギーを下げようとしても,試行関数がうっかり基底状態成分を含んでいれば,いくらでも下がってしまう.励起状態を変分的に扱うには「基底状態と直交する」という拘束条件を課す必要がある(それでも第1励起状態までが限界である).また,密度汎関数理論では交換相関汎関数が近似であるため,得られた全エネルギーが真の値より低くなることがある.「変分的である」と主張できるのは,厳密なハミルトニアンの期待値を厳密に評価している場合に限られる.本書では第9章と第10章でこの点に立ち返る.

例題3.5 基底関数を増やすとエネルギーは必ず下がるか

基底関数系 $B_1=\{\chi_1,\dots,\chi_K\}$ を使って変分計算して得たエネルギーを $E^{(1)}$,$B_1$ に新しい関数 $\chi_{K+1}$ を1本加えた $B_2$ で得たエネルギーを $E^{(2)}$ とする.$E^{(2)}\le E^{(1)}$ が常に成り立つことを説明せよ.

解答 $B_2$ で作れる試行関数の集合は,$B_1$ で作れる試行関数の集合を含んでいる.実際,$B_2$ の展開 $\phi=\sum_{p=1}^{K+1}c_p\chi_p$ において $c_{K+1}=0$ と置けば $B_1$ の任意の試行関数が再現される.最小化は「より広い集合の中での最小値」を取る操作だから,集合が広がれば最小値は下がるか,変わらないかのどちらかである.よって $E^{(2)}\le E^{(1)}$.

この単調性のおかげで,基底関数を系統的に増やしていったときのエネルギーの下がり方を見れば,「基底関数系の収束」を判定できる.実際の計算では,基底を増やしてもエネルギーがもう下がらなくなった時点を「完全基底極限に十分近い」とみなす.ただし,変分原理が保証するのはエネルギーの単調性だけであり,双極子モーメントやスピン密度といった他の物理量が単調に改善する保証はないことに注意してほしい.

3.9 寄り道:期待値があるなら分散もある

なぜ?:期待値の「隣」にあるもの

本章では一貫してエネルギーの期待値を計算してきた.しかし確率密度 $\abs{\psi}^2$ から期待値が出せるのなら,統計学でいう分散や標準偏差も当然計算できるはずである.実際それらには量子力学的な深い意味があり,さらに熱統計力学まで足を伸ばすと,比熱という身近な物理量がエネルギーのゆらぎそのものであるという美しい結果に行き着く.本節はやや寄り道だが,「期待値」という道具がどこまで広がるかを見ておく価値は大きい.

3.9.1 分散と標準偏差の定義

確率論では,$N$ 回の観測で得られた観測量 $x_i$($i=1,\dots,N$)の分散(variance)$\sigma^2$ と標準偏差(standard deviation)$\sigma$ を

$$ \begin{equation} \sigma^2=\frac{1}{N}\sum_{i=1}^{N}\bigl(x_i-\braket{x}\bigr)^2, \qquad \sigma=\sqrt{\frac{1}{N}\sum_{i=1}^{N}\bigl(x_i-\braket{x}\bigr)^2} \label{eq:3-variance-def} \end{equation} $$

と定義する.$\braket{x}$ は期待値(平均)である.「観測値が期待値からどれだけ散らばっているか」を測る量である.

3.9.2 量子力学における分散

量子力学では,観測量は演算子で表され,期待値は波動関数で挟んだ積分で与えられる.したがって,たとえば運動量 $\hat p$ の分散は,式 \eqref{eq:3-variance-def} の $(x_i-\braket{x})^2$ を演算子 $(\hat p-\braket{p})^2$ に置き換えて

$$ (\Delta p)^2=\int\psi^*\bigl(\hat p-\braket{p}\bigr)^2\psi\,\dd x $$

と定義される.$\braket{p}$ はただの数(演算子ではない)であることに注意して,2乗を展開しよう:

$$ \bigl(\hat p-\braket{p}\bigr)^2=\hat p^2-2\braket{p}\hat p+\braket{p}^2 . $$

これを代入すると3つの積分に分かれる:

$$ \begin{align} (\Delta p)^2 &=\int\psi^*\hat p^2\psi\,\dd x -2\braket{p}\int\psi^*\hat p\,\psi\,\dd x +\braket{p}^2\underbrace{\int\psi^*\psi\,\dd x}_{=1\ (\text{規格化条件})} \nonumber\\ &=\braket{p^2}-2\braket{p}\braket{p}+\braket{p}^2 \nonumber\\ &=\braket{p^2}-\braket{p}^2 . \label{eq:3-dp2} \end{align} $$

したがって運動量の標準偏差(=不確定性,uncertainty)は

$$ \begin{equation} \Delta p=\sqrt{\braket{p^2}-\braket{p}^2} \label{eq:3-delta-p} \end{equation} $$

である.位置でもエネルギーでも,任意の物理量について同じ形が成り立つ.「2乗の期待値」から「期待値の2乗」を引く——この形は本節の後半でも,そしてこの先の章でも繰り返し現れる.

例題3.6 水素原子1s軌道の運動量の不確定性

水素原子の1s軌道について,(a) $\braket{p}$ を求めよ.(b) $\braket{p^2}$ を運動エネルギーの期待値から求めよ.(c) $\Delta p$ を数値で求め,$\Delta x\simeq a_0$ と見積もったときの $\Delta x\,\Delta p$ を評価せよ.

解答 (a) 1s軌道は実関数であり,束縛された定常状態である.定常状態では電子は「どこかへ流れていく」ことがないので,運動量の期待値はゼロである:$\braket{p}=0$.(電子は動いていないという意味ではない.右向きの成分と左向きの成分が等量あって打ち消し合っている,という意味である.)

(b) 3.5.5節の導出より,水素($Z'=1$)の1s軌道では

$$ \frac{\braket{p^2}}{2m}=-\epsilon_{1s}=13.606\ \mathrm{eV} \quad\Longrightarrow\quad \braket{p^2}=2m\times13.606\ \mathrm{eV} . $$

(c) $13.606\ \mathrm{eV}=13.606\times1.602\times10^{-19}\ \mathrm{J}=2.180\times10^{-18}$ J,$m=9.109\times10^{-31}$ kg より

$$ \braket{p^2}=2\times9.109\times10^{-31}\times2.180\times10^{-18}=3.971\times10^{-48}\ \mathrm{kg^2\,m^2/s^2}, $$ $$ \Delta p=\sqrt{\braket{p^2}-0}=1.993\times10^{-24}\ \mathrm{kg\,m/s} . $$

$\Delta x\simeq a_0=0.529\times10^{-10}$ m と見積もると

$$ \Delta x\,\Delta p\simeq0.529\times10^{-10}\times1.993\times10^{-24}=1.054\times10^{-34}\ \mathrm{J\,s}=\hbar . $$

Heisenbergの不確定性関係 $\Delta x\,\Delta p\ge\hbar/2$ の下限のすぐ近くにある.水素原子の1s軌道は,不確定性関係が許すぎりぎりまで電子を核に近づけた状態だと読める.もっと縮めれば引力によるエネルギーは下がるが,$\Delta p$ が増えて運動エネルギーが上がる.この綱引きの釣り合いが原子の大きさ $a_0$ を決めているのである.

3.9.3 熱力学:比熱を自由エネルギーで書く

ここから熱統計力学へ足を伸ばす.古典熱力学では,定積比熱(定積熱容量)は

$$ \begin{equation} C_V\equiv\left(\frac{\partial U}{\partial T}\right)_V \label{eq:3-cv-def} \end{equation} $$

と定義される.$U$ は内部エネルギーである.一方,Helmholtzの自由エネルギーは

$$ \begin{equation} F\equiv U-TS \label{eq:3-f-def} \end{equation} $$

で定義される($S$ はエントロピー).式 \eqref{eq:3-f-def} を体積一定の条件で温度微分すると,積の微分法により

$$ \begin{align} \left(\frac{\partial F}{\partial T}\right)_V &=\left(\frac{\partial U}{\partial T}\right)_V-T\left(\frac{\partial S}{\partial T}\right)_V-S\left(\frac{\partial T}{\partial T}\right)_V \nonumber\\ &=C_V-T\left(\frac{\partial S}{\partial T}\right)_V-S \label{eq:3-dfdt-thermo} \end{align} $$

となる.ここで熱力学の基本関係式 $\dd F=-S\,\dd T-p\,\dd V$ から得られる

$$ \left(\frac{\partial F}{\partial T}\right)_V=-S $$

を式 \eqref{eq:3-dfdt-thermo} の左辺に代入すると

$$ -S=C_V-T\left(\frac{\partial S}{\partial T}\right)_V-S $$

となり,両辺の $-S$ が消えて

$$ \begin{equation} C_V=T\left(\frac{\partial S}{\partial T}\right)_V \label{eq:3-cv-entropy} \end{equation} $$

を得る.これは「体積一定の条件では,温度上昇によってエントロピーが増加するから,熱が物質内に蓄積される」ことを意味している.エントロピーが増えない物質は熱を溜められない,というわけである.

さらに $S=-(\partial F/\partial T)_V$ を式 \eqref{eq:3-cv-entropy} に代入すれば,比熱を自由エネルギーだけで書ける:

$$ \begin{equation} C_V=T\left(\frac{\partial S}{\partial T}\right)_V =T\left(\frac{\partial}{\partial T}\right)_V\left(-\frac{\partial F}{\partial T}\right)_V =-T\left(\frac{\partial^2F}{\partial T^2}\right)_V . \label{eq:3-cv-second-derivative} \end{equation} $$

3.9.4 熱統計力学:分配関数から出発する

熱統計力学では,系がとりうるエネルギー準位を $E_1,E_2,\dots,E_N$ として分配関数(partition function)

$$ \begin{equation} Z=\sum_{i=1}^{N}e^{-E_i/k_{\mathrm B}T} \label{eq:3-partition} \end{equation} $$

を定義し,そこから内部エネルギーとHelmholtz自由エネルギーを

$$ \begin{equation} U\equiv\frac{1}{Z}\sum_{i=1}^{N}E_i\,e^{-E_i/k_{\mathrm B}T}=\braket{E}, \qquad F\equiv-k_{\mathrm B}T\ln Z \label{eq:3-u-stat} \end{equation} $$

と定義する.$k_{\mathrm B}$ はBoltzmann定数である.$U$ の式は「各準位のエネルギー $E_i$ を,その出現確率 $e^{-E_i/k_{\mathrm B}T}/Z$ で重みづけして平均したもの」であり,まさにエネルギーの期待値 $\braket{E}$ である.同様に,エネルギーの2乗の期待値は

$$ \braket{E^2}=\frac{1}{Z}\sum_{i=1}^{N}E_i^2\,e^{-E_i/k_{\mathrm B}T} $$

と書ける.これらを式 \eqref{eq:3-cv-second-derivative} に代入して比熱を計算してみよう.

導出:$C_V$ をエネルギーのゆらぎで表す

準備:分配関数の温度微分.指数関数の微分から始める.$u=-E_i/(k_{\mathrm B}T)$ と置くと $\dfrac{\dd u}{\dd T}=+\dfrac{E_i}{k_{\mathrm B}T^2}$ であるから($T$ が分母にあるので符号が反転することに注意),

$$ \frac{\partial}{\partial T}e^{-E_i/k_{\mathrm B}T}=\frac{E_i}{k_{\mathrm B}T^2}\,e^{-E_i/k_{\mathrm B}T} . $$

したがって

$$ \begin{equation} \left(\frac{\partial Z}{\partial T}\right)_V=\frac{1}{k_{\mathrm B}T^2}\sum_{i=1}^{N}E_i\,e^{-E_i/k_{\mathrm B}T} =\frac{Z\braket{E}}{k_{\mathrm B}T^2} . \label{eq:3-dzdt} \end{equation} $$

ステップ1:$F$ の1階微分.$F=-k_{\mathrm B}T\ln Z$ を積の微分法で微分する:

$$ \begin{align} \left(\frac{\partial F}{\partial T}\right)_V &=-k_{\mathrm B}\ln Z-k_{\mathrm B}T\cdot\frac{1}{Z}\left(\frac{\partial Z}{\partial T}\right)_V \nonumber\\ &\overset{\eqref{eq:3-dzdt}}{=}-k_{\mathrm B}\ln Z-k_{\mathrm B}T\cdot\frac{1}{Z}\cdot\frac{Z\braket{E}}{k_{\mathrm B}T^2} \nonumber\\ &=-k_{\mathrm B}\ln Z-\frac{1}{T}\sum_{i=1}^{N}E_i\,e^{-E_i/k_{\mathrm B}T}\Big/Z =-k_{\mathrm B}\ln Z-\frac{U}{T} . \label{eq:3-dfdt-stat} \end{align} $$

ちなみにこの式は $F=U-TS$ と $(\partial F/\partial T)_V=-S$ から $S=k_{\mathrm B}\ln Z+U/T$,すなわち $F=-k_{\mathrm B}T\ln Z$ が確かに $U-TS$ と整合することを示している.

ステップ2:$F$ の2階微分(第1項).式 \eqref{eq:3-dfdt-stat} をもう一度 $T$ で微分する.まず第1項は,\eqref{eq:3-dzdt} を使って

$$ \frac{\partial}{\partial T}\bigl(-k_{\mathrm B}\ln Z\bigr) =-k_{\mathrm B}\frac{1}{Z}\left(\frac{\partial Z}{\partial T}\right)_V =-k_{\mathrm B}\cdot\frac{\braket{E}}{k_{\mathrm B}T^2} =-\frac{U}{T^2} . $$

ステップ3:$F$ の2階微分(第2項).第2項 $-U/T$ の微分は商の微分法により

$$ \frac{\partial}{\partial T}\left(-\frac{U}{T}\right) =-\frac{1}{T}\left(\frac{\partial U}{\partial T}\right)_V+\frac{U}{T^2} . $$

ステップ4:足し合わせる.ステップ2と3を足すと $-U/T^2$ と $+U/T^2$ が打ち消し合い,

$$ \left(\frac{\partial^2F}{\partial T^2}\right)_V=-\frac{1}{T}\left(\frac{\partial U}{\partial T}\right)_V . $$

あとは $\left(\dfrac{\partial U}{\partial T}\right)_V$ を分配関数から具体的に計算すればよい.

ステップ5:$U$ の温度微分.$U=\dfrac{1}{Z}\displaystyle\sum_i E_ie^{-E_i/k_{\mathrm B}T}$ を,商の微分法で微分する.分子の微分と分母の微分の2つの寄与が出る:

$$ \begin{align} \left(\frac{\partial U}{\partial T}\right)_V &=\frac{1}{Z}\sum_{i=1}^{N}E_i\cdot\frac{E_i}{k_{\mathrm B}T^2}e^{-E_i/k_{\mathrm B}T} \ -\ \frac{1}{Z^2}\left(\frac{\partial Z}{\partial T}\right)_V\sum_{i=1}^{N}E_ie^{-E_i/k_{\mathrm B}T} \nonumber\\ &=\frac{1}{k_{\mathrm B}T^2}\cdot\frac{1}{Z}\sum_{i=1}^{N}E_i^2\,e^{-E_i/k_{\mathrm B}T} \ -\ \frac{1}{Z^2}\cdot\frac{Z\braket{E}}{k_{\mathrm B}T^2}\cdot Z\braket{E} \nonumber\\ &=\frac{1}{k_{\mathrm B}T^2}\braket{E^2}-\frac{1}{k_{\mathrm B}T^2}\braket{E}^2 =\frac{1}{k_{\mathrm B}T^2}\Bigl(\braket{E^2}-\braket{E}^2\Bigr) . \label{eq:3-dudt} \end{align} $$

ここで2行目の第2項では,$\sum_iE_ie^{-E_i/k_{\mathrm B}T}=Z\braket{E}$ を2度使った.

ステップ6:仕上げ.ステップ4と5を合わせると

$$ \begin{equation} \left(\frac{\partial^2F}{\partial T^2}\right)_V =-\frac{1}{T}\cdot\frac{1}{k_{\mathrm B}T^2}\Bigl(\braket{E^2}-\braket{E}^2\Bigr) \label{eq:3-d2fdt2} \end{equation} $$

であり,これを式 \eqref{eq:3-cv-second-derivative} に代入すれば

$$ C_V=-T\left(\frac{\partial^2F}{\partial T^2}\right)_V =\frac{1}{k_{\mathrm B}T^2}\Bigl(\braket{E^2}-\braket{E}^2\Bigr) $$

を得る.$-T$ と $-1/T$ が掛かって $+1$ になったところが最後の詰めである.

結論を書き下そう.

$$ \begin{equation} C_V=\frac{\braket{E^2}-\braket{E}^2}{k_{\mathrm B}T^2} =\frac{(\Delta E)^2}{k_{\mathrm B}T^2} . \label{eq:3-cv-fluctuation} \end{equation} $$

分子は,式 \eqref{eq:3-dp2} とまったく同じ形——エネルギーの分散,すなわちエネルギーの不確定性の2乗 $(\Delta E)^2$ である.

物理的意味:比熱はエネルギーのゆらぎである

式 \eqref{eq:3-cv-fluctuation} は驚くべきことを言っている.比熱という,熱量計で測れるごく日常的な量が,系のエネルギーのゆらぎ(不確定性)そのものだというのである.直感的には次のように読める.温度 $T$ の熱浴に接した系は,平均 $\braket{E}$ のまわりでエネルギーが揺らいでいる.この揺らぎが大きい系は,熱をやりとりする「余地」が大きい系である.だから温度を少し上げたとき,たくさんのエネルギーを飲み込むことができる——それが比熱が大きいということである.逆に,エネルギー準位が1つしかない(=ゆらげない)系は,いくら温めてもエネルギーを蓄えられず $C_V=0$ となる.

この「応答係数=ゆらぎ」という構造は熱統計力学の至るところに現れる(ゆらぎ・散逸定理,fluctuation–dissipation theorem).たとえば磁化率は磁化のゆらぎ,圧縮率は密度のゆらぎである.第1章で述べた「計算科学で求まる物理量」の多くが,こうしたゆらぎの計算に帰着することを覚えておいてほしい.分子動力学法で比熱を求めるときは,実際に式 \eqref{eq:3-cv-fluctuation} を使ってエネルギーの時系列の分散を取る.

例題3.7 2準位系の比熱

エネルギー $0$ と $\Delta$($\Delta\gt0$)の2つの準位だけをもつ系について,式 \eqref{eq:3-cv-fluctuation} から $C_V$ を求めよ.また高温極限 $k_{\mathrm B}T\gg\Delta$ と低温極限 $k_{\mathrm B}T\ll\Delta$ での振る舞いを調べよ.

解答 $x\equiv\Delta/k_{\mathrm B}T$ と置く.分配関数は

$$ Z=e^{0}+e^{-x}=1+e^{-x} . $$

期待値は

$$ \braket{E}=\frac{0\cdot1+\Delta e^{-x}}{1+e^{-x}}=\frac{\Delta e^{-x}}{1+e^{-x}}, \qquad \braket{E^2}=\frac{0^2\cdot1+\Delta^2e^{-x}}{1+e^{-x}}=\frac{\Delta^2e^{-x}}{1+e^{-x}} . $$

分散は,共通因子 $\Delta^2\dfrac{e^{-x}}{1+e^{-x}}$ でくくると

$$ \braket{E^2}-\braket{E}^2 =\Delta^2\frac{e^{-x}}{1+e^{-x}}\left[1-\frac{e^{-x}}{1+e^{-x}}\right] =\Delta^2\frac{e^{-x}}{1+e^{-x}}\cdot\frac{1}{1+e^{-x}} =\frac{\Delta^2e^{-x}}{(1+e^{-x})^2} . $$

分子分母に $e^{2x}$ を掛けて整理すると $\dfrac{\Delta^2e^{x}}{(e^{x}+1)^2}$.よって

$$ \begin{equation} C_V=\frac{1}{k_{\mathrm B}T^2}\cdot\frac{\Delta^2e^{x}}{(e^{x}+1)^2} =k_{\mathrm B}\,x^2\,\frac{e^{x}}{(e^{x}+1)^2}, \qquad x=\frac{\Delta}{k_{\mathrm B}T} . \label{eq:3-schottky} \end{equation} $$

これはSchottky比熱(Schottky specific heat)と呼ばれる有名な形である.

高温極限($x\to0$):$e^x\to1$ より $C_V\to k_{\mathrm B}x^2/4=k_{\mathrm B}(\Delta/k_{\mathrm B}T)^2/4\to0$.温度が高すぎると2つの準位がほぼ等確率で占有され,これ以上エネルギーを吸収できなくなるため比熱は減る.

低温極限($x\to\infty$):$e^x/(e^x+1)^2\simeq e^{-x}$ より $C_V\simeq k_{\mathrm B}x^2e^{-x}\to0$.温度が低すぎると励起準位に上がれないので,やはり熱を蓄えられない.

したがって $C_V$ は中間の温度($x\simeq2.4$,すなわち $k_{\mathrm B}T\simeq0.42\Delta$)でピークをもつ山型になる.エネルギーが最もよく揺らぐ温度で比熱が最大になるという,式 \eqref{eq:3-cv-fluctuation} の主張どおりの振る舞いである.この山は結晶場分裂した磁性イオンの比熱などで実際に観測される(文献[7]).

注意:記法の揺れについて

熱統計力学の教科書では,分配関数を $Z$ ではなく $\Xi$ や $Q$ と書いたり,$\beta\equiv1/k_{\mathrm B}T$ を使って $Z=\sum_ie^{-\beta E_i}$ と書いたりする.$\beta$ を使うと $C_V$ の導出はもっと簡潔になり,$\braket{E}=-\partial\ln Z/\partial\beta$,$\braket{E^2}-\braket{E}^2=\partial^2\ln Z/\partial\beta^2$ の2行で済む.本書では講義に合わせて $T$ で微分したが,$\beta$ 表記に慣れておくと後々便利である(文献[8]).また,本章の $Z$(分配関数)と,3.5節までの $Z$(原子番号)はまったく別物である.同じ文字を使うのは伝統的な不便さであり,文脈で判断してほしい.

3.10 まとめと演習

3.10.1 この章のまとめ

  1. ヘリウム原子のハミルトニアンは式 \eqref{eq:3-he-hamiltonian} の5項からなる.運動エネルギー2項,核-電子引力2項,電子間反発1項である.
  2. 厳密には解けない.$\psi=\phi_1\phi_2$ と置いて両辺を $\phi_1\phi_2$ で割ると式 \eqref{eq:3-separation-divided} になるが,電子間反発項だけが2つの座標を同時に含むため定数として括り出せず,変数分離が破綻する.
  3. 電子間反発は無視できない.無視すると $E=2Z^2\epsilon_{1s}=-108.8$ eV(式 \eqref{eq:3-nonint-energy})となり,実験値 $-79.005$ eV から 29.8 eV もずれる.
  4. ブラケット記法は多重積分の記号地獄を避けるための記法である.$\langle a|b\rangle$(式 \eqref{eq:3-braket-def})はスカラー,$\ket{b}\bra{a}$ は行列.規格直交性 $\langle\psi_n|\psi_m\rangle=\delta_{nm}$(式 \eqref{eq:3-orthonormal})と完全性条件 $\sum_n\ket{\psi_n}\bra{\psi_n}=\mathbb{I}$(式 \eqref{eq:3-completeness})が中心的な道具である.
  5. 試行関数 $\psi_T=\phi_{Z'}(\bm r_1)\phi_{Z'}(\bm r_2)$(式 \eqref{eq:3-trial-function})のエネルギー期待値は5つの項に分解でき,$\braket{1/r}=Z'/a_0$(式 \eqref{eq:3-mean-inv-r})と $e_0^2/4\pi\varepsilon_0a_0=-2\epsilon_{1s}$(式 \eqref{eq:3-hartree})を使って $\braket{E}=\left(-2Z'^2+4ZZ'-\tfrac54Z'\right)\epsilon_{1s}$(式 \eqref{eq:3-energy-quadratic})となる.電子間反発の項に現れる係数 $5/8$ の導出は第4章4.7節で行う(式 \eqref{eq:3-term5}).
  6. 平方完成(式 \eqref{eq:3-complete-square})により,$\epsilon_{1s}\lt0$ すなわち $-2\epsilon_{1s}\gt0$ に注意して最小値を読むと $Z'=Z-\tfrac{5}{16}=1.6875$,$\braket{E}=2\epsilon_{1s}\left(Z-\tfrac{5}{16}\right)^2=-77.49$ eV(式 \eqref{eq:3-emin}).誤差は 37.8% から 1.9% に激減する.
  7. 変分原理(式 \eqref{eq:3-variational}):任意の規格化された試行関数について $\bra{\Psi}\Ham\ket{\Psi}\ge E_0$.証明は固有関数の完全系展開 $\ket{\Psi}=\sum_n\ket{\psi_n}c_n$,$c_n=\langle\psi_n|\Psi\rangle$(式 \eqref{eq:3-expansion},\eqref{eq:3-cn})から,$\bra{\Psi}\Ham\ket{\Psi}=\sum_mE_m\abs{\langle\psi_m|\Psi\rangle}^2\ge E_0\sum_m\abs{\langle\psi_m|\Psi\rangle}^2=E_0$(式 \eqref{eq:3-norm-sum},\eqref{eq:3-proof-final})と進む.
  8. 試行関数は基底関数の線形結合(式 \eqref{eq:3-basis-expansion})で作る.局在基底 STO-3G・6-31G の構成は第5・6章で扱う.Gaussian, GAMESS, VASP, Quantum ESPRESSO などのソフトウェアはすべて変分原理に立脚している.
  9. 分散:$(\Delta p)^2=\braket{p^2}-\braket{p}^2$(式 \eqref{eq:3-dp2}).同じ構造が熱統計力学にも現れ,$C_V=\dfrac{\braket{E^2}-\braket{E}^2}{k_{\mathrm B}T^2}$(式 \eqref{eq:3-cv-fluctuation}),すなわち比熱はエネルギーのゆらぎである.

次の第4章では,本章と対をなすもう一つの近似法である摂動論を学ぶ.変分法が「試行関数を用意して最小化する」のに対し,摂動論は「解ける問題からのズレを次数展開する」方針である.ヘリウムに1次摂動を適用すると $\left(2Z^2-\tfrac54Z\right)\epsilon_{1s}=-74.83$ eV が得られ,本章の例題3.3(a) と一致する.そして,本章で結果だけ引用した電子間反発積分の係数 $5/8$ は,そこで完全に計算される.

3.10.2 演習問題

演習3.1 リチウムイオン Li$^+$ の変分計算

Li$^+$ は $Z=3$,電子2個の系であり,ヘリウムと同じ電子配置 $1s^2$ をもつ.(a) 式 \eqref{eq:3-zprime-opt} を使って有効核電荷 $Z'$ と全エネルギーを求めよ.(b) Li の第2イオン化エネルギー $75.64$ eV と第3イオン化エネルギー $122.45$ eV から実験値を求め,誤差を評価せよ.(c) 誤差(%)をヘリウムの場合と比べ,$Z$ が大きくなると変分法の相対誤差がどうなるか論ぜよ.

ヒント:(a) 式 \eqref{eq:3-zprime-opt} は $Z$ を含んだままの一般式である.$Z=3$ を入れれば $Z'=3-0.3125$.エネルギーは $2\epsilon_{1s}Z'^2$.(b) 実験値は2つのイオン化エネルギーの和に負号を付けたもの.(c) $\braket{E}\propto(Z-5/16)^2$ と,正しい値 $\propto Z^2$ の相対差を考えると,$Z$ が大きいほど $5/16$ の重みが相対的に小さくなる.

演習3.2 遮蔽定数とSlater則

(a) 本章で得た遮蔽定数 $s=Z-Z'=5/16$ が $Z$ に依らないことを,式 \eqref{eq:3-zprime-opt} から確認せよ.(b) 無機化学で使うSlater則では,1s軌道にいる相手電子の遮蔽定数を $0.30$ と取る.$5/16$ との差はどれくらいか.(c) He の有効核電荷が実測でおよそ $1.7$ とされることと,本章の $1.6875$ を比較せよ.

ヒント:(b) $5/16=0.3125$.差は 0.0125,すなわち約4%.Slater則は多数の原子のデータに合うよう丸めた経験則であり,1s殻については変分法の結果とほぼ一致する.(c) 実測の「有効核電荷」は,X線散乱因子や軌道半径から逆算した値であり,定義がやや異なることにも注意する.

演習3.3 $\braket{1/r}$ と $\braket{r}$ を自分で積分する

規格化された1s軌道 $\phi_{Z'}(\bm r)=\pi^{-1/2}(Z'/a_0)^{3/2}e^{-Z'r/a_0}$ について,(a) $\braket{1/r}=Z'/a_0$ を積分から示せ.(b) $\braket{r}=\dfrac{3a_0}{2Z'}$ を示せ.(c) $\braket{1/r}\neq1/\braket{r}$ であることを確認し,なぜ一致しないのか説明せよ.

ヒント:球対称なので $\int f(r)\dd^3r=4\pi\int_0^\infty f(r)r^2\dd r$ を使う.(a) では被積分関数が $r^1$,(b) では $r^3$ になる.積分公式 $\int_0^\infty r^ne^{-ar}\dd r=n!/a^{n+1}$ を $a=2Z'/a_0$ で使う.(c) 一般に凸関数 $g$ について $\braket{g(r)}\ge g(\braket{r})$(Jensenの不等式).$1/r$ は凸関数である.

演習3.4 なぜ $Z'=2$ では $-77.49$ eV に届かないのか

試行関数として $Z'=Z=2$ の場合(=遮蔽を考えない場合)のエネルギー期待値を式 \eqref{eq:3-energy-quadratic} から求め,$-74.83$ eV になることを確かめよ.次に,この値が $Z'=1.6875$ の場合より高い理由を,変分原理の観点から述べよ.また,$Z'=1.6875$ のときのエネルギー $-77.49$ eV が,実験値 $-79.005$ eV より高い理由も述べよ.

ヒント:変分原理は「試行関数の族の中でエネルギーが最も低いものが,その族の中では最良の近似」と言っている.$Z'=2$ は族の中の1点にすぎず,最適点ではない.一方,$Z'=1.6875$ が実験値に届かないのは,族そのもの(単純な積の形の関数)が真の波動関数を含んでいないからである.図3.6の左図を思い出すとよい.真の波動関数は $r_{12}$ に陽に依存する(電子どうしが避け合う).

演習3.5 変分原理の等号成立条件

変分原理 $\bra{\Psi}\Ham\ket{\Psi}\ge E_0$ において,等号が成り立つのは $\ket{\Psi}$ が基底状態のときに限ることを,式 \eqref{eq:3-proof-final} の証明をたどって示せ.ただし基底状態には縮退がない($E_0\lt E_1$)とする.また,縮退がある場合はどうなるか述べよ.

ヒント:式 \eqref{eq:3-proof-final} の不等号は項ごとの評価 $E_m\abs{\langle\psi_m|\Psi\rangle}^2\ge E_0\abs{\langle\psi_m|\Psi\rangle}^2$ から来ている.等号が成り立つには,すべての $m$ について $(E_m-E_0)\abs{\langle\psi_m|\Psi\rangle}^2=0$ が必要である.$m\ge1$ では $E_m-E_0\gt0$ だから $\langle\psi_m|\Psi\rangle=0$,したがって $\abs{\langle\psi_0|\Psi\rangle}^2=1$ となる.縮退がある場合($E_0=E_1=\cdots=E_g$)は,その縮退した固有関数どうしの任意の線形結合でも等号が成り立つ.

演習3.6 調和振動子の比熱をゆらぎから求める

エネルギー準位が $E_n=\left(n+\tfrac12\right)h\nu$($n=0,1,2,\dots$)である1次元調和振動子について,(a) 分配関数 $Z$ を等比級数の和として求めよ.(b) 式 \eqref{eq:3-cv-fluctuation} を使って $C_V$ を求め,Einsteinの比熱の式

$$ C_V=k_{\mathrm B}\left(\frac{h\nu}{k_{\mathrm B}T}\right)^2\frac{e^{h\nu/k_{\mathrm B}T}}{\bigl(e^{h\nu/k_{\mathrm B}T}-1\bigr)^2} $$

が得られることを示せ.(c) 高温極限で $C_V\to k_{\mathrm B}$(Dulong–Petitの法則の1自由度分)になることを確かめよ.

ヒント:(a) $x\equiv h\nu/k_{\mathrm B}T$ と置くと $Z=e^{-x/2}\sum_{n=0}^{\infty}e^{-nx}=\dfrac{e^{-x/2}}{1-e^{-x}}$.(b) $\beta=1/k_{\mathrm B}T$ を使い $\braket{E}=-\partial\ln Z/\partial\beta$,$\braket{E^2}-\braket{E}^2=\partial^2\ln Z/\partial\beta^2$ を用いるのが早い(3.9.4節の注意を参照).零点エネルギー $h\nu/2$ は定数なので分散には寄与しない.(c) $x\to0$ で $e^x-1\simeq x$ を使う.この結果は第1章で触れたフォノンの比熱の出発点である.

3.10.3 参考文献

  1. 原田 義也『量子化学(上巻)』裳華房(2007). 第5章に,ヘリウム原子の変分計算と有効核電荷 $Z'=Z-5/16$ の導出が詳しい.
  2. A. Szabo, N. S. Ostlund『新しい量子化学 — 電子構造の理論入門(上)』大野公男・阪井健男・望月祐志 訳, 東京大学出版会(1987)(原著:Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory, Dover, 1996). 第1章に変分原理,第3章に STO-3G・6-31G の説明と H$_2$ のポテンシャル曲線の比較図がある.
  3. P. Atkins, J. de Paula, R. Friedman, Atkins' Physical Chemistry, 12th ed., Oxford University Press(2022). ヘリウム原子・遮蔽・Slater則の標準的な記述.
  4. L. Pauling, E. B. Wilson, Introduction to Quantum Mechanics with Applications to Chemistry, Dover(1985). 第7章が変分法の古典的名解説.
  5. E. A. Hylleraas, "Neue Berechnung der Energie des Heliums im Grundzustande, sowie des tiefsten Terms von Ortho-Helium", Z. Physik 54, 347 (1929). 試行関数に $r_{12}$ を陽に含めることで精度を飛躍的に高めた歴史的論文.
  6. C. L. Pekeris, "$1{}^1S$, $2{}^1S$, and $2{}^3S$ States of H$^-$ and of He", Phys. Rev. 115, 1216 (1959). 変分法を極限まで押し進めたヘリウムの高精度計算.
  7. C. Kittel『固体物理学入門(上)』宇野良清ほか 訳, 丸善(第8版, 2005). Schottky比熱,Einstein模型・Debye模型による格子比熱.
  8. 田崎 晴明『統計力学 I』培風館(2008). 分配関数からの $C_V$ の導出と,ゆらぎと応答の関係.$\beta$ 表記に慣れるのによい.
  9. Q. Sun et al., "PySCF: the Python-based simulations of chemistry framework", WIREs Comput. Mol. Sci. 8, e1340 (2018). 付録Bで使う量子化学計算ライブラリ.本章のヘリウム計算を数行で追試できる.
  10. Y. Mochizuki et al., "Theoretical exploration of mixed-anion antiperovskites as candidate semiconductors", Phys. Rev. Materials 4, 044601 (2020). 変分原理に立脚した第一原理計算(密度汎関数理論)を新規半導体探索に応用した例.
  11. NIST Atomic Spectra Database, National Institute of Standards and Technology. 本章で用いたイオン化エネルギーの実験値(He: 24.587 eV, 54.418 eV / Li: 75.640 eV, 122.454 eV)の出典.