クラマース・クローニッヒの関係(Kramers-kronig relaion)についての直感的理解

はじめに

クラマース・クローニッヒの関係(Kramers-kronig relaion)とは線形応答における周波数応答関数の実部虚部がヒルベルト変で関係づけられることを示した式である。1926年にラルフ・クローニッヒ、1927年にヘンリク・アンソニー・クラマースによって電磁波の分散現象に対して導かれた。(クラマース・クローニッヒの関係式 – Wikipedia)

要は屈折率と吸収(※もちろんこれだけには限らない)がヒルベルト変換の関係にあるのである。この関係は今まであまり直感的な理解はしておらず、なんとなくこんな関係があるんだなと理解していた程度だった。
しかし最近になって、会社の同期におすすめの講義用スライドを教えてもらい腑に落ちたので、今回はその紹介と自分の理解を書いておきたい。他にも光学関連の記事を書いているので是非参考にしてください。

本編

Kramars-Kronigの関係とは

Kramars-Kronigの関係の直感的理解の前に基本的な理解を教科書を借りて行う[1]。
光の電場\(E\)に対して実数部分を\(E^{(r)}(t)(- \infty \leq t \leq \infty)\)とする。ここでこれを自乗可積分と仮定して以下のようなフーリエ変換の形
\begin{equation}
E^{(r)}(t) = \int_0^\infty a(\nu ) \cos [\phi(\nu )-2\pi\nu t]d\nu
\end{equation}
で表すことができ、この実数部分に対応して、複素関数
\begin{equation}
E(t) = \int_0^\infty a(\nu ) \exp \left( i[\phi(\nu )-2\pi\nu t] \right) d\nu
\end{equation}
を考えれば
\begin{equation}
E(t) = E^{(r)}(t) + i E^{(i)}(t)
\end{equation}
が成立する。ただし、\(E^{(i)}(t)\)は
\begin{equation}
E^{(i)}(t) = \int_0^\infty a(\nu ) \sin [\phi(\nu )-2\pi\nu t]d\nu
\end{equation}
で与えられる。これらの式から実数部分と虚数部分はsin関数かcos関数かの違いしかなく\( \phi(\nu) \)を\( \phi(\nu) – 2\pi \)で置き換えることによって得られるので、実数部分\(E^{(r)}(t)\)と虚数部分\(E^{(i)}(t)\)ないしは複素関数\(E(t)\)の一つがわかれば他の部分も一意に定まることがわかる。また実数部分\(E^{(r)}(t)\)と虚数部分\(E^{(i)}(t)\)は双対Fourier積分または共役関数ともよばれ互いにヒルベルト変換
\begin{align}
E^{(i)}(t) &= \frac{1}{\pi} P \int_{-\infty }^\infty \frac{ E^{(r)}(t’)}{t’-t} dt’ \\
E^{(r)}(t) &= \frac{-1}{\pi} P \int_{-\infty }^\infty \frac{ E^{(i)}(t’)}{t’-t} dt’
\end{align}
の関係にある。

直感的理解

今回教えていただいた資料は東京大学の種村拓夫准教授の授業用スライド(OE2019_tanemura_slides1-4.pdf)である。他の資料内容も興味深かったが特に今回はKramers-Kronigの関係周辺を学んでいく。

光と物質は以下の図(Polarization Diagram Images, Photos | Mungfali)のように電場\(E\)内では物質の分子\(P\)が分極を作ることで相互作用を行う。これはドルーデモデルや、ローレンツモデルと呼ばれる方程式を解くことで近似的な数式理解はできるが、ここでは触れない。

ここで電場と分極の相互作用を考える。周波数領域で分極\(P(\omega)\)は
\begin{equation}
P(\omega) = \varepsilon_0 \chi (\omega )E(\omega )
\end{equation}
の関係にあり(\(\chi (t)\)は複素感受率)、時間領域においては
\begin{equation}
P(t) = \varepsilon_0 \int_{-\infty }^\infty \chi (t’)E(t-t’ )dt’
\end{equation}
で表される。ここで周波数領域と時間領域の式の変化は畳み込み積とフーリエ変換の関係を用いており(参考:畳み込み積分（合成積）とは？フーリエ変換したらどうなる？ | 理系大学生の数学駆け込み寺)それぞれの()内の変化はフーリエ変換の関係にある(例: \( \mathscr{F}[E(t)] = E(\omega )\) )。ここで注目されたいのが電場をインパルス関数と仮定して( \( E(t) = \delta (t)\) )みると分極は \( P(t) = \chi (t)\)となることが時間領域の分極の関係からわかる。
さらに、この\(\chi (t)\)をフーリエ変換すると
\begin{align}
\chi (t) &= \frac{1}{2\pi} \int_{-\infty }^\infty \chi (\omega ) \exp (-i\omega t)d\omega\\
&= \frac{1}{2\pi} \int_{-\infty }^\infty [\chi_r (\omega ) \cos \omega t + \chi_i(\omega ) \sin \omega t ] + i[\chi_i (\omega ) \cos \omega t – \chi_r(\omega ) \sin \omega t ] d \omega
\end{align}

特に今、実分極\(P(t)\)の応答を考えているので、\(\chi(t)\)も実数になる。
よって、虚数部分は0になるので
\begin{align}
0 &= \int_{-\infty }^\infty i[\chi_r (\omega ) \cos \omega t – \chi_i(\omega ) \sin \omega t ] d \omega
\end{align}
ここで、それぞれの項に注目すると第一項は\( \cos × \chi_i \)の正負無限の積分、第二項は\( \sin × \chi_r \)の正負無限の積分なので、奇関数の積分が0になることに注目する(参考:偶関数と奇関数の積分 – 微分積分 – 基礎からの数学入門)と
\begin{equation}
\left\{
\begin{array}{rl}
\chi_i(\omega) &= – \chi_i(-\omega )  (奇関数)  \\
\chi_r(\omega) &= \chi_r(-\omega )  (偶関数)
\end{array}
\right.
\end{equation}
よって\(\chi(\omega)^* = \chi(-\omega)\)となる。ここで\(\chi (t)\)を再度整理すると
\begin{equation}
\chi (t) = \frac{1}{2\pi} \int_{-\infty }^\infty [\chi_r (\omega ) \cos \omega t + \chi_i(\omega ) \sin \omega t ] d \omega
\end{equation}
ここでフーリエサイン変換、フーリエコサイン変換(物理とか-偶関数・奇関数のフーリエ変換)を思い出すとフーリエ空間上(周波数空間上)でのcosの項\(\chi_r (\omega )\)、sinの項\(\chi_i (\omega )\)は実空間上では偶関数\(\chi_r (t )\)、奇関数\(\chi_i (t)\)となる。これを数式に落とし込むと
\begin{equation}
\left\{
\begin{array}{rl}
\chi(t) &\equiv \chi_{even}(t) + \chi_{odd} (t)\\
\chi_{even} (t) &\equiv \frac{1}{2\pi} \int_{-\infty }^\infty \chi_r (\omega ) \exp (-i\omega t)d\omega\\
\chi_{odd} (t) &\equiv \frac{1}{2\pi} \int_{-\infty }^\infty i\chi_i (\omega ) \exp (-i\omega t)d\omega
\end{array}
\right.
\end{equation}
と書ける。

これらをグラフで確認する

import numpy as np
import matplotlib.pyplot as plt

# パラメータ
alpha = 1.0   # 減衰率
omega = 10.0   # 振動の周波数
t = np.linspace(-5, 5, 400)  # 時間軸

# 1. インパルス関数（デルタ関数の近似）
t_impulse = [0]
y_impulse = [1]
plt.figure(figsize=(10, 6))

plt.subplot(2, 1, 1)
plt.stem(t_impulse, y_impulse)
plt.title("E(t) = δ(t)")
plt.xlabel("t")
plt.ylabel("Amplitude")
plt.grid()

# 2. t=0からの緩和振動 (t >= 0)
y_relax = np.heaviside(t, 1) * np.exp(-alpha * t) * np.cos(omega * t)
plt.subplot(2, 1, 2)
plt.plot(t, y_relax)
plt.title("P(t)/ε0 = χ(t)")
plt.xlabel("t")
plt.ylabel("Amplitude")
plt.grid()

plt.figure(figsize=(10, 6))
# 3. 偶関数での緩和振動
y_even = np.exp(-alpha * np.abs(t)) * np.cos(omega * t)
plt.subplot(2, 1, 1)
plt.plot(t, y_even)
plt.title("χ_even")
plt.xlabel("t")
plt.ylabel("Amplitude")
plt.grid()

# 4. 奇関数での緩和振動
y_odd = np.sign(t) * np.exp(-alpha * np.abs(t)) * np.cos((omega) * t)
plt.subplot(2, 1, 2)
plt.plot(t, y_odd)
plt.title("χ_odd")
plt.xlabel("t")
plt.ylabel("Amplitude")
plt.grid()

plt.tight_layout()
plt.show()

#確認用
plt.plot(t, y_odd+y_even)
plt.title("χ_odd")
plt.xlabel("t")
plt.ylabel("Amplitude")
plt.grid()

上記の図が、下式のうち１枚目が電場(今回はδ関数)と２枚目が分極の変化
\begin{equation}
P(\omega) = \varepsilon_0 \chi (\omega )E(\omega )
\end{equation}

上記の図が、下式のうち１枚目が\(\chi_r (t )\)と２枚目が\(i\chi_i (t )\)を意味している。
この図からKramers-Kronigの関係の重要な特性が読み取れる。分極の図というのは改めて以下の式で表される。
\begin{equation}
P(t) = \varepsilon_0 \int_{-\infty }^\infty (\chi_r (t’ ) + \chi_i (t’ ))E(t-t’ )dt’
\end{equation}
つまり、\( P(t) ∝ \chi_r (t ) + \chi_i (t)\)となるのである。ここで\(P(t)\)というのは電場が入射した際に初めて誘起されるものであり、\(t<0\)ではあくまで0になっている。(これは当たり前の因果律であり、電場が入る前に分極が誘起されるのがおかしいというのは当たり前であるからだ)

ということは以下の式が成り立つことがわかる。
\begin{align}
P(t<0) &= 0 \\
\chi_r(t<0) + \chi_i(t<0) &=0
\end{align}

これを満たすのにsgn関数(符号関数)という関数を導入することでうまく説明ができる。

\(\chi_{even} (t )\)と\(\chi_{odd} (t )\)は以下の関係が成り立つ。
\begin{equation}
\left\{
\begin{array}{rl}
\chi_{even} (t) = \chi_{odd}(t)\cdot　sgn (t) \\
\chi_{odd} (t) = \chi_{even}(t)\cdot　sgn (t)
\end{array}
\right.
\end{equation}
ここでsgn関数のフーリエ変換は
\begin{equation}
\mathscr{F}[sgn] = \left\{
\begin{array}{rl}
i\frac{2}{\omega} (\omega \ne 0) \\
0 (\omega = 0)
\end{array}
\right.
\end{equation}
よってフーリエ変換の畳み込み定理からKramers-Kronigの関係が導かれる。
\begin{align}
\chi_r(\omega) &= -\frac{1}{\pi} P.V. \int_{-\infty }^\infty \frac{\chi_i (\omega ‘)}{\omega -\omega’} d\omega’ \\
\chi_i(\omega) &= \frac{1}{\pi} P.V. \int_{-\infty }^\infty \frac{\chi_r (\omega ‘)}{\omega -\omega’} d\omega’
\end{align}

実際のKramars-Kronigの関係のイメージ

例えば、実際のレーザー媒質などでは以下のような吸収断面積となる。

Ybファイバーレーザーの吸収/誘導放出断面積(参考:Paschotta_2005_ytterbium_doped_laser_gain_media)

これを踏まえたうえで実際に媒質の吸収付近の屈折率分布がどのようになっているかを計算してみると以下のようなグラフになる。

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad
import japanize_matplotlib

# パラメータ設定
gamma = 1.0  # 減衰係数
omega_0 = 3.0  # 共鳴周波数

# 虚部 (Lorentzian)
def ImF(omega,intensity=1):
    return intensity * gamma / ((omega - omega_0)**2 + gamma**2)

# 実部
def ReF(omega,intensity=1):
    def integrand(omega_prime):
        return ImF(omega_prime,intensity) / (omega_prime - omega)
    
    # 特異点 ω' = ω を避けるために分割して積分
    eps = 1e-6  # 小さな値を設定
    integral1, _ = quad(integrand, -10, omega - eps)  # 左側
    integral2, _ = quad(integrand, omega + eps, 10)  # 右側
    
    return (integral1 + integral2) / np.pi

# 周波数軸
omega_values = np.linspace(-5, 10, 500)

# 実部・虚部を計算
Re_values = np.array([ReF(omega) for omega in omega_values])
Im_values = np.array([ImF(omega) for omega in omega_values])

# 虚部が大きい場合の実部・虚部を計算
Re_values_big = np.array([ReF(omega,intensity=2) for omega in omega_values])
Im_values_big = np.array([ImF(omega,intensity=2) for omega in omega_values])

# 虚部が普通の場合/大きい場合のグラフ描画
plt.figure(figsize=(8, 5))
plt.plot(omega_values, Re_values, label="Re[F(ω)] (実部)", color="b")
plt.plot(omega_values, Im_values, label="Im[F(ω)] (虚部)", color="r")
plt.plot(omega_values, Re_values_big, label="Re[F(ω)_big] (実部)", color="b",linestyle="dashed")
plt.plot(omega_values, Im_values_big, label="Im[F(ω)_big] (虚部)", color="r", linestyle="dashed")
plt.axhline(0, color="black", linewidth=0.5)
plt.axvline(omega_0, color="orange", linestyle="dotted", label="ω_0")
plt.xlabel("ω (周波数)")
plt.ylabel("F(ω)")
plt.title("Kramers-Kronig関係 (ローレンツ関数)")
plt.legend()
plt.grid()
plt.show()

吸収の様子をローレンツ関数で近似してあげた結果が以上のグラフとなる。F(ω)はローレンツ関数のピークが1のときでF(ω)_bigがローレンツ関数のピークが2のときである。このグラフから屈折率(\( n \sim \sqrt{1+\chi_r} \))と吸収断面積(\(\sigma_{abs} \sim \chi_i\))の関係について以下のことがわかる。

参考文献