金属のフェルミ準位の温度依存性の計算

プログラムの動作

EF-T-metal.py は、金属中の自由電子モデルに基づき、フェルミ準位（$E_F$）の温度依存性を数値的に計算し、その結果をソマーフェルト展開による近似式と比較するPythonプログラムです。主な機能は以下の通りです。

所定の温度範囲（0Kから4000K）において、フェルミ・ディラック統計と3次元自由電子ガスの状態密度を用いて、電子数密度からフェルミ準位を算出します。
フェルミ準位の算出には、数値積分とニュートン法を組み合わせた手法を使用します。
計算された厳密なフェルミ準位と、ソマーフェルト展開による近似式の結果を比較し、標準出力に表形式で表示します。
結果は matplotlib を用いてグラフとして可視化され、温度に対するフェルミ準位の変化がプロットされます。
電子密度と有効質量はコマンドライン引数として指定でき、柔軟なシミュレーションが可能です。

このプログラムは、金属の電子物性において重要なフェルミ準位が温度によってどのように変化するかを理解し、その厳密解と近似解の乖離を視覚的に確認することを目的としています。

原理

本プログラムは、金属中の電子がフェルミ・ディラック統計に従う3次元自由電子ガスとして振る舞うというモデルに基づいています。

状態密度 (Density of States, DOS) 3次元自由電子ガスにおける状態密度 $D(E)$ は、エネルギー $E$ に対して以下の式で与えられます。 $$D(E) = D_0 \sqrt{E}$$ ここで $D_0$ は以下の定数です。 $$D_0 = (2S+1) \frac{2\pi (2m^*)^{3/2}}{h^3}$S$ は電子のスピン ($1/2$)、 $m^$ は電子の有効質量、 $h$ はプランク定数です。プログラム中では、`D0` は $m^$ を電子質量 $m_e$ に対する相対値として受け取り、最終的に単位を [cm$^{-3}$ eV$^{-1.5}$] に変換しています。
フェルミ・ディラック分布関数 電子が特定のエネルギー準位 $E$ を占有する確率を表すフェルミ・ディラック分布関数 $f(E, T, E_F)$ は、以下の式で与えられます。 $$f(E, T, E_F) = \frac{1}{\exp\left(\frac{E - E_F}{k_B T}\right) + 1}$$ ここで $T$ は絶対温度、 $E_F$ はフェルミ準位、 $k_B$ はボルツマン定数です。プログラムでは、エネルギー $E$ と $E_F$ を [eV] 単位として扱い、指数部の $\frac{E - E_F}{k_B T}$ を計算するために、分子を素電荷 $e$ で乗じることで [J] 単位に変換しています（exp((E - EF) * e / kB / T)）。
電子数密度 全電子数密度 $N$ は、状態密度とフェルミ・ディラック分布関数の積を全エネルギー範囲で積分することによって得られます。 $$N = \int_0^\infty D(E) f(E, T, E_F) dE$$ プログラムは、この関係式から与えられた電子数密度 $N$ に対応するフェルミ準位 $E_F$ を探します。
フェルミ準位の計算
- 0 K におけるフェルミ準位 ($E_{F0}$): $T=0$ K ではフェルミ・ディラック分布関数はヘビサイド関数となるため、上記の積分は解析的に解くことができ、$E_{F0}$ は以下の式で計算されます。 $$E_{F0} = \left(\frac{1.5 N}{D_0}\right)^{2/3}$$
- $T > 0$ K におけるフェルミ準位: $T > 0$ K では、CalEF 関数がニュートン法（scipy.optimize.newton）を用いてフェルミ準位 $E_F$ を数値的に求めます。ニュートン法は、関数 $g(E_F) = N_e(T, E_F, dE) - N = 0$ の根を見つけます。 Ne(T, EF, dE) 関数では、計算効率を高めるため、フェルミ準位 $E_F$ の近傍 ($E_F - dE$ から $E_F + dE$) は scipy.integrate.quad による数値積分で、それより低いエネルギー領域 ($0$ から $E_F - dE$) は解析的に積分することで計算されています。数値積分範囲 dE は $n_{range} \cdot k_B T / e$ で決定されます。
ソマーフェルト展開による近似式 高温におけるフェルミ準位の近似値は、ソマーフェルド展開を用いて以下の式で与えられます。 $$E_F(T) \approx E_{F0} - \frac{\pi^2}{6} \left(\frac{k_B T}{e}\right)^2 \frac{1}{D(E_{F0})} \left.\frac{dD(E)}{dE}\right|_{E=E_{F0}}$$ プログラムでは、この近似式も計算し、厳密な数値計算結果と比較します。ここで $e$ は素電荷であり、$k_B T / e$ は熱エネルギーを [eV] 単位に変換したものです。

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

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

numpy: 数値計算を効率的に行うための基盤ライブラリ。
scipy: 科学技術計算ライブラリ。特にintegrate（数値積分）とoptimize（最適化、ニュートン法）モジュールを使用します。
matplotlib: グラフ描画ライブラリ。

これらのライブラリは pip を使ってインストールできます。コマンドプロンプトやターミナルで以下のコマンドを実行してください。

pip install numpy scipy matplotlib

必要な入力ファイル

このプログラムは、コマンドライン引数以外の外部入力ファイル（例: テキストファイル、CSVファイル）を必要としません。必要なパラメータはプログラム内にハードコードされているか、コマンドライン引数で指定されます。

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

プログラムの実行により、以下の出力が生成されます。

標準出力: 計算された各温度におけるフェルミ準位（厳密値と近似値）、および電子数密度のチェック値が表形式で表示されます。出力例：
```
 T(K)   	  EF(eV)  	EFapprox(eV)	Ncheck(cm-3)
0.0000	  3.120610	  3.120610	1.000000e+23
200.0000	  3.119290	  3.119290	1.000000e+23
...
```
グラフウィンドウ: matplotlib によって生成されたグラフウィンドウが表示されます。このグラフには、温度（横軸）に対するフェルミ準位（縦軸）の厳密解と近似解が異なる色でプロットされます。ウィンドウは、ユーザーがEnterキーを押すまで開いたままになります。

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

基本的な実行コマンドと引数の説明は以下の通りです。

python EF-T-metal.py [電子密度N] [有効質量m]

EF-T-metal.py: 実行するPythonスクリプトの名前。
[電子密度N]: オプション。電子数密度を [cm$^{-3}$] 単位で指定します。デフォルト値は 5.0e22 cm$^{-3}$ です。
[有効質量m]: オプション。有効質量を電子質量 ($m_e$) に対する相対値として指定します。デフォルト値は 1.0 ($1.0 m_e$) です。

引数を指定しない場合、プログラムはデフォルト値を使用します。引数を指定する場合、電子密度と有効質量は両方、または電子密度のみを指定できます（有効質量のみを指定することはできません）。

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

例1: デフォルト値での実行

電子密度 $N = 5.0 \times 10^{22}$ cm$^{-3}$、有効質量 $m = 1.0 m_e$ を使用してプログラムを実行します。

python EF-T-metal.py

実行結果の説明: 標準出力に、0Kから4000Kまでの温度範囲で200K刻みで計算された、フェルミ準位（厳密値と近似値）および電子数密度の検証値が表形式で出力されます。同時に、matplotlib によってグラフウィンドウがポップアップ表示され、厳密解と近似解のフェルミ準位の温度依存性が視覚的に比較されます。通常、低温では両者はよく一致しますが、高温になるにつれて乖離が見られる場合があります。プログラムはユーザーがEnterキーを押すまでグラフウィンドウを開いたまま待機します。

例2: 特定の電子密度と有効質量を指定して実行

電子密度を $6.0 \times 10^{22}$ cm$^{-3}$、有効質量を $1.2 m_e$ に設定してプログラムを実行します。

python EF-T-metal.py 6.0e22 1.2

実行結果の説明: 上記例1と同様に、標準出力に計算結果の表が出力され、グラフウィンドウが表示されます。ただし、使用される電子密度と有効質量が指定された値に変更されるため、フェルミ準位の絶対値や温度依存性の傾向はデフォルト値の場合とは異なります。例えば、電子密度が高いほどフェルミ準位は高くなり、有効質量が大きいほどフェルミ準位は低くなる傾向があります。グラフはこれらの新しいパラメータに基づいた結果を反映します。