バンド理論の勉強 のバックアップ(No.8)

更新


公開メモ

若林克法先生の授業のノートをもとに勉強しなおす

2015-02-24~26 に、当時 NIMS にいらした若林克法先生 (現在は関西学院大学) が 筑波大学でグラフェンの特異なナノスケール効果について授業をしてくださいました。

その際のノートなどを見ながら、バンド理論について勉強しなおそうと思います。

原子軌道の線形結合により分子軌道を構成(LCAO法の基礎)

wikipedia:LCAO法 などに紹介されている方法です。

水素原子 H

水素原子では1つの陽子の周りを1つの電子が回っている。この電子に対する時間を含まないシュレーディンガー方程式は、

\left\{-\frac{\hbar^2}{2m}\nabla^2-\frac{e^2}{r}\right\}\phi(\bm r)=\varepsilon\phi(\bm r)

これを解いた基底状態は 1s 軌道と呼ばれる。

基底状態のエネルギーを \varepsilon_{1s}

固有関数を \phi_{1s}(\bm r) としよう。

水素分子 H2

陽子が2つあるとき、その周辺での「電子1つに対する」ハミルトニアンは、

\hat H=\frac{-\hbar^2}{2m}\nabla^2-\frac{e^2}{r_1}-\frac{e^2}{r_2}

ただし、 r_1, r_2 はそれぞれの陽子からの距離、となる。

それぞれの陽子を中心とした2つの 1s 軌道の波動関数を元にして、 それらの線形結合によりこの「分子軌道」を近似的に表したい。

\Phi(\bm r)=c_1\phi_1(\bm r)+c_2\phi_2(\bm r)

ものすごく大雑把な近似

\phi_1(\bm r) は陽子 1 の周り以外でほぼゼロになると考えて、

-\frac{e^2}{r_2}\times\phi_1(\bm r)=0

と近似すれば、

\hat H\phi_1(\bm r)=\varepsilon_{1s}\phi_1(\bm r)

同様に、

\hat H\phi_2(\bm r)=\varepsilon_{1s}\phi_2(\bm r)

とすれば、任意の c_1,c_2 に対して

\hat H\Phi(\bm r)=\varepsilon_{1s}\Phi(\bm r)

となって、 \Phi(\bm r) はエネルギー \varepsilon_{1s} を持つ固有状態となる。

これは2つの原子が遠く離れて別々に置かれており、 陽子1の周りに |c_1|^2 の確率で、 陽子2の周りに |c_2|^2 の確率で、電子が見出される状態に対応する。

ヒュッケル法

\phi_1(\bm r) にある電子に対する陽子2からのポテンシャルが完全には無視できない場合にも、 \phi_1,\phi_2 の重なりが小さい場合には c_1,c_2 を適切に選ぶことで \Phi は近似的に \hat H の固有関数となる。

\Phi(\bm r)=c_1\phi_1(\bm r)+c_2\phi_2(\bm r) \ket \Phi=c_1\ket 1+c_2\ket 2 と書くことにする。

\hat H\ket\Phi=E\ket\Phi となるように c_1,c_2 を決めるため、左から \bra 1 をかけてみる。

&math(c_1\underbrace{\bra 1\hat H\ket 1}_{h_{11}}+c_2\underbrace{\bra 1\hat H\ket 2}_{h_{12}}= c_1E\underbrace{\braket{1|1}}_{=\,1}+c_2E\underbrace{\braket{1|2}}_{=\,\delta});

ここで、 \delta\ll 1 つまり \phi_1,\phi_2 の重なり積分 \delta=\braket{1|2} が十分に小さいとして無視すると(これを wikipedia:ヒュッケル法 と呼ぶ)、

&math( \begin{cases} c_1h_{11}+c_2h_{12}=c_1E\\ c_1h_{21}+c_2h_{22}=c_2E\\ \end{cases} ); つまり &math( \begin{pmatrix} h_{11}&h_{12}\\ h_{21}&h_{12}\\ \end{pmatrix}\begin{pmatrix}c_1\\c_2\end{pmatrix}= E\begin{pmatrix}c_1\\c_2\end{pmatrix} );

を得る。(2行目の式は左から \bra 2 をかけて得た)

h_{12}=\bra 1\hat H\ket 2=\big(\bra 2\hat H^\dagger\ket 1\big)^\dagger=\big(\bra 2\hat H\ket 1\big)^*=h_{21}^*

より左辺の行列はエルミートであり、実数固有値と正規直交固有ベクトルを持つ。

$h_{11},h_{22}$ について

&math( h_{11}&=\braket{1|\hat H|1}\\ &=\int d\bm r\phi_1^*(\bm r)\left[-\frac{\hbar^2}{2m}\nabla^2-\frac{e^2}{r_1}-\frac{e^2}{r_2}\right]\phi_1(\bm r)\\ &=\underbrace{\int d\bm r\phi_1^*(\bm r)\left[-\frac{\hbar^2}{2m}\nabla^2-\frac{e^2}{r_1}\right]\phi_1(\bm r)}_{\displaystyle \varepsilon_{1s}}

  • \underbrace{\int d\bm r\phi_1^*(\bm r)\frac{e^2}{r_2}\phi_1(\bm r)}_{\displaystyle \tilde\alpha}\\ &=\varepsilon_{1s}-\tilde\alpha\equiv-\alpha );

同様にして h_{22}=-\alpha となる。

$h_{21},h_{12}$ について

&math( h_{21}&=\braket{2|\hat H|1}\\ &=\int d\bm r\phi_2^*(\bm r)\left[-\frac{\hbar^2}{2m}\nabla^2-\frac{e^2}{r_1}-\frac{e^2}{r_2}\right]\phi_1(\bm r)\\ &=\int d\bm r\phi_2^*(\bm r)\varepsilon_1\phi_1(\bm r)

  • \underbrace{\int d\bm r\phi_2^*(\bm r)\frac{e^2}{r_2}\phi_1(\bm r)}_{\displaystyle \gamma}\\ &=\varepsilon_1\underbrace{\braket{2|1}}_{=\,\delta\ \sim\,0}-\gamma=-\gamma );

時間によらない波動関数は必ず実数関数に取れる。 \phi_1,\phi_2 が実数関数であるとき、 \gamma は実数になるため、以降では \gamma\in\mathbb R とする。

このとき、 h_{12}=h_{21}^*=h_{21}=-\gamma

$c_1,c_2$ を求める

&math( \begin{pmatrix} \alpha-E&-\gamma\\

  • \gamma&\alpha-E \end{pmatrix}\begin{pmatrix}c_1\\c_2\end{pmatrix}=\bm 0 );

より、自明でない解を持つためには

&math( \begin{vmatrix}

  • \alpha-E&-\gamma\\
  • \gamma&-\alpha-E \end{vmatrix}=(\alpha-E)^2-\gamma^2=0 );

したがって、

E=-\alpha\pm\gamma

E=-\alpha+\gamma のとき、 \begin{pmatrix}c_1\\c_2\end{pmatrix}=\frac{1}{\sqrt 2}\begin{pmatrix}1\\-1\end{pmatrix}

E=-\alpha-\gamma のとき、 \begin{pmatrix}c_1\\c_2\end{pmatrix}=\frac{1}{\sqrt 2}\begin{pmatrix}1\\1\end{pmatrix}

つまり、

\Phi\propto\phi_1-\phi_2 のとき、エネルギーは E=\epsilon_{1s}-\tilde\alpha+\gamma  (反結合性軌道)

\Phi\propto\phi_1+\phi_2 のとき、エネルギーは E=\epsilon_{1s}-\tilde\alpha-\gamma  (結合性軌道)

である。

後者は前者に比べ、両方の陽子からの引力相互作用を受けうる分子中心付近での状態密度が高いため、 \braket{1|\hat H|2} の分だけエネルギーが低くなっている。

ヒュッケル近似の意味

波動関数は重ならないが( \braket{1|2}=0
相互作用は取り入れる( \braket{1|\,e^2/r\,|2}\ne 0
という近似が必要になる意味があまり良く理解できていない・・・

重なり積分をゼロとせずにまじめに考えるには \braket{1|2}=\braket{1|2}=\delta\ne 0 として、
(固有関数を実数に取ると内積を入れ替えても値が変らない)

&math( \begin{pmatrix} \varepsilon_{1s}-\tilde\alpha&\varepsilon_{1s}\delta-\gamma\\ \varepsilon_{1s}\delta-\gamma&\varepsilon_{1s}-\tilde\alpha\\ \end{pmatrix}\begin{pmatrix}c_1\\c_2\end{pmatrix}= \begin{pmatrix} E&E\delta\\ E\delta&E\\ \end{pmatrix}\begin{pmatrix}c_1\\c_2\end{pmatrix} );

これが自明でない解を持つには、

&math( \begin{vmatrix} \varepsilon_{1s}-E-\tilde\alpha&(\varepsilon_{1s}-E)\delta-\gamma\\ (\varepsilon_{1s}-E)\delta-\gamma&\varepsilon_{1s}-E-\tilde\alpha\\ \end{vmatrix}=(\varepsilon_{1s}-E-\tilde\alpha)^2-((\varepsilon_{1s}-E)\delta-\gamma)^2=0 );

&math( &\varepsilon_{1s}-E-\tilde\alpha=\mp\{(\varepsilon_{1s}-E)\delta-\gamma\}\\ &(1\pm\delta)(\varepsilon_{1s}-E)=\tilde\alpha\pm\gamma\\ &\varepsilon_{1s}-E=\frac{\tilde\alpha\pm\gamma}{1\pm\delta}\\ );

この時の解は、

&math( &\begin{pmatrix} \varepsilon_{1s}-E-\tilde\alpha&(\varepsilon_{1s}-E)\delta-\gamma\\ (\varepsilon_{1s}-E)\delta-\gamma&\varepsilon_{1s}-E-\tilde\alpha\\ \end{pmatrix}\begin{pmatrix}c_1\\c_2\end{pmatrix}=\\ &\frac{1}{1\pm\delta}\begin{pmatrix} \tilde\alpha\pm\gamma-(1\pm\delta)\tilde\alpha&(\tilde\alpha\pm\gamma)\delta-(1\pm\delta)\gamma\\ (\tilde\alpha\pm\gamma)\delta-(1\pm\delta)\gamma&\tilde\alpha\pm\gamma-(1\pm\delta)\tilde\alpha\\ \end{pmatrix}\begin{pmatrix}c_1\\c_2\end{pmatrix}=\\ &\frac{1}{1\pm\delta}\begin{pmatrix} \pm\gamma\mp\delta\tilde\alpha&\tilde\alpha\delta-\gamma\\ \tilde\alpha\delta-\gamma&\pm\gamma\mp\delta\tilde\alpha\\ \end{pmatrix}\begin{pmatrix}c_1\\c_2\end{pmatrix}=\\ &\frac{\gamma-\tilde\alpha\delta}{1\pm\delta}\begin{pmatrix} \pm 1&-1\\

  • 1&\pm 1 \end{pmatrix}\begin{pmatrix}c_1\\c_2\end{pmatrix}=0\\ );

すなわち、

\begin{pmatrix}c_1\\c_2\end{pmatrix}\propto\begin{pmatrix}1\\\pm1\end{pmatrix} のとき、 E=\varepsilon_{1s}-\frac{\tilde\alpha\pm\gamma}{1\pm\delta}

となる。解の規格化は、

&math( \braket{\Phi|\Phi}&=A^2\big[\braket{1|1}\pm\braket{1|2}\pm\braket{2|1}+\braket{2|2}\big]\\ &=A^2(2\pm 2\delta)=1 );

より、 A=\frac{1}{\sqrt{2(1\pm\delta)}} となるから、

\begin{pmatrix}c_1\\c_2\end{pmatrix}=\frac{1}{\sqrt{2(1\pm\delta)}}\begin{pmatrix}1\\\pm1\end{pmatrix}

である。

ヒュッケル近似によりエネルギーが少しだけ変化したが、解の形は変らず、 定性的には良い近似になっていることが分かる(?)。

以下、特に断らない場合にはヒュッケル法を使う。

4つの原子からなる仮想的な四角い分子

4個の原子軌道(AO : atomic orbital): \ket 1,\ket 2,\ket 3,\ket 4

分子軌道(MO : molecular orbital): \ket\Phi=c_1\ket 1+c_2\ket 2+c_3\ket 3+c_4\ket 4

ヒュッケル近似で

  • 波動関数は重ならない \braket{i|j}=0\ \ (i\ne j)
  • クーロン積分 \braket{i|\hat H|i}=-\alpha
  • 飛び移り積分 \braket{i|\hat H|j}=-\gamma\ \ (i,j は隣り合う)
  • 隣接する原子核からだけ相互作用を受ける \braket{i|\hat H|j}=0\ \ (i,j は離れている)

とすると、1-2, 2-3, 3-4, 4-1 の間にのみ相互作用があるため、

&math( \begin{pmatrix}

  • \alpha & -\gamma & 0 & -\gamma\\
  • \gamma & -\alpha & -\gamma & 0 \\ 0 & -\gamma & -\alpha & -\gamma \\
  • \gamma & 0 & -\gamma & -\alpha \\ \end{pmatrix} \begin{pmatrix}c_1\\c_2\\c_3\\c_4\end{pmatrix} = E\begin{pmatrix}c_1\\c_2\\c_3\\c_4\end{pmatrix} );

これを解くとエネルギーは

E=-a-2g, -a, -a, -a+2g

対応するベクトルは、

\psi=\frac{1}{\sqrt 4}\big(|1\rangle+|2\rangle+|3\rangle+|4\rangle\big)

\psi=\frac{1}{\sqrt 2}\big(-|2\rangle+|4\rangle\big)

\psi=\frac{1}{\sqrt 2}\big(-|1\rangle+|3\rangle\big)

\psi=\frac{1}{\sqrt 4}\big(-|1\rangle+|2\rangle-|3\rangle+|4\rangle\big)

つまり、

 ++  0-  -0  -+
 ++  +0  0+  +-

の形になる。

一番左は4つのボンドすべてが結合性、
一番右は4つのボンドすべてが反結合性に対応する。

真ん中の2つは対角線なのでまったく相互作用のない形。

この2つは縮退しているため、その一時結合、例えば

 --  +-
 ++  +-

なども固有値 -\alpha の固有関数になる。

これら結合性の相互作用が2つと反結合性の相互作用が2つとで 打ち消し合って、ちょうど何の相互作用もない時と同じ値になっている。

結晶を作る

一次元結晶 = $N$個の原子を直線状に並べる

上の4つの時を拡張すれば、解くべき固有値問題は

&math( \begin{pmatrix}

  • \alpha&-\gamma&&&\\
  • \gamma&-\alpha&-\gamma\\ &-\gamma&-\alpha&\ddots\\ &&\ddots&\ddots&-\gamma\\ &&&-\gamma&-\alpha \end{pmatrix} \begin{pmatrix} c_1\\c_2\\\vdots\\\vdots\\c_N \end{pmatrix} = E \begin{pmatrix} c_1\\c_2\\\vdots\\\vdots\\c_N \end{pmatrix} );

あるいは同じ意味だが、

&math( \begin{pmatrix} 0&-\gamma&&&\\

  • \gamma&0&-\gamma\\ &-\gamma&0&\ddots\\ &&\ddots&\ddots&-\gamma\\ &&&-\gamma&0 \end{pmatrix} \begin{pmatrix} c_1\\c_2\\\vdots\\\vdots\\c_N \end{pmatrix} = (E+\alpha) \begin{pmatrix} c_1\\c_2\\\vdots\\\vdots\\c_N \end{pmatrix} );

Mathematica を使えば解くのは簡単。

上記の通り \alpha は固有値をずらすだけなので、 便宜的に \alpha=0,\gamma=1 として解いて、固有値を小さいものから順にプロットすると、

LANG:mathematica
a = Array[If[Abs[#1 - #2] == 1 , -1, 0] &, {50, 50}];
solution = Eigensystem[a];
ListPlot[Sort[solution[[1]] // N]]
ListPlot[Table[SortBy[Transpose[{solution[[1]] // N, solution[[2]]}], First][[n]] // Last, 
   {n, 1, 5}], PlotMarkers -> Automatic, PlotLegends -> Automatic]

n-atoms-energy.png

何やら cos 的な関数が見えてくる。 (縦軸は \gamma を単位としたエネルギー、横軸はエネルギーの低いほうから順に振った量子数)

エネルギーの低いほうから固有ベクトルを5つプロットすると、 固定端に対する弦の振動のような関数が見える。

n-atoms-wave-function.png

エネルギーが低いものから順に、波長が短くなり、波数が大きくなっていく。

ブロッホ波

なぜ正弦波が現れたのか?

上記方程式を整理すると、

c_{n+1}+\frac{E+\alpha}{\gamma}c_n+c_{n-1}=0

となる。
ただし、 n=1\sim N c_0=c_{n+1}=0 とした。

この方程式は固体物理学で1次元格子振動を扱う際に学んだ、 「互いにバネで繋がれたN個の質点」の方程式とそっくりな形をしている。

互いにバネで繋がれたN個の質点 を思い出す

個体物理では1次元格子振動のモデルとして、 各質点の変位を u_i として、

K(u_{n+1}-u_n)-K(u_n-u_{n-1})=M\ddot u_n

という形の方程式を解いた。ただし K は質点を繋ぐバネのバネ定数、 M は質点の質量。

ここでも固定端の条件は u_0=u_{N+1}=0 で与えられる。

質点間の距離を \Delta x と置き、

u_n(t)=\sin(kn\Delta x)\tau(t)

と変数分離できることを仮定する。

ただし k は固定端の条件より、整数 m に対して k(N+1)\Delta x=m\pi を満たすとする。

&math( \frac{M}{K}\frac{\ddot \tau(t)}{\tau(t)}\sin\big(kn\Delta x\big)= &\sin\big(k(n+1)\Delta x\big)-2\sin\big(kn\Delta x\big)+\sin\big(k(n-1)\Delta x\big)\\ );

&math( \frac{M}{K}\frac{\ddot \tau(t)}{\tau(t)}\sin\big(kn\Delta x\big)= &\sin\big(kn\Delta x\big)\cos\big(k\Delta x\big)-\cos\big(kn\Delta x\big)\sin\big(k\Delta x\big)

  • 2\sin\big(kn\Delta x\big)+\\ &\sin\big(kn\Delta x\big)\cos\big(-k\Delta x\big)-\cos\big(kn\Delta x\big)\sin\big(-k\Delta x\big) \\ );

&math( \frac{M}{K}\frac{\ddot \tau(t)}{\tau(t)}\cancel{\sin\big(kn\Delta x\big)}= &\cancel{\sin\big(kn\Delta x\big)}\cos\big(k\Delta x\big)-\cancel{\cos\big(kn\Delta x\big)\sin\big(k\Delta x\big)}

  • 2\cancel{\sin\big(kn\Delta x\big)}+\\ &\cancel{\sin\big(kn\Delta x\big)}\cos\big(k\Delta x\big)+\cancel{\cos\big(kn\Delta x\big)\sin\big(k\Delta x\big)}\\ );

&math( \frac{M}{K}\frac{\ddot \tau(t)}{\tau(t)}=&2\cos\big(k\Delta x\big)-2\\ );

&math( \frac{\ddot \tau(t)}{\tau(t)}=&-2\frac{K}{M}\big\{1-\cos(k\Delta x)\big\}\equiv -\omega^2\\ );

のように右辺を位置にも時刻にもよらない定数 -\omega^2 に置ける。

ここから \tau(t)=\sin(\omega t+B) となって、解を

u_n^m(t)=A\sin(k_mn\Delta x)\sin(\omega_m t+B)

と書ける。ただし、 m=1,\dots,N に対して、

k_m=\frac{m\pi}{(N+1)\Delta x}

かつ、

\omega_m=\sqrt{2\frac{K}{M}\big\{1-\cos(m\pi n/(N+1))\big\}}

このとき分散関係( \omega(k) のグラフ)は \sqrt{K/M} を単位として、

LANG:mathematica
Plot[Sqrt[2 (1 - Cos[Pi x])], {x, 0, 1}]

&ref(): File not found: "1d-chain_dispersion1.png" at page "バンド理論の勉強";

こんな形になる。

閑話休題(ブロッホ波再び)

n 番目の原子位置を R_n=n\Delta R として、

c_n=e^{ikR_n}

と置けば、

e^{ik\Delta R}+\frac{E+\alpha}{\gamma}+e^{-ik\Delta R}=0

E=-\alpha-2\gamma\cos(k\Delta R)

より、例えば \alpha=\gamma=1 とすれば、

LANG:mathematica
Plot[-1 - Cos[Pi x], {x, 0, 1}]

&ref(): File not found: "1d-chain_dispersion2.png" at page "バンド理論の勉強";

そのまま cos 的な、こんなグラフになる。

上記の固体物理で学んだ分散関係とは一見異なって見えるが、それは E\leftrightarrow \omega^2 の関係にあるためであり、固体物理の分散関係で \omega^2 をプロットすればまったく同じ プロットが得られることになる。

このときの解は上で見たように 0 サイトと n+1 サイトを固定端とする正弦波的な形になる。

c_n^m\propto \sin(k_mR_n)

&math(k_m=\frac{m\pi}{(N+1)\Delta R}

周期境界条件

周期境界条件では

&math(k_m=\frac{2m\pi}{(N+1)\Delta R}

に対して、

c_n^{m,s}\propto \sin(k_mR_n)

と、

c_n^{m,c}\propto \cos(k_mR_n)

が縮退するため、固有値問題

&math( \begin{pmatrix} 0&-\gamma&&&\gamma\\

  • \gamma&0&-\gamma\\ &-\gamma&0&\ddots\\ &&\ddots&\ddots&-\gamma\\ \gamma&&&-\gamma&0 \end{pmatrix} \begin{pmatrix} c_1\\c_2\\\vdots\\\vdots\\c_N \end{pmatrix} = (E+\alpha) \begin{pmatrix} c_1\\c_2\\\vdots\\\vdots\\c_N \end{pmatrix} );

に対する固有値を小さいものから並べれば、

LANG:mathematica
solution = Array[ If[Abs[#1 - #2] == 1 || (#1 == 1 && #2 == 50) || (#1 == 50 && #2 == 1), -1, 0] &, {50, 50}] //
  Eigensystem[#] &;
 ListPlot[Sort[solution[[1]] // N]]
 Transpose[{solution[[1]] // N, solution[[2]] // N}] //
   SortBy[#, First] & // Map[Last, #] & // Orthogonalize //
       ListPlot[#[[{1, 2, 3, 4, 5}]], PlotMarkers -> Automatic, PlotLegends -> Automatic] &

cyclic-boundary_energy.png

固有ベクトルをエネルギーの低いものから5つ並べれば、

cyclic-boundary_waves.png

となる。

エネルギー最低は定数関数なので sin, cos の区別が無く、縮退していない。

エネルギー最高は互い違いに入れ替わるのでやはり縮退していない。

原子数が奇数の場合、互い違いに入れ替わるベクトルは周期境界条件を満たさず解にならないため、 エネルギー最高はやはり縮退する。

2次元格子:グラフェンの場合


Counter: 42035 (from 2010/06/03), today: 7, yesterday: 0