バンド理論の勉強 のバックアップソース(No.17)

更新

[[公開メモ]]

#contents
&mathjax();

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

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

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

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

[[wikipedia:LCAO法]] などに紹介されている、広い分野で利用されている方法です。

** 水素原子 H = 原子軌道 [#se7e6e74]

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

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

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

基底状態のエネルギーを&math(\varepsilon_{1s});

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

** 水素分子 H__2__ = 分子軌道[#w25fa0fb]

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

&math(\hat H=\frac{-\hbar^2}{2m}\nabla^2-\frac{e^2}{r_1}-\frac{e^2}{r_2});

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

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

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

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

*** ものすごく大雑把な近似 [#eda06a14]

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

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

と近似すれば、

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

同様に、

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

とすれば、任意の &math(c_1,c_2); に対して

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

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

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

*** ヒュッケル法 [#edd86937]

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

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

&math(\hat H\ket\Phi=E\ket\Phi); となるように &math(c_1,c_2); を決めるため、左から &math(\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});

ここで、 &math(\delta\ll 1); つまり &math(\phi_1,\phi_2); の重なり積分 &math(\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行目の式は左から &math(\bra 2); をかけて得た)

&math(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}$ について [#s359b189]

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

同様にして &math(h_{22}=-\alpha); となる。

*** $h_{21},h_{12}$ について [#s359b189]

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

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

このとき、&math(h_{12}=h_{21}^*=h_{21}=-\gamma);

*** $c_1,c_2$ を求める [#g4b2b66d]

&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(E=-\alpha\pm\gamma);

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

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

つまり、

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

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

である。

&ref(bonding_and_antibonding.png);

結合性軌道は反結合性軌道に比べ、両方の陽子からの引力相互作用を受ける
分子中心付近での状態密度が高いため、
&math(\braket{1|\hat H|2}); の分だけエネルギーが低くなっている。

*** ヒュッケル近似の意味 [#zb50686f]

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

このあたりを理解するために、重なり積分をゼロとせずに、もう少しまじめに考えた結果と比べてみる。

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

すなわち、

&math(\begin{pmatrix}c_1\\c_2\end{pmatrix}\propto\begin{pmatrix}1\\\pm1\end{pmatrix}); のとき、
&math(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
); 

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

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

である。

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

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

*** 電子が複数ある影響 [#y43e2c3a]

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

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

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

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

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

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

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

** 4つの原子からなる仮想的な四角い分子 [#g8136cc3]

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

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

ヒュッケル近似で
- 波動関数は重ならない &math(\braket{i|j}=0\ \ (i\ne j));
- クーロン積分 &math(\braket{i|\hat H|i}=-\alpha);
- 飛び移り積分 &math(\braket{i|\hat H|j}=-\gamma\ \ (i,j は隣り合う));
- 隣接する原子核からだけ相互作用を受ける &math(\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}
);

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

&math(E=-a-2g,\ -a,\ -a,\ -a+2g);

対応するベクトルは、

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

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

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

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

つまり、

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

の形になる。

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

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

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

  --  +-
  ++  +-

なども固有値 &math(-\alpha); の固有関数になる。

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

* 結晶を作る [#x835baba]

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

上の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 を使えば解くのは簡単。

上式の通り &math(\alpha); は固有値をずらすだけなので、
便宜的に &math(\alpha=0,\gamma=1); として解いて、固有値を小さいものから順にプロットすると、
(&math(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]

&ref(n-atoms-energy.png);

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

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

&ref(n-atoms-wave-function.png);

エネルギーが高くなるに従い、波長が短くなり、波数が大きくなっていく。


** ブロッホ波 [#tdf3e943]

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

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

&math(c_{n+1}+\frac{E+\alpha}{\gamma}c_n+c_{n-1}=0);

となる。~
ただし、&math(n=1\sim N); で &math(c_0=c_{n+1}=0); とした。

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

** 互いにバネで繋がれたN個の質点 を思い出す [#a7b30c25]

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

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

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

ここでも固定端の条件は &math(u_0=u_{N+1}=0); で与えられる。少々変形して、

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

より、

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

として方程式が解ける。

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

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

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

置けば上記方程式は

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

となり、両辺を &math(e^{ink\Delta x}); で割って、

&math(e^{ik\Delta x}-2+e^{-ik\Delta x}=-\omega^2\frac{M}{K});

&math(\omega=\pm\sqrt{\frac{2K}{M}}\sqrt{1-\cos ik\Delta x});

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

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

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

&attachref(1d-chain-dispersion1.svg);

こんな形になる。

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

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

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

と置き、&math(k_m=m\pi/\Delta x(N+1)); とすれば、&math(m=1\dots N); に対して
&math(\bm u_m); は線形独立になる。例えば &math(N=5); の場合、&math(\{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]
 ]

&attachref(sine-waves.svg);

グラフにはサイン波を表示したが、実際には格子点上の値しか意味を持たないため、
例えば &math(m=(N+1)-1); と &math(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]]]
 ]

&attachref(1d-chain-dispersion3.svg);

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

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

** 閑話休題(ブロッホ波再び) [#xa17b785]

&math(n); 番目の原子位置を &math(R_n=n\Delta R); として、

&math(c_n=e^{ikR_n});

と置けば、このベクトルは

&math(e^{ik\Delta R}+\frac{E+\alpha}{\gamma}+e^{-ik\Delta R}=0);

より、

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

の固有ベクトルとなっている。例えば &math(\alpha=\gamma=1); とすると、
&math(E(k)); は上式のとおり cos 的な、こんなグラフになる。

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

&attachref(1d-chain-dispersion2.svg);

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

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

&math(c_n^m\propto \sin(k_mR_n));

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

** 周期境界条件 [#jd161846]

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

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

であり、同じ波数に対して

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

と、

&math(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] &

&ref(cyclic-boundary_energy.png);

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

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

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

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

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

&ref(cyclic-boundary_waves.png);

となる。

** 状態密度 [#l71bfbb6]

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

&ref(1d-chain_dos.png,,50%);

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

** 1次元状態密度 [#gc672e7e]

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

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

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

より、

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

である。

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

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

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

としなければならない。

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

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

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

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

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

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

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

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

** 有効質量 [#efee3301]

[[量子力学Ⅰ/群速度と波束の崩壊/メモ#b368cfc6]] において、波数 &math(k); の電子の群速度は分散関係から

&math(v=\frac{\PD \omega}{\PD k});

で求められることを導いた。&math(E=\hbar \omega); とすれば、

&math(v=\frac{1}{\hbar}\frac{\PD E}{\PD k});       (*)

である。(これは &math(\dot x=\frac{\PD H}{\PD p}); に対応する)

外力 &math(f); によるエネルギー変化 

&math(\Delta E=(f v)\Delta t); 

がすべて運動量すなわち &math(k); の変化にあてられた場合、

&math(\Delta E=\frac{\PD E}{\PD k}\Delta k=\hbar v \Delta k);

と表せるから、

&math(f \Delta t=\hbar \Delta k); すなわち &math(\frac{dk}{dt}=\frac{f}{\hbar});

の関係が導かれる。このとき (*) 式より、

&math(
\frac{dv}{dt}
&=\frac{1}{\hbar}\frac{d}{dt}\frac{\PD E}{\PD k}\\
&=\frac{1}{\hbar}\frac{\PD}{\PD k}\frac{d E}{dt}\\
&=\frac{1}{\hbar}\frac{\PD}{\PD k}\frac{\PD E}{\PD k}\frac{dk}{dt}\\
&=\frac{1}{\hbar}\frac{\PD^2 E}{\PD k^2}\frac{dk}{dt}\\
&=\frac{1}{\hbar^2}\frac{\PD^2 E}{\PD k^2}f\\
);

これとニュートン方程式 &math(f=m^*\frac{dv}{dt}); とを比べれば、

&math(\frac{1}{m^*}=\frac{1}{\hbar^2}\frac{\PD^2 E}{\PD k^2});

を得る。((この部分は http://www.sci.osaka-cu.ac.jp/~matrsc/DENSHI/yoshino/basictxt_old.htm を参考にさせていただきました))

というか、以上は古典論の話なので、正準方程式から

&math(\dot p=-\frac{\PD H}{\PD x}=-\frac{\PD V}{\PD x}=f);

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

より、&math(\frac{1}{m}=\frac{\PD^2 H}{\PD p^2}); とした方が見通しがいいのかもしれない。
(たぶんここではハミルトニアンが運動エネルギーと保存力によるポテンシャルのみを含んでいることを仮定している)

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

&math(\dot{\bm p}=-\bm\nabla_x H=-\bm\nabla_x V=\bm f);

&math(\dot{\bm x}=\bm \nabla_p H);

&math(\ddot{\bm x}&=\frac{d}{dt}\bm\nabla_p H\\
&=\bm\nabla_p \frac{dH}{dt}\\
&=\bm\nabla_p (\bm\nabla_p H\cdot \dot{\bm p})\\
&=\Bigg(\frac{\PD^2 H}{\PD p_i\PD p_j}\Bigg) \dot{\bm p}\\
&=\Bigg(\frac{\PD^2 H}{\PD p_i\PD p_j}\Bigg) \bm f\\
);

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

*** 今の場合 [#xac50074]

&math(\frac{\PD E}{\PD k}=2\Delta R\gamma\sin(k\Delta R));

&math(\frac{\PD^2 E}{\PD k^2}=2\Delta R^2\gamma\cos(k\Delta R));

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

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

&attachref(effective-mass1.svg);

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

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

** 運動量について [#fe9354ea]


* 2次元三角格子 [#p818be3f]

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

~

&ref(triangular-lattice.png,,20%);

~

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

&math(|\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個あり、基本ベクトルを &math(\bm a,\bm b); と置けば

&math(\bm a, -\bm a, \bm b, -\bm b,-\bm a+\bm b, \bm a-\bm b);

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

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

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

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

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

今の場合
&math(\bm a=\begin{pmatrix}a\\0\end{pmatrix});、
&math(\bm b=\begin{pmatrix}a/2\\\sqrt{3}\,a/2\end{pmatrix});、
&math(\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}]

&ref(triangular-band1.jpg,,33%);
&ref(triangular-band2.jpg,,33%);
&ref(triangular-band3.jpg,,33%);

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

** 逆格子 [#ha2be66f]

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

&math(\bm a=(a,0));

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

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

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

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

となる。

&math(\bm a\cdot \bm a^*=\bm b\cdot \bm b^*=2\pi);

&math(\bm a\cdot \bm b^*=\bm b\cdot \bm a^*=0);

に注意せよ。
A

&ref(triangular-brillouin-zone.png,,75%);

* グラフェンの場合 [#e0b64894]

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

&ref(graphene-lattice.png,,20%);

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

&math(\ket \Phi=C_A\ket{\Phi_A}+C_B\ket{\Phi_B}); 

と置き、&math(\braket{\Phi_A|\Phi_B}=0); を仮定すると、

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

を得る。

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

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

&math(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行目が &math(\sim); で始まっているのは、
境界上の原子では周囲に3原子存在しない場合があり、
そのあたりの境界条件がうやむやになっているため。

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

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

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

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

となるから、これらを代入すれば &math(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}
);

より、

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

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

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

** 分散関係 [#k06ca563]

グラフェンの特異的な分散関係を見るために、まずは &math(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]

&ref(graphene-dispersion1.jpg,,25%);
&ref(graphene-dispersion2.jpg,,33%);
&ref(graphene-dispersion3.jpg,,33%);
&ref(graphene-dispersion4.jpg,,25%);

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

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

中でもK点では、

&math(f_k=\sum_{l=1}^3 e^{i\bm k\cdot \Delta R_l});

に対して &math(\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
);

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

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

** 状態密度 [#je78b47f]

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

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

&ref(graphene-dos.png,,50%);

ここでは &math(\gamma=1); に取っている。

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

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

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

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

** 2次元状態密度 [#we71c6a7]

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

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

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

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

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

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

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

ただし &math(dE=|\nabla_k E|dk_1);。したがって、

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

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

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

*** ある点 $\bm k_0$ を中心に回転対称性があるとき [#cb3bb0c5]

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

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

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

&math(\frac{\PD \delta k}{\PD E}=\mathit{const.});

であるから、

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

のように &math(\delta k); に対して1次で立ち上がる。

さらに

&math(\frac{\PD \delta k}{\PD E}=\mathit{const.});

であるから、

&math(\rho(E)\propto E);

となる。

~

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

&math(E\propto \delta k^2);

より、

&math(\frac{\PD \delta k}{\PD E}\propto \frac{1}{\delta k});

であるから、

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

のように定数値になる。

*** 鞍点では [#wa53d645]

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

-&math(E=E_0); では2つの直線が交わる形
-&math(E\ne E_0); では双曲線

 LANG:mathematica
 Plot3D[x^2 - y^2, {x, -1, 1}, {y, -1, 1}, MeshFunctions -> {#3 &}, 
 ImageSize -> Large, PlotPoints -> 101]

&ref(suddle-point.jpg,,50%);

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

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

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

のように発散する。

#clear


** モンテカルロじゃなくても計算できるか? [#h618c568]

電子状態密度を上記のデルタ関数を用いた形で 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 はうまく扱えないみたいだった。

** 三角格子の分散関係を加える [#m915707f]

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

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

単純三角格子の時の方が原子間距離が &math(\sqrt{3}); 倍遠いので、
単純計算ではクーロン相互作用が &math(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 &}]

&ref(graphene-dispersion5.jpg,,66%);

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

&ref(graphene-dispersion6.jpg,,50%);

一方、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}]

&ref(graphene-cone+.jpg,,75%);
&ref(graphene-cone-.jpg,,75%);

左が正側、右が負側

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

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

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

&ref(graphene-cone0.jpg);

* 質量のない Dirac 方程式 [#y1612907]

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

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

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