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

更新


公開メモ

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

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

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

である。

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 とすると、

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

である。

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

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

電子が複数ある影響

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

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

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

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

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

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

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

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 として解いて、固有値を小さいものから順にプロットすると、 ( 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

を行列で表せば、

&math( \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} );

したがって、

&math( \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} );

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

&math( \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)

が縮退する。

したがって、固有値問題

&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

このように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} を集めて回ることができるのである。

運動量

&math( \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... を参考にした

まず、

&math( \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)\\ );

&math( \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)\\ );

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

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

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

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

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

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

&math( \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 であり、

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

を得る。一方で、

&math( \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\\ );

より、

&math( \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') と表したうえで、

&math( \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 が間に挟まる形を導ける。

&math( \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)\\ );

&math( \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)\\ );

結局これは、

&math( \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)\\ );

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

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

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

&math( 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'} に対するシュレーディンガー方程式が次のように得られる。

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

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

&math( \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 より、

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

一方で、

&math( \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\\ );

より、

&math( \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 で書いた結果を信用すれば、

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

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

より、

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

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 で再び静止する。

有効質量

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

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

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

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

の関係を導ける。

同様に、上で得た式から

\langle v\rangle=\frac{\PD E}{\PD \hbar k}

\frac{d}{dt}\langle v\rangle=\frac{\PD }{\PD \hbar k}\frac{d}{dt}E

\frac{d}{dt}\langle v\rangle=\frac{\PD }{\PD \hbar k}\frac{\PD E}{\PD \hbar k}\frac{d\hbar k}{dt}

として、

\frac{1}{m^*}=\frac{\PD^2 H}{\PD \hbar k^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

&math(\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} のように負の質量が現れる。

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つの原子 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} より、

&math( C_A\underbrace{\braket{\Phi_A|\hat H|\Phi_A}}_{h_{AA}}+ C_BA\underbrace{\braket{\Phi_A|\hat H|\Phi_B}}_{h_{AB}}=EC_A);

などとなり、

&math( \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原子波動関数の間に片方の原子核からのポテンシャルを挟んだものを 原子数分足したものになる。

&math( h_{AB}&=\braket{\Phi_A|H|\Phi_B}\\ &\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}\\ &= -\gamma\sum_{R_A}\sum_{R_{B(隣接)}} \exp\big[i\bm k\cdot(\bm R_B-\bm R_A)\big]\\ &= -\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 N \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)

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

このとき、

&math( \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} );

&math( \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|f_k|+E_\mathrm{tri}(k)-\tilde\alpha

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

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

&math( \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|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) より、

&math( 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) ガンマ点の周りの円錐型の部分では

\frac{\PD \delta k}{\PD E}=\mathit{const.}

であるから、

\rho(E)\propto\frac{1}{\Delta k_x\Delta k_y}2\pi\delta k

のように \delta k に対して1次で立ち上がる。

さらに

\frac{\PD \delta k}{\PD E}=\mathit{const.}

であるから、

\rho(E)\propto E

となる。


(2)ガンマ点の周りの回転放物面的な部分では [#x6fa440b]

E\propto \delta k^2

より、

\frac{\PD \delta k}{\PD E}\propto \frac{1}{\delta k}

であるから、

\rho(E)\propto\frac{2\pi}{\Delta k_x\Delta k_y}

のように定数値になる。

鞍点では

E-E_0=x^2-y^2 の形の理想的な鞍点において、 特定のエネルギーでの切断面は

  • E=E_0 では2つの直線が交わる形
  • E\ne E_0 では双曲線
LANG:mathematica
Plot3D[x^2 - y^2, {x, -1, 1}, {y, -1, 1}, MeshFunctions -> {#3 &}, 
ImageSize -> Large, PlotPoints -> 101]

suddle-point.jpg

\bm \nabla E=\begin{pmatrix}\PD E/\PD x\\\PD E/\PD y\end{pmatrix}=\begin{pmatrix}2x\\-2y\end{pmatrix} だから、

E=E_0 の直線上 y=\pm x の上の 原点から距離 t の点 \frac{1}{\sqrt{2}}\begin{pmatrix}t\\\pm t\end{pmatrix} において \bm \nabla E=\begin{pmatrix}\sqrt{2}t\\\mp\sqrt{2}t\end{pmatrix} となって、

\rho(E_0)\propto\int_{-\infty}^{\infty}|\nabla E|^{-1}dt=\int_{-\infty}^{\infty}\frac{dt}{2|t|}=+\infty

のように発散する。

モンテカルロじゃなくても計算できるか?

電子状態密度を上記のデルタ関数を用いた形で Mathematica に与えても、

LANG:mathematica
rho[e_] := 
  Integrate[Integrate[
    DiracDelta[Abs[Fk[{x, y}]] - e], 
    {x, 0, y/Sqrt[3]}], 
    {y, 0, 2 Pi/Sqrt[3]}
  ]

Mathematica はうまく扱えないみたいだった。

三角格子の分散関係を加える

正確には、上記のまさにグラフェンらしい分散関係に、単純三角格子の分散関係を足す必要がある。 ( \tilde\alpha は定数なので無視する)

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

単純三角格子の時の方が原子間距離が \sqrt{3} 倍遠いので、 単純計算ではクーロン相互作用が 1/3 になる。 原子が複数いる場合には遮蔽効果によりクーロン相互作用はより急峻に減衰するため。 実際にはさらに影響は小さいと考えられる。

実際に相互作用がどのくらい違うかよく分からないまま、 とりあえず 1/10 だけ、単純三角格子の効果を入れたところ、次のようになった。

Plot3D[{honeycomb[{x, y}]/10 + Abs[Fk[{x, y}]], 
  honeycomb[{x, y}]/10 - Abs[Fk[{x, y}]]}, {x, 2 Pi 2/3 - 1, 
  2 Pi 2/3 + 1}, {y, -1, 1}, ViewPoint -> Left, ImageSize -> Large, 
 PlotPoints -> 101, BoxRatios -> Automatic, MeshFunctions -> {#3 &}]

graphene-dispersion5.jpg

ガンマ点での高さや有効質量の絶対値が正側と負側とで異なる

graphene-dispersion6.jpg

一方、K点付近では正側・負側の非対称性は小さく見える。

ContourPlot[
  honeycomb[{x, y}]/10 + Abs[Fk[{x, y}]] - 
  honeycomb[{2 Pi 2/3, 0}]/10, 
  {x, 2 Pi 2/3 - 1, 2 Pi 2/3 + 1}, {y, -1, 1}]

ContourPlot[
  honeycomb[{x, y}]/10 - Abs[Fk[{x, y}]] - 
  honeycomb[{2 Pi 2/3, 0}]/10, 
  {x, 2 Pi 2/3 - 1, 2 Pi 2/3 + 1}, {y, -1, 1}]

graphene-cone+.jpg graphene-cone-.jpg

左が正側、右が負側

このようにプロットしてみると、コーンの裾の部分は完全な丸じゃなくて丸まった三角形になることがよく分かる。

コーンの中心部分は円錐になっている。

ContourPlot[
  honeycomb[{x, y}]/10 + Abs[Fk[{x, y}]] - 
  honeycomb[{2 Pi 2/3, 0}]/10, 
  {x, 2 Pi 2/3 - 0.1, 2 Pi 2/3 + 0.1}, {y, -0.1, 0.1}]

graphene-cone0.jpg

質量のない Dirac 方程式

電子の有効質量はバンドの湾曲に対応するため、 直線的な円錐状のバンド構造は有効質量ゼロの電子系を与える。

以下このことを見てみる。


Counter: 42026 (from 2010/06/03), today: 4, yesterday: 0