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

更新


公開メモ

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

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 を取った場合、

&math( \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 で微分すると

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

より、

&math( \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 増えたとすると、

&math( \Delta E=\frac{\PD E}{\PD k}\Delta k=\braket v \Delta t f );

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

&math( \frac{\PD E}{\PD k}\Delta k=\frac{1}{\hbar}\frac{\PD E}{\PD k} \Delta t f );

より、

&math( \frac{dk}{dt}=f/\hbar );

を得る。

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

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

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

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

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

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

負の有効質量

負の有効質量をもつ電子は、電場と同じ方向へ加速される。(通常の電子は電場と逆の方向へ加速される)

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

電子が一つのとき

バンドの下端に電子が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点は鞍点で、曲率および有効質量は一方向に正、一方向に負の値を持つ。

グラフェンの場合

三角格子の単位胞に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_B\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}\\ &\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}^* が求まる。

このとき、

&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\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} と置けば、

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

&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 方程式

K 点周りの解

理想的なグラフェンでは2枚のバンドのうち下のバンドのみが電子で埋まっていて、 ちょうどディラックコーンの先端部分に「フェルミ点」が来る。

このディラックコーン付近の波動関数を記述する方程式は質量のない Dirac 方程式になることを以下確かめる。

上記の、

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

に現れる f_k について、K 点周りの様子を知るため、 \bm k-\bm K=\delta \bm k と書くことにして、 \delta \bm k の一次の項までを残そう。

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

より、

&math( f_k &=\sum_{l=1}^3 e^{-i\bm k\cdot\Delta\bm R_l}\\ &=\sum_{l=1}^3 e^{i\bm K\cdot\Delta\bm R_l}e^{-i\delta\bm k\cdot\Delta\bm R_l}\\ &=\sum_{l=1}^3 e^{i\bm K\cdot\Delta\bm R_l}(\cancel{1}-i\delta\bm k\cdot\Delta\bm R_l+\dots)\\ &\sim-i\delta\bm k\cdot\sum_{l=1}^3 \Delta\bm R_l e^{i\bm K\cdot\Delta\bm R_l}\\ &=-i\delta\bm k\cdot\frac{\sqrt 3}{2}a\begin{pmatrix} i\\1 \end{pmatrix}\\ &=\frac{\sqrt 3}{2}a(\delta k_x-i\delta k_y) );

ただしここで a は格子定数であり、 \sum_{l=1}^3 e^{i\bm K\cdot\Delta\bm R}=0 および \sum_{l=1}^3 \Delta\bm R e^{i\bm K\cdot\Delta\bm R}=\frac{\sqrt 3}{2}a\begin{pmatrix}i\\1\end{pmatrix} の関係を用いた。

LANG:mathematica
In[1]:= tau={0,1/Sqrt[3]}
Out[1]= {0,1/Sqrt[3]}
In[2]:= r=RotationMatrix[2Pi/3]
Out[2]= {{-(1/2),-(Sqrt[3]/2)},{Sqrt[3]/2,-(1/2)}}
In[3]:= tauList={tau, r.tau, r.r.tau}
Out[3]= {{0,1/Sqrt[3]},{-1/2,-1/(2 Sqrt[3])},{1/2,-1/(2 Sqrt[3])}}
In[4]:= k=2 Pi {2/3,0}
Out[4]= {(4 \[Pi])/3,0}
In[5]:= Sum[Exp[I k.tau], {tau, tauList}] // Simplify
Out[5]= 0
In[6]:= Sum[tau Exp[I k.tau],{tau,tauList}] // Simplify
Out[6]= {(I Sqrt[3])/2,Sqrt[3]/2}
In[7]:= Sum[-tau Exp[-I k.tau], {tau, tauList}] // Simplify
Out[7]= {(I Sqrt[3])/2,-Sqrt[3]/2}

さらに、 E-E_\mathrm{tri}(k)+\tilde\alpha\to E と置きなおすと、

&math( \frac{\sqrt 3}{2}\gamma a\begin{pmatrix} 0&\delta k_x+i\delta k_y\\ \delta k_x-i\delta k_y&0\\ \end{pmatrix} \begin{pmatrix}C_A\\C_B\end{pmatrix} =E\begin{pmatrix}C_A\\C_B\end{pmatrix} );

ただし、

E=\pm\gamma |f_k|=\pm\gamma \frac{\sqrt 3}{2}a|\delta \bm k|

\delta k_x+i\delta k_y=|\delta \bm k|e^{i\phi_{\delta k}} と置けば、

&math( \begin{pmatrix}C_A\\C_B\end{pmatrix}\propto \begin{pmatrix} \pm 1\\ e^{i\phi_{\delta k}} \end{pmatrix} );

が上記方程式の解となる。

Dirac 方程式が出てくる

質量のない Dirac 方程式が出てくる気配は、 上記の方程式を次のようにパウリ行列 \sigma_x,\sigma_y,\sigma_z で書き直すことで感じられる。

あーと、どこかで i -i とが入れ替わってしまったみたい。 本来別にどちらでもよいのだけれど、一般的に扱われるパウリ行列に合わせて、 上記の方程式の i -i とを入れ替えておく。

&math( &\frac{\sqrt 3}{2}\gamma a\begin{pmatrix} 0&\delta k_x-i\delta k_y\\ \delta k_x+i\delta k_y&0\\ \end{pmatrix} \begin{pmatrix}C_A\\C_B\end{pmatrix}\\ &=\frac{\sqrt 3}{2}\gamma a\left\{ \begin{pmatrix}0&1\\1&0\end{pmatrix}\delta k_x+ \begin{pmatrix}0&-i\\i&0\end{pmatrix}\delta k_x+ \begin{pmatrix}1&0\\0&-1\end{pmatrix}\cdot 0 \right\} \begin{pmatrix}C_A\\C_B\end{pmatrix}\\ &=\frac{\sqrt 3}{2}\gamma a\big( \sigma_x\delta k_x+ \sigma_y\delta k_x+ \sigma_z\cdot 0 \big) \begin{pmatrix}C_A\\C_B\end{pmatrix}\\ &=\frac{\sqrt 3}{2}\gamma a\big(\bm\sigma\cdot\delta \bm k\big) \begin{pmatrix}C_A\\C_B\end{pmatrix}=E\begin{pmatrix}C_A\\C_B\end{pmatrix}\\ );

そこで、両辺に e^{i\delta\bm k\cdot\bm R} を掛けると、

&math( \hbar\delta \bm ke^{i\delta\bm k\cdot\bm R}=\frac{\hbar}{i}\bm\nabla_Re^{i\delta\bm k\cdot\bm R} );

より、 \delta \bm k=\frac{1}{\hbar}\hat{\bm p} におきかえられて、

&math( \frac{\sqrt 3}{2}\frac{\gamma a}{\hbar}\big(\bm\sigma\cdot\hat{\bm p}\big) \begin{pmatrix}C_A\\C_B\end{pmatrix}e^{i\delta\bm k\cdot\bm R} =E\begin{pmatrix}C_A\\C_B\end{pmatrix}e^{i\delta\bm k\cdot\bm R} );

これは自由な電子に対する一般的な Dirac 方程式、

&math( \Big(i\hbar\frac{\PD}{\PD t}+i\hbar c\sum_{i=1}^3 \alpha_i\frac{\PD}{\PD x^i}-\beta mc^2\Big)\psi=0 );

m=0 と置いた式と瓜二つである。*3http://members3.jcom.home.ne.j... を参考にした

というのも、Dirac 方程式では

\{\alpha_i,\alpha_j\}=0\ \ \ i\ne j

\{\alpha_i,\beta\}=0

\alpha_i^2=\beta^2=1

を満たすために \alpha,\beta は4次行列になるが、ここでは \beta に対する条件が必要ないため、 \alpha_i にパウリ行列を入れるだけで上記の条件を満たせるためである。

このとき、

&math( \frac{\hbar}{-i}\frac{\PD}{\PD t}\psi=c\frac{\hbar}{i}\bm\sigma\cdot\bm\nabla\psi=E\psi );

として変数分離した

&math( c\frac{\hbar}{i}\bm\sigma\cdot\bm\nabla\phi=E\phi );

と上記 \begin{pmatrix}C_A\\C_B\end{pmatrix}e^{i\delta\bm k\cdot\bm R} の方程式とがまったく同じ形をしている。

つまり \begin{pmatrix}C_A\\C_B\end{pmatrix}e^{i\delta\bm k\cdot\bm R} は質量のない事由粒子に対する Dirac 方程式の定常状態の解になっている。

以下、

&math( \begin{pmatrix}C_A\\C_B\end{pmatrix}e^{i\delta\bm k\cdot\bm R}= \begin{pmatrix}F_A\\F_B\end{pmatrix} );

と書く。

$F_A,F_B$ の意味

&math( \ket{\Phi} &=\frac{C_A}{\sqrt N}\sum_{\bm R_A}e^{i\bm k\cdot\bm R_A}\ket{\bm R_A}+ \frac{C_B}{\sqrt N}\sum_{\bm R_B}e^{i\bm k\cdot\bm R_B}\ket{\bm R_B}\\ &=\frac{C_Ae^{i\bm k\cdot\bm R}}{\sqrt N}\sum_{\bm R_A}e^{-i\bm k\cdot(\bm R-\bm R_A)}\ket{\phi(\bm R-\bm R_A)}+ \frac{C_Be^{i\bm k\cdot\bm R}}{\sqrt N}\sum_{\bm R_B}e^{-i\bm k\cdot(\bm R-\bm R_B)}\ket{\phi(\bm R-\bm R_B)}\\ &=e^{i\bm k\cdot\bm R}\Bigg(\frac{C_A}{\sqrt N}\sum_{\bm R_A}e^{-i\bm k\cdot(\bm R-\bm R_A)}\ket{\phi(\bm R-\bm R_A)}+ \frac{C_B}{\sqrt N}\sum_{\bm R_B}e^{-i\bm k\cdot(\bm R-\bm R_B)}\ket{\phi(\bm R-\bm R_B)}\Bigg)\\ );

の括弧内が周期関数になるから、 \Phi がブロッホ関数である、というのが基本の話。

\bm k\sim\bm K のとき、 e^{i\bm k\cdot\bm R} は隣り合う基本単位胞でほぼ位相が \pi だけ回るから、 隣り合う基本単位胞で波動関数の符号が反転する。

\bm K を別にしてやると、 e^{i(\bm k-\bm K)\cdot\bm R}=e^{i\delta \bm k\cdot\bm R} はゆっくりと変化する関数となるが、残りの部分は基本単位胞2つで元に戻る、周期が二倍の関数になる。

&math( \ket{\Phi} &=e^{i\delta\bm k\cdot\bm R}\Bigg(\frac{C_A}{\sqrt N}e^{i\bm K\cdot\bm R}\sum_{\bm R_A}e^{-i\bm k\cdot(\bm R-\bm R_A)}\ket{\phi(\bm R-\bm R_A)}+ \frac{C_B}{\sqrt N}e^{i\bm K\cdot\bm R}\sum_{\bm R_B}e^{-i\bm k\cdot(\bm R-\bm R_B)}\ket{\phi(\bm R-\bm R_B)}\Bigg)\\ );

この e^{i\delta \bm k\cdot\bm R} C_A,C_B をまとめたのが F_A,F_B なのだけれど、、、だから何?という感じも残る???

まだちょっとしっくり来ていない。

ワイル方程式

非相対論的な、普通の Dirac 方程式では

&math( \alpha_i=\begin{pmatrix} 0&\sigma_i\\ \sigma_i&0 \end{pmatrix} );, &math( \beta=\begin{pmatrix} 0&\bm 1\\

  • \bm 1&0 \end{pmatrix} );

と取ることが多いが、これは |\bm p|\ll m のときに方程式が対角形になるようにする配慮だった。 *4http://osksn2.hep.sci.osaka-u....

上記の条件を満たす \alpha_i,\beta はこの他にもいろいろ考えられる。

質量のない Dirac 方程式は逆に |\bm p|\gg m の極限と考えられ、 その場合にはカイラル表示、あるいはワイル表示と呼ばれる表示法が用いられる。

&math(\alpha_i=\begin{pmatrix} \sigma_i&0\\ 0&-\sigma_i\\ \end{pmatrix} );, &math( \beta=\begin{pmatrix} \bm 1&0\\ 0&\bm 1\\ \end{pmatrix} );, &math( \gamma_i=\beta\bm\alpha_i\begin{pmatrix} 0&-\sigma_i\\ \sigma_i&0\\ \end{pmatrix} );

たとえばスピン 1/2 で質量がごくわずかしかなく、ほぼ光速で飛び回るニュートリノなどがこの形式で扱われる。

この4成分表示では、

&math( \bm\alpha\cdot\hat{\bm p}\begin{pmatrix} F_A\\F_B\\F'_A\\F'_B \end{pmatrix}=E\begin{pmatrix} F_A\\F_B\\F'_A\\F'_B \end{pmatrix} );

と書けば F_A,F_B の2成分固有値は E F'_A,F'_B の2成分固有値は -E となって、 4成分合わせて上下バンドの A, B サイトに関する波動関数を表せることになる。

ここまで来るともはや完全に質量のない Dirac 方程式と言える。

磁場との相互作用

E\to E+e\phi , \bm p\to \bm p+e\bm A の置き換えをすればよい。

&math( &\pm\bm\alpha\cdot(\hat{\bm p}+e\bm A)\begin{pmatrix} F^\pm_A\\F^\pm_B \end{pmatrix}=\\ &\pm \begin{pmatrix} 0&p_x+eA_x-ip_y-ieA_y\\ p_x+eA_x+ip_y+ieA_y&0\\ \end{pmatrix}\begin{pmatrix} F^\pm_A\\F^\pm_B \end{pmatrix}=(E+e\phi)\begin{pmatrix} F^\pm_A\\F^\pm_B \end{pmatrix} );

すなわち、

&math( \begin{cases} \pm(p_x+eA_x-ip_y-ieA_y)F^\pm_B=(E+e\phi)F^\pm_A\\ \pm(p_x+eA_x+ip_y+ieA_y)F^\pm_A=(E+e\phi)F^\pm_B \end{cases} );

F^\pm_B を消去すると、

&math( &(p_x+eA_x-ip_y-ieA_y)(p_x+eA_x+ip_y+ieA_y)F^\pm_A\\ &=\Big\{(p_x+eA_x)^2+(p_y+eA_y)^2+i(p_x+eA_x)(p_y+eA_y)-i(p_y+eA_y)(p_x+eA_x)\Big\}F^\pm_A\\ &=\Big\{(\bm p_\parallel+e\bm A_\parallel)^2+ie(p_xA_y-p_yA_x-A_yp_x+A_xp_y)\Big\}F^\pm_A\\ &=\Big\{(\bm p_\parallel+e\bm A_\parallel)^2+e(\bm\nabla\times\bm A)_z\Big\}F^\pm_A\\ &=\Big\{(\bm p_\parallel+e\bm A_\parallel)^2+e\bm B_z\Big\}F^\pm_A\\ &=(E+e\phi)^2F^\pm_A );

F^\pm_A を消去すると、

&math( &(p_x+eA_x+ip_y+ieA_y)(p_x+eA_x-ip_y-ieA_y)F^\pm_B\\ &=\Big\{(p_x+eA_x)^2+(p_y+eA_y)^2-i(p_x+eA_x)(p_y+eA_y)+i(p_y+eA_y)(p_x+eA_x)\Big\}F^\pm_B\\ &=\Big\{(\bm p_\parallel+e\bm A_\parallel)^2-ie(p_xA_y-p_yA_x-A_yp_x+A_xp_y)\Big\}F^\pm_B\\ &=\Big\{(\bm p_\parallel+e\bm A_\parallel)^2-e(\bm\nabla\times\bm A)_z\Big\}F^\pm_B\\ &=\Big\{(\bm p_\parallel+e\bm A_\parallel)^2-eB_z\Big\}F^\pm_B\\ &=(E+e\phi)^2F^\pm_B );

\bm\nabla\times\bm A=\bm B は磁束密度。 B_z はその z 成分。

F^\pm_A F^\pm_B とで磁場に対する応答が逆になることが分かる。

固有値がエネルギーの2乗になっているところが良くわからないけれど、 磁場をかけるとエネルギーが分裂したり・・・するのかしら???

ヘリシティ

F_A,F_B z 方向の疑スピン +1, -1 の状態に対応するが、 上で解いた通り、実際には F_A,F_B の状態に同じだけ電子が入るため、 z 方向の疑スピンは消えてなくなる

一方、 x,y 方向の疑スピンは生き残る。その値がいくつになるかを見ておく。

F_A,F_B 成分の位相は \bm k の向きに依存するのだった。

&math( \begin{pmatrix}F_A\\F_B\end{pmatrix}= \begin{pmatrix}1\\\pm e^{i\phi_{\delta\bm k}}\end{pmatrix} );

疑スピンの一の期待値は、

&math( &\frac{1}{\sqrt 2}\begin{pmatrix}1&\pm e^{-i\phi}\end{pmatrix}\ \bm\sigma\ \frac{1}{\sqrt 2}\begin{pmatrix}1\\\pm e^{i\phi}\end{pmatrix}\\ &= \frac{1}{2}\begin{pmatrix}1&\pm e^{-i\phi}\end{pmatrix}\Bigg( \begin{pmatrix}0&1\\1&0\end{pmatrix}\ \begin{pmatrix}0&-i\\i&0\end{pmatrix}\ \begin{pmatrix}1&0\\0&-1\end{pmatrix} \Bigg)\begin{pmatrix}1\\\pm e^{i\phi}\end{pmatrix}\\ &= \frac{1}{2}\begin{pmatrix}1&\pm e^{-i\phi}\end{pmatrix}\Bigg( \begin{pmatrix}\pm e^{i\phi}\\1\end{pmatrix}\ \begin{pmatrix}\mp ie^{i\phi}\\i\end{pmatrix}\ \begin{pmatrix}1\\\mp e^{i\phi}\end{pmatrix} \Bigg)\\ &=\begin{pmatrix}\pm \cos\phi&\pm \sin\phi&0\end{pmatrix}\\ );

より、 x,y 平面内の疑スピンの方向は上バンドでは \bm k 方向、 下バンドでは -\bm k 方向を向く。

ヘリシティ

F_A,F_B z 方向の疑スピン +1, -1 の状態に対応するが、 上で解いた通り、実際には F_A,F_B の状態に同じだけ電子が入るため、 z 方向の疑スピンは消えてなくなる

一方、 x,y 方向の疑スピンは生き残る。

F_A,F_B 成分の位相は \bm k の向きに依存するのだった。

&math( \begin{pmatrix}F_A\\F_B\end{pmatrix}= \begin{pmatrix}1\\\pm e^{i\phi_{\delta\bm k}}\end{pmatrix} );

疑スピンの一の期待値を出すにはヘリシティ演算子 \bm\sigma\cdot\bm k の期待値を求めればよい。

Berry 位相

\phi_k\to\phi_k+2\pi のように \bm k がK点周りを一周したとき、 上記の表式では

&math( \begin{pmatrix}se^{i\phi_{\delta k}+2\pi}\\1\end{pmatrix} =\begin{pmatrix}se^{i\phi_{\delta k}}\\1\end{pmatrix} );

となって、 \begin{pmatrix}C_A\\C_B\end{pmatrix} は元の値に戻る。

ところが、 \begin{pmatrix}C_A\\C_B\end{pmatrix}

&math( \begin{pmatrix} C_A\\C_B \end{pmatrix}=\frac{1}{\sqrt 2}\begin{pmatrix} se^{i\phi_{\delta k}/2}\\e^{i\phi_{\delta k}/2} \end{pmatrix} );

と書くこともできて、このとき

&math( \frac{1}{\sqrt 2}\begin{pmatrix} se^{i(\phi_{\delta k}+2\pi)/2}\\e^{i(\phi_{\delta k}+2\pi)/2} \end{pmatrix}=\frac{1}{\sqrt 2}\begin{pmatrix}

  • se^{i\phi_{\delta k}/2}\\-e^{i\phi_{\delta k}/2} \end{pmatrix} );

のように一周回ると符号が反転してしまう。

Klein tunneling

質量のない Dirac 方程式からは 自身のエネルギーよりも高い障壁を、透過率1でトンネルする解が得られる。

障壁高さを V として、 0<E<V の電子を考える。

この電子は、障壁の外では上側のコーンを、障壁の中では下側のコーンを通過する。

2次元空間で考えるので、障壁は 0<x<D として、入射角度を \phi_k 、 障壁内部での進行方向を \theta_k とする。

入射、出射界面で屈折するので \theta_k\ne \phi_k であることが予想される。

また、それぞれ反射波は x 軸方向のみが折り返されるため、 \pi-\phi_k,\pi-\theta_k となる。

I) x<0

&math( \begin{pmatrix}C_A\\C_B\end{pmatrix}_\mathrm{I} &=\underbrace{\begin{pmatrix}se^{i\phi_k}\\1\end{pmatrix}e^{ik_x x}e^{ik_y y}}_{入射波}\ +\ \underbrace{r\begin{pmatrix}se^{i(\pi-\phi_k)}\\1\end{pmatrix}e^{-ik_x x}e^{ik_y y}}_{反射波}\\ &=\underbrace{\begin{pmatrix}se^{i\phi_k}\\1\end{pmatrix}e^{ik_x x}e^{ik_y y}}_{入射波}\ +\ \underbrace{r\begin{pmatrix}-se^{i\phi_k}\\1\end{pmatrix}e^{-ik_x x}e^{ik_y y}}_{反射波}\\ );

II) 0<x<D

&math( \begin{pmatrix}C_A\\C_B\end{pmatrix}_\mathrm{II} =\underbrace{a\begin{pmatrix}se^{i\theta_k}\\1\end{pmatrix}e^{ik'_x x}e^{ik'_y y}}_{進行波}\ +\ \underbrace{b\begin{pmatrix}-se^{i\theta_k}\\1\end{pmatrix}e^{-ik'_x x}e^{ik'_y y}}_{後退波} );

III) D<x

&math( \begin{pmatrix}C_A\\C_B\end{pmatrix}_\mathrm{III} &=\underbrace{t\begin{pmatrix}se^{i\phi_k}\\1\end{pmatrix}e^{ik_x x}e^{ik_y y}}_{透過波} );

境界条件は界面で連続になること。

&math( \begin{cases} \Phi_\mathrm{I}(0,y)=\Phi_\mathrm{II}(0,y)\\ \Phi_\mathrm{II}(0,y)=\Phi_\mathrm{III}(0,y) \end{cases} );

さて、一見当然のように思えるこの条件はどこからくるのだろう?

C_A,C_B は各原子サイトに置いた原子の波動関数にかける係数だから、 これらの条件がなくても波動関数の連続性は保たれる。当然、その微係数も。 つまり、これらは波動関数の連続性からくる要件ではない。

これは界面に隣接する原子に対応するシュレーディンガー方程式が成り立つようにする条件なのである。

界面の向こう側の原子が、こちら側と異なる係数を持っていたら、界面に隣接する原子に対するシュレーディンガー方程式が変化してしまうでしょ?ということ。

だから、微係数の連続性はここでは必要にならない。

ここの話は http://www.nature.com/nphys/journal/v2/n9/full/nphys384.html の論文でも取り扱われていて、そこに出てくる答えは

r=2ie^{i\phi_k}\sin(qD)\frac{\sin\phi_k-ss'\sin\theta}{ss'[e^{iqD}\cos(\theta-\phi)+e^{-iqD}\cos(\theta+\phi)]-2i\sin(qD)}

\theta=\tan^{-1}\left(\frac{k_y}{q}\right)

上記の境界条件から自分でも解いてみる。まずそれぞれの領域の分散関係から、

\varepsilon=\gamma\sqrt{k_x^2+k_y^2}

に対して、

\varepsilon=\gamma\sqrt{q^2+k_y^2}+V

より、

q=\sqrt{(\varepsilon-V)^2/\gamma^2-k_y^2}

条件式は4つ出てくる。

&math(\left\{\begin{array}{l} 1+r=a+b\\ se^{i\phi}-sre^{-i\phi}=as'e^{i\theta}-bs'e^{-i\theta}\\ ae^{iqD}+be^{-iqD}=te^{ikD}\\ as'e^{iqD}e^{i\theta}-bs'e^{-iqD}e^{-i\theta}=ste^{ikD}e^{i\phi} \end{array}\right. );

ここから a,b,t を消去すれば r が求まる。

第1式より b=-a+(1+r)

第2式に代入して、

 &math( &2\cos\phi-(1+r)e^{-i\phi}=2ass'\cos\theta-(1+r)ss'e^{-i\theta}\\ &(1+r)(e^{-i\phi}-ss'e^{-i\theta})=2\cos\phi-2ss'a\cos\theta );

第3式を第4式に代入し、さらに第1式を使って、

 &math( &as'e^{i(qD+\theta)}-bs'e^{-iqD}e^{-i\theta}=se^{i\phi}(ae^{iqD}+be^{-iqD})\\ &2ass'\cos(qD+\theta)-(1+r)ss'e^{-i(qD+\theta)}= 2iae^{i\phi}\sin(qD)+(1+r)e^{i\phi}e^{-iqD}\\ &(1+r)(e^{i\phi}e^{-iqD}+ss'e^{-i(qD+\theta)})= 2a\big(ss'\cos(qD+\theta)-ie^{i\phi}\sin(qD)\big)\\ );

両式から a を消去すると、

&math( &(1+r)(e^{-i\phi}-ss'e^{-i\theta})=2\cos\phi-ss'\cos\theta (1+r)\frac{e^{i\phi}e^{-iqD}+ss'e^{-i(qD+\theta)}}{ss'\cos(qD+\theta)-ie^{i\phi}\sin(qD)}\\ &(1+r)\left[e^{-i\phi}-ss'e^{-i\theta}+\cos\theta \frac{ss'e^{i\phi}e^{-iqD}+e^{-i(qD+\theta)}}{ss'\cos(qD+\theta)-ie^{i\phi}\sin(qD)}\right]=2\cos\phi\\ &r=2\cos\phi\left[e^{-i\phi}-ss'e^{-i\theta}+e^{-iqD}\cos\theta \frac{ss'e^{i\phi}+e^{-i\theta}}{ss'\cos(qD+\theta)-ie^{i\phi}\sin(qD)}\right]^{-1}-1\\ );

似ても似つかないように見えるが、 値を入れて計算してみるとこれで論文の式と等しいらしい。

他の値は、

&math( &a=\frac{2\cos\phi-(1+r)(e^{-i\phi}-ss'e^{-i\theta})}{2ss'\cos\theta}\\ &b=-a+(1+r)\\ &t=e^{-ikD}(ae^{iqD}+be^{-iqD})\\ );

となる。

反射率を計算してみる

あまり適当な値を入れると r>1 になってしまうので注意が必要。

(E-V)^2-k_y^2<0 となる条件では q が虚数になって、 エネルギー障壁中で指数関数的に減衰する解になるから、 そこは別に計算しなければならない。恐らくこれは通常のトンネル条件。

ここでは (E-V)^2-k_y^2<0 の領域はプロットしないことにする。

\gamma=1,s=1,D=1 で、 V=10,20,30,40,50,100,10000,1000000 について |r|^2 をグラフ化してみた。

r2-v10.jpg r2-v20.jpg r2-v30.jpg r2-v40.jpg

r2-v50.jpg r2-v100.jpg r2-v10000.jpg r2-v1000000.jpg

障壁高さが高いほど反射率の低い範囲が広がるという不可思議な解が得られる。

論文には方向によって反射率が大きく異なる不思議な形のグラフが書かれているのだけれど、 実はこれは特殊な状況でしか起きない。 というか、近似の条件である V\gg E が成り立つ場合にはあんな風にならない。

下は V=50 の時の反射率を濃淡表示して、同心円を重ねた物。

円の半径が小さいとき、すなわち E\ll V の時は反射率の振動はほぼ同心円状に起きているため、反射率におかしな非等方性は現れない。

一方、 E\sim V となる、「ふち」ギリギリを通るような円ではふちのあたりで同心円からずれるために非等方性の強い反射率が得られる。

右は |k| の値を 0 から 30 まで増やしながら、透過率 |t|^2=1-|r|^2 を方向に対してプロットした。

r2-v50-circle.jpg  r2-v50-direction.gif

これが V=1000 になると、 |r|^2 の分布がほぼ同心円になることを反映して、 透過率の特異的な非等方性も失われる。

r2-v1000000-circle.jpg  r2-v1000000-direction.gif

「電子のエネルギーとディラックポイントのエネルギーとの差のちょうど二倍のエネルギー障壁では」 透過確率が1になる条件がある。ということ?

フェルミレベル位置が決まると透過確率が低くなる障壁高さが決まってしまうので、 これはかなり厳しい条件に思える。異なる障壁高さではどうなるのだろうか???


Counter: 42030 (from 2010/06/03), today: 2, yesterday: 0