3.1. 水素原子を並べると電子バンドは生まれるのか

固体物理の教科書では、原子がたくさん集まると、原子の離散的なエネルギー準位が密集して「バンド」になると説明されます。この過程を、実際の電子状態計算で確かめてみましょう。

3.1.1. 計算対象と計算条件

計算するのは、等間隔に一直線に並んだ中性の水素原子鎖 \(\mathrm{H}_N\) です（下図）。隣り合う水素原子の間隔はすべて1.8 Åとし、原子核の位置は最適化せずに固定します。鎖の両端は開いており、周期境界条件は課しません。これを表 3.1にあるような設定で、量子化学計算ソフト Psi4 [1] を用いて計算します。

../_images/model_chains.png — 図 3.1 計算モデルの模式図。全原子間隔を1.8 Åに保ち、原子核を固定したまま原子数 \(N\) を増やします。原子を結ぶ線は並び方を示すためのもので、化学結合を表すものではありません。

表 3.1 計算条件
項目	条件
計算コード	Psi4 1.10
モデル	中性の有限水素原子鎖。開放端、原子核位置固定
原子間距離	1.8 Å
原子数	主に \(N\) = 2, 4, 8, 16, 24, 32, 48, 64
基底関数	STO-3G（H原子あたり1個の空間基底関数）
スピン制限計算	RHF、RKS-PBE
スピン非制限計算	UHF、UKS-PBE（\(N\) = 16, 32）
電子数とスピン	電子数 \(N\)、偶数原子数、\(N_\alpha = N_\beta = N/2\)
SCF収束条件	エネルギー・密度ともに \(10^{-9}\)、最大500反復

原子数依存性の計算には密度フィッティング（DF）SCFを、後半のスピン制限解・非制限解の比較と安定性解析には近似のない4中心積分（PK法）を使いました。PBEの積分グリッドはPsi4の既定値です。以下に示す数値はすべて収束したSCF解から取得しています。

3.1.2. Tight-Binding 模型による予想

孤立した水素原子の基底状態では、電子が1s軌道を占めます。水素原子を2個近づけると、2つの1s軌道から結合性軌道と反結合性軌道ができます。さらに原子を増やすと、 \(N\) 個の1s基底関数から \(N\) 個の分子軌道ができます。

中性の \(\mathrm{H}_N\) 鎖には \(N\) 個の電子があり、スピン制限の閉殻計算では、エネルギーの低い \(N/2\) 個の分子軌道を2電子ずつ占有します。1原子あたり1電子なので、1軌道のバンド模型でいえば 半充填 の状態です。

この状況を説明する最も簡単なモデルが、隣接原子間のホッピングだけを考える強束縛（tight-binding）模型です。各原子の軌道エネルギーを \(\varepsilon_0\) 、隣接原子間のホッピングを \(t\) とすると、開放端の \(N\) 原子鎖の固有値は

\[E_j = \varepsilon_0 - 2t\cos\left(\frac{j\pi}{N+1}\right), \qquad j = 1, \ldots, N\]

となります。一方、原子間隔 \(a\) の無限周期鎖のバンド分散は

(3.1)\[E(k) = \varepsilon_0 - 2t\cos(ka)\]

です。有限鎖の \(j\) 番目の準位に形式的に \(k_j a = j\pi/(N+1)\) という 擬似波数 を割り当てると、原子数を増やすほど、このコサインカーブを細かくサンプリングしていくことが分かります。

3.1.3. 分子軌道エネルギーがバンドになるか

RKS-PBE計算で得た分子軌道エネルギーを、上記の疑似波数に対して並べたのが図 3.2 です。

../_images/pbe_levels.png — 図 3.2 RKS-PBEで求めた分子軌道エネルギー（原子数ごとに軌道エネルギーの平均を差し引いています）。黒い破線は、式 (3.1) で \(t \simeq 1.87\) eV とした場合（後述）のバンドです。横軸は擬似波数 \(ka = j\pi/(N+1)\) 、中央にある垂直な点線は半充填位置（電子占有／非占有の境界） \(ka = \pi/2\) を示しています。

原子数が少ないうちは準位がまばらですが、 \(N\) = 32、64になると点が1本の滑らかな曲線を描きます。たしかに、分子軌道の集まりからバンド状のエネルギースペクトルが生まれる様子が、数値実験として直接確認できました。

ただし、ここで計算しているのは、あくまでも有限の孤立分子であり、周期結晶のバンド構造そのものではありません。擬似波数が実際の波動関数の空間振動と対応しているかは、軌道係数を解析して検証する必要があります。RKS-PBEの準位は黒い破線に近いものの、特に上下端付近ではずれています。

例えば \(N\) = 64の全軌道エネルギーに \(A+B\cos(ka)\) を最小二乗フィットすると、二乗平均平方根誤差は約 0.312 eV になります。\(C\cos(2ka)\) を加えて、フィッティングの自由度を上げると、約 0.041 eV まで減らせます。これは、最近接以遠の相互作用、軌道の重なり、自己無撞着ポテンシャルなどが、単純な模型との差に寄与しているためだと考えられますが、今は、これ以上の深堀はしないことにします。

3.1.4. RHFの場合

同じ水素鎖の計算を、今度はRHF法で行い、先ほどと同様に、原子数を変えて軌道エネルギーをプロットしてみると図 3.3 のようになります。RKS-PBEの場合とは大きく様子が異なり、はっきりとバンドギャップが現れているのがわかります。

../_images/rhf_levels.png — 図 3.3 RHFで求めた分子軌道エネルギー。軸の定義と黒い破線は図 3.2 と同じです。破線はPBEのギャップから求めた \(t \simeq 1.87\) eVの強束縛模型を比較の基準として示しており、RHFへのフィットではありません。縦軸の表示範囲はRHFの準位に合わせています。

表 3.2 および図 3.4 に、RHFおよびRKS-PBEによる、最高被占有軌道（HOMO）と最低空軌道（LUMO）のエネルギー差（HOMO-LUMOギャップ）の、水素鎖長に対する変化を示しています。これをみると、RKS-PBEの場合は、水素鎖長とともに HOMO-LUMO ギャップがどんどん小さくなっていくのに対し、RHFの場合は、水素鎖長によらず、明らかに有意な差が残ることが分かります。その結果、RHF法で水素鎖のバンドを描くと、図 3.3 のようにバンドギャップが現れることになります。

表 3.2 HOMO–LUMO軌道ギャップ（DF-SCF）
原子数 \(N\)	RHF (eV)	RKS-PBE (eV)
2	12.124	3.806
4	9.325	2.327
8	7.089	1.302
16	5.679	0.691
24	5.214	0.470
32	5.007	0.356
48	4.833	0.240
64	4.764	0.181

../_images/orbital_gaps.png — 図 3.4 原子数に対するHOMO–LUMO軌道ギャップ（縦軸は対数）。破線は有限サイズ則 \(11.75/(N+1)\) eVです。

RKS-PBEのギャップは原子数とともに速やかに小さくなり、原子数を2倍にするとおおむね半分になります。実際、次のような式

(3.2)\[\Delta_N^{\mathrm{PBE}} \simeq \frac{11.75\ \mathrm{eV}}{N+1}\]

を考えると、図 3.4 の点線で示されるように、HOMO-LUMO ギャップの水素鎖長依存性のプロットによくフィットすることが分かります。上式 (3.2) は、半充填の最近接強束縛模型での有限サイズギャップ

\[\Delta_N^{\mathrm{TB}} = 4t\sin\left(\frac{\pi}{2(N+1)}\right) \simeq \frac{2\pi t}{N+1}\]

と同じ形をしています。偶数 \(N\) の半充填鎖では左の式が厳密に成り立ち、右の近似は大きな \(N\) に対するものです。係数を比較し、フェルミ準位付近の準位間隔を対応させると、有効ホッピングは

\[t \approx 11.75/(2\pi) \approx 1.87 {\rm eV}\]

になります。この模型のバンド幅 \(4t \approx 7.48\) eV は、図 3.2 から読み取れる、RKS-PBEの \(N\) = 64での軌道エネルギーの最大値と最小値の差、約 8.01 eVに近い値です。

一方、RHFでは \(N\) = 64まで増やしても約4.8 eVのギャップが残ります。

3.1.5. 波動関数（分子軌道）の解析

あらためて、２つの手法で計算した \(N\) = 64の軌道エネルギーを並べてプロットしたものを図 3.5 に示します。PBEの準位は半充填位置 \(ka = \pi/2\) をはさんで途切れずにつながるのに対し、RHFでは占有軌道と空軌道の間に大きな段差があります。同じ原子配置と電子数なのに、なぜこれほど違う結果になるのでしょうか？

../_images/rhf_pbe_levels_n64.png — 図 3.5 \(N\) = 64でのRHFとRKS-PBEの分子軌道エネルギー。塗りつぶしは占有軌道、白抜きは空軌道です。軸の定義は図 3.2 と同じです。

原因を探るため、隣り合う原子の間の電子的な結びつきを調べます。重なり行列 \(S\) と（スピンについて和をとった）密度行列 \(D\) から、Löwdin直交化した密度行列

\[P^{\mathrm{L}} = S^{1/2} D S^{1/2}\]

をつくり、その隣接原子間の要素を 結合次数に相当する量

(3.3)\[B_i = P^{\mathrm{L}}_{i,i+1}\]

として使います（これはこの解析のための指標で、結合次数の唯一の定義ではありません）。例えば、最近接強束縛模型の場合であれば、波動関数

\[\phi_j(i) = \sqrt{\frac{2}{N+1}}\sin\left(\frac{ij\pi}{N+1}\right)\]

を使って以下のような評価を行うことになります。

(3.4)\[ B_i = 2\sum_{j=1}^{N/2}\phi_j(i)\phi_j(i+1).\]

結合次数の解析結果を図 3.6 に示します。

../_images/hydrogen_chain_bond_order_n32_from_csv.png — 図 3.6 N=32の各結合におけるLöwdin密度行列要素（式 (3.3) ）。Uniform TBは相互作用のない一様な最近接強束縛模型の場合（式 (3.4) ）。

RHFでは、鎖の中央でも \(B_i\) がおおよそ0.86、0.38、0.86、0.38、……と大きく交替しています。PBEでは0.67と0.61程度の小さな交替にとどまります。これは、水素原子間の結合の強さが、図 3.7 のようになっていることを示しています。

../_images/bond-order-image.png — 図 3.7 \(N\) = 32の水素鎖中央における結合次数の交替を表す模式図。数値は隣接原子間のLöwdin密度行列要素 \(B_i\) の代表値です。原子間距離は一定で、青と橙はそれぞれ大きい値と小さい値を示します。

このように、 原子核は等間隔のまま、電子による結合の強さだけが「強・弱・強・弱」と交替する状態 は、結合次数波あるいは結合秩序波（bond-order wave、BOW）と呼ばれています（原子そのものが交互に変位する Peierls 歪みとは区別が必要です）。

BOW状態が無限鎖でも保たれるとすると、電子状態の周期は \(a\) から \(2a\) に倍化します。これは、強弱2種類のホッピング \(t_1, t_2\) を交互に並べた Su-Schrieffer-Heeger (SSH) 模型 [5] と呼ばれる強束縛模型で記述できる状態となっており、簡単な計算から、折り返されたバンドの境界に

\[\Delta_{\mathrm{BOW}} = 2|t_1 - t_2|\]

のギャップが開くことが示されます。図 3.5 で、RHFの準位が半充填位置 \(ka = \pi/2\) で段差を持つことも、この描像と整合しています。

3.1.6. スピン分極を考慮した場合

ここまでの計算は、上向きスピンと下向きスピンの電子に同じ空間軌道を使わせるスピン制限計算でした。そこで、両スピンの軌道を独立に最適化できるUHFとUKS-PBEでも計算してみました。バンド風プロットの結果を図 3.8 および図 3.9 に示します。

../_images/hydrogen_chain_uhf_orbital_spectrum.png — 図 3.8 UHFの軌道エネルギー。横軸は軌道番号から割り当てた擬似波数で、塗りつぶしは占有軌道、白抜きは空軌道を示します。

../_images/hydrogen_chain_uks_pbe_orbital_spectrum.png — 図 3.9 UKS-PBEの軌道エネルギー。横軸と記号の意味は左図と同じですが、縦軸の範囲は異なります。

この結果から、まず \(\alpha\) スピン、 \(\beta\) スピンの軌道エネルギーが全く同じであることがわかります。そして今度は UHF と UKS-PBE の両方で大きなバンドギャップが現れています。表 3.3 に、バンドギャップ（HOMO-LUMOギャップ）とLöwdin のスピン電荷密度解析の結果を示します。この結果から、各サイトには平均して0.6 \(\sim\) 0.9 のスピン差が生じていることがわかります。

表 3.3 \(N\) = 32でのスピン制限解と非制限解の比較（PK積分）
手法	HOMO–LUMO軌道ギャップ (eV)	Löwdinスピン密度の平均絶対値
RHF	5.005	0
RKS-PBE	0.356	0
UHF	19.493	0.898
UKS-PBE	3.611	0.626

続いて図 3.10 に、UHF および UKS-PBE の結合次数解析の結果と、スピン密度分布の結果を示します。

../_images/spin_and_bonds_n32.png — 図 3.10 \(N\) = 32の鎖のLöwdin結合次数（上）とLöwdinスピン密度 \(\alpha-\beta\) （下）。RHFの結合交替に対し、UHFとUKS-PBEでは結合次数がほぼ一様になり、スピン密度が原子ごとに反転します。

RHFで大きかったBOWはほぼ消え、代わりに 隣り合う原子ごとにスピン密度の符号が反転する（反強磁性的）な状態 が得られていることがわかります。

非制限計算では、もう一つ、強磁性的な状態も実現する可能性がありますが、計算してみたところ、原子あたり 0.6 \(\sim\) 1 eV 近くエネルギーの高い不安定な状態でした。もちろん、RHF および RKS-PBE も、よりエネルギーの高い（\(N\) = 32では、UHF→RHFで1原子あたり約1.49 eV、UKS-PBE→RKS-PBEで1原子あたり約0.125 eVだけ不安定化する）状態です。

したがって、今回の計算で得られた中では、反強磁性的な解が最も安定な電子状態ということになります。そしてこの「スピンの向きが異なるサイトが隣同士になっている」という状況は、とりもなおさず、系の周期を2倍にするという効果につながることになります。

3.1.7. どの結果が「正しい」のか？

ここで注意したいのは、軌道ギャップが開いているのが正しいわけではない、ということです。反対に、スピン制限のPBE計算で軌道ギャップがほぼ \(1/N\) で小さくなったことも、真の多電子系が金属であることの証明ではありません。今回行った計算は全て「それぞれの理論手法の特徴を見た」という程度の話だと思っておくのが安全です。

誤解のないように言っておくと、DFTによる第一原理計算は、通常の状況であれば非常に信頼できる結果が得られる手法です。しかしその一方で「強相関電子系」と呼ばれる特別な系の記述を不得手とすることも知られています。

実は、水素原子鎖は、電子相関を扱う手法のベンチマークとして広く使われている系です。多体計算による研究では、原子間距離に応じてMott絶縁体的な振る舞い、準長距離の反強磁性相関、絶縁体–金属転移などの性質が報告されています [2] 。有限の水素鎖でHF解の安定性と対称性の破れを調べた古典的な研究もあります [3] 。また別の研究では、原子間距離1.4 Åの水素鎖について、PBEなどがほぼゼロのギャップを与える一方、密度行列繰り込み群（DMRG）による多電子計算では絶縁体ギャップが得られています [4] 。

今回の計算での教訓は、 同じ原子配置と電子数であっても、計算手法が違えば、得られる電子状態の解釈そのものが変わるケースもある ということです。