EF-bisection.py 技術ドキュメント

プログラムの動作

EF-bisection.py は、半導体中のフェルミエネルギー $E_F$ を2分法（bisection method）を用いて計算するPythonプログラムです。このプログラムの主な目的は、半導体の伝導帯、価電子帯、ドナー準位、アクセプター準位、およびそれぞれの有効状態密度やキャリア密度が与えられた条件下で、電荷中性条件を満たすフェルミエネルギーを数値的に探索することです。

プログラムは以下の主要な機能を提供します。

半導体物理定数とパラメータの定義: 電気素量、ボルツマン定数、温度、伝導帯下端・価電子帯上端エネルギー、有効状態密度、ドナー・アクセプター準位とその密度などの物理定数および半導体パラメータをコード内に設定します。
キャリア密度とイオン化準位密度の計算: フェルミ・ディラック分布関数に基づき、電子密度、正孔密度、イオン化ドナー密度、イオン化アクセプター密度を計算する関数を実装しています。
電荷中性条件の適用: これらの密度を用いて電荷中性条件（$n + N_A^- = p + N_D^+$）から得られる残余電荷 $\Delta Q(E_F)$ を定義します。
2分法による数値解法: $\Delta Q(E_F) = 0$ となる $E_F$ を見つけるために、初期範囲内で2分法を適用し、指定された精度または最大繰り返し数に達するまで $E_F$ を探索します。
途中経過と結果の出力: 各繰り返しステップでの $E_F$ の値、$\Delta Q$ の値、および各キャリア密度の値を標準出力に表示し、収束した $E_F$ を報告します。

原理

このプログラムは、半導体における電荷中性条件を満たすフェルミエネルギー $E_F$ を見つけるために、以下の物理式と2分法アルゴリズムを使用します。

1. 電荷中性条件

半導体中の電荷中性条件は、電子、正孔、イオン化ドナー、およびイオン化アクセプターの電荷バランスによって与えられます。

\[n + N_A^- = p + N_D^+\]

ここで、

$n$: 電子密度
$N_A^-$: イオン化アクセプター密度
$p$: 正孔密度
$N_D^+$: イオン化ドナー密度

この式を並べ替えると、フェルミエネルギー $E_F$ の関数として残余電荷 $\Delta Q(E_F)$ を定義できます。

\[\Delta Q(E_F) = n(E_F) + N_A^-(E_F) - p(E_F) - N_D^+(E_F) = 0\]

この $\Delta Q(E_F) = 0$ となる $E_F$ を探索します。

2. 各項の物理式

各密度は以下の式で計算されます。

フェルミ・ディラック分布関数: 特定のエネルギー準位 $E$ が電子によって占有される確率を表します。 $$f(E, E_F, T) = \frac{1}{\exp\left(\frac{E - E_F}{k_B T}\right) + 1}$$ ここで、$k_B$ はボルツマン定数、$T$ は絶対温度です。プログラムでは $ekBT = e / (k_B T)$ という定数を利用して $\frac{1}{k_B T}$ を表現しています（$e$ は電子の素電荷ですが、ここでは単位変換の役割を果たしています）。
電子密度: 伝導帯中の電子密度は、伝導帯下端 $E_c$ と伝導帯有効状態密度 $N_c$ を用いて計算されます。 $$n = N_c \exp\left(-\frac{E_c - E_F}{k_B T}\right)$$
正孔密度: 価電子帯中の正孔密度は、価電子帯上端 $E_v$ と価電子帯有効状態密度 $N_v$ を用いて計算されます。 $$p = N_v \exp\left(-\frac{E_F - E_v}{k_B T}\right)$$
イオン化ドナー密度: ドナー準位 $E_D$ とドナー密度 $N_D$ を用いて計算されます。ドナー準位の電子が励起されて伝導帯へ移動した場合、ドナーは正にイオン化します。 $$N_D^+ = N_D (1 - f(E_D, E_F, T))$$
イオン化アクセプター密度: アクセプター準位 $E_A$ とアクセプター密度 $N_A$ を用いて計算されます。アクセプター準位が価電子帯から電子を受け取った場合、アクセプターは負にイオン化します。 $$N_A^- = N_A f(E_A, E_F, T)$$

3. 2分法（Bisection Method）

2分法は、連続関数 $g(x)=0$ の根を見つけるための数値解析アルゴリズムです。このプログラムでは $g(E_F) = \Delta Q(E_F)$ となります。

初期区間の設定: 根が存在すると思われる区間 $[E_{F, \min}, E_{F, \max}]$ を設定します。この区間では、$g(E_{F, \min})$ と $g(E_{F, \max})$ の符号が異なる必要があります（すなわち、$g(E_{F, \min}) \cdot g(E_{F, \max}) < 0$）。
中点の計算: 区間の中点 $E_{F, \text{half}} = (E_{F, \min} + E_{F, \max}) / 2$ を計算します。
区間の更新: $g(E_{F, \text{half}})$ の符号を評価します。
- もし $g(E_{F, \min}) \cdot g(E_{F, \text{half}}) < 0$ なら、根は $[E_{F, \min}, E_{F, \text{half}}]$ の間にあるため、$E_{F, \max}$ を $E_{F, \text{half}}$ に更新します。
- そうでなければ（$g(E_{F, \text{half}}) \cdot g(E_{F, \max}) < 0$ の場合）、根は $[E_{F, \text{half}}, E_{F, \max}]$ の間にあるため、$E_{F, \min}$ を $E_{F, \text{half}}$ に更新します。
収束判定: 区間の幅 $|E_{F, \max} - E_{F, \min}|$ が許容誤差 eps より小さくなるまで、または最大繰り返し回数 nmaxiter に達するまで手順2と3を繰り返します。

この方法により、指定された精度でフェルミエネルギー $E_F$ の近似解を得ることができます。

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

このプログラムは、数値計算のために numpy ライブラリを使用しています。csv および sys ライブラリもインポートされていますが、現在のコードでは直接使用されていません。

numpy がインストールされていない場合は、以下の pip コマンドでインストールできます。

pip install numpy

必要な入力ファイル

このプログラムは、外部の入力ファイルを必要としません。必要なすべての物理定数、半導体パラメータ、および2分法の設定値は、EF-bisection.py のソースコード内に直接ハードコードされています。

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

このプログラムは、外部ファイルを出力として生成しません。計算結果（各繰り返しステップでのフェルミエネルギーの推定値、残余電荷、各キャリア密度、および最終的な収束結果）はすべて標準出力（コンソール）に表示されます。

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

EF-bisection.py はコマンドライン引数を取らずに直接実行されます。

python EF-bisection.py

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

以下のコマンドを実行すると、プログラムが起動し、標準出力にフェルミエネルギーの計算過程と結果が表示されます。

python EF-bisection.py

実行結果の例

（上記のプログラムを実行した場合の出力例。実際の出力値は設定されたパラメータによって異なります。）

Solution of EF by bisection method

  EFmin =   0.00000000  dQmin =  1.0e+15
  EFmax =   1.00000000  dQmax = -1.0e+16
  Iter 0: EFhalf =   0.50000000  dQhalf = -1.0e+16
     Ne=1.7589e+04  Nh=2.2618e+10  NA-=1.5034e+10  ND+=1.0000e+16  dQ=-1.0000e+16
  Iter 1: EFhalf =   0.25000000  dQhalf = -1.0e+16
     Ne=4.1487e+07  Nh=1.1784e+14  NA-=1.4720e+10  ND+=1.0000e+16  dQ=-1.0000e+16
  Iter 2: EFhalf =   0.12500000  dQhalf = -1.0e+16
     Ne=9.7891e+09  Nh=6.1362e+15  NA-=1.0183e+10  ND+=1.0000e+16  dQ=-1.0000e+16
  Iter 3: EFhalf =   0.06250000  dQhalf = -1.0e+16
     Ne=2.3090e+12  Nh=3.1979e+16  NA-=1.0182e+10  ND+=1.0000e+16  dQ=-1.0000e+16
  Iter 4: EFhalf =   0.03125000  dQhalf =  1.0e+15
     Ne=5.4468e+14  Nh=1.6669e+17  NA-=1.0182e+10  ND+=1.0000e+16  dQ=1.0182e+15
  Iter 5: EFhalf =   0.04687500  dQhalf =  4.0e+14
     Ne=7.3982e+13  Nh=6.4485e+16  NA-=1.0182e+10  ND+=1.0000e+16  dQ=4.0903e+14
  Iter 6: EFhalf =   0.05468750  dQhalf =  1.0e+14
     Ne=2.0487e+13  Nh=4.5126e+16  NA-=1.0182e+10  ND+=1.0000e+16  dQ=1.3544e+14
  Iter 7: EFhalf =   0.05859375  dQhalf =  5.0e+13
     Ne=1.1667e+13  Nh=3.8290e+16  NA-=1.0182e+10  ND+=1.0000e+16  dQ=5.5601e+13
  Iter 8: EFhalf =   0.06054688  dQhalf =  3.0e+13
     Ne=7.9719e+12  Nh=3.5042e+16  NA-=1.0182e+10  ND+=1.0000e+16  dQ=2.9649e+13
  Iter 9: EFhalf =   0.06152344  dQhalf =  1.0e+13
     Ne=6.5921e+12  Nh=3.3547e+16  NA-=1.0182e+10  ND+=1.0000e+16  dQ=1.5034e+13
  Iter 10: EFhalf =   0.06201172  dQhalf =  6.0e+12
     Ne=5.9922e+12  Nh=3.2842e+16  NA-=1.0182e+10  ND+=1.0000e+16  dQ=6.2323e+12
  Iter 11: EFhalf =   0.06225586  dQhalf =  2.0e+12
     Ne=5.7001e+12  Nh=3.2490e+16  NA-=1.0182e+10  ND+=1.0000e+16  dQ=2.6828e+12
  Iter 12: EFhalf =   0.06237793  dQhalf =  1.0e+12
     Ne=5.5501e+12  Nh=3.2317e+16  NA-=1.0182e+10  ND+=1.0000e+16  dQ=1.1578e+12
  Iter 13: EFhalf =   0.06243896  dQhalf =  5.0e+11
     Ne=5.4754e+12  Nh=3.2230e+16  NA-=1.0182e+10  ND+=1.0000e+16  dQ=5.0069e+11
  Iter 14: EFhalf =   0.06246948  dQhalf =  2.0e+11
     Ne=5.4382e+12  Nh=3.2186e+16  NA-=1.0182e+10  ND+=1.0000e+16  dQ=2.1643e+11
  Iter 15: EFhalf =   0.06248474  dQhalf =  9.0e+10
     Ne=5.4196e+12  Nh=3.2165e+16  NA-=1.0182e+10  ND+=1.0000e+16  dQ=9.3562e+10
  Iter 16: EFhalf =   0.06249237  dQhalf =  4.0e+10
     Ne=5.4103e+12  Nh=3.2154e+16  NA-=1.0182e+10  ND+=1.0000e+16  dQ=4.0430e+10
  Iter 17: EFhalf =   0.06249619  dQhalf =  1.0e+10
     Ne=5.4057e+12  Nh=3.2149e+16  NA-=1.0182e+10  ND+=1.0000e+16  dQ=1.7483e+10
  Iter 18: EFhalf =   0.06249809  dQhalf =  8.0e+09
     Ne=5.4034e+12  Nh=3.2146e+16  NA-=1.0182e+10  ND+=1.0000e+16  dQ=7.5683e+09
  Iter 19: EFhalf =   0.06249904  dQhalf =  3.0e+09
     Ne=5.4022e+12  Nh=3.2145e+16  NA-=1.0182e+10  ND+=1.0000e+16  dQ=3.2770e+09
  Success: Convergence reached at EF = 0.06249904

上記の出力例では、プログラムが2分法を繰り返し実行し、各ステップでのフェルミエネルギーの候補 EFhalf と残余電荷 dQhalf、および各キャリア密度を表示しています。最終的に、設定された収束条件 eps を満たした EF = 0.06249904 で計算が成功裏に終了したことが示されています。