Lagrange の未定乗数法はなぜ効くのか

― 拘束条件つきの極値から，永年方程式・共有結合まで（マテリアル計算科学第 11 回）

「ある条件を守りながら，何かを最小にする」．物理はこの形の問題ばかりです． 与えられた電子数のもとでエネルギーを最小に， 与えられたエネルギーのもとでエントロピーを最大に， 与えられた時間で作用を最小に．これを解く道具が Lagrange の未定乗数法です．ここでは，まず「なぜあれで解けるのか」を絵で納得したうえで，分子 AB のエネルギー期待値に適用して永年方程式を導きます．

結論を先に拘束条件つきの極値では，目的関数の等高線と，拘束の曲線がちょうど接します．接するとは「2 つの勾配が平行」ということで，その比例定数が未定乗数です．そしてこの未定乗数は，ただの計算の道具ではありません ―― 分子軌道の問題では，エネルギー準位の半分（E = 2ε）という物理量そのものになります．

1. 絵で見る ―― なぜ「接する」のか

2 変数の関数 f(x, y) を，条件 g(x, y) = 0 のもとで最小にしたいとします． g = 0 は平面上の 1 本の曲線なので，その曲線の上だけを歩いて，いちばん低いところを探す問題です．

歩いている途中で，f がまだ変化しているなら ―― まだ下げる余地があります． 下げようがなくなるのは，曲線に沿って動いても f が変わらなくなったときです．曲線に沿った方向を t とすると，これは

(1)

\nabla f \cdot t = 0

ということです．ところが，拘束曲線に沿って歩いても g は 0 のままなので， $\nabla g \cdot t = 0$ も いつでも成り立っています．つまり ∇f も ∇g も，どちらも t に垂直． 2 次元で同じ方向に垂直なベクトルは平行しかありません（∇g ≠ 0，つまり拘束曲線がその点でなめらかで，つぶれていないことが前提です）．だから

(2)

\nabla f = λ \nabla g

この λ が未定乗数です．「未定」なのは，値が最初から分かっていないから ―― 解いてはじめて決まる量だからです．

図 1．等高線と拘束曲線が接するところ． 背景の色が f の値（青いほど小さく，赤いほど大きい），その上の細い線が等高線，緑の太線が拘束条件 g = 0．赤い点をスライダーで曲線に沿って動かすと， 接した瞬間に 2 本の矢印が平行になり，そこが曲線に沿って f が停留する点 ―― 最小（最大）の候補になります．

目的関数 f 拘束条件 g = 0 曲線の上の位置 t = 0.60

2. 3 番目の式が，拘束条件を呼び戻す

(2) は 2 本の式ですが，未知数は x, y, λ の 3 つ．1 本足りません．足りない 1 本は拘束条件そのもの g = 0 です．この 3 本をまとめて自動的に出す仕掛けが，次の関数です．

(3)

F (x, y, λ) = f (x, y) - λ g (x, y) \frac{\partial F}{\partial x} = \frac{\partial F}{\partial y} = \frac{\partial F}{\partial λ} = 0

ここが仕掛けです

\partial F / \partial λ = - g = 0

―― λ で微分すると，拘束条件そのものが戻ってきます．だから「拘束つきの 2 変数問題」を「拘束なしの 3 変数問題」に丸ごと置き換えられる．これが未定乗数法の全部です．λ を「もう 1 つの変数」として扱う，というだけの発想が効きます．

符号は好みです

F = f + λ g

と書く教科書も多くあります． λ の符号が変わるだけで，答えは同じです．講義ノートでは

F = E - ε (⟨ Ψ | Ψ ⟩ - 2)

と引く形を使っているので，以下それに合わせます．

3. 例題 ―― 分子 AB のエネルギーを最小にする

ここからが本題です．原子 A, B の原子軌道 ψ_A, ψ_B を混ぜて，分子軌道をつくります（LCAO 近似）．

(4)

Ψ = c_{A} ψ_{A} + c_{B} ψ_{B}

決めたいのは 2 つの係数だけです．決め方は「エネルギーがいちばん低くなるように」 ―― 変分原理です．使う積分は 3 種類．

記号	定義	意味
S	$∭ ψ_{A}^{*} ψ_{B} d V$	重なり積分．2 つの軌道がどれだけ空間的に重なっているか（0 ≦ S ≦ 1）
α_A, α_B	$∭ ψ_{A}^{*} \hat{H} ψ_{A} d V$ など	クーロン積分．結合する前の，原子の電子準位にあたる量（負）． $\hat{H}$ には相手の核による引力も入っています
β	$∭ ψ_{A}^{*} \hat{H} ψ_{B} d V$	共鳴積分．電子が A と B を行き来する効果（ふつう負）

係数が実数（ $c_{A} = c_{A}^{*}$ ）とすると，展開はすぐに終わります．

(5)

⟨ Ψ | \hat{H} | Ψ ⟩ = c_{A}^{2} α_{A} + c_{B}^{2} α_{B} + 2 c_{A} c_{B} β

(6)

⟨ Ψ | Ψ ⟩ = c_{A}^{2} + c_{B}^{2} + 2 c_{A} c_{B} S = 2

なぜ「= 2」なのか結合性軌道に入る電子が 2 個だからです．Ψ を「電子 2 個ぶんの密度」に規格化しています．したがって

E = ⟨ Ψ | \hat{H} | Ψ ⟩ / ⟨ Ψ | Ψ ⟩

は電子 1 個あたりのエネルギー ―― つまり分子軌道の準位で， 2 電子ぶんの合計はその 2 倍になります．
そして (6) を使うと，割り算が消えて

E = \frac{1}{2} (c_{A}^{2} α_{A} + c_{B}^{2} α_{B} + 2 c_{A} c_{B} β)

と，ただの 2 次式になります．

4. 未定乗数法を当てる

(7)

F = E - ε (⟨ Ψ | Ψ ⟩ - 2) = \frac{1}{2} (c_{A}^{2} α_{A} + c_{B}^{2} α_{B} + 2 c_{A} c_{B} β) - ε (c_{A}^{2} + c_{B}^{2} + 2 c_{A} c_{B} S - 2)

3 つの偏微分を 0 と置きます．

(8)

\frac{\partial F}{\partial c_{A}} = c_{A} α_{A} + c_{B} β - ε (2 c_{A} + 2 c_{B} S) = 0

(9)

\frac{\partial F}{\partial c_{B}} = c_{B} α_{B} + c_{A} β - ε (2 c_{B} + 2 c_{A} S) = 0

(10)

\frac{\partial F}{\partial ε} = - (c_{A}^{2} + c_{B}^{2} + 2 c_{A} c_{B} S - 2) = 0 （拘束条件そのもの）

(8)(9) を c_A, c_B について整理して行列で書くと，

(11)

[\begin{matrix} α_{A} - 2 ε & β - 2 ε S \\ β - 2 ε S & α_{B} - 2 ε \end{matrix}] [\begin{matrix} c_{A} \\ c_{B} \end{matrix}] = [\begin{matrix} 0 \\ 0 \end{matrix}]

5. 永年方程式 ―― なぜ行列式を 0 と置くのか

(11) には，いつでも $c_{A} = c_{B} = 0$ という解があります．でもそれは拘束条件 (6) の $⟨ Ψ | Ψ ⟩ = 2$ を満たしません ―― Ψ が消えてしまい，電子がどこにもいないことになる． 欲しいのは「自明でない解」です．

逆行列があると，自明な解しか残らないもし行列 M に逆行列があれば，Mc = 0 の両辺に M⁻¹ を掛けて c = 0 が確定してしまいます．だから逆行列があってはいけない．逆行列が無い条件は det M = 0 ―― これが永年方程式です．

幾何で言えば：(8) と (9) はそれぞれ原点を通る直線です．ふつう 2 直線は原点でしか交わりません． 行列式が 0 になる瞬間，2 本が重なって 1 本になり，交点が「直線まるごと」に化けます．そこではじめて，原点以外の解が生まれます．

(12)

| \begin{matrix} α_{A} - 2 ε & β - 2 ε S \\ β - 2 ε S & α_{B} - 2 ε \end{matrix} | = 0

展開すると ε の 2 次方程式になります．

(13)

4 (1 - S^{2}) ε^{2} + (4 β S - 2 α_{A} - 2 α_{B}) ε + (α_{A} α_{B} - β^{2}) = 0

(14)

ε = \frac{1}{4 (1 - S^{2})} [α_{A} + α_{B} - 2 β S \pm \sqrt{{(α_{A} + α_{B} - 2 β S)}^{2} - 4 (1 - S^{2}) (α_{A} α_{B} - β^{2})}]

これは一般化固有値問題です (11) は Hc = 2ε Sc と書けます（H はハミルトニアン行列，S は重なり行列）．原子軌道を n 個並べれば n×n になり，分子軌道が n 本出てきます． 第一原理計算のプログラムが実際に解いているのは，まさにこの式です．

6. 水素分子の場合 ―― 結合性と反結合性

同種の原子なら $α_{A} = α_{B} = α$ ． (14) の根号の中は $(2 α - 2 β S)^{2} - 4 (1 - S^{2}) (α^{2} - β^{2}) = (2 α S - 2 β)^{2}$ ときれいに平方になり，

(15)

ε_{b} = \frac{α + β}{2 (1 + S)} ε_{a} = \frac{α - β}{2 (1 - S)}

係数は (11) の 1 行目に戻して求めます（1 行目は $S α - β$ でくくれる形になるので， $S α \neq β$ ，つまり 2 つの根が重なっていないことが前提です）． $ε = ε_{b}$ のとき $c_{A} - c_{B} = 0$ ， $ε = ε_{a}$ のとき $c_{A} + c_{B} = 0$ ．拘束条件 (6) で大きさをそろえると，

(16)

Ψ_{b} = \frac{1}{\sqrt{1 + S}} (ψ_{A} + ψ_{B}) Ψ_{a} = \frac{1}{\sqrt{1 - S}} (ψ_{A} - ψ_{B})

未定乗数 ε の正体ここで

E = 2 ε

になります．つまり ε は分子軌道のエネルギー準位（の半分）そのものでした．「拘束条件の強さ」として導入した補助変数が，物理量そのものだったのです．解析力学の一般化力（拘束力），熱力学の化学ポテンシャルも，まったく同じ構図をしています．

図 2．永年方程式の根と，エネルギー準位． 左が行列式 D(ε) のグラフで，0 を横切る 2 点が根．右がそこから決まる準位図です． α_A と α_B を離していくと共有結合からイオン結合へ移り，共有結合度が下がっていくのが見えます．ただし読み取り値の「安定化エネルギー」は中性の原子 2 個（α_A + α_B）を基準に測った量なので，準位差を開くほど絶対値はむしろ大きくなります ―― 基準の取り方にご注意ください．

α_A = −13.6 eV α_B = −13.6 eV β = −9.0 eV S = 0.50

7. ギャップの式から読めること

結合性と反結合性の隔たりを計算します．(15) の差を取ると

(17)

Δ ε = ε_{a} - ε_{b} = \frac{S α - β}{1 - S^{2}} \underset{S \to 0}{⟶} - β

講義スライドとの違い（1 か所）スライド 14 では

Δ ε = (- β + S α) / 2 (1 - S^{2})

となっていますが，通分を最後まで進めると分母の 2 が約分され，(17) のように分母は $1 - S^{2}$ になります（分母が半分になるぶん，Δε の値はスライドの 2 倍です）．
数値で確かめると，α = −13.6, β = −9, S = 0.5 のとき ε_a − ε_b = −4.600 − (−7.533) = 2.933， (17) の値も 2.933，スライドの式では 1.467 になります．
スライド 16 の $Δ ε = \sqrt{{((α_{A} - α_{B}) / 2)}^{2} + β^{2}}$ や，Δε_c = −β の方は (17) と整合していますので，結論はそのままで大丈夫です．

ここから 3 つのことが読めます．

|β| が大きいほどギャップが開く． β は電子が A ↔ B を行き来する効果なので， 原子が近づいて軌道がよく重なるほど，結合は強くなる．
内殻の電子は，ふつう結合に使えない． 深い準位ほど軌道が原子核に強く縛られて広がらず，S も β も 0 に近づきます．するとギャップがほぼ 0 になり，結合性・反結合性の両方に電子が入って 下がった分と上がった分が打ち消し合う ―― 得がありません．
準位差はイオン結合ぶん． S = 0 とすると (14) から $Δ ε = \sqrt{{(\frac{α_{A} - α_{B}}{2})}^{2} + β^{2}}$ ． 直角三角形の斜辺の形です（Δε は根の差 ε_a − ε_b なので，図2の読み取り値の「ギャップ 2ε_a − 2ε_b」は，この 2 倍です）．

共有結合度とイオン結合度は「足して 100 %」ではありません

Δ ε_{c} = - β

（共有結合ぶん），

Δ ε_{i} = | α_{A} - α_{B} | / 2

（イオン結合ぶん）と置くと，

Δ ε^{2} = Δ ε_{c}^{2} + Δ ε_{i}^{2}

（S = 0 のとき）．共有結合度

Δ ε_{c} / Δ ε

とイオン結合度

Δ ε_{i} / Δ ε

は， 2 乗して足すと 1 になります．だから「共有 71 %・イオン 71 %」と出ても矛盾ではありません（どちらも 1/√2 ＝ちょうど半々のとき）．

準位差を開くと，共有結合ぶんの得はどう痩せるか S = 0 として (14) の低い方の根を書くと

2 ε_{b} = \frac{1}{2} [α_{A} + α_{B} - \sqrt{{(α_{A} - α_{B})}^{2} + 4 β^{2}}]

です．準位差が

| β |

に比べて大きいとき，根号を展開すると，深い方（ここでは B）の準位を基準にして

2 ε_{b} ≒ α_{B} - \frac{β^{2}}{| α_{A} - α_{B} |}

―― 電子 2 個ぶんでは

2 β^{2} / | α_{A} - α_{B} |

しか下がりません．準位差を開くほど，共鳴でかせげる得は痩せていくということで，次の §8 の「100 % のイオン結合に共有結合ぶんの得はない」は，その行き着く先です．
（図 2 の読み取り値「安定化エネルギー」は中性の原子 2 個が基準なので，ここの「深い方の準位を基準」とは基準が違います．）

8. 100 % イオン結合に，なぜ共有結合の得がないのか

電子が全部 B に行ってしまった状態を入れてみます． $c_{A} = 0$ ， $c_{B} = \sqrt{2}$ なら拘束条件 (6) を満たします．このとき

(18)

E = \frac{1}{2} (0 + 2 α_{B} + 0) = α_{B}

この模型では，孤立した原子 B の準位とまったく同じ β の項に

c_{A} c_{B}

が掛かっているので， 片方の係数が 0 になった瞬間，共鳴積分の効果が丸ごと消えます． 共鳴による安定化はゼロです ―― この 2 軌道の模型では， 100 % のイオン結合に共有結合ぶんの得はありません．
ただし，現実の NaCl のようなイオン結合が安定なのは，主にイオンどうしの静電引力 （結晶なら Madelung エネルギー）によるもので，その寄与はこの模型には入っていません．ここで言えるのは「共鳴積分による得が消える」ところまでです．

逆に言えば ―― 「電子がどちらの原子から来たのか分からない」ことそのものが，結合のエネルギー源です．講義ノートの言葉を借りれば，粒子は「かくれんぼ」をしているときにいちばん安定している，ということになります．

9. 数値で確かめる

ここまでの式を，シミュレーターに載っているコードそのもので確かめました（35 項目，すべて通過）．

確かめたこと	結果
(14)：2 次方程式の解が行列式 (12) の根になっているか	\|D(ε)\| < 8×10⁻¹⁵
(15)：同種原子の閉じた式 (α±β)/2(1±S) と一致するか	相対差 < 10⁻¹²
(2)：解の点で ∇(½⟨Ψ\|H\|Ψ⟩) と ∇⟨Ψ\|Ψ⟩ が平行か（外積 = 0）	\|∇F\| < 2×10⁻¹⁴，外積 < 5×10⁻¹⁴
(6)：解が拘束条件 ⟨Ψ\|Ψ⟩ = 2 を満たし，E = 2ε になるか	相対差 < 10⁻¹⁰
(16)：係数が 1/√(1±S) になるか	相対差 < 10⁻¹⁰
楕円上を 2000 点＋三分探索して数値的に見つけた最小点が，解析解と一致するか（4 例）	位置・エネルギーとも 10⁻⁸ 以下で一致
(18)：c_A = 0, c_B = √2 で E = α_B になるか	厳密に一致（本当の最小はさらに 0.26 eV 低い）
(17)：Δε = (Sα−β)/(1−S²)（スライド 14 の式は 2 倍ずれ）	実測 2.9333 ＝ (17) の値
共有結合度・イオン結合度が √(Δε_c² + Δε_i²) = Δε（S = 0，Δε = ε_a − ε_b）を満たすか	相対差 < 10⁻¹²
S(R) = e^−R(1+R+R²/3) の逆解き（波動関数の絵に使う核間距離）	往復誤差 < 2×10⁻¹⁶，S = 0.7529 → R = 1.400 a₀
β > Sα にすると，低い方の解が ψ_A − ψ_B（反結合の形）に入れ替わるか	入れ替わることを確認（β = −2, Sα = −6.8 eV）

10. まとめ

拘束条件つきの極値では，目的関数の等高線と拘束の曲線が接します．接する＝勾配が平行＝ $\nabla f = λ \nabla g$ (2)．
$F = f - λ g$ と置くと， $\partial F / \partial λ = 0$ が 拘束条件を自動的に呼び戻す．これが未定乗数法の全部です (3)．
分子軌道に当てると，2 つの係数についての連立方程式が出ます (11)．
自明でない解が要るので逆行列があってはならない ⇒ 行列式 = 0 ⇒ 永年方程式 (12)．
根は 2 つ． $β < S α$ のふつうの場合は，低い方が結合性軌道，高い方が反結合性軌道 (15)(16)．（ $β > S α$ になると，低い方の解が ψ_A − ψ_B の形に入れ替わります ―― §9 の最後の行．）
未定乗数 ε は，分子軌道のエネルギー準位の半分でした（E = 2ε）．
ギャップは 共有結合ぶん (−β) とイオン結合ぶん (|α_A−α_B|/2) の直角三角形 （(14) で S = 0 とした形）．片方の係数が 0 になった瞬間 β の効果が消えるので， 100 % のイオン結合には共有結合ぶんの得がありません (18)．

東京理科大学望月研究室／講義「マテリアル計算科学」第 11 回「永年方程式と共有結合」より． §1・§2 の一般論と，§7 のギャップの整理は講義スライドの外を補ったものです．数値検証は倍精度の直接計算，楕円上の三分探索，1s 重なり積分の二分法による．グラフで確かめたい方はシミュレーターへ．数式はブラウザ標準の MathML で組んでいます（外部ライブラリ・通信は一切ありません）．