tkfermi_integral.py ライブラリの技術ドキュメント

tkfermi_integral.py ドキュメント

1. ライブラリの機能や目的

tkfermi_integral.py は、高速かつ堅牢なフェルミ・ディラック積分 $F_r(\eta)$ を計算するためのPythonライブラリです。半導体物理学、固体物理学、統計力学など、キャリア統計や電子状態の分析が不可欠な分野で広く利用されます。

主な機能:

一般化されたフェルミ・ディラック積分 $F_r(\eta)$ の計算: 任意の次数 $r$ と化学ポテンシャル $\eta$ に対する積分値を効率的に計算します。
高精度と広範囲な対応:
- 退化極限 ($\eta$ が非常に大きい場合): Sommerfeld展開を用いて高精度かつ高速に計算します。
- 非退化極限 ($\eta$ が非常に小さい場合): 交代級数展開を用いて効率的に計算します。
- 中間領域: 適応型数値積分（scipy.integrate.quad）を内部的に利用し、精度と速度のバランスを取ります。
特定の整数次数での高速化: $r = -1, 0, 1$ の場合、解析解または特殊関数（scipy.special.expit, numpy.log1p, scipy.special.spence）を用いて非常に高速に計算します。
フェルミ・ディラック分布関数: 電子、ホール、ドーパントの電離分数など、フェルミ・ディラック分布関数自体も提供します。
キャリア密度計算のヘルパー関数: $F_{1/2}(\eta)$ を用いた電子・正孔密度の計算を簡略化する関数を提供します。
キャッシュ機構: functools.lru_cache を用いて、同じ引数での再計算を避け、性能を向上させます。

解決する課題:

フェルミ・ディラック積分は、閉じた形式での解析解が存在しない場合が多く、特に広いパラメータ範囲で高精度かつ効率的に計算することは困難です。本ライブラリは、異なる物理的極限（退化、非退化、中間）に合わせた最適なアルゴリズムを自動的に選択することで、この課題を解決し、信頼性の高い結果を提供します。

フェルミ・ディラック積分 $F_r(\eta)$ の定義は以下の通りです。 $$ F_r(\eta) = \frac{1}{\Gamma(r+1)} \int_0^\infty \frac{x^r}{1 + \exp(x-\eta)} dx $$ ただし、このライブラリでは `FermiIntegral_fast` など主要な関数は**ガンマ関数による正規化を行わない**積分 $$ I_r(\eta) = \int_0^\infty \frac{x^r}{1 + \exp(x-\eta)} dx $$ を計算します。これは、関数 `electron_density` や `hole_density` で利用される定数 $kF12 = 2/\sqrt{\pi}$ と組み合わせて、物理的な濃度計算に直接適用できるように設計されています。例えば、$I_{1/2}(\eta)$ は $\Gamma(3/2) F_{1/2}(\eta)$ に相当します。

2. importする方法

このライブラリを他のPythonプログラムからインポートするには、以下のいずれかの方法を使用できます。

ライブラリ全体をインポートする場合:

import tkfermi_integral

# フェルミ・ディラック積分F_1/2(eta)を計算
result = tkfermi_integral.FermiIntegral_half(5.0)
print(result)

特定の関数のみをインポートする場合:

from tkfermi_integral import FermiIntegral_fast, electron_density

# 一般的なフェルミ・ディラック積分を計算
result_F_1_5 = FermiIntegral_fast(2.0, 1.5)
print(result_F_1_5)

# 電子密度を計算
Nc = 1e19
eta_c = 3.0
n = electron_density(Nc, eta_c)
print(n)

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

tkfermi_integral.py は以下の非標準ライブラリに依存しています。

numpy: 数値計算、特に配列操作に必要です。
scipy: 科学技術計算、特に数値積分 (scipy.integrate.quad) や特殊関数 (scipy.special.expit, scipy.special.log_expit, scipy.special.zeta, scipy.special.spence) に必要です。

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

pip install numpy scipy

4. importできる変数と関数

変数

xlim_exp (float): フェルミ・ディラック分布関数を計算する際に、exp() 関数の引数をクリップするための閾値です。この値を超える引数は、数値的な安定性を保つためにこの閾値にクリップされます。デフォルト値は 40.0 です。
kF12 (float): $2/\sqrt{\pi}$ の定数。電子・正孔密度計算の正規化係数として使用されます。

関数

フェルミ・ディラック分布関数

fermi_dirac_e_npvector(x)
- 動作: フェルミ・ディラック分布関数 $f(x) = 1 / (1 + \exp(x))$ を、NumPy配列 x に対してベクトル化された形式で計算します。内部的に xlim_exp を用いて数値的な安定性を確保します。
- 引数:
  - x (ndarray): エネルギーを $kT$ で割った無次元量。
- 戻り値:
  - (ndarray): x に対応するフェルミ・ディラック分布関数の値。
fermi_dirac_e(x, eta, g=1.0)
- 動作: 一般的なフェルミ・ディラック分布関数 $f(x) = 1 / (1 + g \cdot \exp(x - \eta))$ を計算します。scipy.special.expit を使用して数値的な安定性を高めています。 $$ f(x, \eta, g) = \frac{1}{1 + g \cdot \exp(x - \eta)} = \text{expit}((\eta - x) - \ln(g)) $$
- 引数:
  - x (float or ndarray): エネルギーを $kT$ で割った無次元量。
  - eta (float or ndarray): フェルミ準位を $kT$ で割った無次元量（化学ポテンシャル）。
  - g (float, optional): スピン縮退度などの係数。デフォルトは 1.0。
- 戻り値:
  - (float or ndarray): フェルミ・ディラック分布関数の値。
log_fermi_dirac_e(x, eta, g=1.0)
- 動作: fermi_dirac_e と同じ定義のフェルミ・ディラック分布関数の自然対数 $\ln(f(x))$ を、scipy.special.log_expit を用いて安定に計算します。
- 引数: fermi_dirac_e と同じ。
- 戻り値:
  - (float or ndarray): フェルミ・ディラック分布関数の自然対数の値。
fermi_dirac_h(x, eta, g=1.0)
- 動作: 正孔（空孔）の存在確率 $1 - f(x)$ を計算します。これも scipy.special.expit を使用して安定に計算されます。 $$ 1 - f(x, \eta, g) = \frac{1}{1 + g \cdot \exp(\eta - x)} = \text{expit}((x - \eta) - \ln(g)) $$
- 引数: fermi_dirac_e と同じ。
- 戻り値:
  - (float or ndarray): 正孔の存在確率。
log_fermi_dirac_h(x, eta, g=1.0)
- 動作: fermi_dirac_h と同じ定義の正孔存在確率の自然対数 $\ln(1 - f(x))$ を、scipy.special.log_expit を用いて安定に計算します。
- 引数: fermi_dirac_e と同じ。
- 戻り値:
  - (float or ndarray): 正孔の存在確率の自然対数の値。
ionized_acceptor_frac(eta, xA, g=1.0)
- 動作: アクセプタが電離している割合を計算します。fermi_dirac_e(xA, eta, g=g) と等価です。
- 引数:
  - eta (float or ndarray): フェルミ準位を $kT$ で割った無次元量。
  - xA (float or ndarray): アクセプタ準位を $kT$ で割った無次元量。
  - g (float, optional): アクセプタの縮退度。デフォルトは 1.0。
- 戻り値:
  - (float or ndarray): 電離しているアクセプタの割合。
ionized_donor_frac(eta, xD, g=1.0)
- 動作: ドナーが電離している割合を計算します。fermi_dirac_h(xD, eta, g=g) と等価です。
- 引数:
  - eta (float or ndarray): フェルミ準位を $kT$ で割った無次元量。
  - xD (float or ndarray): ドナー準位を $kT$ で割った無次元量。
  - g (float, optional): ドナーの縮退度。デフォルトは 1.0。
- 戻り値:
  - (float or ndarray): 電離しているドナーの割合。

フェルミ・ディラック積分

FermiIntegral_fast(eta, r, *, epsabs=1.0e-10, epsrel=1.0e-8, limit=50, use_cache=True)
- 動作: 一般化されたフェルミ・ディラック積分 $I_r(\eta)$ を計算します。 $$ I_r(\eta) = \int_0^\infty \frac{x^r}{1 + \exp(x - \eta)} dx $$ この関数は、引数 `eta` と `r` に応じて、Sommerfeld展開、交代級数展開、または数値積分を自動的に選択します。特に、整数次数 $r = -1, 0, 1$ の場合は、解析解または特殊関数を利用して非常に高速に計算されます。
  - $r = -1$: $I_{-1}(\eta) = \frac{1}{1 + \exp(-\eta)} = \text{expit}(\eta)$
  - $r = 0$: $I_0(\eta) = \ln(1 + \exp(\eta))$
  - $r = 1$: $I_1(\eta) = -Li_2(-\exp(\eta))$ (ポリ対数関数 $Li_2$)
    - 注記: 実装では scipy.special.spence(1.0 + np.exp(eta)) を用いています。$\eta > 100$ の場合は漸近展開 $0.5 \eta^2 + \pi^2/6$ を使用します。
- 引数:
  - eta (float): 化学ポテンシャルを $kT$ で割った無次元量。
  - r (float): 積分の次数。
  - epsabs (float, optional): 数値積分における絶対誤差許容値。デフォルトは 1.0e-10。
  - epsrel (float, optional): 数値積分における相対誤差許容値。デフォルトは 1.0e-8。
  - limit (int, optional): 数値積分における区分数上限。デフォルトは 50。
  - use_cache (bool, optional): 計算結果をキャッシュするかどうか。デフォルトは True。
- 戻り値:
  - (float): フェルミ・ディラック積分の値 $I_r(\eta)$。
FermiIntegral_half(eta, *, epsabs=1.0e-10, epsrel=1.0e-8, limit=50)
- 動作: $r=0.5$ のフェルミ・ディラック積分 $I_{1/2}(\eta)$ を計算するラッパー関数です。
- 引数: FermiIntegral_fast と同じ（ただし r は固定）。
- 戻り値:
  - (float): $I_{1/2}(\eta)$ の値。
FermiIntegral_3half(eta, *, epsabs=1.0e-10, epsrel=1.0e-8, limit=50)
- 動作: $r=1.5$ のフェルミ・ディラック積分 $I_{3/2}(\eta)$ を計算するラッパー関数です。
- 引数: FermiIntegral_fast と同じ（ただし r は固定）。
- 戻り値:
  - (float): $I_{3/2}(\eta)$ の値。

キャリア密度計算ヘルパー

electron_density(Nc, eta_c, *, epsabs=1.0e-10, epsrel=1.0e-8, limit=50)
- 動作: 伝導帯内の電子密度 $n$ を計算します。標準的な放物線バンドを仮定し、伝導帯有効状態密度 $Nc$ と無次元化学ポテンシャル $\eta_c$ を用います。ライブラリの定義する非正規化積分 $I_{1/2}(\eta)$ を使用すると、電子密度は以下の式で与えられます。 $$ n = N_c \cdot \frac{2}{\sqrt{\pi}} \cdot I_{1/2}(\eta_c) $$ ここで $\eta_c = (E_f - E_c) / (kT)$ です。
- 引数:
  - Nc (float): 伝導帯有効状態密度（例: [cm^-3] または [m^-3])。
  - eta_c (float or array-like): 電子の無次元化学ポテンシャル。
  - epsabs, epsrel, limit は FermiIntegral_half に渡されます。
- 戻り値:
  - (float or ndarray): 電子密度。Nc と同じ単位。
hole_density(Nv, eta_v, *, epsabs=1.0e-10, epsrel=1.0e-8, limit=50)
- 動作: 価電子帯内の正孔密度 $p$ を計算します。標準的な放物線バンドを仮定し、価電子帯有効状態密度 $Nv$ と無次元化学ポテンシャル $\eta_v$ を用います。ライブラリの定義する非正規化積分 $I_{1/2}(\eta)$ を使用すると、正孔密度は以下の式で与えられます。 $$ p = N_v \cdot \frac{2}{\sqrt{\pi}} \cdot I_{1/2}(-\eta_v) $$ ここで $\eta_v = (E_f - E_v) / (kT)$ です。-eta_v は、正孔のフェルミ準位が価電子帯上限から下向きに測られることを示しています。
- 引数:
  - Nv (float): 価電子帯有効状態密度（例: [cm^-3] または [m^-3])。
  - eta_v (float or array-like): 正孔の無次元化学ポテンシャル。
  - epsabs, epsrel, limit は FermiIntegral_half に渡されます。
- 戻り値:
  - (float or ndarray): 正孔密度。Nv と同じ単位。

5. main scriptとして実行したときの動作

tkfermi_integral.py を直接実行 (python tkfermi_integral.py) すると、ライブラリの簡単な自己テストが実行されます。これは、$F_{1/2}(\eta)$ と $F_{3/2}(\eta)$ のいくつかの eta 値での計算結果、およびそれらの積分を用いた電子・正孔密度の計算例を表示します。また、キャッシュの利用状況も表示されます。

実行時の出力例:

Fermi–Dirac integral test

r = 1/2
  eta=-10.00  F1/2=3.09709033e-05
  eta=-3.00  F1/2=6.55931086e-02
  eta=-1.00  F1/2=3.13289069e-01
  eta= 0.00  F1/2=6.78216773e-01
  eta= 2.00  F1/2=2.16431989e+00
  eta= 5.00  F1/2=7.11475753e+00
  eta=10.00  F1/2=2.45717327e+01

r = 3/2
  eta=-10.00  F3/2=3.09709033e-05  n=3.48625299e+14  p=3.48625299e+23
  eta=-3.00  F3/2=1.31186217e-01  n=7.39174549e+17  p=1.65215017e+22
  eta=-1.00  F3/2=4.70775537e-01  n=3.53503102e+18  p=3.70588691e+20
  eta= 0.00  F3/2=9.02058993e-01  n=7.62002302e+18  p=7.62002302e+18
  eta= 2.00  F3/2=3.13627702e+00  n=2.64893706e+19  p=4.70775537e+17
  eta= 5.00  F3/2=1.23354728e+01  n=1.04253160e+20  p=1.31186217e+16
  eta=10.00  F3/2=4.85192131e+01  n=4.09919018e+20  p=3.09709033e+14

cache info:
Hits=14, Misses=14, Maxsize=1024, Currsize=14