このオッサンがすごい2025～o3-miniを利用した光学応答のPython実装

 【イントロ】

前回のアドベントカレンダー記事では、固体の光学応答を解析する際に利用されるクラマース・クローニッヒ変換(KK変換)をAIを用いて実装に取り組む記事を報告しました。

その際は、Anthoropic社の提供するチャット型AI、Claudeを利用することでコード生成を行いました。

生成においては、Claudeが出力したPythonコードをGoogle Colab上で実行し、エラーをClaudeに再質問することで、期待される振る舞いを出力するコードを作成しました。

つまり、コード実行にはある程度の修正対応が必要でした。

さて、2025年1月31日、OpenAI社が新型のAI、o3-miniをリリースしました。

「小規模で高速ながら、特に科学、数学、コーディングに強みを発揮」するとの触れ込みで、ますますAIの活躍の幅が広がる可能性を期待させてくれます。(Impress watch参考)

今回の記事では、このオッサンミニ、ではなくo3-miniを利用することで、よりよいPythonコードの作成体験が実現できるのではないかと考え、前回に引き続き、KK変換の実装に取り組みました。

結論、このオッサン・・・すごいです。

元気なおじさんの画像(いらすとや)

【方法】

o3-miniはChatGPTの無料ユーザーにも提供されていたため、その機能を利用しました。

以下のプロンプトでコード生成を指示しました

”””

あなたは物理学の専門家で、Pythonのプロフェッショナルです。波数と反射率のデータをクラマース・クローニッヒ変換し、光学伝導度を出力するPythonコードを書いてください。ただし、反射率はDrudeの式に基づく金属の場合を例としてください。

”””

ChatGPTの指示と応答

【結果】

プロンプトを送信すると、o3-miniにより推論が始まり、その次にKK変換の実行手順の解説とPythonコードの実装が始まります。そして出力されたPythonコードをGoogle Colabで実行すると・・・なんと！一発で動作しました。

なんということでしょう。これは驚き屋が湧き上がるのも納得です。

ドルーデモデルの誘電率から直接求めた光学伝導度と、KK変換を挟んで出力した光学伝導度の比較

ここで得られた結果を利用して、KK変換の振る舞い、そしてモデルのパラメータ依存性を可視化することを行いました。

【基本の説明】

一般に固体、特に金属の光学応答は、自由電子による応答を反映するドルーデモデルと、フォノンやバンド間遷移による応答を反映するローレンツモデルによりモデル化することができます。

反射率と光学伝導度に対して、それぞれのモデルは以下のような振る舞いを示します。光学伝導度に関しては特に断りがなければ実部の振る舞いを議論します。また、波数と周波数を区別せず説明している部分がありますが、相互に変換可能な単位のため、必要に応じて読み替えるものとします。揃えるの大変なので。

参考文献：田島節子 著、高温超伝導体の電荷応答: 強い電子相互作用がもたらすエキゾチックな物性 (物質・材料テキストシリーズ)

• 反射率
• ドルーデ成分：ゼロ波数に向かって反射率が１に向かって上昇する振る舞い。特にプラズマ周波数と呼ばれるエネルギー以下で反射率が急上昇する。これは、入射電磁波に応じる自由電子により電磁波が放出されるが、これ以上の周波数では自由電子が電磁波の振動に追随できず、電磁波を反射できないと解釈できる。
• ローレンツ成分：有限周波数にピークを持つ構造。バンド間遷移やフォノン振動に対応するエネルギーで電磁波が反射されると解釈できる
• 光学伝導度
• ドルーデ成分：ゼロ波数に向かってピークをもつ構造
• ローレンツ成分：有限周波数にピークをもつ構造

【KK変換の振る舞い】

よくある実験としては、フーリエ変換型赤外分光法(FT-IR)で得られた光学反射率の波数依存性をKK変換して光学伝導度を求めることが行われます。

KK変換の難しい点は、その計算中に発散のある主値積分が含まれる点と、積分範囲として0から無限大までの波数成分が必要な点です。

実験で得られる光学反射率は有限成分までに限られるため、可能な限り高エネルギーまで測定を行い、更にそれ以上の波数成分は外挿して補間する処理が必要となります。

ここで気になるのは、以下の2点です。

①高波数成分が少ないとどうなるのか

②波数の分解能の影響はどの程度あるか

まず、①の点について比較を行いました。

測定点数を5万に固定し、積分範囲を１から10万cm-1までにした場合と、1から50万cm-1にした場合で、KK変換の結果を比較しました。結果として、積分範囲を広く取らないと、KK変換後の光学伝導度の振る舞いがマイナスの値を示すという非物理的な結果になることが確認できました。KK変換のために高周波数側までのデータ必要となる所以がここにあります。

次に②の点について比較を行いました。

積分範囲を１から10万cm-1に固定し、測定点数を1万から2万に増やした場合で、KK変換の結果を比較しました。結果として、測定点数を増やしても非物理的な振る舞いを示すKK変換後の結果は変わらないことが確認できました。

上記比較から、KK変換の妥当性を確保するためには、測定分解能をむやみに上げるよりも、高周波数側（と低周波数側）のデータの取得が重要であると考えられます。

【モデルのパラメータ依存性】

このパートでは、モデルのパラメータ依存性、特にドルーデモデルの散乱率とプラズマ周波数依存性について可視化を行いました。

まず、散乱率の依存性です。散乱率を100から1000cm-1まで変化させた場合の振る舞いです。

散乱率は、文字通り電子の散乱される確率を意味しています。一般に金属の抵抗率は低温ほど小さくなります。これは低温では電子の格子による散乱が抑制され、さらに電子間の散乱の抑制されるため、低温ほど散乱率が小さくなることが期待されます。

以下図では、散乱率が小さくなるほど、反射率が高くなっており、低温ほど反射率が高くなることが確認できます。この反射率をKK変換したあとの光学伝導度も、散乱率が小さいゼロ波数のピークの幅が狭くなりゼロ波数の伝導度が高くなっており、抵抗率の振る舞いと一致しています（ゼロ波数の光学伝導度、つまり直流伝導度の逆数が抵抗率のため）。

次にプラズマ周波数の依存性です。プラズマ周波数(の自乗)は自由電子のキャリア密度に比例し、有効質量に反比例します。軽い電子ほど電磁波の振動に追随しやすいのでプラズマ周波数は高くなりそうなので、直感にも一致しています。尻軽ですね、電子ちゃん。

実際プラズマ周波数を1万から2万cm-1まで変化させたとき、反射率の立ち上がりが高波数側に移動していることがわかります。また同時に、光学伝導度のピークの高さが大きくなっていることがみてとれます。キャリア密度が多く電子が軽いほど電流が流れやすい、つまり伝導度が増加するという直感とも一致する振る舞いになっています。

金属はこの反射率の１への立ち上がり周波数の存在と、光学伝導度のゼロ周波数ピークの存在が特徴となっています。

【まとめ】

オッサンミニ、じゃなくてo3-miniすごい！

簡単な質問を投げかけるだけで、実行可能なPythonコードを出力してくれました。この他にも、電気抵抗の久保公式やMuSRのKubo-Toyabe関数、比熱のBCS的温度依存性も簡単な質問だけで実装できることが確認できました。AIの数学、物理学の知識とコード実装能力の進化を感じます。この能力を利用すれば、論文に記載されている式の手元での再現や、パラメータ変化させたときの比較、そして実験結果とのフィッティングなども容易になることが期待できます。

もちろん、出力される結果が妥当な内容か検証できる能力は重要ですが、最初のたたき台をすごいスピードで作ってくれるのは助かります。

研究において基本的にコーディングは本来悩むべき目的ではなく、手段に過ぎません。AIと上手に付き合って、研究のスピードを加速させていきたいですね。

走り続けるおっさん（いらすとや）

【実行コード】

以下は、実際に実行した際のGoogle ColabのPythonコードです。

In [6]:

!pip install nbconvert

Requirement already satisfied: nbconvert in /usr/local/lib/python3.11/dist-packages (7.16.6)
Requirement already satisfied: beautifulsoup4 in /usr/local/lib/python3.11/dist-packages (from nbconvert) (4.12.3)
Requirement already satisfied: bleach!=5.0.0 in /usr/local/lib/python3.11/dist-packages (from bleach[css]!=5.0.0->nbconvert) (6.2.0)
Requirement already satisfied: defusedxml in /usr/local/lib/python3.11/dist-packages (from nbconvert) (0.7.1)
Requirement already satisfied: jinja2>=3.0 in /usr/local/lib/python3.11/dist-packages (from nbconvert) (3.1.5)
Requirement already satisfied: jupyter-core>=4.7 in /usr/local/lib/python3.11/dist-packages (from nbconvert) (5.7.2)
Requirement already satisfied: jupyterlab-pygments in /usr/local/lib/python3.11/dist-packages (from nbconvert) (0.3.0)
Requirement already satisfied: markupsafe>=2.0 in /usr/local/lib/python3.11/dist-packages (from nbconvert) (3.0.2)
Requirement already satisfied: mistune<4,>=2.0.3 in /usr/local/lib/python3.11/dist-packages (from nbconvert) (3.1.1)
Requirement already satisfied: nbclient>=0.5.0 in /usr/local/lib/python3.11/dist-packages (from nbconvert) (0.10.2)
Requirement already satisfied: nbformat>=5.7 in /usr/local/lib/python3.11/dist-packages (from nbconvert) (5.10.4)
Requirement already satisfied: packaging in /usr/local/lib/python3.11/dist-packages (from nbconvert) (24.2)
Requirement already satisfied: pandocfilters>=1.4.1 in /usr/local/lib/python3.11/dist-packages (from nbconvert) (1.5.1)
Requirement already satisfied: pygments>=2.4.1 in /usr/local/lib/python3.11/dist-packages (from nbconvert) (2.18.0)
Requirement already satisfied: traitlets>=5.1 in /usr/local/lib/python3.11/dist-packages (from nbconvert) (5.7.1)
Requirement already satisfied: webencodings in /usr/local/lib/python3.11/dist-packages (from bleach!=5.0.0->bleach[css]!=5.0.0->nbconvert) (0.5.1)
Requirement already satisfied: tinycss2<1.5,>=1.1.0 in /usr/local/lib/python3.11/dist-packages (from bleach[css]!=5.0.0->nbconvert) (1.4.0)
Requirement already satisfied: platformdirs>=2.5 in /usr/local/lib/python3.11/dist-packages (from jupyter-core>=4.7->nbconvert) (4.3.6)
Requirement already satisfied: jupyter-client>=6.1.12 in /usr/local/lib/python3.11/dist-packages (from nbclient>=0.5.0->nbconvert) (6.1.12)
Requirement already satisfied: fastjsonschema>=2.15 in /usr/local/lib/python3.11/dist-packages (from nbformat>=5.7->nbconvert) (2.21.1)
Requirement already satisfied: jsonschema>=2.6 in /usr/local/lib/python3.11/dist-packages (from nbformat>=5.7->nbconvert) (4.23.0)
Requirement already satisfied: soupsieve>1.2 in /usr/local/lib/python3.11/dist-packages (from beautifulsoup4->nbconvert) (2.6)
Requirement already satisfied: attrs>=22.2.0 in /usr/local/lib/python3.11/dist-packages (from jsonschema>=2.6->nbformat>=5.7->nbconvert) (25.1.0)
Requirement already satisfied: jsonschema-specifications>=2023.03.6 in /usr/local/lib/python3.11/dist-packages (from jsonschema>=2.6->nbformat>=5.7->nbconvert) (2024.10.1)
Requirement already satisfied: referencing>=0.28.4 in /usr/local/lib/python3.11/dist-packages (from jsonschema>=2.6->nbformat>=5.7->nbconvert) (0.36.2)
Requirement already satisfied: rpds-py>=0.7.1 in /usr/local/lib/python3.11/dist-packages (from jsonschema>=2.6->nbformat>=5.7->nbconvert) (0.22.3)
Requirement already satisfied: pyzmq>=13 in /usr/local/lib/python3.11/dist-packages (from jupyter-client>=6.1.12->nbclient>=0.5.0->nbconvert) (24.0.1)
Requirement already satisfied: python-dateutil>=2.1 in /usr/local/lib/python3.11/dist-packages (from jupyter-client>=6.1.12->nbclient>=0.5.0->nbconvert) (2.8.2)
Requirement already satisfied: tornado>=4.1 in /usr/local/lib/python3.11/dist-packages (from jupyter-client>=6.1.12->nbclient>=0.5.0->nbconvert) (6.4.2)
Requirement already satisfied: six>=1.5 in /usr/local/lib/python3.11/dist-packages (from python-dateutil>=2.1->jupyter-client>=6.1.12->nbclient>=0.5.0->nbconvert) (1.17.0)
Requirement already satisfied: typing-extensions>=4.4.0 in /usr/local/lib/python3.11/dist-packages (from referencing>=0.28.4->jsonschema>=2.6->nbformat>=5.7->nbconvert) (4.12.2)

In [ ]:

from google.colab import drive
drive.mount('/content/drive')

In [ ]:

!jupyter nbconvert --to html /content/drive/MyDrive/Colab_Notebooks/20250201_KK変換_ChatGPT.ipynb

o3-miniの出力結果¶

以下は、Drude‐モデルに基づく金属の反射率データから、クラマース・クローニッヒ変換（KK変換）を用いて位相情報を復元し、そこから複素屈折率・誘電率を求め、最終的に光学伝導度（ここでは ( \sigma(\omega)=\frac{\omega,\mathrm{Im},\epsilon(\omega)}{4\pi} ) として計算）を得る例です。

以下のコードは、

1. 与えたパラメータ（プラズマ振動数 (\omega_p)、散逸率 (\gamma)、背景定数 (\epsilon_\infty)）からDrudeモデルの誘電率を計算し、そこから通常のFresnelの式（正入射の場合）で反射率 (R(\omega)) を算出します。
2. 数値的にKK変換を行い、反射率から反射振幅の位相 (\theta(\omega)) を復元します。
3. 反射振幅 (r(\omega)=\sqrt{R(\omega)}e^{i\theta(\omega)}) から逆Fresnel式 ( n(\omega)=(1+r)/(1-r) ) により複素屈折率（複素リフラクティブ指数）を求め、そこから複素誘電率 (\epsilon(\omega)=n(\omega)^2) を計算します。
4. 誘電率から、上記の式により光学伝導度 (\sigma(\omega)) を求め、グラフにプロットします。

なお、数値積分には特に分割点付近での特異点処理（主値積分）を意識しており、簡易的な例として実装しています。実際の解析ではデータの補完や積分区間の工夫が必要となります。

以下にサンプルコードを示します。

---

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import simps

# --- パラメータの設定 ---
# 単位系はここでは波数(cm^-1)を用いると仮定します
omega_min = 1.0        # 最小周波数 [cm^-1]
omega_max = 30000.0    # 最大周波数 [cm^-1]
N = 2000               # 周波数の分割数
omega = np.linspace(omega_min, omega_max, N)  # 周波数配列

# Drudeモデルのパラメータ
omega_p = 15000.0      # プラズマ振動数 [cm^-1]
gamma   = 500.0        # 散逸率 [cm^-1]
eps_inf = 1.0          # 高周波の背景誘電率

# --- Drudeモデルに基づく誘電率、屈折率、反射率の計算 ---
# Drudeの誘電率: eps(omega)= eps_inf - omega_p^2/(omega^2 + i*omega*gamma)
eps = eps_inf - omega_p**2/(omega**2 + 1j*omega*gamma)
# 複素屈折率 n = sqrt(eps) (注意：複素平方根)
n_complex = np.sqrt(eps)
# 正入射でのFresnelの反射係数: r = (n - 1)/(n + 1)
r_complex_true = (n_complex - 1)/(n_complex + 1)
# 反射率
R = np.abs(r_complex_true)**2

# --- クラマース・クローニッヒ変換により位相を求める ---
def kk_phase(omega, R):
    """
    与えられた周波数配列 omega と反射率 R から、
    クラマース・クローニッヒ変換により反射振幅の位相 theta(omega) を計算する。
    
    theta(omega) = - (omega/π) P ∫_0^∞ ln R(omega')/(omega'^2 - omega^2) d omega'
    ここでは数値積分により近似する。
    """
    theta = np.zeros_like(omega)
    d_omega = omega[1] - omega[0]
    lnR = np.log(R)
    
    for i, w in enumerate(omega):
        # 積分区間：全周波数（ただし、分母に特異点があるので近傰点を除く）
        integrand = lnR / (omega**2 - w**2)
        # 主値積分：omega[i] の近傾部分は除外（絶対値が非常に大きくなる点）
        # ここでは、該当インデックス i を飛ばす
        integrand[i] = 0.0  # この点はゼロとする
        theta[i] = - w/np.pi * simps(integrand, omega)
    return theta

theta = kk_phase(omega, R)

# --- KK変換から反射振幅 r を復元し、屈折率を求める ---
r_complex_kk = np.sqrt(R) * np.exp(1j * theta)

# 逆Fresnel式より、反射係数 r = (n - 1)/(n + 1) なので、
# 解いて n = (1 + r)/(1 - r)
n_complex_kk = (1 + r_complex_kk) / (1 - r_complex_kk)
eps_kk = n_complex_kk**2

# --- 光学伝導度の計算 ---
# ここでは、CGS単位系における一般的な関係
#   σ(ω) = [ω * Im(ε(ω))] / (4π)
sigma = omega * np.imag(eps_kk) / (4*np.pi)

# --- プロット ---
plt.figure(figsize=(10,6))
plt.plot(omega, sigma, label='Optical Conductivity (KK derived)')
plt.xlabel("Wave number ω (cm$^{-1}$)")
plt.ylabel("σ(ω) [arb. units]")
plt.title("Optical Conductivity from KK Analysis of Drude Reflectance")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

# --- 参考: Drudeの理論式による光学伝導度 ---
# Drudeの理論式 (CGS) では、実部伝導度は
#   σ(ω) = (ω_p^2/(4π)) * (γ/(ω^2 + γ^2))
sigma_drude = (omega_p**2/(4*np.pi)) * (gamma/(omega**2 + gamma**2))

plt.figure(figsize=(10,6))
plt.plot(omega, sigma, label='KK derived σ(ω)', lw=2)
plt.plot(omega, sigma_drude, '--', label='Drude model σ(ω)', lw=2)
plt.xlabel("Wave number ω (cm$^{-1}$)")
plt.ylabel("σ(ω) [arb. units]")
plt.title("Comparison of KK derived and Drude model Optical Conductivity")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

---

コードのポイント解説¶

1. Drude反射率の生成

   • (\epsilon(\omega)=\epsilon_\infty - \frac{\omega_p^2}{\omega^2 + i,\omega,\gamma}) として誘電率を計算し，
   • ( n(\omega)=\sqrt{\epsilon(\omega)}) を用い，
   • 正入射での反射係数 ( r(\omega)=\frac{n-1}{n+1} ) から反射率 ( R=|r|^2 ) を得ています。

2. クラマース・クローニッヒ変換

   • 関数 kk_phase 内で、数値積分（Simpson法）を用いて
     (\theta(\omega) = -\frac{\omega}{\pi} , P \int_0^\infty \frac{\ln R(\omega')}{\omega'^2-\omega^2},d\omega')
     を計算しています。
   • 積分の特異点（(\omega'=\omega)）は該当点の値をゼロとすることで主値積分を近似しています。

3. 反射振幅から屈折率への変換

   • KK変換で得た位相と元の反射率から ( r(\omega)=\sqrt{R},e^{i\theta(\omega)} ) を復元し，
   • 逆Fresnel式 ( n=(1+r)/(1-r) ) により複素屈折率を求めています。

4. 光学伝導度の計算

   • 誘電率から ( \sigma(\omega)=\frac{\omega,\mathrm{Im}(\epsilon(\omega))}{4\pi} ) を算出しています（単位系により定数は変わります）。

5. 比較プロット

   • 最後に、Drudeモデルの既知の理論式
     (\sigma(\omega)=\frac{\omega_p^2}{4\pi},\frac{\gamma}{\omega^2+\gamma^2})
     とKK変換から得た結果を比較しています。

このコードを実行することで、Drudeの反射率からクラマース・クローニッヒ変換を用いて位相を復元し、光学伝導度を数値的に求める一連の流れを体験できます。

※ 注意：実際の実験データに対しては、周波数範囲の拡大（低周波・高周波補完）やノイズ除去などの前処理が必要です。また、積分の精度向上のための手法（例えば積分区間の分割の工夫など）も検討してください。

In [ ]:

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import simps

# --- パラメータの設定 ---
# 単位系はここでは波数(cm^-1)を用いると仮定します
omega_min = 1.0        # 最小周波数 [cm^-1]
omega_max = 3000000.0    # 最大周波数 [cm^-1]
N = 200000               # 周波数の分割数
omega = np.linspace(omega_min, omega_max, N)  # 周波数配列

# Drudeモデルのパラメータ
omega_p = 15000.0      # プラズマ振動数 [cm^-1]
gamma   = 500.0        # 散逸率 [cm^-1]
eps_inf = 1.0          # 高周波の背景誘電率

# --- Drudeモデルに基づく誘電率、屈折率、反射率の計算 ---
# Drudeの誘電率: eps(omega)= eps_inf - omega_p^2/(omega^2 + i*omega*gamma)
eps = eps_inf - omega_p**2/(omega**2 + 1j*omega*gamma)
# 複素屈折率 n = sqrt(eps) (注意：複素平方根)
n_complex = np.sqrt(eps)
# 正入射でのFresnelの反射係数: r = (n - 1)/(n + 1)
r_complex_true = (n_complex - 1)/(n_complex + 1)
# 反射率
R = np.abs(r_complex_true)**2

# --- クラマース・クローニッヒ変換により位相を求める ---
def kk_phase(omega, R):
    """
    与えられた周波数配列 omega と反射率 R から、
    クラマース・クローニッヒ変換により反射振幅の位相 theta(omega) を計算する。

    theta(omega) = - (omega/π) P ∫_0^∞ ln R(omega')/(omega'^2 - omega^2) d omega'
    ここでは数値積分により近似する。
    """
    theta = np.zeros_like(omega)
    d_omega = omega[1] - omega[0]
    lnR = np.log(R)

    for i, w in enumerate(omega):
        # 積分区間：全周波数（ただし、分母に特異点があるので近傰点を除く）
        integrand = lnR / (omega**2 - w**2)
        # 主値積分：omega[i] の近傾部分は除外（絶対値が非常に大きくなる点）
        # ここでは、該当インデックス i を飛ばす
        integrand[i] = 0.0  # この点はゼロとする
        theta[i] = - w/np.pi * simps(integrand, omega)
    return theta

theta = kk_phase(omega, R)

# --- KK変換から反射振幅 r を復元し、屈折率を求める ---
r_complex_kk = np.sqrt(R) * np.exp(1j * theta)

# 逆Fresnel式より、反射係数 r = (n - 1)/(n + 1) なので、
# 解いて n = (1 + r)/(1 - r)
n_complex_kk = (1 + r_complex_kk) / (1 - r_complex_kk)
eps_kk = n_complex_kk**2

# --- 光学伝導度の計算 ---
# ここでは、CGS単位系における一般的な関係
#   σ(ω) = [ω * Im(ε(ω))] / (4π)
sigma = omega * np.imag(eps_kk) / (4*np.pi)

# --- プロット ---
plt.figure(figsize=(10,6))
plt.plot(omega, sigma, label='Optical Conductivity (KK derived)')
plt.xlabel("Wave number ω (cm$^{-1}$)")
plt.ylabel("σ(ω) [arb. units]")
plt.title("Optical Conductivity from KK Analysis of Drude Reflectance")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

# --- 参考: Drudeの理論式による光学伝導度 ---
# Drudeの理論式 (CGS) では、実部伝導度は
#   σ(ω) = (ω_p^2/(4π)) * (γ/(ω^2 + γ^2))
sigma_drude = (omega_p**2/(4*np.pi)) * (gamma/(omega**2 + gamma**2))

plt.figure(figsize=(10,6))
plt.plot(omega, sigma, label='KK derived σ(ω)', lw=2)
plt.plot(omega, sigma_drude, '--', label='Drude model σ(ω)', lw=2)
plt.xlabel("Wave number ω (cm$^{-1}$)")
plt.ylabel("σ(ω) [arb. units]")
plt.title("Comparison of KK derived and Drude model Optical Conductivity")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

<ipython-input-23-3fcfb173ccbd>:42: RuntimeWarning: divide by zero encountered in divide
  integrand = lnR / (omega**2 - w**2)
<ipython-input-23-3fcfb173ccbd>:46: DeprecationWarning: 'scipy.integrate.simps' is deprecated in favour of 'scipy.integrate.simpson' and will be removed in SciPy 1.14.0
  theta[i] = - w/np.pi * simps(integrand, omega)

In [ ]:

# --- プロット ---
plt.figure(figsize=(10,6))
plt.plot(omega[0:2000], sigma[0:2000], label='Optical Conductivity (KK derived)')
plt.xlabel("Wave number ω (cm$^{-1}$)")
plt.ylabel("σ(ω) [arb. units]")
plt.title("Optical Conductivity from KK Analysis of Drude Reflectance")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

# --- 参考: Drudeの理論式による光学伝導度 ---
# Drudeの理論式 (CGS) では、実部伝導度は
#   σ(ω) = (ω_p^2/(4π)) * (γ/(ω^2 + γ^2))
sigma_drude = (omega_p**2/(4*np.pi)) * (gamma/(omega**2 + gamma**2))

plt.figure(figsize=(10,6))
plt.plot(omega[0:2000], sigma[0:2000], label='KK derived σ(ω)', lw=2)
plt.plot(omega[0:2000], sigma_drude[0:2000], '--', label='Drude model σ(ω)', lw=2)
plt.xlabel("Wave number ω (cm$^{-1}$)")
plt.ylabel("σ(ω) [arb. units]")
plt.title("Comparison of KK derived and Drude model Optical Conductivity")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

高速化¶

以下の方法などで、計算時間の短縮が期待できます。コード全体の中でボトルネックとなっているのは、各周波数 ( \omega_i ) に対して主値積分を行うループ部分です。以下、いくつかの高速化手法をご紹介します。

---

1. 並列化による高速化（joblib を用いる）¶

各 ( \omega_i ) に対する積分は互いに独立であるため、Python の並列処理ライブラリ joblib を用いて並列化することが可能です。以下は joblib を用いた例です。

from joblib import Parallel, delayed

def compute_theta(i, omega, lnR, d_omega):
    w = omega[i]
    integrand = lnR / (omega**2 - w**2)
    # 特異点（i 番目）は主値積分のため 0 に置換
    integrand[i] = 0.0
    # Simpson ルールで積分
    return - w/np.pi * simps(integrand, omega)

# 並列処理を用いて各 theta[i] を計算
theta_list = Parallel(n_jobs=-1)(
    delayed(compute_theta)(i, omega, np.log(R), omega[1]-omega[0]) for i in range(len(omega))
)
theta = np.array(theta_list)

上記コードでは n_jobs=-1 としているため、利用可能な全ての CPU コアを使って計算します。

---

2. ベクトル化による高速化¶

ループを完全に排除するため、2次元配列を用いて積分対象の行列全体を一度に計算する方法もあります。各 ( \omega_i ) に対する積分は、以下のように行列を作成することでベクトル化できます。

# lnR の準備
lnR = np.log(R)

# omega^2 のベクトル
omega2 = omega**2

# 行と列でブロードキャストして分母行列を作成
# 行: 固定した積分変数 w_i, 列: 積分変数 ω'
denom = omega2[None, :] - omega2[:, None]

# ゼロ除算となる対角成分を一旦適当な値にしておく（後で置換）
# ※ここでは np.eye を用いて後で対角成分だけ 0 に置換する方法をとる
integrand = lnR[None, :] / denom

# 対角成分（特異点）の処理：主値積分なので、対角成分は 0 にする
np.fill_diagonal(integrand, 0.0)

# 各行について積分を行う
theta = - omega/np.pi * simps(integrand, omega, axis=1)

この方法では (N \times N) の行列を作成します。N=2000 であれば約400万要素となり、現代のPCでは十分扱えますが、より大きなデータの場合はメモリ使用量に注意してください。

---

3. Numba を用いた JIT コンパイル¶

Python のループ部分に Numba を適用して JIT コンパイルすることで、ループ内の数値計算を高速化する方法もあります。たとえば、以下のように関数をデコレータ @njit を付与して高速化できます。

from numba import njit

@njit
def kk_phase_numba(omega, lnR):
    N = omega.shape[0]
    theta = np.empty(N, dtype=np.float64)
    d_omega = omega[1] - omega[0]
    # シンプソン則はループ内で単純な重み付き和として実装可能
    for i in range(N):
        w = omega[i]
        s = 0.0
        for j in range(N):
            if j == i:
                continue  # 特異点はスキップ
            # 積分の中間値
            s += lnR[j] / (omega[j]**2 - w**2)
        theta[i] = - w/np.pi * s * d_omega
    return theta

# 使用例
theta = kk_phase_numba(omega, np.log(R))

※ 上記はあくまで単純化した例です。シンプソン則は厳密には各区間ごとの重みを考慮する必要があるため、積分精度を確保するためには重みの計算を追加するなどの改良が必要です。

---

4. 高速フーリエ変換 (FFT) を用いたヒルベルト変換¶

クラマース・クローニッヒ変換はヒルベルト変換の一種であるため、理論上は FFT を用いて高速に計算できる方法もあります。ただし、FFT を適用するためにはデータの補間や端点処理（周期性の仮定など）が必要になるため、厳密性や計算範囲に注意が必要です。実際の応用例としては、以下のような手法が検討できます。

1. ( \ln R(\omega) ) のデータを適切に拡張し、対称性を持たせる。
2. FFT を用いてヒルベルト変換を計算する（SciPy の hilbert 関数も利用可能）。
3. 得られた位相情報から元の反射率振幅を再構成する。

ただし、この方法は実装の複雑性が増すため、精度とのトレードオフを考慮して検討してください。

---

まとめ¶

• 並列化: 各周波数に対する積分計算は互いに独立なので、joblib.Parallel などを用いると大幅な高速化が期待できます。
• ベクトル化: (N \times N) の行列を用いたベクトル化によって、Python のループを回避し、NumPy の高速な内部ルーチンを利用できます。
• Numba の利用: JIT コンパイルにより、純粋な Python ループを高速化可能です。
• FFT を用いたヒルベルト変換: 理論的には最も高速な手法ですが、前処理や端点処理の検討が必要です。

用途や求める精度、データサイズに応じて上記の手法を組み合わせることで、現在の処理の高速化が期待できます。

In [ ]:

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import simps
from joblib import Parallel, delayed

# --- パラメータの設定 ---
omega_min = 1.0        # 最小周波数 [cm^-1]
omega_max = 100000    # 最大周波数 [cm^-1]
N = 10000               # 周波数の分割数
omega = np.linspace(omega_min, omega_max, N)  # 周波数配列

# Drudeモデルのパラメータ
omega_p = 15000.0      # プラズマ振動数 [cm^-1]
gamma   = 500.0        # 散逸率 [cm^-1]
eps_inf = 1.0          # 高周波の背景誘電率

# --- Drudeモデルに基づく誘電率、屈折率、反射率の計算 ---
# Drudeの誘電率: eps(omega)= eps_inf - omega_p^2/(omega^2 + i*omega*gamma)
eps = eps_inf - omega_p**2 / (omega**2 + 1j * omega * gamma)
# 複素屈折率 n = sqrt(eps) (注意：複素平方根)
n_complex = np.sqrt(eps)
# 正入射でのFresnelの反射係数: r = (n - 1)/(n + 1)
r_complex_true = (n_complex - 1) / (n_complex + 1)
# 反射率
R = np.abs(r_complex_true)**2

# --- クラマース・クローニッヒ変換による位相の計算 (並列化版) ---
def compute_theta(i, omega, lnR):
    """
    指定されたインデックス i に対して、積分を実行して
    theta(omega_i) を計算する。
    """
    w = omega[i]
    # 各周波数での積分対象となる関数 lnR/(omega^2 - w^2)
    integrand = lnR / (omega**2 - w**2)
    # 分母がゼロとなる特異点 (j == i) は主値積分のため 0 にする
    integrand[i] = 0.0
    # Simpsonのルールで数値積分
    return - w/np.pi * simps(integrand, omega)

# lnR の事前計算
lnR = np.log(R)
# 並列化処理: 利用可能な全てのCPUコアを使用（n_jobs=-1）
theta_list = Parallel(n_jobs=-1)(
    delayed(compute_theta)(i, omega, lnR) for i in range(len(omega))
)
theta = np.array(theta_list)

# --- KK変換から反射振幅 r を復元し、屈折率および誘電率を計算 ---
r_complex_kk = np.sqrt(R) * np.exp(1j * theta)
n_complex_kk = (1 + r_complex_kk) / (1 - r_complex_kk)
eps_kk = n_complex_kk**2

# --- 光学伝導度の計算 ---
# CGS単位系における関係: σ(ω) = [ω * Im(ε(ω))] / (4π)
sigma = omega * np.imag(eps_kk) / (4 * np.pi)

# --- プロット: KK変換による光学伝導度 ---
plt.figure(figsize=(10, 6))
plt.plot(omega, sigma, label='Optical Conductivity (KK derived)')
plt.xlabel("Wave number ω (cm$^{-1}$)")
plt.ylabel("σ(ω) [arb. units]")
plt.title("Optical Conductivity from KK Analysis of Drude Reflectance")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

# --- 参考: Drudeの理論式による光学伝導度 ---
# Drudeモデル (CGS) による実部伝導度:
#   σ(ω) = (ω_p^2/(4π)) * (γ/(ω^2 + γ^2))
sigma_drude = (omega_p**2 / (4 * np.pi)) * (gamma / (omega**2 + gamma**2))

plt.figure(figsize=(10, 6))
plt.plot(omega, sigma, label='KK derived σ(ω)', lw=2)
plt.plot(omega, sigma_drude, '--', label='Drude model σ(ω)', lw=2)
plt.xlabel("Wave number ω (cm$^{-1}$)")
plt.ylabel("σ(ω) [arb. units]")
plt.title("Comparison of KK derived and Drude model Optical Conductivity")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

ローレンツモデルの追加¶

In [ ]:

以下は、Drude項に加えて Lorentz 項（単一オシレーター）を導入し、全体の誘電率から反射率を計算、KK変換で位相を求めた上で光学伝導度を計算する完全なコード例です。
なお、Lorentz項は
[ \varepsilon_{\mathrm{Lorentz}}(\omega)= \frac{S}{\omega_{0}^{2}-\omega^{2}-i\gamma_{L}\omega} ] という形で実装しており、ここではパラメータとして以下の例を採用しています。パラメータはお好みや実験データに合わせて調整してください。

• Drude項
  • プラズマ振動数: (\omega_p = 15000,\mathrm{cm^{-1}})
  • 散逸率: (\gamma = 500,\mathrm{cm^{-1}})
• Lorentz項
  • 共鳴周波数: (\omega_{0} = 20000,\mathrm{cm^{-1}})
  • オシレーター強度: (S = f,\omega_p^2) （ここでは (f=0.3) としています）
  • Lorentz の減衰定数: (\gamma_{L} = 1000,\mathrm{cm^{-1}})

以下に、joblib を用いた並列化処理を含む完全なコードを示します。

---

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import simps
from joblib import Parallel, delayed

# --- パラメータの設定 ---
omega_min = 1.0         # 最小周波数 [cm^-1]
omega_max = 30000.0     # 最大周波数 [cm^-1]
N = 2000                # 周波数の分割数
omega = np.linspace(omega_min, omega_max, N)  # 周波数配列

# --- Drudeモデルのパラメータ ---
omega_p = 15000.0       # Drudeのプラズマ振動数 [cm^-1]
gamma   = 500.0         # Drudeの散逸率 [cm^-1]
eps_inf = 1.0           # 高周波の背景誘電率

# --- Lorentzモデルのパラメータ ---
omega0_L = 20000.0      # Lorentz共鳴周波数 [cm^-1]
f = 0.3                 # オシレーター強度の比率
S = f * omega_p**2      # オシレーター強度
gamma_L = 1000.0        # Lorentz項の減衰定数 [cm^-1]

# --- 誘電率の計算 ---
# Drude項 (ω₀=0 の項)
eps_drude = eps_inf - omega_p**2 / (omega**2 + 1j * omega * gamma)

# Lorentz項
eps_lorentz = S / (omega0_L**2 - omega**2 - 1j * omega * gamma_L)

# 全体の誘電率 (Drude + Lorentz)
eps_total = eps_drude + eps_lorentz

# --- 複素屈折率と反射率の計算 ---
# 複素屈折率 n = sqrt(ε)
n_complex = np.sqrt(eps_total)
# 正入射での Fresnel の反射係数: r = (n-1)/(n+1)
r_complex_true = (n_complex - 1) / (n_complex + 1)
# 反射率 R = |r|^2
R = np.abs(r_complex_true)**2

# --- クラマース・クローニッヒ変換による位相の計算（並列化） ---
def compute_theta(i, omega, lnR):
    """
    インデックス i に対応する ω_i に対して、主値積分により
    位相 θ(ω_i) を計算する。
    
    θ(ω_i) = - (ω_i/π) P ∫_0∞  ln R(ω')/(ω'^2 - ω_i^2) dω'
    """
    w = omega[i]
    # 積分対象：全周波数で lnR/(ω'^2 - ω_i^2)
    integrand = lnR / (omega**2 - w**2)
    # 特異点 (j == i) は主値積分として 0 に置換
    integrand[i] = 0.0
    # Simpsonの公式による数値積分
    return - w/np.pi * simps(integrand, omega)

# lnR の事前計算
lnR = np.log(R)

# 並列処理: 各 ω_i での積分を全CPUコアで計算 (n_jobs=-1)
theta_list = Parallel(n_jobs=-1)(
    delayed(compute_theta)(i, omega, lnR) for i in range(len(omega))
)
theta = np.array(theta_list)

# --- KK変換から反射振幅の復元 ---
# 反射振幅: r(ω) = √R * exp(iθ)
r_complex_kk = np.sqrt(R) * np.exp(1j * theta)
# 逆Fresnel式より、r = (n-1)/(n+1) なので、解いて n = (1+r)/(1-r)
n_complex_kk = (1 + r_complex_kk) / (1 - r_complex_kk)
# 誘電率: ε = n^2
eps_kk = n_complex_kk**2

# --- 光学伝導度の計算 ---
# CGS単位系における関係: σ(ω) = [ω * Im(ε(ω))] / (4π)
sigma_KK = omega * np.imag(eps_kk) / (4 * np.pi)
# 理論値: 誘電率 eps_total から直接計算
sigma_theory = omega * np.imag(eps_total) / (4 * np.pi)

# --- プロット ---
# (1) KK変換から求めた光学伝導度
plt.figure(figsize=(10, 6))
plt.plot(omega, sigma_KK, label='Optical Conductivity (KK derived)', lw=2)
plt.xlabel("Wave number ω (cm$^{-1}$)")
plt.ylabel("σ(ω) [arb. units]")
plt.title("Optical Conductivity from KK Analysis (Drude + Lorentz)")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

# (2) 理論値 (eps_total から直接計算) との比較
plt.figure(figsize=(10, 6))
plt.plot(omega, sigma_KK, label='KK derived σ(ω)', lw=2)
plt.plot(omega, sigma_theory, '--', label='Direct σ(ω) from ε_total', lw=2)
plt.xlabel("Wave number ω (cm$^{-1}$)")
plt.ylabel("σ(ω) [arb. units]")
plt.title("Comparison of KK derived and Direct Optical Conductivity")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

---

解説¶

1. 誘電率の計算

   • Drude 項は従来通り eps_drude = eps_inf - omega_p**2/(omega**2 + iωγ) として計算し、
   • Lorentz 項は ( \varepsilon_{\mathrm{Lorentz}}(\omega)= \frac{S}{\omega_{0}^{2}-\omega^{2}-i\gamma_{L}\omega} ) として計算します。
   • これらを足し合わせたものを全体の誘電率 eps_total としています。

2. 反射率の計算

   • 複素屈折率は (\sqrt{\varepsilon_{\text{total}}}) で求め、Fresnel の式から反射振幅・反射率を算出しています。

3. クラマース・クローニッヒ変換

   • 並列化（joblib）を用いて、各周波数点での主値積分による位相計算を行っています。

4. 光学伝導度の計算と比較

   • KK変換により復元した ( \varepsilon(\omega) ) から計算した光学伝導度と、もともと eps_total から直接計算した伝導度を比較プロットしています。

このコードを実行することで、Drude項に Lorentz項を加えた場合の反射率および光学伝導度の計算と、その KK 変換による位相復元がどのように行われるかを確認できます。

In [ ]:

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import simps
from joblib import Parallel, delayed

# --- パラメータの設定 ---
omega_min = 1.0         # 最小周波数 [cm^-1]
omega_max = 100000.0     # 最大周波数 [cm^-1]
N = 20000               # 周波数の分割数
omega = np.linspace(omega_min, omega_max, N)  # 周波数配列

# --- Drudeモデルのパラメータ ---
omega_p = 15000.0       # Drudeのプラズマ振動数 [cm^-1]
gamma   = 500.0         # Drudeの散逸率 [cm^-1]
eps_inf = 1.0           # 高周波の背景誘電率

# --- Lorentzモデルのパラメータ ---
omega0_L = 20000.0      # Lorentz共鳴周波数 [cm^-1]
f = 0.3                 # オシレーター強度の比率
S = f * omega_p**2      # オシレーター強度
gamma_L = 1000.0        # Lorentz項の減衰定数 [cm^-1]

# --- 誘電率の計算 ---
# Drude項 (ω₀=0 の項)
eps_drude = eps_inf - omega_p**2 / (omega**2 + 1j * omega * gamma)

# Lorentz項
eps_lorentz = S / (omega0_L**2 - omega**2 - 1j * omega * gamma_L)

# 全体の誘電率 (Drude + Lorentz)
eps_total = eps_drude + eps_lorentz

# --- 複素屈折率と反射率の計算 ---
# 複素屈折率 n = sqrt(ε)
n_complex = np.sqrt(eps_total)
# 正入射での Fresnel の反射係数: r = (n-1)/(n+1)
r_complex_true = (n_complex - 1) / (n_complex + 1)
# 反射率 R = |r|^2
R = np.abs(r_complex_true)**2

# --- クラマース・クローニッヒ変換による位相の計算（並列化） ---
def compute_theta(i, omega, lnR):
    """
    インデックス i に対応する ω_i に対して、主値積分により
    位相 θ(ω_i) を計算する。

    θ(ω_i) = - (ω_i/π) P ∫_0∞  ln R(ω')/(ω'^2 - ω_i^2) dω'
    """
    w = omega[i]
    # 積分対象：全周波数で lnR/(ω'^2 - ω_i^2)
    integrand = lnR / (omega**2 - w**2)
    # 特異点 (j == i) は主値積分として 0 に置換
    integrand[i] = 0.0
    # Simpsonの公式による数値積分
    return - w/np.pi * simps(integrand, omega)

# lnR の事前計算
lnR = np.log(R)

# 並列処理: 各 ω_i での積分を全CPUコアで計算 (n_jobs=-1)
theta_list = Parallel(n_jobs=-1)(
    delayed(compute_theta)(i, omega, lnR) for i in range(len(omega))
)
theta = np.array(theta_list)

# --- KK変換から反射振幅の復元 ---
# 反射振幅: r(ω) = √R * exp(iθ)
r_complex_kk = np.sqrt(R) * np.exp(1j * theta)
# 逆Fresnel式より、r = (n-1)/(n+1) なので、解いて n = (1+r)/(1-r)
n_complex_kk = (1 + r_complex_kk) / (1 - r_complex_kk)
# 誘電率: ε = n^2
eps_kk = n_complex_kk**2

# --- 光学伝導度の計算 ---
# CGS単位系における関係: σ(ω) = [ω * Im(ε(ω))] / (4π)
sigma_KK = omega * np.imag(eps_kk) / (4 * np.pi)
# 理論値: 誘電率 eps_total から直接計算
sigma_theory = omega * np.imag(eps_total) / (4 * np.pi)

In [ ]:

# --- プロット ---
range_max = 20000
# (1) KK変換から求めた光学伝導度
plt.figure(figsize=(10, 6))
plt.plot(omega[0:range_max], sigma_KK[0:range_max], label='Optical Conductivity (KK derived)', lw=2)
plt.xlabel("Wave number ω (cm$^{-1}$)")
plt.ylabel("σ(ω) [arb. units]")
plt.title("Optical Conductivity from KK Analysis (Drude + Lorentz)")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

# (2) 理論値 (eps_total から直接計算) との比較
plt.figure(figsize=(10, 6))
plt.plot(omega[0:range_max], sigma_KK[0:range_max], label='KK derived σ(ω)', lw=2)
plt.plot(omega[0:range_max], sigma_theory[0:range_max], '--', label='Direct σ(ω) from ε_total', lw=2)
plt.xlabel("Wave number ω (cm$^{-1}$)")
plt.ylabel("σ(ω) [arb. units]")
plt.title("Comparison of KK derived and Direct Optical Conductivity")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

反射率のプロット¶

反射率 ( R(\omega) ) の波数依存性をプロットする方法¶

Drude成分とLorentz成分を組み合わせた誘電率 ( \varepsilon(\omega) ) から反射率 ( R(\omega) ) の波数依存性をプロットするには、以下の手順で行います。

---

手順¶

1. 誘電率を計算する

   • Drude項:
     [ \varepsilon_{\mathrm{Drude}}(\omega) = \varepsilon_{\infty} - \frac{\omega_p^2}{\omega^2 + i \omega \gamma} ]
   • Lorentz項:
     [ \varepsilon_{\mathrm{Lorentz}}(\omega) = \frac{S}{\omega_0^2 - \omega^2 - i \gamma_L \omega} ]
   • 総合的な誘電率: [ \varepsilon_{\mathrm{total}}(\omega) = \varepsilon_{\mathrm{Drude}}(\omega) + \varepsilon_{\mathrm{Lorentz}}(\omega) ]

2. 複素屈折率を求める

   • ( n_{\mathrm{complex}}(\omega) = \sqrt{\varepsilon_{\mathrm{total}}(\omega)} )

3. Fresnelの公式から反射係数を求める

   • 正入射におけるFresnel反射係数:
     [ r(\omega) = \frac{n(\omega) - 1}{n(\omega) + 1} ]
   • 反射率は ( R(\omega) = |r(\omega)|^2 )

4. 波数（周波数）依存性をプロットする

   • ( \omega )（波数）を横軸、( R(\omega) ) を縦軸にしてプロットする。

---

実装コード¶

以下のPythonコードを実行することで、Drude + Lorentz モデルに基づく反射率の波数依存性をプロットできます。

import numpy as np
import matplotlib.pyplot as plt

# --- パラメータの設定 ---
omega_min = 1.0        # 最小周波数 [cm^-1]
omega_max = 30000.0    # 最大周波数 [cm^-1]
N = 2000               # 周波数の分割数
omega = np.linspace(omega_min, omega_max, N)  # 周波数配列

# --- Drudeモデルのパラメータ ---
omega_p = 15000.0      # プラズマ振動数 [cm^-1]
gamma   = 500.0        # 散逸率 [cm^-1]
eps_inf = 1.0          # 高周波の背景誘電率

# --- Lorentzモデルのパラメータ ---
omega0_L = 20000.0     # Lorentz共鳴周波数 [cm^-1]
f = 0.3                # オシレーター強度の比率
S = f * omega_p**2     # オシレーター強度
gamma_L = 1000.0       # Lorentzの減衰定数 [cm^-1]

# --- 誘電率の計算 ---
# Drude項 (ω₀=0 の項)
eps_drude = eps_inf - omega_p**2 / (omega**2 + 1j * omega * gamma)

# Lorentz項
eps_lorentz = S / (omega0_L**2 - omega**2 - 1j * omega * gamma_L)

# 全体の誘電率 (Drude + Lorentz)
eps_total = eps_drude + eps_lorentz

# --- 複素屈折率と反射率の計算 ---
# 複素屈折率 n = sqrt(ε)
n_complex = np.sqrt(eps_total)

# Fresnelの反射係数: r = (n - 1) / (n + 1)
r_complex = (n_complex - 1) / (n_complex + 1)

# 反射率 R = |r|^2
R = np.abs(r_complex)**2

# --- 反射率の波数依存性をプロット ---
plt.figure(figsize=(10, 6))
plt.plot(omega, R, label='Reflectance $R(\omega)$', color='b', lw=2)
plt.xlabel("Wave number ω (cm$^{-1}$)")
plt.ylabel("Reflectance $R(\omega)$")
plt.title("Reflectance Spectrum (Drude + Lorentz Model)")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

---

プロットの解釈¶

• 低波数領域では Drude の影響が強く、反射率が高くなる傾向がある。
• 高波数領域では Lorentz 共鳴（( \omega_0 = 20000 ) cm(^{-1})）に対応する特徴的なピークが現れる。
• 散逸率 ( \gamma ) や Lorentz の ( \gamma_L ) を調整すると、共鳴幅が変化する。

---

パラメータ調整による変化¶

• Drude 項の影響を強くしたい場合

  • ( \omega_p )（プラズマ振動数）を増やすと、低周波数側の反射率が上がる。
  • ( \gamma )（散逸率）を増やすと、Drude成分の影響範囲が広がる。

• Lorentz 共鳴を調整したい場合

  • ( \omega_0 ) を変えると、共鳴の位置が変わる。
  • ( S )（オシレーター強度）を変えると、共鳴の強さが変わる。
  • ( \gamma_L ) を増やすと、共鳴の幅が広がる。

---

まとめ¶

• Drude + Lorentz モデルに基づく反射率 ( R(\omega) ) の波数依存性をプロットする方法を示した。
• Fresnel の公式を用いて反射率を計算し、波数の関数としてプロット。
• パラメータ調整により、低波数での Drude の寄与や、高波数での Lorentz 共鳴の寄与を調整可能。

この方法を使えば、実験データと比較しながら物性パラメータを最適化することができます！ 🚀

In [ ]:

import numpy as np
import matplotlib.pyplot as plt

# --- パラメータの設定 ---
omega_min = 1.0        # 最小周波数 [cm^-1]
omega_max = 30000.0    # 最大周波数 [cm^-1]
N = 2000               # 周波数の分割数
omega = np.linspace(omega_min, omega_max, N)  # 周波数配列

# --- Drudeモデルのパラメータ ---
omega_p = 15000.0      # プラズマ振動数 [cm^-1]
gamma   = 500.0        # 散逸率 [cm^-1]
eps_inf = 1.0          # 高周波の背景誘電率

# --- Lorentzモデルのパラメータ ---
omega0_L = 20000.0     # Lorentz共鳴周波数 [cm^-1]
f = 0.3                # オシレーター強度の比率
S = f * omega_p**2     # オシレーター強度
gamma_L = 1000.0       # Lorentzの減衰定数 [cm^-1]

# --- 誘電率の計算 ---
# Drude項 (ω₀=0 の項)
eps_drude = eps_inf - omega_p**2 / (omega**2 + 1j * omega * gamma)

# Lorentz項
eps_lorentz = S / (omega0_L**2 - omega**2 - 1j * omega * gamma_L)

# 全体の誘電率 (Drude + Lorentz)
eps_total = eps_drude + eps_lorentz

# --- 複素屈折率と反射率の計算 ---
# 複素屈折率 n = sqrt(ε)
n_complex = np.sqrt(eps_total)

# Fresnelの反射係数: r = (n - 1) / (n + 1)
r_complex = (n_complex - 1) / (n_complex + 1)

# 反射率 R = |r|^2
R = np.abs(r_complex)**2

# --- 反射率の波数依存性をプロット ---
plt.figure(figsize=(10, 6))
plt.plot(omega, R, label='Reflectance $R(\omega)$', color='b', lw=2)
plt.xlabel("Wave number ω (cm$^{-1}$)")
plt.ylabel("Reflectance $R(\omega)$")
plt.title("Reflectance Spectrum (Drude + Lorentz Model)")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

In [ ]:

import numpy as np
import matplotlib.pyplot as plt

# --- パラメータの設定 ---
omega_min = 1.0        # 最小周波数 [cm^-1]
omega_max = 30000.0    # 最大周波数 [cm^-1]
N = 2000               # 周波数の分割数
omega = np.linspace(omega_min, omega_max, N)  # 周波数配列

# --- 共通パラメータ ---
gamma   = 500.0        # Drudeモデルの散逸率 [cm^-1]
eps_inf = 1.0          # 高周波の背景誘電率

# Lorentzモデルのパラメータ（固定）
omega0_L = 20000.0     # Lorentz共鳴周波数 [cm^-1]
f = 0.3                # オシレーター強度の比率
gamma_L = 1000.0       # Lorentzの減衰定数 [cm^-1]

# --- プラズマ振動数のリスト ---
omega_p_list = [10000.0, 15000.0, 20000.0]

# プロットの準備
plt.figure(figsize=(10, 6))

# 各omega_pについて計算とプロットを実施
for omega_p in omega_p_list:
    # Drudeモデルの項
    eps_drude = eps_inf - omega_p**2 / (omega**2 + 1j * omega * gamma)

    # Lorentzモデルの項
    S = f * omega_p**2  # オシレーター強度
    eps_lorentz = S / (omega0_L**2 - omega**2 - 1j * omega * gamma_L)

    # 全体の誘電率 (Drude + Lorentz)
    eps_total = eps_drude + eps_lorentz

    # 複素屈折率の計算: n = sqrt(ε)
    n_complex = np.sqrt(eps_total)

    # Fresnelの反射係数: r = (n - 1) / (n + 1)
    r_complex = (n_complex - 1) / (n_complex + 1)

    # 反射率 R = |r|^2
    R = np.abs(r_complex)**2

    # omega_p の値をラベルに含めてプロット
    plt.plot(omega, R, lw=2, label=f"$\omega_p$ = {omega_p} cm$^{{-1}}$")

# グラフの装飾
plt.xlabel("Wave number ω (cm$^{-1}$)")
plt.ylabel("Reflectance $R(\omega)$")
plt.title("Reflectance Spectrum (Drude + Lorentz Model)")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

In [ ]:

import numpy as np
import matplotlib.pyplot as plt

# --- パラメータの設定 ---
omega_min = 1.0         # 最小周波数 [cm^-1]
omega_max = 30000.0     # 最大周波数 [cm^-1]
N = 2000                # 周波数の分割数
omega = np.linspace(omega_min, omega_max, N)  # 周波数配列

# --- Drudeモデルのパラメータ ---
omega_p = 15000.0       # プラズマ振動数 [cm^-1]
gamma   = 500.0         # 散逸率 [cm^-1]
eps_inf = 1.0           # 高周波の背景誘電率

# --- Lorentzモデルのパラメータ ---
omega0_L = 20000.0      # Lorentz共鳴周波数 [cm^-1]
f = 0.3                 # オシレーター強度の比率
S = f * omega_p**2      # オシレーター強度
gamma_L = 1000.0        # Lorentzの減衰定数 [cm^-1]

# --- 誘電率の計算 ---
# Drude項 (ω₀=0 の項)
eps_drude = eps_inf - omega_p**2 / (omega**2 + 1j * omega * gamma)

# Lorentz項
eps_lorentz = S / (omega0_L**2 - omega**2 - 1j * omega * gamma_L)

# 全体の誘電率 (Drude + Lorentz)
eps_total = eps_drude + eps_lorentz

# --- 反射率の計算 ---
# 複素屈折率 n = sqrt(ε)
n_complex = np.sqrt(eps_total)
# Fresnelの反射係数: r = (n - 1) / (n + 1)
r_complex = (n_complex - 1) / (n_complex + 1)
# 反射率 R = |r|^2
R = np.abs(r_complex)**2

# --- 光学伝導度の計算 ---
# ε(ω) = ε_inf + iσ(ω)/(ω)  (ここでは単位や ε0 を省略している)
# よって、σ(ω) = -i ω [ε(ω) - ε_inf]
sigma = -1j * omega * (eps_total - eps_inf)
sigma_real = np.real(sigma)
sigma_imag = np.imag(sigma)

# --- グラフのプロット (反射率と光学伝導度を横に並べる) ---
fig, axs = plt.subplots(1, 2, figsize=(14, 6))

# 左側: 反射率のプロット
axs[0].plot(omega, R, lw=2, color='blue')
axs[0].set_xlabel("Wave number ω (cm$^{-1}$)")
axs[0].set_ylabel("Reflectance $R(ω)$")
axs[0].set_title("Reflectance Spectrum (Drude + Lorentz)")
axs[0].grid(True)

# 右側: 光学伝導度のプロット（実部と虚部）
axs[1].plot(omega, sigma_real, lw=2, label='Re σ', color='red')
axs[1].plot(omega, sigma_imag, lw=2, label='Im σ', color='green')
axs[1].set_xlabel("Wave number ω (cm$^{-1}$)")
axs[1].set_ylabel("Optical Conductivity σ (arb. units)")
axs[1].set_title("Optical Conductivity Spectrum")
axs[1].legend()
axs[1].grid(True)

plt.tight_layout()
plt.show()

In [ ]:

import numpy as np
import matplotlib.pyplot as plt

# --- パラメータの設定 ---
omega_min = 1.0         # 最小周波数 [cm^-1]
omega_max = 30000.0     # 最大周波数 [cm^-1]
N = 2000                # 周波数の分割数
omega = np.linspace(omega_min, omega_max, N)  # 周波数配列

# --- 共通パラメータ ---
omega_p = 15000.0       # プラズマ振動数 [cm^-1]
eps_inf = 1.0           # 高周波の背景誘電率

# Lorentzモデルのパラメータ（固定）
omega0_L = 20000.0      # Lorentz共鳴周波数 [cm^-1]
f = 0.3                 # オシレーター強度の比率
S = f * omega_p**2      # オシレーター強度
gamma_L = 1000.0        # Lorentzの減衰定数 [cm^-1]

# --- Drude成分の散逸率γのリスト ---
gamma_list = [100.0, 500.0, 1000.0]

# --- プロットの準備 ---
fig, axs = plt.subplots(1, 2, figsize=(14, 6))
colors = ['red', 'blue', 'green']  # 各γに対応する色

# --- 各γでの計算とプロット ---
for i, gamma in enumerate(gamma_list):
    # Drudeモデルの項（Lorentz項とは独立にγのみ変更）
    eps_drude = eps_inf - omega_p**2 / (omega**2 + 1j * omega * gamma)

    # Lorentzモデルの項
    eps_lorentz = S / (omega0_L**2 - omega**2 - 1j * omega * gamma_L)

    # 全体の誘電率 (Drude + Lorentz)
    eps_total = eps_drude + eps_lorentz

    # --- 反射率の計算 ---
    # 複素屈折率 n = √ε
    n_complex = np.sqrt(eps_total)
    # Fresnelの反射係数: r = (n - 1) / (n + 1)
    r_complex = (n_complex - 1) / (n_complex + 1)
    # 反射率 R = |r|^2
    R = np.abs(r_complex)**2

    # --- 光学伝導度の計算 ---
    # σ(ω) = -i ω [ε(ω) - ε_inf]
    sigma = -1j * omega * (eps_total - eps_inf)
    sigma_real = np.real(sigma)
    # sigma_imag = np.imag(sigma)

    # --- 反射率のプロット ---
    axs[0].plot(omega, R, lw=2, color=colors[i],
                label=f'γ = {gamma:.0f} cm⁻¹')

    # --- 光学伝導度のプロット ---
    # 実部は実線、虚部は破線でプロット
    axs[1].plot(omega, sigma_real, lw=2, color=colors[i],
                linestyle='-', label=f'Re σ, γ = {gamma:.0f}')
    # axs[1].plot(omega, sigma_imag, lw=2, color=colors[i],
    #             linestyle='--', label=f'Im σ, γ = {gamma:.0f}')

# --- グラフの装飾 ---
# 反射率プロット
axs[0].set_xlabel("Wave number ω (cm⁻¹)")
axs[0].set_ylabel("Reflectance R(ω)")
axs[0].set_title("Reflectance Spectrum (varying Drude γ)")
axs[0].grid(True)
axs[0].legend(fontsize=9)

# 光学伝導度プロット
axs[1].set_xlabel("Wave number ω (cm⁻¹)")
axs[1].set_ylabel("Optical Conductivity σ (arb. units)")
axs[1].set_title("Optical Conductivity Spectrum (varying Drude γ)")
axs[1].grid(True)
axs[1].legend(fontsize=9, ncol=2)

plt.tight_layout()
plt.show()

In [ ]:

import numpy as np
import matplotlib.pyplot as plt

# --- パラメータの設定 ---
omega_min = 1.0         # 最小周波数 [cm^-1]
omega_max = 30000.0     # 最大周波数 [cm^-1]
N = 2000                # 周波数の分割数
omega = np.linspace(omega_min, omega_max, N)  # 周波数配列

# --- 共通パラメータ ---
gamma = 500.0           # Drudeモデルの散逸率 [cm^-1]
eps_inf = 1.0           # 高周波の背景誘電率

# Lorentzモデルのパラメータ（固定）
omega0_L = 20000.0      # Lorentz共鳴周波数 [cm^-1]
f = 0.3                 # オシレーター強度の比率
gamma_L = 1000.0        # Lorentzの減衰定数 [cm^-1]

# --- プラズマ周波数のリスト ---
omega_p_list = [10000.0, 15000.0, 20000.0]

# --- プロットの準備 ---
fig, axs = plt.subplots(1, 2, figsize=(14, 6))
colors = ['red', 'blue', 'green']  # 各ωₚに対応する色

# --- 各プラズマ周波数での計算とプロット ---
for i, omega_p in enumerate(omega_p_list):
    # Lorentz項のオシレーター強度 S = f * ωₚ²
    S = f * omega_p**2

    # Drude項 (ω₀=0 の項)
    eps_drude = eps_inf - omega_p**2 / (omega**2 + 1j * omega * gamma)

    # Lorentz項
    eps_lorentz = S / (omega0_L**2 - omega**2 - 1j * omega * gamma_L)

    # 全体の誘電率 (Drude + Lorentz)
    eps_total = eps_drude + eps_lorentz

    # --- 反射率の計算 ---
    # 複素屈折率 n = √ε
    n_complex = np.sqrt(eps_total)
    # Fresnelの反射係数: r = (n - 1) / (n + 1)
    r_complex = (n_complex - 1) / (n_complex + 1)
    # 反射率 R = |r|²
    R = np.abs(r_complex)**2

    # --- 光学伝導度の計算 ---
    # σ(ω) = -i ω [ε(ω) - ε_inf]
    sigma = -1j * omega * (eps_total - eps_inf)
    sigma_real = np.real(sigma)
    # sigma_imag = np.imag(sigma)

    # --- 反射率のプロット ---
    axs[0].plot(omega, R, lw=2, color=colors[i],
                label=f'$\omega_p$ = {omega_p:.0f} cm$^{{-1}}$')

    # --- 光学伝導度のプロット ---
    # 実部は実線、虚部は破線でプロット
    axs[1].plot(omega, sigma_real, lw=2, color=colors[i],
                linestyle='-', label=f'Re $\sigma$, $\omega_p$ = {omega_p:.0f}')
    # axs[1].plot(omega, sigma_imag, lw=2, color=colors[i],
    #             linestyle='--', label=f'Im $\sigma$, $\omega_p$ = {omega_p:.0f}')

# --- グラフの装飾 ---
# 反射率プロット
axs[0].set_xlabel("Wave number ω (cm$^{-1}$)")
axs[0].set_ylabel("Reflectance R(ω)")
axs[0].set_title("Reflectance Spectrum (varying plasma frequency)")
axs[0].grid(True)
axs[0].legend(fontsize=9)

# 光学伝導度プロット
axs[1].set_xlabel("Wave number ω (cm$^{-1}$)")
axs[1].set_ylabel("Optical Conductivity σ (arb. units)")
axs[1].set_title("Optical Conductivity Spectrum (varying plasma frequency)")
axs[1].grid(True)
axs[1].legend(fontsize=9, ncol=2)

plt.tight_layout()
plt.show()