バンド理論の勉強 の履歴(No.22)
更新
若林克法先生の授業のノートをもとに勉強しなおす†
2015-02-24~26 に、当時 NIMS にいらした若林克法先生 (現在は関西学院大学) が 筑波大学でグラフェンの特異なナノスケール効果について授業をしてくださいました。
その際のノートなどを見ながら、バンド理論について勉強しなおそうと思います。
原子軌道の線形結合により分子軌道を構成(LCAO法の基礎)†
wikipedia:LCAO法 などに紹介されている、広い分野で利用されている方法です。
水素原子 H = 原子軌道†
水素原子では1つの陽子の周りを1つの電子が回っている。この電子に対する時間を含まないシュレーディンガー方程式は、
これを解いた基底状態は 軌道と呼ばれる。
基底状態のエネルギーを
固有関数を としよう。
水素分子 H2 = 分子軌道†
陽子が2つあるとき、その周辺での「電子1つに対する」ハミルトニアンは、
ただし、 はそれぞれの陽子からの距離、となる。
それぞれの陽子を中心とした2つの 1s 軌道の波動関数を元にして、 それらの線形結合によりこの「分子軌道」を近似的に表したい。
陽子2つに電子1つなので、正確には水素分子ではなく水素分子イオン H2+ を考えていることになる。
ものすごく大雑把な近似†
は陽子 1 の周り以外でほぼゼロになると考えて、
と近似すれば、
同様に、
とすれば、任意の に対して
となって、 はエネルギー を持つ固有状態となる。
これは2つの原子が遠く離れて別々に置かれており、 陽子1の周りに の確率で、 陽子2の周りに の確率で、電子が見出される状態に対応する。
ヒュッケル法†
にある電子に対する陽子2からのポテンシャルが完全には無視できない場合にも、 の重なりが小さい場合には を適切に選ぶことで は近似的に の固有関数となる。
を と書くことにする。
となるように を決めるため、左から をかけてみる。
&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});
ここで、 つまり の重なり積分 が十分に小さいとして無視すると(これを 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行目の式は左から をかけて得た)
より左辺の行列はエルミートであり、実数固有値と正規直交固有ベクトルを持つ。
$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_{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 );
時間によらない波動関数は必ず実数関数に取れる。 が実数関数であるとき、 は実数になるため、以降では とする。
このとき、
$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 );
したがって、
のとき、
のとき、
つまり、
のとき、エネルギーは (反結合性軌道)
のとき、エネルギーは (結合性軌道)
である。
結合性軌道は反結合性軌道に比べ、両方の陽子からの引力相互作用を受ける 分子中心付近での状態密度が高いため、 の分だけエネルギーが低くなっている。
ヒュッケル近似の意味†
波動関数は重ならないが(
)
相互作用は取り入れる(
)
という近似が必要になる意味があまり良く理解できていない。。。
このあたりを理解するために、重なり積分をゼロとせずに、もう少しまじめに考えた結果と比べてみる。
とすると、
&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\\ );
すなわち、
のとき、
となる。解の規格化は、
&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 );
より、 となるから、
である。
ヒュッケル近似によりエネルギーが少しだけ変化したが、解の形は変らず、 定性的には良い近似になっていることが分かる(?)。
以下、特に断らない場合にはヒュッケル法を使う。
電子が複数ある影響†
上で考えたのは陽子2つの周りを1つに電子がいる状況。
すなわち、水素分子 H2 ではなく水素分子イオン H2+ である。
水素分子では電子が2つになるため、電子・電子相互作用が入る。
通常、電子・電子相互作用は平均場近似で取り入れる。 すなわち他の電子により部分的に遮蔽された原子核からのポテンシャル中での 1電子の運動を考えることになる。
電子の空間分布が偏りを持つ場合、それぞれの原子核からのポテンシャルは等しくならないが、 ヒュッケル近似ではその影響は無視される。
最終的に1電子軌道から多電子軌道を構成する際には交換相互作用が現れる。
このようなことから、話はそれほど単純ではないとはいえ、 ヒュッケル近似から導かれる以下で見るような性質は、 実際の系でも強く表れる。
4つの原子からなる仮想的な四角い分子†
4個の原子軌道(AO : atomic orbital):
分子軌道(MO : molecular orbital):
ヒュッケル近似で
- 波動関数は重ならない
- クーロン積分
- 飛び移り積分
- 隣接する原子核からだけ相互作用を受ける
とすると、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} );
これを解くとエネルギーは
対応するベクトルは、
つまり、
++ 0- -0 -+ ++ +0 0+ +-
の形になる。
一番左は4つのボンドすべてが結合性、
一番右は4つのボンドすべてが反結合性に対応する。
真ん中の2つは対角線なのでまったく相互作用のない形。
この2つは縮退しているため、その一時結合、例えば
-- +- ++ +-
なども固有値 の固有関数になる。
これら結合性の相互作用が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 を使えば解くのは簡単。
上式の通り は固有値をずらすだけなので、 便宜的に として解いて、固有値を小さいものから順にプロットすると、 ( とした)
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]
何やら cos 的な関数が見えてくる。 (縦軸は を単位としたエネルギー、横軸はエネルギーの低いほうから順に振った量子数)
エネルギーの低いほうから固有ベクトルを5つプロットすると、 固定端に対する弦の振動のような関数が見える。
エネルギーが高くなるに従い、波長が短くなり、波数が大きくなっていく。
ブロッホ波†
なぜ正弦波が現れたのか?
上記方程式を整理すると、
となる。
ただし、
で
とした。
この方程式は固体物理学で1次元格子振動を扱う際に学んだ、 「互いにバネで繋がれたN個の質点」の方程式とそっくりな形をしている。
互いにバネで繋がれたN個の質点 を思い出す†
個体物理では1次元格子振動のモデルとして、 各質点の変位を として、
という形の方程式を解いた。ただし は質点を繋ぐバネのバネ定数、 は質点の質量。
ここでも固定端の条件は で与えられる。少々変形して、
を行列で表せば、
&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} );
より、
として方程式が解ける。
この固有値問題は 次元対称行列に対するものであるから、 重複度を含めて 個の実数固有値と、対応する一時独立な 個の固有ベクトルが存在する。
良く知られるように、質点間の距離を として、
置けば上記方程式は
となり、両辺を で割って、
を満たす が固有値、固有ベクトルの候補となる。
このとき分散関係( のグラフ)は を単位として、
LANG:mathematica Plot[Sqrt[2 (1 - Cos[Pi x])], {x, 0, 1}]
こんな形になる。
ただし、上でも述べた通り や の値は このグラフ上の任意の値をとれるわけではなく、 の条件から許される値はとびとびになる。
の形でこの「境界条件」を満たすことはできないが、 ある の値に足して が条件を満たすとき、同時に も条件を満たすことに注意して( )、
と置き、 とすれば、 に対して は線形独立になる。例えば の場合、 は、
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] ]
グラフにはサイン波を表示したが、実際には格子点上の値しか意味を持たないため、 例えば と とは符号を除いて同一の値を取り、 独立な解を与えない。
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]]] ]
これが の範囲に限られる理由となる。
これらに対応する のそれぞれに対して の2つの が対応するが、これらの線形結合で実数解を求めれば の形に限られることになる。
閑話休題(ブロッホ波再び)†
番目の原子位置を として、
と置けば、このベクトルは
より、
の固有ベクトルとなっている。例えば とすると、 は上式のとおり cos 的な、こんなグラフになる。
LANG:mathematica Plot[-1 - Cos[Pi x], {x, 0, 1}]
上記の固体物理で学んだ分散関係とは一見異なって見えるが、それは の関係にあるためであり、固体物理の分散関係で をプロットすれば同じ形の プロットが得られる。
固有ベクトルの形は上で見たように、 サイトと サイトを固定端とする、正弦波的な形になる。
周期境界条件†
上記の例では固定端となっていたが、これを周期的境界条件にしてみると、 分散関係は固定端の場合と同じだが、周期境界条件で許される波数は
であり、同じ波数に対して
と、
が縮退する。
したがって、固有値問題
&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] &
このように2つずつが等しいエネルギーを持つことになる。
ただし基底状態は定数関数なので sin, cos の区別が無く、縮退しない。
また上のように原子数が偶数の場合、エネルギー最大の状態は隣通しで符号が入れ替わる形になるため、やはり sin, cos の区別はなくなり、縮退しない。
原子数が奇数の場合は、互い違いに入れ替わるベクトルは周期境界条件を満たさず解にならないため、 エネルギー最大の状態はやはり縮退する。
固有ベクトルをエネルギーの低いものから5つ並べれば、
となる。
状態密度†
今の場合 は等間隔に並ぶため、状態密度は に比例する。 ここでは となって、
cos の上下端で状態密度は発散するのは1次元結晶の特徴である。
1次元状態密度†
空間で状態が等間隔に存在しているとき、 つまり例えば に状態があるとき、 空間での状態密度(単位波数当たりの状態密度)は一定値 である。
分散関係 が単調であれば、 空間での状態密度は
より、
である。
したがって、上のように となる点においては、 一次元状態密度は発散する。
一般には となるような は複数あるので、
としなければならない。
この式は、形式的には次のように書くこともできる。
の因子がなくなったのは、 それが暗に に含まれているためだ。
つまり、 であるとき、その周辺で
が成り立つ。すなわち、 で となるとき、
となる。 で書かれたデルタ関数を で積分すると を集めて回ることができるのである。
運動量†
&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 );
すなわち、 に対するシュレーディンガー方程式が次のように得られる。
&math( \underbrace{\left\{\frac{(\hbar k+\hat p)^2}{2m}+V\right\}}_{\hat H(k)}u_k=\varepsilon(k)u_k );
ここに現れるハミルトニアンは に依存するが、 このような場合 ファインマンの定理 により、
&math( \int u_k^*\frac{\PD \hat H}{\PD k}u_k\,dx=\frac{\PD \varepsilon}{\PD k} );
が成り立つ。ここでは であり、
&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} );
を得る。
これと正準方程式 とを比べると、 が に対応する正準運動量になっていそうな気がしてくる。
ただ、ここまで良く見返してみると、上の議論では や の周期性など物理的な条件を一切使っていない。
それでは、任意のハミルトニアンの任意の固有関数に対して、 任意に決めた を使って同様の議論ができてしまうのではないか?
おかしな波数をとってみる†
以下、ちょっと納得のいかない結果になってしまうため、慎重に。
あるポテンシャル に対する時間を含まないシュレーディンガー方程式の エネルギー固有値及び固有関数が、 ほとんど連続な量子数 を用いて と書けるとする。
ポテンシャルの周期性や、波動関数がブロッホ状態であることなどは仮定していないことに注意せよ。
を別の量子数 を使って と表したうえで、
&math( \Phi_k(x) &=e^{ik'x}\underbrace{e^{-ik'x}\Phi_k(x)}_{u_{k'}(x)}\\ );
のように無理やり をくくりだした形から始めて、速度の期待値を求める。
この関数に を作用させると、上記と同様にして が間に挟まる形を導ける。
&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'} );
すなわち、 に対するシュレーディンガー方程式が次のように得られる。
&math( \underbrace{\left\{\frac{(\hbar k'+\hat p)^2}{2m}+V\right\}}_{\hat H(k')}u_{k'}=\varepsilon(k(k'))u_{k'} );
ここに現れるハミルトニアンは に依存するので、 「ファインマンの定理が成り立つならば」
&math( \int u_{k'}^*\frac{\PD \hat H}{\PD k'}u_{k'}\,dx=\frac{\PD \varepsilon}{\PD k'} );
となる。
で、あとは一直線。 より、
&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'} );
となって、右辺は の取り方によって値が変わってしまう。
どうしてこうなった?
怪しいのはファインマンの定理のところだけのように思われる。
英語版の Wikipedia にも書かれているように、 ファインマンの定理が成立するためには関数空間やハミルトニアンに対して一定の条件が必要になる。 Carfì, David (2010). "The pointwise Hellmann–Feynman theorem". AAPP: Physical, Mathematical, and Natural Sciences. 88 (1). no. C1A1001004. doi:10.1478/C1A1001004
たぶん、正しい を取ったときのみ上記の演算が許される?
ファインマンの定理再考†
正しく を取った場合、
&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)}\\ );
この は結晶の端を除いて周期 の周期関数になっている。
を で微分すると
&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) );
となって、これも周期 の周期関数。
おそらくこの、「周期関数になっている」というのが重要な条件になっているのだと思われる。
で気づいたのが、そもそも周期的境界条件以外の束縛状態の定常状態を考える限り、 運動量あるいは速度の期待値が非ゼロになるような関数は解空間に含まれていないことだ。
箱型ポテンシャルや調和振動子では波数 の成分は必ず 波数 の成分とある決まった係数で線形結合を取り、 全体での運動量の期待値をゼロにしないと境界条件を満たさない。
定常的に非ゼロの運動量を持った解は束縛されていないのだからこれは当然。
ファインマンの定理では や を含む 関数ベクトル空間でハミルトニアンがエルミートであることを仮定しているが、 適当な を取ってしまうとその空間でハミルトニアンがエルミートでは なくなってしまうのが問題。*2演算子のエルミート性については 量子力学Ⅰ/固有値と期待値#qf83b503(./?%E9%87%8F%E5%AD%90%E5%8A%9B...) を参照のこと
周期的境界条件では領域の端でちょうど関数値が等しくなることから 内積を部分積分した際に現れる積分値がゼロになり、 運動量演算子のエルミート性が保たれる。
ただ実は、 に対して を周期関数にするには 必ずしも としなければならないわけではなくて、 となっていても良い。ただしこの は結晶格子の逆格子ベクトル。
このとき なので、 を使っても、 を使っても、速度の期待値は変わらない。
これが結晶の分散関係がブリルアンゾーンをまたいで周期的に繰り返す理由になっている。
今の場合†
より、
&math( \langle v\rangle_k=\frac{2\gamma\Delta x}{\hbar}\sin(k\Delta x) );
確かに、 としても の値は変わらない形になっている。
LANG:Mathematica Plot[Sin[Pi x], {x, -1, 1}]
では だが、 ある程度 が大きくなると速度は徐々に上がらなくなり、 で最大値 を取り、 それ以上では が増えると は減少し、 で再び静止する。
有効質量†
外力 が短時間 加わって、 速度の期待値が 増えたとする。 このとき粒子は 進んでいるので、 外力から得たエネルギーは である。
この過程で波数が 増えたとすると、
&math( \Delta E=\frac{\PD E}{\PD k}\Delta k=\braket v \Delta t f );
ここに上で得た の表式を代入すれば、
&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 );
を得る。
すなわち、一定の外力の下では は時間とともに線形に増加する。
一方、このとき粒子の加速度は、
&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 );
より、ニュートン方程式 と比較すると、 を有効質量として
の関係が得られる。
エネルギー対波数の分散関係を基準とすれば、
その傾きが速度であり、
波数の時間変化が力そのものであるため、
分散関係の二次微分が質量の逆数に対応する。
二次微分が大きければ同じ力に対する速度の増加率が大きい = 質量は小さい
二次微分が小さければ同じ力に対する速度の増加率が小さい = 質量は大きい
古典論との対応を見ておこう。
ハミルトニアンが運動エネルギーと保存力によるポテンシャルのみを含んでいる場合、 古典論において正準方程式から
&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\\ );
のようにして、
の関係を導けることに対応している。
3次元空間での一般形は、
&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}\\ );
のように、逆有効質量はテンソルになる。
今の場合†
LANG:mathematica Plot[1/Cos[Pi k], {k, 0, 1}]
では より となる。
一方、 では より のように負の質量が現れる。
負の有効質量†
負の有効質量をもつ電子は、電場と同じ方向へ加速される。(通常の電子は電場と逆の方向へ加速される)
ここからどのような結論が導かれるだろうか?
電子が一つのとき†
バンドの下端に電子が1つある状況から、外場 によりこの電子を加速する。
散乱がなければ電子の波数は時間とともに一様に増加して、
より、
速度および位置は、
のように振動する。
力が強くなれば速度変化が速くなり、振幅が小さくなる。
この振動はブロッホ振動と呼ばれる。
外場がゼロになると なので、等速直線運動を続けることになる。
実際には散乱があるので徐々に減衰して際安定状態に落ち込み、静止する。
散乱が強ければ電子が最高速度を超えることはなく、 平均的な速度が外場に比例するオーミックな電流が流れる。
電子が底の方に詰まっているとき†
最安定状態では 側と 側に同じだけ電子がいるため、 電流は相殺されてゼロになる。
外場 がかかると、 すべての電子の が で変化して、 側の電子が減り、 側の電子が増える。
その差の分だけ電流が流れる。
電子が半分以上詰まっているとき†
外場 がかかると、 すべての電子の が で変化して、 側の電子が減り、 側の電子が増える。
の電子は負の有効質量のために正方向の速度が減速され、
の電子は負の有効質量のために負方向の速度が加速されるが、
全体として見るとやはり負方向へ動く電子が減って、正方向へ動く電子が増えるため、
電流は力の加わった方向へ流れることになる。
詰まり具合と電流応答†
電子が底の方に少しだけ詰まっている時、 外場によって が少しだけ変わっても、 正側、負側とも速度が遅いので、電流の応答は小さい、と思われる。
電子が半分まで詰まっている時、 外場によって が少し変わるだけでも、 最高速で負側に動く電子が減って、最高速で正側に動く電子が増えるので、 電流の応答は大きい、と思われる。
電子がかなり上まで詰まっている時、 やはり正側、負側とも速度が遅いので、電流の応答は小さい、と思われる。
1個だけ残して残りが埋まっているとき†
外場 がかかると、 すべての電子の が で変化するため、 空いた電子状態の も同じように変化する。
この電子状態は負の有効質量をもっているため、電子に加わる力とは反対方向に加速される。
電子の空状態なので、相対的に正の電荷を帯びている。
この電子の空準位がホールと呼ばれる。
2次元三角格子†
1原子状態をつなげてブロッホ波的な三角格子を作ると、これが実際の三角格子の解になっていることを確かめる。
個の原子からなる三角格子中の原子位置を で表し、 これらの原子を中心とする1原子状態を元にブロッホ状態を作る。
3角格子では最隣接原子が6個あり、基本ベクトルを と置けば
に位置する。したがって係数間の方程式は
であり、 と置いてみるとこれは、
を満たす固有値 の固有ベクトルとなることが分かる。
今の場合 、 、 であるから、その分散関係は次のようになる。
honeycomb[k_] = Cos[k.{1, 0}] + Cos[k.{1/2, Sqrt[3]/2}] + Cos[k.{-1/2, Sqrt[3]/2}]
これらは、 に対して、対応する をプロットしたものである。
逆格子†
実空間での基本ベクトルが
であるのに対応して、逆格子空間における基本ベクトルは、
となる。
に注意せよ。 A
有効質量†
2次元なので、方向によって有効質量が異なる場合がある。
より、
三角格子の場合†
上で計算した分散関係を再掲すると、
ガンマ点、K'点、M点で速度ゼロ。その他の点では傾きに応じた方向の速度を持つ。
ガンマ点でエネルギー最大、曲率および有効質量は等方的で負の値を持つ。
K' 点でエネルギー最小、曲率および有効質量は当方的で正の値を持つ。
M点は鞍点で、曲率および有効質量は一方向に正、一方向に負の値を持つ。
グラフェンの場合†
三角格子の単位胞に2つの原子 A, B がある構造。図ではそれぞれ白と赤。
系全体の波動関数を、A に対する三角格子波動関数と、B に対する三角格子波動巻子との線形結合で
と置き、 を仮定すると、
より、
&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} );
を得る。
は等しく、三角格子のエネルギーと隣接原子からのクーロン相互作用とを加えたものになる。前者は上で解いたように波数に依存して変化するが、後者はブロッホ波の位相によらず、状態密度のみに依存する量であるため定数である。
であり、 隣り合う原子上の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行目が で始まっているのは、 境界上の原子では周囲に3原子存在しない場合があり、 そのあたりの境界条件がうやむやになっているため。
上図の白丸を A サイトとすると、隣接 B サイトは
となるから、これらを代入すれば が求まる。
このとき、
&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} );
より、
より、 と置けば、
&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} );
分散関係†
グラフェンの特異的な分散関係を見るために、まずは のグラフを見る。
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]
分散関係はゼロを境に上下に対象であり、「卵のパック」のような形になる。
大きく膨らんでいる場所がガンマ点、 鋭くへこんでいる場所がK点、 鞍点となっているのがM点である。
中でもK点では、
に対して より、
&math( f_{\bm k}&=\exp(0)+\exp(i2\pi/3)+\exp(-i2\pi/3)\\ &=1+w+w^{-1}\\ &=0 );
ただし、 は 1 の3乗根。
となって、 すなわちギャップは閉じている。
状態密度†
このとき、状態密度をモンテカルロ法により求めると、
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}]
ここでは に取っている。
では状態密度はステップ関数的に立ち上がる。
これはΓ点のパラボリックな分散関係を反映している。
では発散する。
これはM点に存在する鞍点を反映している。
では線形に立ち上がる。
これはK点、K’点に存在する円錐状の分散関係を反映している。
以下、少し詳しく見てみる。
2次元状態密度†
1次元の時と同様に、2次元状態密度は形式的に
と書ける。この形は解析的に考えることができる場合には有用であるが、 数値計算可能な形ではないので、計算できる形に変形したい。
と書ければ良いのだが、 右辺の範囲が必ずしも と一致しないため難しい。
方向ではなく、 方向と、それに垂直な 方向に軸を取れば、
ただし 。したがって、
ここで、 は となる等高線に沿って計算する。
とすると、格段にイメージはしやすくなるが、 これも数値計算可能な形は限られている。。。
ある点 $\bm k_0$ を中心に回転対称性があるとき†
を用いて、 と書けるから、
は円周方向を向くノルム一定のベクトルで、 は円周に沿った長さ ( は半径) に渡る積分になり、
となる。
(1) ガンマ点の周りの円錐型の部分では
であるから、
のように に対して1次で立ち上がる。
さらに
であるから、
となる。
(2)ガンマ点の周りの回転放物面的な部分では [#x6fa440b]
より、
であるから、
のように定数値になる。
鞍点では†
の形の理想的な鞍点において、 特定のエネルギーでの切断面は
- では2つの直線が交わる形
- では双曲線
LANG:mathematica Plot3D[x^2 - y^2, {x, -1, 1}, {y, -1, 1}, MeshFunctions -> {#3 &}, ImageSize -> Large, PlotPoints -> 101]
だから、
の直線上 の上の 原点から距離 の点 において となって、
のように発散する。
モンテカルロじゃなくても計算できるか?†
電子状態密度を上記のデルタ関数を用いた形で 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 はうまく扱えないみたいだった。
三角格子の分散関係を加える†
正確には、上記のまさにグラフェンらしい分散関係に、単純三角格子の分散関係を足す必要がある。 ( は定数なので無視する)
単純三角格子の時の方が原子間距離が 倍遠いので、 単純計算ではクーロン相互作用が になる。 原子が複数いる場合には遮蔽効果によりクーロン相互作用はより急峻に減衰するため。 実際にはさらに影響は小さいと考えられる。
実際に相互作用がどのくらい違うかよく分からないまま、 とりあえず 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 &}]
ガンマ点での高さや有効質量の絶対値が正側と負側とで異なる
一方、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}]
左が正側、右が負側
このようにプロットしてみると、コーンの裾の部分は完全な丸じゃなくて丸まった三角形になることがよく分かる。
コーンの中心部分は円錐になっている。
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}]
Dirac 方程式†
理想的なグラフェンでは2枚のバンドのうち下のバンドのみが電子で埋まっていて、 ちょうどディラックコーンの先端部分に「フェルミ点」が来る。
電場を印可すればこのコーンの近くで電子の分布が変化するが、 電子の有効質量はバンドの湾曲に対応するため、 直線的な円錐状のバンド構造では、円錐の中心から見て放射方向への有効質量が発散している。
そのため、円錐を回る方向には慣性が働くが、放射方向には慣性が働かない?
以下、そのような運動を記述する Dirac 方程式を導く。
上では2つの三角格子を重ね合わせる形でグラフェンの電子状態を構成したが、 ここでは原子を1つずつ足すことで再構成しよう。
&math(\ket \phi=\sum_{n'}\left( C_Ae^{i\bm k\cdot \bm R_{An'}}\ket{\bm R_{An'}}+ C_Be^{i\bm k\cdot \bm R_{Bn'}}\ket{\bm R_{Bn'}} \right));
を求めるため、シュレーディンガー方程式に左から を掛ける。 上では3角格子内部の相互作用、すなわち A-A や B-B の相互作用を考えたが、 ここでは本当の再隣接原子間の相互作用以外を無視する。 また、 も に含めることにして、便宜上ゼロに置く。
&math( \braket{\bm R_{An}|H|\phi}=\braket{\bm R_{An}|E|\phi} );
&math(
- \gamma C_B\sum_{l=1}^3 e^{i\bm k\cdot \delta \bm R_{BA,l}}=EC_A );
、 と書くことにすると、 ( は K 点の座標) &math( &-\gamma C_B\sum_{l=1}^3 e^{i\bm K\cdot\bm\tau_l} e^{i\delta\bm k\cdot\bm\tau_l}=EC_A\\ );
今考えている範囲で は小さいためその1次までを残し、
&math( &-\gamma C_B\sum_{l=1}^3 e^{i\bm K\cdot\bm\tau_l} (1+i\delta\bm k\cdot\bm\tau_l+\dots)=EC_A\\ &-i\gamma \left(\sum_{l=1}^3 \bm\tau_l e^{i\bm K\cdot\bm\tau_l}\right)\cdot\delta\bm k C_B \sim EC_A\\ &-i\gamma \frac{\sqrt 3}{2}\begin{pmatrix}i\\-1\end{pmatrix}\cdot\delta\bm k C_B \sim EC_A\\ );
ここで、
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}
より、 および、 を用いた。
以上より、
&math( \gamma \frac{\sqrt 3}{2}(\delta k_x+i\delta k_y) C_B \sim EC_A );
左から を掛けた式からは、 を使って、
&math( \gamma \frac{\sqrt 3}{2}(\delta k_x-i\delta k_y) C_A \sim EC_B );
を得る。
まとめると、
&math( \gamma\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} );
同様のことを K' 点で行えば、
&math( \gamma\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} );
となり、これらは質量のない Dirac 方程式 = Weyl 方程式と呼ばれる。
解いてみるとどうなるか。
より、
確かに円錐形の分散関係が得られる。
代入して固有ベクトルを求める。
ただし、 で、 の方位を表す。
&math( \begin{pmatrix} C_A\\C_B \end{pmatrix}=\frac{1}{\sqrt 2}\begin{pmatrix} se^{i\phi_{\delta k}}\\1 \end{pmatrix} );
以下、 を と書く。
ええと、、、この結果と、上で得られた結果で と置いた
&math( \begin{pmatrix}C_A\\C_B\end{pmatrix}\propto \begin{pmatrix} \pm 1\\
- e^{i\phi_k} \end{pmatrix} );
との整合性はどうなっているのだろう???
何か、ほとんどそのままのような???
Berry 位相†
のように がK点周りを一周したとき、 上記の表式では
&math( \begin{pmatrix}se^{i\phi_{\delta k}+2\pi}\\1\end{pmatrix} =\begin{pmatrix}se^{i\phi_{\delta k}}\\1\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でトンネルする解が得られる。
障壁高さを として、 の電子を考える。
この電子は、障壁の外では上側のコーンを、障壁の中では下側のコーンを通過する。
2次元空間で考えるので、障壁は として、入射角度を 、 障壁内部での進行方向を とする。
入射、出射界面で屈折するので であることが予想される。
また、それぞれ反射波は 軸方向のみが折り返されるため、 となる。
I)
&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)
&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)
&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} );
さて、一見当然のように思えるこの条件はどこからくるのだろう?
は各原子サイトに置いた原子の波動関数にかける係数だから、 これらの条件がなくても波動関数の連続性は保たれる。当然、その微係数も。 つまり、これらは波動関数の連続性からくる要件ではない。
これは界面に隣接する原子に対応するシュレーディンガー方程式が成り立つようにする条件なのである。
界面の向こう側の原子が、こちら側と異なる係数を持っていたら、界面に隣接する原子に対するシュレーディンガー方程式が変化してしまうでしょ?ということ。
だから、微係数の連続性はここでは必要にならない。
ここの話は http://www.nature.com/nphys/journal/v2/n9/full/nphys384.html の論文でも取り扱われていて、そこに出てくる答えは
上記の境界条件から自分でも解いてみる。まずそれぞれの領域の分散関係から、
に対して、
より、
条件式は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. );
ここから を消去すれば が求まる。
第1式より
第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)\\ );
両式から を消去すると、
&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})\\ );
となる。
反射率を計算してみる†
あまり適当な値を入れると になってしまうので注意が必要。
となる条件では が虚数になって、 エネルギー障壁中で指数関数的に減衰する解になるから、 そこは別に計算しなければならない。恐らくこれは通常のトンネル条件。
ここでは の領域はプロットしないことにする。
で、 について をグラフ化してみた。
障壁高さが高いほど反射率の低い範囲が広がるという不可思議な解が得られる。
論文には方向によって反射率が大きく異なる不思議な形のグラフが書かれているのだけれど、 実はこれは特殊な状況でしか起きない。 というか、近似の条件である が成り立つ場合にはあんな風にならない。
下は の時の反射率を濃淡表示して、同心円を重ねた物。
円の半径が小さいとき、すなわち の時は反射率の振動はほぼ同心円状に起きているため、反射率におかしな非等方性は現れない。
一方、 となる、「ふち」ギリギリを通るような円ではふちのあたりで同心円からずれるために非等方性の強い反射率が得られる。
右は の値を 0 から 30 まで増やしながら、透過率 を方向に対してプロットした。
これが になると、 の分布がほぼ同心円になることを反映して、 透過率の特異的な非等方性も失われる。
*1 http://www1.doshisha.ac.jp/~makato/lecture/IC_S1/eff_mass.pdf を参考にした
*2 演算子のエルミート性については 量子力学Ⅰ/固有値と期待値#qf83b503 を参照のこと