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

プログラムの動作

bose_condensation.py は、ボーズ凝縮現象をシミュレートし、ボーズ粒子系の巨視的な特性を計算・可視化するためのPythonプログラムです。特に、化学ポテンシャル \(\mu\)、粒子濃度 \(N\)、およびボーズ・アインシュタイン積分 \(F_\sigma(\alpha)\) の温度依存性や引数依存性を解析します。

このプログラムの主な機能は以下の通りです。

物理定数の定義: 物理学で一般的に使用される主要な定数（プランク定数 \(h\)、ボルツマン定数 \(k_B\)、電子電荷 \(e\) など）がプログラム冒頭で定義されています。
ガンマ関数の計算: \(\Gamma(\sigma)\) を再帰的に計算する関数 Gamma を提供します。
ボーズ・アインシュタイン積分の計算: \(F_{\sigma}(\alpha)\) を数値積分により計算する関数 Fsalpha を提供します。この関数は、ボーズ粒子系の粒子数密度を決定する上で中心的な役割を果たします。
粒子濃度計算: 化学ポテンシャルと温度が与えられたときの粒子濃度を計算する関数 Ne を提供します。
化学ポテンシャル決定: 与えられた粒子濃度と温度に対して、化学ポテンシャル（プログラム中では EF として扱われ、eV 単位で表現されます）を二分法で求める関数 CalEF を提供します。
ボーズ凝縮臨界温度の計算: 指定された粒子数密度に基づき、ボーズ凝縮が起こる臨界温度 \(T_c\) を計算します。
モードに応じた計算とプロット:
- Fs モードでは、様々な温度における粒子濃度 \(N\) を化学ポテンシャルに関連するパラメータ \(\alpha\) の関数として計算し、プロットします。
- mu モードでは、温度範囲にわたって化学ポテンシャル、熱的波長、臨界温度、凝縮粒子数などを計算し、化学ポテンシャルに関連するパラメータ \(\alpha\) を温度の関数としてプロットします。
コマンドライン引数による柔軟な設定: 実行時に計算モードや温度範囲、ステップサイズなどを指定できます。

このプログラムは、ボーズ凝縮の基本的な物理現象を理解し、その挙動を定量的に解析・可視化したい研究者や学生にとって有用です。

原理

本プログラムは、ボーズ・アインシュタイン統計に従う自由粒子系の理論に基づいています。

粒子数密度 \(N\): 3次元の自由ボーズガスにおいて、粒子数密度 \(N\) は温度 \(T\) と化学ポテンシャル \(\mu\) の関数として、以下のボーズ・アインシュタイン積分 \(F_\sigma(\alpha)\) を用いて表されます。

\[N = \frac{1}{\lambda_T^3} F_{3/2}(\alpha)\]

ここで、\(\lambda_T\) は熱的波長であり、次式で定義されます。

\[\lambda_T = \frac{h}{\sqrt{2\pi m k_B T}}\]

また、\(\alpha\) は無次元のパラメータで、\(\alpha = -\mu / (k_B T)\) と定義されます。bose_condensation.py では EF という変数が mu / e （化学ポテンシャルを電子ボルト単位で表したもの）に対応し、alpha = -EF / (kB * T / e) の関係が用いられています。
ボーズ・アインシュタイン積分 \(F_\sigma(\alpha)\): この関数は、一般的に以下のように定義されます。

\[F_{\sigma}(\alpha) = \frac{1}{\Gamma(\sigma)} \int_0^\infty \frac{y^{\sigma - 1}}{e^{y + \alpha} - 1} dy\]

プログラムでは、Fsalpha(sigma, alpha) 関数がこれに対応し、scipy.integrate.quad を用いて数値的に計算されます。ボーズ粒子系では \(\sigma = 3/2\) がよく使われます。
ガンマ関数 \(\Gamma(\sigma)\): 階乗関数の一般化であるガンマ関数は、Gamma(sigma) 関数によって再帰的に計算されます。特定の引数 (\(1.0\), \(0.5\)) に対しては特殊な値が返されます。
臨界温度 \(T_c\): ボーズ凝縮は、温度が臨界温度 \(T_c\) 以下になると、巨視的な数の粒子が最低エネルギー状態（基底状態）に凝縮する現象です。\(T_c\) は、化学ポテンシャル \(\mu=0\) （すなわち \(\alpha=0\)）のときに、すべての粒子が基底状態以外の励起状態に存在できる最大粒子数密度から決定されます。このとき、\(F_{3/2}(0) = \zeta(3/2)\)（\(\zeta\) はリーマンのゼータ関数）となることを利用し、\(T_c\) は次のように計算されます。

\[T_c = \frac{h^2}{2\pi m k_B} \left( \frac{N}{\zeta(3/2)} \right)^{2/3}\]

プログラムでは、zeta32 = Fsalpha(1.5, 0.0) として \(\zeta(3/2)\) の値が数値的に計算されています。
化学ポテンシャルの決定: CalEF 関数は、与えられた粒子数密度 \(N\) と温度 \(T\) に対して、適切な化学ポテンシャル（EF）を見つけるために二分法を使用します。これは、\(Ne(EF, T) - N = 0\) となる EF を探索するものです。
ボーズ凝縮下の粒子数: \(T < T_c\) の場合、化学ポテンシャルは \(\mu = 0\) となり、基底状態以外の励起状態にある粒子数 \(N'\) は \(N' = N (T/T_c)^{3/2}\) となります。したがって、基底状態に凝縮した粒子数 \(N_0\) は \(N_0 = N - N'\) と計算されます。
近似式: ExecMu 関数内では、\(T > T_c\) の場合に \(\alpha\) の近似式が用いられています。これは、\(T \approx T_c\) 付近での化学ポテンシャルの振る舞いを近似的に記述するものです。

これらの原理に基づいて、プログラムはボーズ粒子系の状態を計算し、温度や化学ポテンシャルの変化に対する応答を解析します。

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

bose_condensation.py は以下の非標準Pythonライブラリを使用しています。

numpy: 数値計算を効率的に行うための基本的なライブラリです。
scipy: 科学技術計算のためのライブラリで、このプログラムでは特に数値積分 (scipy.integrate) の機能を使用しています。optimize モジュールもインポートされていますが、newton 関数などの具体的な機能は使われておらず、二分法はプログラム内で独自に実装されています。
matplotlib: グラフ描画のためのライブラリです。

これらのライブラリは、Pythonのパッケージマネージャーである pip を使ってインストールできます。

pip install numpy scipy matplotlib

必要な入力ファイル

bose_condensation.py は、実行に必要な外部入力ファイルを一切使用しません。すべての設定パラメータ（粒子の質量 m、粒子濃度 N、温度範囲、\(\alpha\) 範囲など）は、プログラムのソースコード内に直接定義されているか、コマンドライン引数を通じて実行時に与えられます。

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

bose_condensation.py は、計算結果をファイルとして保存しません。プログラムの実行結果は、以下の形式で出力されます。

標準出力 (コンソール): 各計算ステップでの温度、熱的波長、化学ポテンシャル、粒子濃度、凝縮粒子数などの数値データが整形された表形式で表示されます。
インタラクティブなグラフウィンドウ: matplotlib を使用して、計算されたデータに基づいたグラフ（例えば、\(\alpha\) 対 \(N_e\) のプロットや、温度対 \(\alpha\) のプロット）がリアルタイムで表示されます。このウィンドウはユーザーが Enter キーを押すまで開いたままになります。

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

bose_condensation.py は、実行モードと計算範囲をコマンドライン引数で指定します。

基本的な書式は以下の通りです。

mu モード（化学ポテンシャルと温度の関係を計算・プロット）:
```
python bose_condensation.py mu Tmin Tmax Tstep
```
- Tmin: 最小温度 (K)
- Tmax: 最大温度 (K)
- Tstep: 温度のステップサイズ (K)
Fs モード（\(F_\sigma(\alpha)\) と粒子濃度の関係を計算・プロット）:
```
python bose_condensation.py Fs Tmin Tmax Tstep alphamin alphamax alphastep
```
- Tmin: 最小温度 (K)
- Tmax: 最大温度 (K)
- Tstep: 温度のステップサイズ (K)
- alphamin: \(\alpha\) の最小値
- alphamax: \(\alpha\) の最大値
- alphastep: \(\alpha\) のステップサイズ

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

例1: 化学ポテンシャルの温度依存性を計算・プロットする (`mu` モード)

この例では、温度が 2.5 K から 4.5 K まで 0.01 K 刻みで変化するときの化学ポテンシャルに関連するパラメータ \(\alpha\) の挙動を計算し、プロットします。

python bose_condensation.py mu 2.5 4.5 0.01

実行結果の説明:

プログラムはまず、粒子の質量 m (kg)、粒子濃度 N (cm^-3)、およびガンマ関数と \(\zeta(3/2)\) の値を出力します。次に、計算された臨界温度 Tc を表示します。

その後、指定された温度範囲で、各温度における熱的波長 (\(\lambda_T\))、化学ポテンシャル（EF、eV単位）、近似化学ポテンシャル（EFapprox）、励起状態の粒子数密度 (N')、最大粒子数密度 (Nmax)、そしてボーズ凝縮した粒子数密度 (N0) が標準出力に表形式で表示されます。

T < Tc の場合、化学ポテンシャルは 0 になり、基底状態以外の粒子数 N' は N * (T/Tc)^1.5 で計算され、凝縮粒子数 N0 が表示されます。
T >= Tc の場合、化学ポテンシャルは非ゼロとなり、CalEF 関数によって二分法で計算されます。

最後に、温度 \(T\) に対する \(\alpha = -\mu/(k_B T)\) のプロットがインタラクティブなウィンドウに表示されます。プロットには、厳密な計算値と近似値の曲線が含まれ、温度が \(T_c\) を下回ると \(\alpha\) がゼロになる様子が視覚的に確認できます。グラフが表示された後、ユーザーが Enter キーを押すまでプログラムは一時停止します。

例2: 粒子濃度 \(N_e\) の \(\alpha\) 依存性を特定の温度で計算・プロットする (`Fs` モード)

この例では、温度を 3 K から 4.5 K まで 0.2 K 刻みで変化させながら、\(\alpha\) の値を 0 から 0.5 まで 0.002 刻みで変化させたときの粒子濃度 \(N_e\) を計算し、プロットします。

python bose_condensation.py Fs 3 4.5 0.2 0 0.5 0.002

実行結果の説明:

プログラムは、例1と同様に定数と臨界温度を出力します。次に、指定された各温度 (T) ごとに、その温度での熱的波長 (lambda_T) と最大粒子数 (Nmax) を表示します。

その後、各温度に対して \(\alpha\) を変化させながら、計算された化学ポテンシャル (mu(eV))、ボーズ・アインシュタイン積分 Fs の値、および粒子濃度 Ne(cm-3) が計算されます（ただし、このモードでは標準出力への詳細な数値はコメントアウトされており、主にグラフに反映されます）。

最終的に、指定された各温度における粒子濃度 \(N_e\) の \(\alpha\) に対する依存性がプロットされたグラフがインタラクティブなウィンドウに表示されます。グラフには、システム全体の総粒子数 N を示す水平な破線も表示され、各温度での \(N_e\) が \(N\) に達する点を確認できます。このグラフが表示された後、ユーザーが Enter キーを押すまでプログラムは一時停止します。