`transfer_matrix.py` 技術ドキュメント

プログラムの動作

transfer_matrix.py は、1次元のポテンシャル中における粒子の波動関数と透過確率を伝達行列法（Transfer Matrix Method）を用いて計算し、その結果を可視化するPythonプログラムです。主に以下の機能を提供します。

波動関数計算モード (wf): 特定のエネルギー $E_z$ を持つ粒子に対する波動関数 $\Psi(z)$、波数ベクトル $k_z(z)$、およびポテンシャル $U(z)$ と有効質量 $m^*(z)$ の空間分布を計算し、グラフィカルに表示します。
透過確率計算モード (tr): 粒子エネルギー $E$ の関数として透過確率 $T(E)$ を計算し、グラフィカルに表示します。また、任意のエネルギー $E_z$ における波動関数も計算し表示します。
ポテンシャル定義の柔軟性: ポテンシャルは、Excelファイルから読み込むか、プログラム内で定義された多重量子井戸 (Multi-Quantum Well, MQW) モデルを使用できます。
物理的洞察の提供: 伝達行列法を適用することで、複雑なポテンシャル構造における量子粒子の挙動（透過、反射、波動関数の挙動など）を視覚的に理解する手助けをします。

本プログラムは、半導体ヘテロ構造における電子の輸送現象や、量子井戸・障壁構造の設計・解析など、量子力学の基礎的な問題を解決するために利用できます。

原理

本プログラムは、1次元の定常状態シュレーディンガー方程式 $$ -\frac{\hbar^2}{2m^*(z)}\frac{d^2\Psi(z)}{dz^2} + U(z)\Psi(z) = E_z\Psi(z) $$ を**伝達行列法**を用いて数値的に解きます。ここで、$m^*(z)$ は有効質量、$U(z)$ はポテンシャルエネルギー、$E_z$ は粒子の全エネルギー、$z$ は空間座標、$\hbar$ はディラック定数です。

伝達行列法のアルゴリズム:

空間の分割: 計算対象の空間 $z \in [z_{\min}, z_{\max}]$ を微小な $N$ 個の領域（ステップ）に分割します。各領域 $i$ （$i=0, \dots, N-1$）内では、ポテンシャル $U(z)$ と有効質量 $m^*(z)$ は一定であると仮定します。
波数ベクトル: 各領域 $i$ における波数ベクトル $k_i$ は、粒子のエネルギー $E_z$ とポテンシャル $U_i$、有効質量 $m_i^*$ から以下のように計算されます。 $$ k_i = \sqrt{\frac{2 m_i^* (E_z - U_i)e}{\hbar^2}} \cdot 1.0 \times 10^{-10} $$ ここで $e$ は素電荷（プログラム中の物理定数 `e`）、$1.0 \times 10^{-10}$ は単位変換（オングストロームからメートル）のための係数です。プログラムでは複素数 `0.0j` を加えることで、$E_z < U_i$ のポテンシャル障壁下では $k_i$ が純虚数（エバネッセント波）となるように扱います。
波動関数: 各領域 $i$ における波動関数 $\Psi_i(z)$ は、前進波と後進波の線形結合で表されます。 $$ \Psi_i(z) = A_i e^{i k_i z} + B_i e^{-i k_i z} $$ ここで $A_i$ と $B_i$ は振幅係数です。
境界条件と伝達行列: 隣接する領域 $i-1$ と $i$ の境界 $z_i$ において、波動関数と確率電流密度（または波動関数の微分に $1/m^*$ を掛けたもの）の連続性の条件を課します。これにより、領域 $i-1$ の係数 $(A_{i-1}, B_{i-1})$ と領域 $i$ の係数 $(A_i, B_i)$ の間に以下の関係が導かれます。プログラムでは、cal_wf 関数において、 $mk = \frac{m_i^*}{m_{i-1}^*} \frac{k_{i-1}}{k_i}$ $ap = 0.5 \left(1 + mk\right)$ $am = 0.5 \left(1 - mk\right)$ $P = \exp\left(i (k_{i-1} - k_i) z_i\right)$ $Q = \exp\left(i (k_{i-1} + k_i) z_i\right)$ として、 $$ A_i = ap \cdot P \cdot A_{i-1} + am \cdot Q^{-1} \cdot B_{i-1} \\ B_i = am \cdot Q \cdot A_{i-1} + ap \cdot P^{-1} \cdot B_{i-1} $$ と計算を進めます。
透過確率の計算: 最も左端の領域（$i=0$）で $A_0 = 1.0$ (入射波の振幅) および $B_0 = 0.0$ (左からの反射波なし) と初期条件を設定します。その後、上記の式を再帰的に適用して、全領域の係数 $A_i, B_i$ を計算します。最終的に、最も右端の領域（$i=N-1$）の係数 $A_{N-1}$ と $B_{N-1}$ が得られます。透過確率 $T$ は、プログラムでは以下の式で計算しています。 $$ T = \left| \frac{m_{N-1}}{m_0} \frac{k_0}{k_{N-1}} \frac{A_0}{A_{N-1}} \right|^2 $$ ここで、添え字 $0$ は入射側、 $N-1$ は透過側の物理量を示します。

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

本プログラムを実行するには、以下の非標準Pythonライブラリが必要です。

numpy: 数値計算を効率的に行うための基盤ライブラリ。
pandas: データ構造とデータ解析ツールを提供し、Excelファイルの読み書きに使用されます。
matplotlib: グラフ描画ライブラリ。計算結果の可視化に使用されます。

これらのライブラリは、Pythonのパッケージインストーラ pip を用いてインストールできます。

pip install numpy pandas matplotlib

必要な入力ファイル

プログラムの動作モード (mode) が 'tr' または 'wf' で、かつポテンシャルの種類 (pottype) が 'mqw' ではない場合、外部のExcelファイルからポテンシャルデータと有効質量データを読み込みます。

ファイル形式: Microsoft Excel (.xlsx) ファイル
ファイル名: デフォルトでは potential_defect.xlsx ですが、コマンドライン引数で変更可能です。
データ構造: 以下の3列を持つデータが必要です。
- 1列目: z (位置) - 単位はオングストローム (Å) を想定。
- 2列目: U (ポテンシャルエネルギー) - 単位は電子ボルト (eV) を想定。
- 3列目: m* (有効質量) - 単位は電子の静止質量 $m_e$ に対する相対値（$m_e$ を1とした値）を想定。

例: potential_defect.xlsx の内容

z (A)	U (eV)	m* (me)
0.0	0.0	1.0
1.0	0.0	1.0
2.0	10.0	1.0
...	...	...

pottype が 'mqw' に設定されている場合、入力ファイルは不要です。プログラム内の IsInBarrier(), U(), meff() 関数によって多重量子井戸のポテンシャルと有効質量が生成されます。

生成される出力ファイル

プログラムの実行モードによって、以下のファイルが生成される可能性があります。

tr (透過確率計算) モードの場合:
- potential.xlsx: 入力として使用された（または内部で生成された）ポテンシャルと有効質量の空間分布データが保存されます。
  - 1列目: z (A) (位置, オングストローム)
  - 2列目: U (eV) (ポテンシャルエネルギー, 電子ボルト)
  - 3列目: m* (me) (有効質量, 電子の静止質量に対する相対値)
- transmission.xlsx: 粒子エネルギーに対する透過確率のデータが保存されます。
  - 1列目: E (eV) (粒子エネルギー, 電子ボルト)
  - 2列目: T (透過確率, 無次元)
両モード共通:
- プログラムは matplotlib を使用して計算結果をリアルタイムでグラフィカルに表示します。このプロットは、デフォルトではファイルとして保存されません。表示されたウィンドウを手動で閉じるまでプログラムは終了しません。

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

transfer_matrix.py は、コマンドライン引数によって動作モード、ポテンシャルの種類、計算パラメータを指定できます。

基本形式:

python transfer_matrix.py (mode) (pottype) (nz) (Ez0) [Emin] [Emax] [nE]

括弧 () で囲まれた引数は必須、角括弧 [] で囲まれた引数はオプションです。オプション引数は、モードによって必須となる場合があります。

mode: プログラムの動作モードを指定します。
- wf: 波動関数計算モード。
- tr: 透過確率計算モード。
pottype: ポテンシャルの種類を指定します。
- Excelファイル名 (例: potential_defect.xlsx)
- mqw: プログラム内部で定義された多重量子井戸ポテンシャルを使用します。
nz: 空間グリッドの点数 (整数)。$z_{\min}$ から $z_{\max}$ までの区間を nz 個の点で分割します。
Ez0: 初期エネルギー (浮動小数点数)。波動関数の計算や、透過確率モードでの特定のエネルギーにおける波動関数プロットに使用されます。単位はeV。

tr モードの場合の追加引数:

Emin: エネルギー範囲の最小値 (浮動小数点数)。透過確率を計算するエネルギー範囲の下限。単位はeV。
Emax: エネルギー範囲の最大値 (浮動小数点数)。透過確率を計算するエネルギー範囲の上限。単位はeV。
nE: エネルギーグリッドの点数 (整数)。Emin から Emax までの区間を nE 個の点で分割します。

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

例1: 波動関数計算モード (`wf`)

Excelファイルで定義されたポテンシャルに対して、特定のエネルギーにおける波動関数を計算しプロットします。

python transfer_matrix.py wf potential_defect.xlsx 201 0.1

引数の説明:

mode: wf (波動関数計算モード)
pottype: potential_defect.xlsx (Excelファイルからポテンシャルを読み込む)
nz: 201 (z軸上のグリッド点数)
Ez0: 0.1 (波動関数を計算するエネルギー, 0.1 eV)

実行結果の説明: プログラムは、potential_defect.xlsx からポテンシャルと有効質量を読み込み、エネルギー $E_z = 0.1 \, \text{eV}$ における波動関数 $\Psi(z)$ を伝達行列法で計算します。最終的に、以下の4つのサブプロットを含むウィンドウが表示されます。

ポテンシャル $U(z)$ と有効質量 $m^*(z)$ の空間分布。
波数ベクトル $k_z(z)$ の実部と虚部の空間分布。
振幅係数 $A_i$ と $B_i$ の絶対値の空間分布。
波動関数 $\Psi(z)$ の実部、虚部、および確率密度 $|\Psi(z)|^2$ の空間分布。

ターミナルには、zminからzmaxまでの透過確率（Ez0における値）が表示されます。ウィンドウを閉じるには、ターミナルでEnterキーを押します。

例2: 透過確率計算モード (`tr`)

プログラム内部で定義された多重量子井戸ポテンシャルに対して、エネルギー依存の透過確率を計算し、波動関数も合わせてプロットします。

python transfer_matrix.py tr mqw 201 0.1 0.01 9.5 1001

引数の説明:

mode: tr (透過確率計算モード)
pottype: mqw (多重量子井戸ポテンシャルを内部生成)
nz: 201 (z軸上のグリッド点数)
Ez0: 0.1 (波動関数を計算するエネルギー, 0.1 eV)
Emin: 0.01 (透過確率計算の最小エネルギー, 0.01 eV)
Emax: 9.5 (透過確率計算の最大エネルギー, 9.5 eV)
nE: 1001 (エネルギーグリッドの点数)

実行結果の説明: プログラムは、内部定義された多重量子井戸ポテンシャルと有効質量を用いて、エネルギー範囲 $0.01 \, \text{eV}$ から $9.5 \, \text{eV}$ までを $1001$ 点で分割し、各エネルギーでの透過確率を計算します。また、$E_z = 0.1 \, \text{eV}$ における波動関数も計算します。以下の2つのファイルが生成されます。

potential.xlsx: 生成されたポテンシャルと有効質量のデータ。
transmission.xlsx: エネルギーと透過確率のデータ。

最終的に、以下の4つのサブプロットを含むウィンドウが表示されます。

ポテンシャル $U(z)$ と有効質量 $m^*(z)$ の空間分布。
エネルギー $E_z = 0.1 \, \text{eV}$ における波数ベクトル $k_z(z)$ の実部と虚部の空間分布。
エネルギー $E$ に対する透過確率 $T(E)$ のグラフ。
エネルギー $E_z = 0.1 \, \text{eV}$ における波動関数 $\Psi(z)$ の実部、虚部、および確率密度 $|\Psi(z)|^2$ の空間分布。

ターミナルには、各エネルギーにおける透過確率の一部が表形式で出力されます。ウィンドウを閉じるには、ターミナルでEnterキーを押します。

transfer_matrix.py 技術ドキュメント