Kronig-Penneyモデルのバンド計算

プログラムの動作

kronig_penney.py は、一次元結晶における電子のエネルギーバンド構造と波動関数をKronig-Penneyモデルに基づいて計算・可視化するためのPythonプログラムです。このモデルは、周期的な矩形ポテンシャルを仮定することで、固体中の電子がどのように振る舞うかを理解するための基礎的な枠組みを提供します。

プログラムは以下の3つの主要なモードで動作します。

  1. graph モード:

    • 特定のブロッホ波数 \(k\) におけるKronig-Penney方程式の左辺（プログラムでは delta と呼ばれる）をエネルギー \(E\) の関数としてプロットします。

    • このプロットは、電子の許容エネルギー準位がどのようにして決定されるか（すなわち、delta が \(\cos(ka)\) に等しくなる点を見つけること）を視覚的に示します。

    • エネルギー準位の探索にはSecant法が用いられます。

  2. band モード:

    • ブロッホ波数 \(k\) の範囲にわたるエネルギーバンド構造を計算しプロットします。

    • これは、結晶中の電子が取りうるエネルギーが連続的なバンドと禁止帯（バンドギャップ）に分かれる様子を示します。

    • 各 \(k\) 点でエネルギー準位を繰り返し探索し、結果をプロットします。

  3. wf モード:

    • 特定のブロッホ波数 \(k\) と準位番号 \(iLevel\) に対応する波動関数 \(\Psi(x)\) およびその電荷密度 \(|\Psi(x)|^2\) を計算しプロットします。

    • 同時に、周期的な矩形ポテンシャルの形状も表示されます。

    • 波動関数の係数はKronig-Penneyモデルの境界条件から導出される連立一次方程式を解くことで求められ、その後正規化されます。

原理

Kronig-Penneyモデルは、自由電子のシュレーディンガー方程式に周期的な矩形ポテンシャルを導入することで、固体中の電子の運動を記述します。ポテンシャルは、幅 \(w\) のポテンシャル井戸と幅 \(b\) のポテンシャル障壁（高さ \(V_0\)）が交互に配置された構造を持ちます。単位胞の格子定数は \(a = w + b\) です。

Kronig-Penney方程式: ブロッホの定理により、周期的なポテンシャル中の電子の波動関数は \(\Psi_k(x) = e^{ikx}u_k(x)\) の形をとります。ここで \(u_k(x)\) はポテンシャルと同じ周期 \(a\) を持つ関数です。シュレーディンガー方程式の解と境界条件（波動関数とその導関数が連続であること）を適用すると、エネルギー \(E\) とブロッホ波数 \(k\) の関係を示す以下の超越方程式（Kronig-Penney方程式）が得られます。

\[ \frac{\beta^2 - \alpha^2}{2\alpha\beta} \sin(\alpha w) \sinh(\beta b) + \cos(\alpha w) \cosh(\beta b) = \cos(ka) \]

ここで、\(E < V_0\) の場合、

  • \(\alpha = \frac{\sqrt{2m_e E}}{\hbar}\)

  • \(\beta = \frac{\sqrt{2m_e (V_0 - E)}}{\hbar}\) そして、

  • \(m_e\): 電子の質量

  • \(\hbar\): 換算プランク定数

  • \(e\): 素電荷

  • \(w\): ポテンシャル井戸の幅（a - bwidth）

  • \(b\): ポテンシャル障壁の幅（bwidth）

  • \(V_0\): ポテンシャル障壁の高さ（bpot）

  • \(k\): ブロッホ波数

プログラム内の cal_delta 関数は、この方程式の左辺 \(\frac{\beta^2 - \alpha^2}{2\alpha\beta} \sin(\alpha w) \sinh(\beta b) + \cos(\alpha w) \cosh(\beta b)\) から右辺 \(\cos(ka)\) を引いた値を返します。したがって、許容されるエネルギー準位 \(E\) は cal_delta(E, ...) がゼロとなる点（あるいは delta = cos(ka) となる点）として見つけられます。

エネルギー準位の探索: find_Elist 関数は、指定されたエネルギー範囲内で cal_delta(E, k, ...) がゼロになる（または符号を反転する）エネルギー \(E\) を探索します。発見されたおおよその根は、refine_E 関数が実装するSecant法によって高精度に精錬されます。Secant法は、関数の根を見つけるための数値反復法の一つで、導関数を計算することなく、2つの初期点から関数の値を線形近似して次の推定値を求めます。

波動関数の計算: cal_wavefunction 関数は、特定のエネルギー \(E\) と波数 \(k\) における波動関数を計算します。これは、それぞれの領域（井戸と障壁）におけるシュレーディンガー方程式の一般解と、境界条件およびブロッホの定理を適用して得られる4つの係数 \(c_0, c_1, c_2, c_3\) （プログラムでは ci リストとして扱われる）を用いて構築されます。これらの係数は、find_Elist 内で連立一次方程式を解くことで決定されます。

必要な非標準ライブラリとインストール方法

このプログラムは以下の非標準ライブラリを使用します。

  • numpy: 数値計算（配列操作、数学関数など）

  • matplotlib: グラフの描画

これらのライブラリは pip コマンドでインストールできます。

pip install numpy matplotlib

以下の標準ライブラリも使用していますが、通常はPythonにプリインストールされています。

  • sys: システム固有のパラメータと機能

  • pprint: データ構造を整形して出力

  • csv: CSVファイルの読み書き

必要な入力ファイル

このプログラムは、外部の入力ファイルを必要としません。 すべての設定パラメータは、プログラム内のグローバル変数として定義されているか、コマンドライン引数として与えられます。

生成される出力ファイル

このプログラムは、ファイルにデータを保存する機能を持っていません。 すべての出力は matplotlib を使用したグラフとして表示されます。

コマンドラインでの使用例 (Usage)

プログラムは、モードに応じて異なる引数を取ります。

python kronig_penney.py
python kronig_penney.py graph (a) (bwidth) (bpot) (k) (Emin) (Emax) (nE)
python kronig_penney.py band (a) (bwidth) (bpot) (kmin) (kmax) (nk)
python kronig_penney.py wf (a) (bwidth) (bpot) (kw) (iLevel) (xwmin) (xwmax) (nxw)

引数の説明:

  • (a): 格子定数 (A, オングストローム)

  • (bwidth): ポテンシャル障壁の幅 (A)

  • (bpot): ポテンシャル障壁の高さ (eV)

  • (k): graph モードでプロットするブロッホ波数 \(k\) (\(\pi/a\) 単位)

  • (Emin): graph モードで走査するエネルギーの下限 (eV)

  • (Emax): graph モードで走査するエネルギーの上限 (eV)

  • (nE): graph モードで走査するエネルギー点数

  • (kmin): band モードでプロットする \(k\) の下限 (\(\pi/a\) 単位)

  • (kmax): band モードでプロットする \(k\) の上限 (\(\pi/a\) 単位)

  • (nk): band モードで走査する \(k\) の点数

  • (kw): wf モードで描画する波動関数のブロッホ波数 \(k\) (\(\pi/a\) 単位)

  • (iLevel): wf モードで描画する波動関数の準位番号（0から始まる整数）

  • (xwmin): wf モードで波動関数を描画する \(x\) 範囲の下限 (A)

  • (xwmax): wf モードで波動関数を描画する \(x\) 範囲の上限 (A)

  • (nxw): wf モードで波動関数を描画する \(x\) の点数

引数を省略した場合、プログラム内のデフォルト値が使用されます。

コマンドラインでの具体的な使用例

1. graph モード

特定の \(k\) 点 (kg = 0.0) におけるKronig-Penney方程式の delta 関数をエネルギー \(E\) の関数としてプロットします。delta = 0 となる点がエネルギー準位を示します。

python kronig_penney.py graph 5.4064 0.5 10.0 0.0 0.0 9.5 51

引数の内訳:

  • graph: モード指定

  • 5.4064: 格子定数 a (A)

  • 0.5: 障壁幅 bwidth (A)

  • 10.0: 障壁高さ bpot (eV)

  • 0.0: \(k\) 点 (kg)

  • 0.0: エネルギー下限 Emin (eV)

  • 9.5: エネルギー上限 Emax (eV)

  • 51: エネルギー点数 nE

実行結果: delta が \(E\) に対してどのように変化するかのグラフが表示されます。グラフ上で delta = 0 のライン（赤の点線）と曲線が交差する点が許容されるエネルギー準位となります。コンソールには、Secant法によって精錬された各エネルギー準位の値と収束情報が表示されます。

2. band モード

ブロッホ波数 \(k\) の範囲にわたるエネルギーバンド構造を計算しプロットします。

python kronig_penney.py band 5.4064 0.5 10.0 -0.5 0.5 21

引数の内訳:

  • band: モード指定

  • 5.4064: 格子定数 a (A)

  • 0.5: 障壁幅 bwidth (A)

  • 10.0: 障壁高さ bpot (eV)

  • -0.5: \(k\) 下限 (kmin, \(\pi/a\) 単位)

  • 0.5: \(k\) 上限 (kmax, \(\pi/a\) 単位)

  • 21: \(k\) 点数 (nk)

実行結果: \(k\) (\(-0.5 \pi/a\) から \(0.5 \pi/a\)) を横軸、エネルギー (eV) を縦軸とするグラフが表示されます。このグラフは、結晶の電子バンド構造を示し、エネルギーバンド（許容エネルギー領域）とバンドギャップ（禁止エネルギー領域）が視覚化されます。コンソールには、各 \(k\) 点で計算されたエネルギー準位が表示されます。

3. wf モード

特定の \(k\) と準位番号における波動関数 \(\Psi(x)\)、その実部・虚部、および電荷密度 \(|\Psi(x)|^2\) をプロットします。同時にポテンシャルの形状も表示されます。

python kronig_penney.py wf 5.4064 0.5 10.0 0.0 0 0.0 16.2192 101

引数の内訳:

  • wf: モード指定

  • 5.4064: 格子定数 a (A)

  • 0.5: 障壁幅 bwidth (A)

  • 10.0: 障壁高さ bpot (eV)

  • 0.0: 波動関数を描画する \(k\) (kw, \(\pi/a\) 単位)

  • 0: 波動関数を描画する準位番号 (iLevel, 0番目の準位)

  • 0.0: \(x\) 範囲下限 (xwmin, A)

  • 16.2192: \(x\) 範囲上限 (xwmax, A, この値は \(3a\) に相当)

  • 101: \(x\) 点数 (nxw)

実行結果: \(x\) 軸に沿ったポテンシャル \(U(x)\) と、対応する波動関数の実部、虚部、および電荷密度 \(|\Psi(x)|^2\) のグラフが表示されます。コンソールには、波動関数の計算に用いられたエネルギー準位、係数、正規化係数などが表示されます。