Kronig-Penney Model Simulation

プログラムの動作

kronig_penney.py は、クロン・ペニー（Kronig-Penney）モデルを用いて、1次元結晶における電子のエネルギーバンド構造と波動関数を計算し、その結果を視覚的に表示するPythonプログラムです。このプログラムは以下の主要な機能を提供します。

1. エネルギー対デルタ関数プロット (graph モード): 特定の波数 \(k\) におけるクロン・ペニー方程式の左辺（プログラム中の delta 関数）を電子エネルギー \(E\) の関数としてプロットします。このグラフで delta = 0 となる点が、許容されるエネルギー準位に対応します。

2. 電子バンド構造プロット (band モード): 逆格子空間の波数 \(k\) を走査し、各 \(k\) における許容エネルギー準位を計算して、電子のエネルギーバンド構造（\(E\)-\(k\) 関係）をプロットします。これにより、エネルギーの許可帯と禁止帯が視覚化されます。

3. 波動関数プロット (wf モード): 特定の波数 \(k\) と特定のエネルギー準位におけるBloch波動関数 \(\Psi(x)\) を実部、虚部、および電荷密度 \(|\Psi(x)|^2\) としてプロットします。同時に、周期的なポテンシャル障壁の形状も表示されます。

このプログラムは、周期的なポテンシャル中で移動する電子の量子力学的振る舞いを理解するための基本的なツールとして機能します。

原理

kronig_penney.py は、クロン・ペニーモデルの物理原理に基づいています。このモデルは、周期的に配置された矩形ポテンシャル障壁を持つ1次元格子における電子の挙動を記述します。

モデルでは、周期 \(a\) の単位格子が、幅 \(w\) のポテンシャル井戸（ポテンシャルエネルギー \(0\)）と、幅 \(b\) のポテンシャル障壁（ポテンシャルエネルギー \(V_0\)）で構成されると仮定します（\(a = w+b\)）。

電子の波動関数は、それぞれの領域（井戸内と障壁内）で異なる形式を持ち、シュレーディンガー方程式の解として与えられます。また、Blochの定理により、全波動関数は \(\Psi(x) = e^{iKx} u(x)\) の形をとります。ここで \(K\) はBloch波の波数、\(u(x)\) は格子の周期 \(a\) を持つ周期関数です。

境界条件（波動関数とその導関数が領域の境界で連続であること）とBlochの定理を適用することで、許容されるエネルギー準位 \(E\) を決定する超越方程式が得られます。プログラム中では、この方程式の変形が cal_delta 関数によって計算されます。

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

ここで、\(\alpha\) と \(\beta\) は次のように定義されます。

\[ \alpha = \frac{\sqrt{2m_e E}}{\hbar}, \quad \beta = \frac{\sqrt{2m_e (V_0 - E)}}{\hbar} \]

\(m_e\) は電子の質量、\(\hbar\) はディラック定数です。プログラムでは、エネルギー \(E\) と障壁高さ \(V_0\) は電子ボルト (eV) 単位で入力されるため、計算時に電子の電荷 \(e\) を掛けてジュール (J) に変換しています。また、井戸幅 \(w\) と障壁幅 \(b\) はオングストローム (A) 単位で入力されるため、計算時に \(1.0 \times 10^{-10}\) を掛けてメートル (m) に変換しています。

プログラム内の変数 k は \(Ka/(2\pi)\) に対応し、ka は \(Ka\) を表します。find_Elist 関数では、Secant法 (refine_E 関数) を用いて、cal_delta の値がゼロとなる \(E\) を見つけることで、エネルギー固有値を探索します。

波動関数の計算では、上記の方程式を満たすエネルギー \(E\) に対し、境界条件から得られる係数（プログラム中の ci リスト）を用いて波動関数を構築します。波動関数は最終的に正規化されます。

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

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

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

• matplotlib: グラフ描画

pprint および csv もインポートされていますが、現在のプログラムの実行パスでは直接使用されていません。

これらのライブラリは、Pythonのパッケージ管理システム pip を使ってインストールできます。

pip install numpy matplotlib

必要な入力ファイル

このプログラムは、外部の入力ファイルを必要としません。 全てのパラメータはプログラム内のデフォルト値として設定されており、コマンドライン引数によって上書きすることができます。

生成される出力ファイル

このプログラムは、計算結果をファイルとして保存しません。 全ての出力は matplotlib を利用したグラフとして、実行時に画面に表示されます。

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

プログラムの基本的な使用方法は以下の通りです。 引数はスペースで区切って指定します。括弧 () で囲まれた変数はオプションです。

Usage: Variables in () are optional
  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)
     ex: python kronig_penney.py graph 5.4064 0.5 10.0 0.0 0.0 9.5 51
     ex: python kronig_penney.py band 5.4064 0.5 10.0 -0.5 0.5 21
     ex: python kronig_penney.py wf 5.4064 0.5 10.0 0.0 0 0.0 16.2192 101

引数の詳細:

• mode: 実行モード (graph, band, wf のいずれか)。デフォルトは graph。

• a: 単位格子の周期 [A]

• bwidth: ポテンシャル障壁の幅 [A]

• bpot: ポテンシャル障壁の高さ [eV]

graph モード固有の引数:

• k: プロットする波数 \(K a / (2\pi)\) [単位なし, -0.5 から 0.5 の範囲]

• Emin: エネルギー走査開始点 [eV]

• Emax: エネルギー走査終了点 [eV]

• nE: エネルギー走査点数

band モード固有の引数:

• kmin: 波数 \(K a / (2\pi)\) 走査開始点 [単位なし, -0.5 から 0.5 の範囲]

• kmax: 波数 \(K a / (2\pi)\) 走査終了点 [単位なし, -0.5 から 0.5 の範囲]

• nk: 波数走査点数

wf モード固有の引数:

• kw: 波動関数を描画する波数 \(K a / (2\pi)\) [単位なし, -0.5 から 0.5 の範囲]

• iLevel: 描画するエネルギー準位のインデックス（0から始まる）

• xwmin: 波動関数を描画するX軸の最小値 [A]

• xwmax: 波動関数を描画するX軸の最大値 [A]

• nxw: 波動関数を描画するX軸の点数

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

1. エネルギー対デルタ関数のプロット (graph モード)

波数 \(K a / (2\pi) = 0.0\) の場合の、クロン・ペニー方程式の左辺 (delta) をエネルギー \(E\) の関数としてプロットします。

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

• a: 5.4064 A (単位格子の周期)

• bwidth: 0.5 A (障壁の幅)

• bpot: 10.0 eV (障壁の高さ)

• k: 0.0 (プロットする波数 \(K a / (2\pi)\))

• Emin: 0.0 eV (エネルギー走査開始点)

• Emax: 9.5 eV (エネルギー走査終了点)

• nE: 51 (エネルギー走査点数)

実行結果の説明: このコマンドを実行すると、エネルギー \(E\) を横軸に、クロン・ペニー方程式の左辺 (delta) を縦軸にしたグラフが表示されます。グラフが横軸 (\(delta = 0\)) と交差する点が、許容されるエネルギー準位（固有値）に対応します。これにより、特定の波数における電子のエネルギー準位を視覚的に確認できます。

2. 電子バンド構造のプロット (band モード)

波数 \(K a / (2\pi)\) を -0.5 から 0.5 の範囲で走査し、各波数における許容エネルギー準位をプロットすることで、エネルギーバンド構造図を生成します。

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

• a: 5.4064 A

• bwidth: 0.5 A

• bpot: 10.0 eV

• kmin: -0.5 (波数 \(K a / (2\pi)\) 走査開始点)

• kmax: 0.5 (波数 \(K a / (2\pi)\) 走査終了点)

• nk: 21 (波数走査点数)

実行結果の説明: このコマンドを実行すると、波数 \(K a / (2\pi)\) を横軸に、エネルギー \(E\) を縦軸にしたバンド図が表示されます。この図は、電子が取りうるエネルギー状態が特定の帯域（許可帯）に集中し、その間に電子が占めることのできないエネルギー領域（禁止帯）が存在することを示します。

3. 波動関数のプロット (wf モード)

波数 \(K a / (2\pi) = 0.0\) における、基底状態（0番目の準位）の波動関数 \(\Psi(x)\) をプロットします。波動関数の実部、虚部、および電荷密度 \(|\Psi(x)|^2\) が表示されます。

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

• a: 5.4064 A

• bwidth: 0.5 A

• bpot: 10.0 eV

• kw: 0.0 (波動関数を描画する波数 \(K a / (2\pi)\))

• iLevel: 0 (描画するエネルギー準位のインデックス、0が基底状態)

• xwmin: 0.0 A (X軸の最小値)

• xwmax: 16.2192 A (X軸の最大値、約 \(3a\) 程度)

• nxw: 101 (X軸の点数)

実行結果の説明: このコマンドを実行すると、X軸座標を横軸に、波動関数の値（実部、虚部）と電荷密度を縦軸にしたグラフが表示されます。背景には周期的なポテンシャル障壁の形状も表示され、波動関数がポテンシャルの各領域でどのように変化するか、そしてBlochの定理に従って周期的な振る舞いを示すかが視覚的に確認できます。電荷密度 \(|\Psi(x)|^2\) は、電子が存在する確率が高い領域を示します。