バンド理論の勉強

(420d) 更新

公開メモ

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

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)

陽子2つに電子1つなので、正確には水素分子ではなく水素分子イオン H2+ を考えていることになる。

ものすごく大雑把な近似

\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 をかけてみる。

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:ヒュッケル法 と呼ぶ)、

\begin{cases} c_1h_{11}+c_2h_{12}=c_1E\\ c_1h_{21}+c_2h_{22}=c_2E\\ \end{cases} つまり \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}$ について

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}$ について

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$ を求める

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

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

\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  (結合性軌道)

である。

bonding_and_antibonding.png

結合性軌道は反結合性軌道に比べ、両方の陽子からの引力相互作用を受ける 分子中心付近での状態密度が高いため、 \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 とすると、

\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}

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

\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

&\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}\\

この時の解は、

&\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}

となる。解の規格化は、

\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}

である。

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

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

電子が複数ある影響

上で考えたのは陽子2つの周りを1つに電子がいる状況。

すなわち、水素分子 H2 ではなく水素分子イオン H2+ である。

水素分子では電子が2つになるため、電子・電子相互作用が入る。

通常、電子・電子相互作用は平均場近似で取り入れる。 すなわち他の電子により部分的に遮蔽された原子核からのポテンシャル中での 1電子の運動を考えることになる。

電子の空間分布が偏りを持つ場合、それぞれの原子核からのポテンシャルは等しくならないが、 ヒュッケル近似ではその影響は無視される。

最終的に1電子軌道から多電子軌道を構成する際には交換相互作用が現れる。

このようなことから、話はそれほど単純ではないとはいえ、 ヒュッケル近似から導かれる以下で見るような性質は、 実際の系でも強く表れる。

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

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}=0\ \ (i,j が離れているとき)
    • 飛び移り積分 \braket{i|\hat H|j}=-\gamma\ \ (i,j が隣り合うとき)

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

\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つの時を拡張すれば、解くべき固有値問題は

\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}

あるいは同じ意味だが、

\begin{pmatrix} 0&-1&&&\\ -1&0&-1\\ &-1&0&\ddots\\ &&\ddots&\ddots&-1\\ &&&-1&0 \end{pmatrix} \begin{pmatrix} c_1\\c_2\\\vdots\\\vdots\\c_N \end{pmatrix} = \frac{E+\alpha}{\gamma} \begin{pmatrix} c_1\\c_2\\\vdots\\\vdots\\c_N \end{pmatrix}

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

上式の通り \alpha は固有値をずらすだけ、 \gamma もスケールするだけなので、 便宜的に \alpha=0,\gamma=1 として解いて、固有値を小さいものから順にプロットすると、 ( N=50 とした)

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 で与えられる。少々変形して、

u_{n+1}-2u_n+u_{n-1}=\frac{M}{K}\ddot u_n

を行列で表せば、

\begin{pmatrix} -2&1\\ 1&-2&1\\ &1&\ddots&\ddots\\ &&\ddots&\ddots&1\\ &&&1&-2 \end{pmatrix} \begin{pmatrix} u_1\\u_2\\\vdots\\\vdots\\u_N \end{pmatrix} =\frac{M}{K}\frac{d^2}{dt^2} \begin{pmatrix} u_1\\u_2\\\vdots\\\vdots\\u_N \end{pmatrix}

したがって、

\begin{pmatrix} -2&1\\ 1&-2&1\\ &1&\ddots&\ddots\\ &&\ddots&\ddots&1\\ &&&1&-2 \end{pmatrix} \begin{pmatrix} u_1\\u_2\\\vdots\\\vdots\\u_N \end{pmatrix} =-\omega^2\frac{M}{K} \begin{pmatrix} u_1\\u_2\\\vdots\\\vdots\\u_N \end{pmatrix}

となる固有値・固有ベクトルが求まれば、

\frac{d^2}{dt^2} \begin{pmatrix} u_1\\u_2\\\vdots\\\vdots\\u_N \end{pmatrix}=-\omega^2 \begin{pmatrix} u_1\\u_2\\\vdots\\\vdots\\u_N \end{pmatrix}

より、

\bm u(t)=e^{\pm i\omega t}\bm u(0)

として方程式が解ける。

この固有値問題は N 次元対称行列に対するものであるから、 重複度を含めて N 個の実数固有値と、対応する一時独立な N 個の固有ベクトルが存在する。

良く知られるように、質点間の距離を \Delta x として、

u_n(t)=e^{ink\Delta x}

と置けば上記方程式は

\Big(e^{ik\Delta x}-2+e^{-ik\Delta x}\Big)e^{ink\Delta x}=-\omega^2\frac{M}{K}e^{ink\Delta x}

となり、両辺を e^{ink\Delta x} で割って、

e^{ik\Delta x}-2+e^{-ik\Delta x}=-\omega^2\frac{M}{K}

\omega=\pm\sqrt{\frac{2K}{M}}\sqrt{1-\cos ik\Delta x}

を満たす \omega,k が固有値、固有ベクトルの候補となる。

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

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

1d-chain-dispersion1.svg

こんな形になる。

ただし、上でも述べた通り k \omega の値は このグラフ上の任意の値をとれるわけではなく、 u_0=u_{N+1}=0 の条件から許される値はとびとびになる。

u_n(t)=e^{ink\Delta x} の形でこの「境界条件」を満たすことはできないが、 ある \omega の値に足して k が条件を満たすとき、同時に -k も条件を満たすことに注意して( \because\cos(\theta)=\cos(-\theta) )、

u_n(t)\propto e^{ink\Delta x}-e^{-ink\Delta x} \propto \sin nk\Delta x

と置き、 k_m=m\pi/\Delta x(N+1) とすれば、 m=1\dots N に対して \bm u_m は線形独立になる。例えば N=5 の場合、 \{u_m\}_{m=1\dots N} は、

LANG:mathematica
Module[{nn}, nn = 5;
  Plot[Table[Sin[n m Pi/(nn + 1)], {m, 1, nn}] // Evaluate, {n, 0, nn + 1},
   Mesh -> 5, MeshStyle -> Directive[PointSize[Large]], ImageSize -> Large]
]

sine-waves.svg

グラフにはサイン波を表示したが、実際には格子点上の値しか意味を持たないため、 例えば m=(N+1)-1 m=(N+1)+1 とは符号を除いて同一の値を取り、 独立な解を与えない。

LANG:mathematica
Module[{nn}, nn = 5;
  Plot[Table[
    If[m > nn, -1, 1] Sin[n m Pi/(nn + 1)], {m, {nn, nn + 2}}] // Evaluate,
  {n, 0, nn + 1}, Mesh -> 5, MeshStyle -> Directive[PointSize[Large]]]
]

1d-chain-dispersion3.svg

これが m=1\dots N の範囲に限られる理由となる。

これらに対応する k_m のそれぞれに対して \omega=\pm\sqrt{\frac{2K}{M}}\sqrt{1-\cos ik\Delta x} の2つの \omega が対応するが、これらの線形結合で実数解を求めれば \cos(\omega t+\delta) の形に限られることになる。

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

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 とすると、 E(k) は上式のとおり cos 的な、こんなグラフになる。

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

1d-chain-dispersion2.svg

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

固有ベクトルの形は上で見たように、 R_0 サイトと R_{N+1} サイトを固定端とする、正弦波的な形になる。

c_n^m\propto \sin(k_mR_n)

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

周期境界条件

上記の例では固定端となっていたが、これを周期的境界条件にしてみると、 分散関係は固定端の場合と同じだが、周期境界条件で許される波数は

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)

が縮退する。

したがって、固有値問題

\begin{pmatrix} 0&-1&&&\red{-1}\\ -1&0&-1\\ &-1&0&\ddots\\ &&\ddots&\ddots&-1\\ \red{-1}&&&-1&0 \end{pmatrix} \begin{pmatrix} c_1\\c_2\\\vdots\\\vdots\\c_N \end{pmatrix} = \frac{E+\alpha}{\gamma} \begin{pmatrix} c_1\\c_2\\\vdots\\\vdots\\c_N \end{pmatrix}

に対する固有値を小さいものから並べれば、 (赤で示した \red{-1} が周期性を担っている)

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

このように2つずつが等しいエネルギーを持つことになる。

ただし基底状態は定数関数なので sin, cos の区別が無く、縮退しない。

また上のように原子数が偶数の場合、エネルギー最大の状態は隣同士で符号が入れ替わる形になるため、やはり sin, cos の区別はなくなり、縮退しない。

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

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

cyclic-boundary_waves.png

となる。

状態密度

今の場合 k_m は等間隔に並ぶため、状態密度は \frac{\PD k}{\PD E} に比例する。 ここでは 1/\sin k\Delta R となって、

1d-chain_dos.png

cos の上下端で状態密度が発散するのは1次元結晶の特徴である。

1次元状態密度

k 空間で状態が等間隔に存在しているとき、 つまり例えば k=n\Delta k に状態があるとき、 k 空間での状態密度(単位波数当たりの状態密度)は一定値 1/\Delta k である。

分散関係 E(k) が単調であれば、 E 空間での状態密度は

\rho(E)dE=\frac{dk}{\Delta k}

より、

\rho(E)=\frac{1}{\Delta k}\left.\frac{\PD k}{\PD E}\right|_E

である。

したがって、上のように \frac{dE}{dk}=0 となる点においては、 一次元状態密度は発散する。

一般には E=E(k) となるような k は複数あるので、

\rho(E)=\frac{1}{\Delta k}\sum_n\left.\frac{\PD k}{\PD E}\right|_{k=k_n}

としなければならない。

この式は、形式的には次のように書くこともできる。

\rho(E)=\frac{1}{\Delta k}\int \delta(E(k)-E)\,dk

\frac{\PD k}{\PD E} の因子がなくなったのは、 それが暗に \delta(E(k)-E) に含まれているためだ。

つまり、 E(k_0)=E_0 であるとき、その周辺で

\delta(E(k_0)-E_0)=\delta(k-k_0)\left.\frac{\PD E}{\PD k}\right|_{k=k_0}

が成り立つ。すなわち、 k=k_1,k_2,k_3,\dots E(k)=k となるとき、

\rho(E)=\frac{1}{\Delta k}\int \delta(E(k)-E)\,dk=\frac{1}{\Delta k}\sum_n\left.\frac{\PD k}{\PD E}\right|_{k=k_n}

となる。 E で書かれたデルタ関数を k で積分すると \frac{\PD E}{\PD k} を集めて回ることができるのである。

運動量

\Phi_k(x) &=\sum_{n=1}^N e^{ikn\Delta x}\phi(x-n\Delta x)\\ &=e^{ikx}\underbrace{\sum_{n=1}^N e^{-ik(x-n\Delta x)}\phi(x-n\Delta x)}_{u_k(x)}\\

で表されるブロッホ状態に対して、速度の期待値を求めてみる。*1http://www1.doshisha.ac.jp/~ma... を参考にした

まず、

\hat p\Phi_k &= \frac{\hbar}{i}\frac{\PD}{\PD x}e^{ikx}u_k(x)\\ &= \hbar ke^{ikx}u_k(x)+e^{ikx}\frac{\hbar}{i}\frac{\PD}{\PD x}u_k(x)\\ &= \hbar ke^{ikx}u_k(x)+e^{ikx}\hat p u_k(x)\\ &= e^{ikx}(\hbar k+\hat p) u_k(x)\\

\hat p^2\Phi_k &= \hat p\,e^{ikx}(\hbar k+\hat p) u_k(x)\\ &= \hbar k e^{ikx}(\hbar k+\hat p) u_k(x)+e^{ikx}(\hbar k+\hat p) \hat p\,u_k(x)\\ &= e^{ikx}(\hbar k+\hat p)^2 u_k(x)\\

を用いてシュレーディンガー方程式を書き直すと、

\left(\frac{\hat p^2}{2m}+V\right)\Phi_k=\varepsilon(k)\Phi_k

e^{ikx}\left\{\frac{(\hbar k+\hat p)^2}{2m}+V\right\}u_k=e^{ikx}\varepsilon(k)u_k

すなわち、 u_k に対するシュレーディンガー方程式が次のように得られる。

\underbrace{\left\{\frac{(\hbar k+\hat p)^2}{2m}+V\right\}}_{\hat H(k)}u_k=\varepsilon(k)u_k

ここに現れるハミルトニアンは k に依存するが、 このような場合 ファインマンの定理 により、

\int u_k^*\frac{\PD \hat H}{\PD k}u_k\,dx=\frac{\PD \varepsilon}{\PD k}

が成り立つ。ここでは \PD \hat H/\PD k=\hbar (\hat p+\hbar k)/m であり、

\int u_k^*\frac{\hbar}{m}(\hat p+\hbar k)u_k\,dx=\frac{\PD \varepsilon}{\PD k}

を得る。一方で、

\langle v\rangle &=\int\Phi_k^*\frac{\hat p}{m}\Phi_k\,dx\\ &=\int \cancel{e^{-ikx}}u_k^*\ \cancel{e^{ikx}}\frac{(\hat p+\hbar k)}{m}u_k\,dx\\ &=\int u_k^*\frac{(\hat p+\hbar k)}{m}u_k\,dx\\

より、

\langle v\rangle=\frac{1}{\hbar}\frac{\PD \varepsilon}{\PD k}

を得る。

これと正準方程式 \dot x=\frac{\PD H}{\PD p} とを比べると、 \hbar k x に対応する正準運動量になっていそうな気がしてくる。

ただ、ここまで良く見返してみると、上の議論では V u_k の周期性など物理的な条件を一切使っていない。

それでは、任意のハミルトニアンの任意の固有関数に対して、 任意に決めた k を使って同様の議論ができてしまうのではないか?

おかしな波数をとってみる

以下、ちょっと納得のいかない結果になってしまうため、慎重に。

あるポテンシャル V に対する時間を含まないシュレーディンガー方程式の エネルギー固有値及び固有関数が、 ほとんど連続な量子数 k を用いて \varepsilon(k), \Phi_k(x) と書けるとする。

ポテンシャルの周期性や、波動関数がブロッホ状態であることなどは仮定していないことに注意せよ。

k を別の量子数 k' を使って k=k(k') と表したうえで、

\Phi_k(x) &=e^{ik'x}\underbrace{e^{-ik'x}\Phi_k(x)}_{u_{k'}(x)}\\

のように無理やり e^{ik'x} をくくりだした形から始めて、速度の期待値を求める。

この関数に \hat p を作用させると、上記と同様にして \hbar k'+\hat p が間に挟まる形を導ける。

\hat p\Phi_k &= \frac{\hbar}{i}\frac{\PD}{\PD x}e^{ik'x}u_{k'}(x)\\ &= \hbar k'e^{ik'x}u_{k'}(x)+e^{ik'x}\frac{\hbar}{i}\frac{\PD}{\PD x}u_{k'}(x)\\ &= e^{ik'x}(\hbar k'+\hat p) u_{k'}(x)\\

\hat p^2\Phi_k &= \hat p\,e^{ik'x}(\hbar k'+\hat p) u_{k'}(x)\\ &= \hbar k' e^{ik'x}(\hbar k'+\hat p) u_{k'}(x)+e^{ik'x}(\hbar k'+\hat p) \hat p\,u_{k'}(x)\\ &= e^{ik'x}(\hbar k'+\hat p)^2 u_{k'}(x)\\

結局これは、

\hat p\,f(x)&=\hat p\,e^{ikx}e^{-ikx}f(x)\\ &=e^{ikx}(\hbar k +\hat p)e^{-ikx}f(x)\\ &=e^{ikx}e^{-ikx}(\hbar k-\hbar k+\hat p)f(x)\\ &=\hat p\,f(x)\\

という恒等式を表しているに過ぎない。

上記を用いてシュレーディンガー方程式を書き直す。

\left(\frac{\hat p^2}{2m}+V\right)\Phi_k=\varepsilon(k(k'))\Phi_k

e^{ik'x}\left\{\frac{(\hbar {k'}+\hat p)^2}{2m}+V\right\}u_{k'}=e^{ik'x}\varepsilon(k(k'))u_{k'}

すなわち、 u_{k'} に対するシュレーディンガー方程式が次のように得られる。

\underbrace{\left\{\frac{(\hbar k'+\hat p)^2}{2m}+V\right\}}_{\hat H(k')}u_{k'}=\varepsilon(k(k'))u_{k'}

ここに現れるハミルトニアンは k' に依存するので、 「ファインマンの定理が成り立つならば」

\int u_{k'}^*\frac{\PD \hat H}{\PD k'}u_{k'}\,dx=\frac{\PD \varepsilon}{\PD k'}

となる。

で、あとは一直線。 \PD \hat H/\PD k'=\hbar (\hat p+\hbar k')/m より、

\int u_{k'}^*\frac{\hbar}{m}(\hat p+\hbar k')u_{k'}\,dx=\frac{\PD \varepsilon}{\PD k'}

一方で、

\langle v\rangle &=\int\Phi_k^*\frac{\hat p}{m}\Phi_k\,dx\\ &=\int \cancel{e^{-ik'x}}u_{k'}^*\ \cancel{e^{ik'x}}\frac{(\hat p+\hbar k')}{m}u_{k'}\,dx\\ &=\int u_{k'}^*\frac{(\hat p+\hbar k')}{m}u_{k'}\,dx\\

より、

\langle v\rangle &=\frac{1}{\hbar}\frac{\PD \varepsilon}{\PD {k'}}\\ &=\frac{1}{\hbar}\frac{\PD \varepsilon}{\PD k}\frac{\PD k}{\PD k'}

となって、右辺は k' の取り方によって値が変わってしまう。

どうしてこうなった?

怪しいのはファインマンの定理のところだけのように思われる。

英語版の Wikipedia にも書かれているように、 ファインマンの定理が成立するためには関数空間やハミルトニアンに対して一定の条件が必要になる。 Carfì, David (2010). "The pointwise Hellmann–Feynman theorem". AAPP: Physical, Mathematical, and Natural Sciences. 88 (1). no. C1A1001004. doi:10.1478/C1A1001004

たぶん、正しい k を取ったときのみ上記の演算が許される?

ファインマンの定理再考

正しく k を取った場合、

\Phi_k(x) &=e^{ikx}\underbrace{\sum_{n=1}^N e^{-ik(x-n\Delta x)}\phi(x-n\Delta x)}_{u_k(x)}\\

この u_k は結晶の端を除いて周期 \Delta x の周期関数になっている。

u_k k で微分すると

\frac{\PD}{\PD k}u_k(x) &=-ie^{ikx}\sum_{n=1}^N (x-n\Delta x)e^{-ik(x-n\Delta x)}\phi(x-n\Delta x)

となって、これも周期 \Delta x の周期関数。

おそらくこの、「周期関数になっている」というのが重要な条件になっているのだと思われる。

で気づいたのが、そもそも周期的境界条件以外の束縛状態の定常状態を考える限り、 運動量あるいは速度の期待値が非ゼロになるような関数は解空間に含まれていないことだ。

箱型ポテンシャルや調和振動子では波数 k の成分は必ず 波数 -k の成分とある決まった係数で線形結合を取り、 全体での運動量の期待値をゼロにしないと境界条件を満たさない。

定常的に非ゼロの運動量を持った解は束縛されていないのだからこれは当然。

ファインマンの定理では u_k \PD u_k/\PD k を含む 関数ベクトル空間でハミルトニアンがエルミートであることを仮定しているが、 適当な k' を取ってしまうとその空間でハミルトニアンがエルミートでは なくなってしまうのが問題。*2演算子のエルミート性については 量子力学Ⅰ/固有値と期待値#qf83b503(http://dora.bk.tsukuba.ac.jp/~...) を参照のこと

周期的境界条件では領域の端でちょうど関数値が等しくなることから 内積を部分積分した際に現れる積分値がゼロになり、 運動量演算子のエルミート性が保たれる。

ただ実は、 \Phi_k に対して u_k を周期関数にするには 必ずしも k'=k としなければならないわけではなくて、 k'=k+G となっていても良い。ただしこの G は結晶格子の逆格子ベクトル。

このとき \PD k'/\PD k=1 なので、 k を使っても、 k' を使っても、速度の期待値は変わらない。

これが結晶の分散関係がブリルアンゾーンをまたいで周期的に繰り返す理由になっている。

今の場合

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

\frac{\PD E}{\PD k}=2\Delta x\gamma\sin(k\Delta x)

より、

\langle v\rangle_k=\frac{2\gamma\Delta x}{\hbar}\sin(k\Delta x)

確かに、 k\to k+G としても \braket{v} の値は変わらない形になっている。

LANG:Mathematica
Plot[Sin[Pi x], {x, -1, 1}]

sine-wave.svg

k\sim 0 では \langle v\rangle\propto k だが、 ある程度 k が大きくなると速度は徐々に上がらなくなり、 k\Delta x=\pi/2 で最大値 2\gamma\Delta x/\hbar を取り、 それ以上では k が増えると \langle v\rangle は減少し、 k\Delta x=\pi で再び静止する。

有効質量

外力 f が短時間 \Delta t 加わって、 速度の期待値が \Delta v 増えたとする。 このとき粒子は \braket v \Delta t 進んでいるので、 外力から得たエネルギーは \Delta E=\braket v \Delta t f である。

この過程で波数が \Delta k 増えたとすると、

\Delta E=\frac{\PD E}{\PD k}\Delta k=\braket v \Delta t f

ここに上で得た \braket v の表式を代入すれば、

\frac{\PD E}{\PD k}\Delta k=\frac{1}{\hbar}\frac{\PD E}{\PD k} \Delta t f

より、

\frac{dk}{dt}=f/\hbar

を得る。

すなわち、一定の外力の下では k は時間とともに線形に増加する。

一方、このとき粒子の加速度は、

\frac{d}{dt}\braket v=\frac{d\braket v}{dk}\frac{dk}{dt}=\frac{1}{\hbar^2}\frac{\PD E}{\PD k^2}f

より、ニュートン方程式 f=m\frac{d}{dt}\braket v と比較すると、 m^* を有効質量として

\frac{1}{m^*}=\frac{1}{\hbar^2}\frac{\PD E}{\PD k^2}

の関係が得られる。

エネルギー対波数の分散関係を基準とすれば、
その傾きが速度であり、
波数の時間変化が力そのものであるため、
分散関係の二次微分が質量の逆数に対応する。

二次微分が大きければ同じ力に対する速度の増加率が大きい = 質量は小さい
二次微分が小さければ同じ力に対する速度の増加率が小さい = 質量は大きい

古典論との対応を見ておこう。

ハミルトニアンが運動エネルギーと保存力によるポテンシャルのみを含んでいる場合、 古典論において正準方程式から

\dot p=-\frac{\PD H}{\PD x}=-\frac{\PD V}{\PD x}=f

\dot x=\frac{\PD H}{\PD p}

\ddot x&=\frac{d}{dt}\frac{\PD H}{\PD p}\\ &=\frac{\PD}{\PD p}\frac{dH}{dt}\\ &=\frac{\PD}{\PD p}\frac{\PD H}{\PD p}\dot p\\ &=\frac{\PD^2 H}{\PD p^2}f\\

のようにして、

\frac{1}{m}=\frac{\PD^2 H}{\PD p^2}

の関係を導けることに対応している。

3次元空間での一般形は、

\hbar\dot{\bm k}=-\bm\nabla_x H=-\bm\nabla_x V=\bm f

\dot{\bm x}=\frac{1}{\hbar}\bm \nabla_k H

\ddot{\bm x}&=\hbar\frac{d}{dt}\bm\nabla_k H\\ &=\frac{1}{\hbar}\bm\nabla_k \frac{dH}{dt}\\ &=\frac{1}{\hbar^2}\bm\nabla_k (\bm\nabla_k H\cdot \hbar\dot{\bm k})\\ &=\frac{1}{\hbar^2}\Bigg(\frac{\PD^2 H}{\PD k_i\PD k_j}\Bigg) \hbar\dot{\bm k}\\

のように、逆有効質量はテンソルになる。

今の場合

m^*=\frac{\hbar^2}{2\Delta R^2\gamma\cos(k\Delta R)}

LANG:mathematica
Plot[1/Cos[Pi k], {k, 0, 1}]

effective-mass1.svg

k\sim 0 では \cos(k\Delta R)\sim 1 より m^*=\frac{\hbar^2}{2\Delta R^2\gamma} となる。

一方、 k\sim \pi/\Delta R では \cos(k\Delta R)\sim -1 より m^*=-\frac{\hbar^2}{2\Delta R^2\gamma} のように負の質量が現れる。

負の有効質量

負の有効質量をもつ電子は、電場と同じ方向へ加速される。というか、もともと電場と逆方向に走っているので電場と逆方向の速度が減速される、というの方が分かりやすいか。通常の電子は電場と逆の方向へ加速されることに注意。

ここからどのような結論が導かれるだろうか?

電子が一つのとき

バンドの下端に電子が1つある状況から、外場 f によりこの電子を加速する。

散乱がなければ電子の波数は時間とともに一様に増加して、

\frac{dk}{dt}=\frac{f}{\hbar} より、 k=\frac{f}{\hbar}t

速度および位置は、

\braket v=\frac{2\gamma\Delta x}{\hbar}\sin(\Delta xft/\hbar)

\braket x=\int\braket v dt=\frac{2\gamma}{f}\cos(\Delta xft/\hbar)

のように振動する。

力が強くなれば速度変化が速くなり、振幅が小さくなる。

この振動はブロッホ振動と呼ばれる。

外場がゼロになると \frac{dk}{dt}=0 なので、等速直線運動を続けることになる。

実際には散乱があるので徐々に減衰して際安定状態に落ち込み、静止する。

散乱が強ければ電子が最高速度を超えることはなく、 平均的な速度が外場に比例するオーミックな電流が流れる。

電子が底の方に詰まっているとき

最安定状態では k 側と -k 側に同じだけ電子がいるため、 電流は相殺されてゼロになる。

外場 f がかかると、 すべての電子の k \frac{dk}{dt}=\frac{f}{\hbar} で変化して、 -k 側の電子が減り、 k 側の電子が増える。

その差の分だけ電流が流れる。

電子が半分以上詰まっているとき

外場 f がかかると、 すべての電子の k \frac{dk}{dt}=\frac{f}{\hbar} で変化して、 -k 側の電子が減り、 k 側の電子が増える。

k\Delta x>\pi/2 の電子は負の有効質量のために正方向の速度が減速され、
k\Delta x<-\pi/2 の電子は負の有効質量のために負方向の速度が加速されるが、
全体として見るとやはり負方向へ動く電子が減って、正方向へ動く電子が増えるため、 電流は力の加わった方向へ流れることになる。

詰まり具合と電流応答

電子が底の方に少しだけ詰まっている時、 外場によって k が少しだけ変わっても、 正側、負側とも速度が遅いので、電流の応答は小さい、と思われる。

電子が半分まで詰まっている時、 外場によって k が少し変わるだけでも、 最高速で負側に動く電子が減って、最高速で正側に動く電子が増えるので、 電流の応答は大きい、と思われる。

電子がかなり上まで詰まっている時、 やはり正側、負側とも速度が遅いので、電流の応答は小さい、と思われる。

1個だけ残して残りが埋まっているとき

外場 f がかかると、 すべての電子の k \frac{dk}{dt}=\frac{f}{\hbar} で変化するため、 空いた電子状態の k も同じように変化する。

この電子状態は負の有効質量をもっているため、電子に加わる力とは反対方向に加速される。

電子の空状態なので、相対的に正の電荷を帯びている。

この電子の空準位がホールと呼ばれる。

2次元三角格子

1原子状態をつなげてブロッホ波的な三角格子を作ると、これが実際の三角格子の解になっていることを確かめる。


triangular-lattice.png


N 個の原子からなる三角格子中の原子位置を \bm R_A で表し、 これらの原子を中心とする1原子状態を元にブロッホ状態を作る。

|\phi_A\rangle=\frac{1}{\sqrt{N}}\sum_{\bm R_A}e^{i\bm k\cdot\bm R_A}|\bm R_A\rangle

3角格子では最隣接原子が6個あり、基本ベクトルを \bm a,\bm b と置けば

\bm a, -\bm a, \bm b, -\bm b,-\bm a+\bm b, \bm a-\bm b

に位置する。したがって係数間の方程式は

-\alpha C(\bm R_i)-\gamma\Big[C(\bm R_i+\bm a)+C(\bm R_i-\bm a)+C(\bm R_i+\bm b)+C(\bm R_i-\bm b)+C(\bm R_i-\bm a+\bm b)+C(\bm R_i+\bm a-\bm b)\Big]=EC(\bm R_i)

であり、 C(\bm R)=C(\bm 0)e^{i\bm k\cdot\bm R} と置いてみるとこれは、

&-\alpha-\gamma\Big[e^{i\bm k\cdot\bm a}+e^{-i\bm k\cdot\bm a}+e^{i\bm k\cdot\bm b}+e^{-i\bm k\cdot\bm b}+e^{i\bm k\cdot(-\bm a+\bm b)}+e^{-i\bm k\cdot(-\bm a+\bm b)}\Big]=\\&-\alpha-2\gamma\Big[\cos \bm k\cdot\bm a+\cos \bm k\cdot\bm b+\cos \bm k\cdot(\bm a-\bm b)\Big]=E\\

を満たす固有値 E の固有ベクトルとなることが分かる。

今の場合 \bm a=\begin{pmatrix}a\\0\end{pmatrix} \bm b=\begin{pmatrix}a/2\\\sqrt{3}\,a/2\end{pmatrix} \bm a-\bm b=\begin{pmatrix}-a/2\\\sqrt{3}\,a/2\end{pmatrix} であるから、その分散関係は次のようになる。

honeycomb[k_] =
  Cos[k.{1, 0}] + Cos[k.{1/2, Sqrt[3]/2}] + Cos[k.{-1/2, Sqrt[3]/2}]

triangular-band1.jpg triangular-band2.jpg triangular-band3.jpg

これらは、 \bm k=(k_x,k_y) に対して、対応する E(\bm k) をプロットしたものである。

逆格子

実空間での基本ベクトルが

\bm a=(a,0)

\bm b=(a/2,\sqrt 3a/2)

であるのに対応して、逆格子空間における基本ベクトルは、

\bm a^*=\frac{2\pi}{a}(1,-1/\sqrt 3)

\bm b^*=\frac{2\pi}{a}(0,2/\sqrt 3)

となる。

\bm a\cdot \bm a^*=\bm b\cdot \bm b^*=2\pi

\bm a\cdot \bm b^*=\bm b\cdot \bm a^*=0

に注意せよ。 A

triangular-brillouin-zone.png

有効質量

2次元なので、方向によって有効質量が異なる場合がある。

\frac{d\bm k}{dt}=\frac{1}{\hbar}\bm f

\braket{\bm v}=\frac{1}{\hbar}\bm \nabla_k H

\frac{d\braket{\bm v}}{dt}=\frac{1}{\hbar^2}\bm \nabla_k(\bm \nabla_k H\cdot \bm f)

より、

\frac{1}{m^*}=\frac{1}{\hbar^2}\Bigg(\frac{\PD^2 H}{\PD k_i\PD k_j}\Bigg)

三角格子の場合

上で計算した分散関係を再掲すると、

triangular-band1.jpg triangular-band4.png

ガンマ点、K'点、M点で速度ゼロ。その他の点では傾きに応じた方向の速度を持つ。

ガンマ点でエネルギー最大、曲率および有効質量は等方的で負の値を持つ。

K' 点でエネルギー最小、曲率および有効質量は等方的で正の値を持つ。

M点は鞍点で、曲率および有効質量は一方向に正、一方向に負の値を持つ。

原子の電子状態の等方性

この手の議論では S 状態を出発点に仮定している場合が多い。

にもかかわらず、例えばグラフェンを扱う際にはπ電子を考えていたりする。

このような仮定は電子状態の等方性が問題にならないようにするためである。

電子状態が強い異方性を持っていると、 \alpha,\gamma の値が向きによって異なってしまう。

π電子は3次元空間で見ると非等方的だが、グラフェンシート内では等方的であるため、 その非等方性は問題にならない。

グラフェンの場合

三角格子の単位胞に2つの原子 A, B がある構造。図ではそれぞれ白と赤。

graphene-lattice.png

系全体の波動関数を、A に対する三角格子波動関数と、B に対する三角格子波動巻子との線形結合で

\ket \Phi=C_A\ket{\Phi_A}+C_B\ket{\Phi_B}

と置き、 \braket{\Phi_A|\Phi_B}=0 を仮定すると、

\braket{\Phi_A|\hat H|\Phi}=E\braket{\Phi_A|\Phi} より、

C_A\underbrace{\braket{\Phi_A|\hat H|\Phi_A}}_{h_{AA}}+ C_B\underbrace{\braket{\Phi_A|\hat H|\Phi_B}}_{h_{AB}}=EC_A

などとなり、

\begin{pmatrix} h_{AA}&h_{AB}\\ h_{BA}&h_{BB}\\ \end{pmatrix} \begin{pmatrix}C_A\\C_B\end{pmatrix} =E\begin{pmatrix}C_A\\C_B\end{pmatrix}

を得る。

h_{AA},h_{BB} は等しく、三角格子のエネルギーと隣接原子からのクーロン相互作用とを加えたものになる。前者は上で解いたように波数に依存して変化するが、後者はブロッホ波の位相によらず、状態密度のみに依存する量であるため定数である。

h_{AA}=E_\mathrm{tri}(\bm k)-\tilde\alpha

h_{AB}=h_{BA}^* であり、 隣り合う原子上の1原子波動関数の間に片方の原子核からのポテンシャルを挟んだものを 原子数分足したものになる。

h_{AB}&=\braket{\Phi_A|H|\Phi_B}\\ &\frac{1}{N}\sim \sum_{R_A}\sum_{R_{B}} \exp\big[-i\bm k\cdot(\bm R_A-\bm R_B)\big] \underbrace{\braket{\bm R_A|H|\bm R_B}}_{-\gamma}\\ &= -\frac{1}{N}\gamma\sum_{R_A}\sum_{R_{B(隣接)}} \exp\big[i\bm k\cdot(\bm R_B-\bm R_A)\big]\\ &= -\frac{1}{N}\gamma\sum_{R_A} \Big( \exp\big[i\bm k\cdot\Delta\bm R_{BA,1}\big]+ \exp\big[i\bm k\cdot\Delta\bm R_{BA,2}\big]+ \exp\big[i\bm k\cdot\Delta\bm R_{BA,3}\big] \Big)\\ &= -\gamma \Big( \exp\big[i\bm k\cdot\Delta\bm R_{BA,1}\big]+ \exp\big[i\bm k\cdot\Delta\bm R_{BA,2}\big]+ \exp\big[i\bm k\cdot\Delta\bm R_{BA,3}\big] \Big)\equiv -\gamma f_k^*\\

2行目が \sim で始まっているのは、 境界上の原子では周囲に3原子存在しない場合があり、 そのあたりの境界条件がうやむやになっているため。

上図の白丸を A サイトとすると、隣接 B サイトは

\Delta\bm R_{BA,1}=a(0,1/\sqrt 3)

\Delta\bm R_{BA,2}=a(-1/2,-1/2\sqrt 3)

\Delta\bm R_{BA,2}=a(1/2,-1/2\sqrt 3)

となるから、これらを代入すれば \bm k の関数として f_{\bm k}^* が求まる。

このとき、

\begin{pmatrix} E_\mathrm{tri}(k)-\tilde\alpha&-\gamma f_{\bm k}^*\\ -\gamma f_{\bm k}&E_\mathrm{tri}(k)-\tilde\alpha\\ \end{pmatrix} \begin{pmatrix}C_A\\C_B\end{pmatrix} =E\begin{pmatrix}C_A\\C_B\end{pmatrix}

\begin{pmatrix} 0&-\gamma f_{\bm k}^*\\ -\gamma f_{\bm k}&0\\ \end{pmatrix} \begin{pmatrix}C_A\\C_B\end{pmatrix} =\big(E-E_\mathrm{tri}(k)+\tilde\alpha\big)\begin{pmatrix}C_A\\C_B\end{pmatrix}

より、

E=\pm\gamma |f_k|+E_\mathrm{tri}(k)-\tilde\alpha

-\gamma f_k C_A=E C_B=\big(\pm\gamma N|f_k|+E_\mathrm{tri}(k)-\tilde\alpha\big)C_B

より、 f_k=|f_k|e^{i\phi_k} と置けば、

\begin{pmatrix}C_A\\C_B\end{pmatrix}\propto \begin{pmatrix} \pm 1+(E_\mathrm{tri}(k)-\tilde\alpha)/\gamma |f_k|\\ -e^{i\phi_k} \end{pmatrix}

分散関係

グラフェンの特異的な分散関係を見るために、まずは E=\pm\gamma|f_k| のグラフを見る。

LANG: mathematica
Fk[k_] =
  Exp[I k.{0, 1/Sqrt[3]}] +
  Exp[I k.{-1/2, -1/(2 Sqrt[3])}] +
  Exp[I k.{1/2, -1/(2 Sqrt[3])}]
Plot3D[{Abs[Fk[{x, y}]], -Abs[Fk[{x , y}]]},
  {x, -10, 10}, {y, -10, 10}, PlotPoints -> 200]

graphene-dispersion1.jpg graphene-dispersion2.jpg graphene-dispersion3.jpg graphene-dispersion4.jpg

分散関係はゼロを境に上下に対象であり、「卵のパック」のような形になる。

大きく膨らんでいる場所がガンマ点、 鋭くへこんでいる場所がK点、 鞍点となっているのがM点である。

中でもK点では、

f_k=\sum_{l=1}^3 e^{i\bm k\cdot \Delta R_l}

に対して \bm k=\frac{2\pi}{a}(2/3,0) より、

f_{\bm k}&=\exp(0)+\exp(i2\pi/3)+\exp(-i2\pi/3)\\ &=1+w+w^{-1}\\ &=0

ただし、 w=e^{i2\pi/3} は 1 の3乗根。

となって、 E_+(\bm k)=E_-(\bm k)=0 すなわちギャップは閉じている。

状態密度

このとき、状態密度をモンテカルロ法により求めると、

LANG: mathematica
points = Table[{
     Abs[Fk[{RandomReal[{-100, 100}], RandomReal[{-100, 100}]}]],
    -Abs[Fk[{RandomReal[{-100, 100}], RandomReal[{-100, 100}]}]]},
     {i, 1, 1000000}] // Flatten;
Histogram[points, {-4, 4, 0.05}]

graphene-dos.png

ここでは \gamma=1 に取っている。

E=\pm 3\gamma では状態密度はステップ関数的に立ち上がる。
これはΓ点のパラボリックな分散関係を反映している。

E=\pm \gamma では発散する。
これはM点に存在する鞍点を反映している。

E=0 では線形に立ち上がる。
これはK点、K’点に存在する円錐状の分散関係を反映している。

以下、少し詳しく見てみる。

2次元状態密度

1次元の時と同様に、2次元状態密度は形式的に

\rho(E)=\frac{1}{\Delta k_x\Delta k_y}\iint \delta(E(\bm k)-E)\,d\bm k

と書ける。この形は解析的に考えることができる場合には有用であるが、 数値計算可能な形ではないので、計算できる形に変形したい。

\rho(E)dE=\frac{dk_x}{\Delta k_x}\frac{dk_y}{\Delta k_y}

と書ければ良いのだが、 右辺の範囲が必ずしも E\sim E+dE と一致しないため難しい。

x,y 方向ではなく、 \bm e_1=\frac{\bm\nabla_k E}{|\bm\nabla_k E|} 方向と、それに垂直な \bm e_2 方向に軸を取れば、

\delta\rho(E)dE=\oint\frac{dk_2}{\Delta k_2}\,\frac{dk_1}{\Delta k_1}

ただし dE=|\nabla_k E|dk_1 。したがって、

\rho(E)=\frac{1}{\Delta k_1\Delta k_2}\oint |\nabla_k E|^{-1}dk_2

ここで、 \oint dk_2 E(\bm k)=E となる等高線に沿って計算する。

とすると、格段にイメージはしやすくなるが、 これも数値計算可能な形は限られている。。。

ある点 $\bm k_0$ を中心に回転対称性があるとき

\delta \bm k\equiv \bm k-\bm k_0 を用いて、 E(\bm k)=E(|\bm k-\bm k_0|)=E(\delta k) と書けるから、

\bm\nabla_{\bm k} E は円周方向を向くノルム一定のベクトルで、 \int dk_2 は円周に沿った長さ 2\pi r ( r は半径) に渡る積分になり、

\rho(E)&=\frac{1}{\Delta k_x\Delta k_y}\int 2\pi\delta k\,\delta(E(\delta k)-E)\,d\delta k\\&=\frac{1}{\Delta k_x\Delta k_y}2\pi\delta k\left.\frac{\PD \delta k}{\PD E}\right|_E\\

となる。


(1) ガンマ点の周り