はじめに
光の勉強をしているとラマン散乱やブルリアン散乱といった言葉に出会うことがある。特に光ファイバーなどの光の閉じ込め作用や作用長が長いときに、この現象は顕著になる。このとき、これらのラマン散乱とは光学フォノンによるもの、ブルリアン散乱とは音響フォノンによるものといった説明がされている文献が見受けられる[1](この説明がどこまで正しいかは保証しかねるところだが…)。
しかし、自分の中で大学の授業で受けっぱなしになっていたところもあって、上記説明はよく咀嚼しきれておらず、特に光学フォノン、音響フォノンといった違いはよくわからなかった。そこで本記事はその説明を行っていくとともにJuliaを用いた具体的な計算例も含めた検証も行っていく。
NaCl型結晶の1次元格子振動
今回は理工学基礎 物性化学といった教科書[2]を参考に説明を行う。こちらの教科書は無機化学初心者や他分野の人間にとっても非常に読みやすく、入門書に向いているので勧めたい。他にも固体化学や無機化学というとキッテルや無機化学なども読みやすいが、入門に近い位置づけにありながら少し難しくも感じた。
また、今回の内容については似た内容の授業スライドが上がっており、東京農工大の佐藤先生のスライドや東大の塩見先生の補助資料などがわかりやすかった。
今回紹介する教科書では「光学的振動」「音響的振動」として触れられていたがJAEAの記事等によるとどうやら「光学フォノン」「音響フォノン」に相当するようである。今回は教科書に倣って「光学的振動」「音響的振動」に統一する。
今回はNaCl型の結晶での1次元格子の振動を考える。図のような同一線上にあり、波数ベクトル\(\rm{q}\)の方向を正としたとき、\(u_n\),\(v_n\)をCl,Naの\(n\)個めの原子の変位として、\(M,m\)をCl,Naの質量とすると以下のような運動方程式が成り立つ。
\begin{equation} M \frac{d^2 u_n}{dt^2} = b(v_n + v_{n-1} -2u_n)\\ m \frac{d^2 v_n}{dt^2} = b(u_n + u_{n-1} -2v_n) \end{equation} ここで結合方向への伸縮であることから\(b\)は結合伸縮力定数(bond-stretching force constant)またはばね定数と呼んだりする。この微分方程式の解の形は、格子定数\(a\)を用いて(~\(na\)が\(x\)に相当) \begin{equation} u_n(na,t) = A e^{i(\omega t-qna)}\\ v_n(na,t) = B e^{i(\omega t-qna)} \end{equation} とすると微分方程式に代入して \begin{equation} -M\omega^2A = bB(1 + e^{iqa})-2bA\\ -m\omega^2B = bA(1 + e^{-iqa})-2bB \end{equation} ここで\(A,B\)がともに0にならないため(変位が0ではない時を考えている)、\(A,B\)の係数についての次の行列式が0になることを確認する。 \begin{equation} \begin{vmatrix} 2b-M\omega^2 && -b(1+e^{iqa})\\ -b(1+e^{-iqa}) && 2b-m\omega^2 \end{vmatrix} ＝ 0 \end{equation} となる。\(b,M,m,a\)は定数なため、上式は以下のような\(\omega^2\)の解を持つ。 \begin{equation} \omega^2(q) = b\left( \frac{1}{M} + \frac{1}{m} \right) \pm b\left\{ \left( \frac{1}{M} + \frac{1}{m} \right)^2 – \frac{4}{Mm} \sin\left( \frac{qa}{2}\right)^2 \right\}^{1/2} \end{equation} この\(q\)と\(\omega\)の関係を分散関係という。右辺の±について+をとったほうを\(\omega_+\)、-のほうを\(\omega_-\)とする。ここでさらにこの数式を深堀していくと \(qa = 0\)の場合： \begin{equation} \omega_+ = \left\{2b\left( \frac{1}{M} + \frac{1}{m} \right)\right\}^{1/2}\\ \omega_- = 0 \end{equation} \(qa \approx 0\)の場合[\( \sin (qa/2) \approx qa/2\)]： \begin{equation} \omega_+ \approx \left\{2b\left( \frac{1}{M} + \frac{1}{m} \right)\right\}^{1/2}\\ \omega_- \approx \left( \frac{b/2}{M+m} \right) qa \end{equation} \(qa = \pi\)の場合： \begin{equation} \omega_+ = \left(\frac{2b}{m} \right)^{1/2}\\ \omega_- = \left(\frac{2b}{M} \right)^{1/2} \end{equation}Juliaを使った分散関係の計算
計算で求めた分散関係についてさらに深堀していこうと思う。実際のNaCl結晶使って実際に計算を行ってみる。
今回の計算で使う定数は以下のようになる。
| 項目 | 値 | 単位 | 出典 |
|---|---|---|---|
| NaClの結合伸縮力定数 \(b\) | 40.6 | N/m | Spec. Chim. Acta A (2016) |
| アボガドロ数 \(N_A\) | 6.022×10²³ | mol⁻¹ | 定数 |
| Naの質量 \(m_{\rm Na}\) | 22.99×10⁻³ / \(N_A\) | kg | 文部科学省 |
| Clの質量 \(m_{\rm Cl}\) | 35.45×10⁻³ / \(N_A\) | kg | 文部科学省 |
| NaClの格子定数 \(a\) | 5.67×10⁻¹⁰ | m | J-STAGE |
using PyPlot
using PyCall
PyCall.PyDict(matplotlib["rcParams"])["font.family"] = "Meiryo";
# define fuction:
ω²₊(b,M,m,qa) = b*(1/M + 1/m) + b*( (1/M + 1/m)^2 - 4/(M*m)*sin(qa/2)^2 )^(1/2)
ω²₋(b,M,m,qa) = b*(1/M + 1/m) - b*( (1/M + 1/m)^2 - 4/(M*m)*sin(qa/2)^2 )^(1/2)
# setting parameters:
# NaCl bond-stretching force constant: 43.3 or 40.6 or 47.5 N/m
# https://diverdi.colostate.edu/C372/experiments/quantum%20computation%20of%20vibrational%20spectra/references/spec_chim_acta_a_2016_v154_p103.pdf?utm_source=chatgpt.com
b = 40.6 # [N/m]
#avogadro's number
NA = 6.022e23 # [mol⁻¹]
#Na's mass: 22.99 g/mol
M_Na = 22.99e-3 / NA # [kg]
#Cl's mass: 35.45 g/mol
M_Cl = 35.45e-3 / NA # [kg]
#http://www.famic.go.jp/ffis/fert/sub6_data/bunatomtable.html
# NaCl lattice constant: 5.67 Å
# https://www.jstage.jst.go.jp/article/jssep/36/0/36_480/_pdf
a = 5.67e-10 # [m]
q_list = [π/a * (i/100) for i in 0:100] # [1/m]
# finished parameter setting
Optical_phonon = map(q -> √(ω²₊(b,M_Na,M_Cl,q*a))/(2π), q_list) # [1/s]
Acoustic_phonon = map(q -> √(ω²₋(b,M_Na,M_Cl,q*a))/(2π), q_list) # [1/s]
# Plot:
fig, ax = subplots()
ax.plot(q_list .* a, Optical_phonon .* 1e-12, label=L"$\omega_+$:Transverse Optical phonon (TO)", color="blue")
ax.plot(q_list .* a, Acoustic_phonon .* 1e-12, label=L"$\omega_-$:Transverse Acoustic phonon (TA)", color="red")
ax.set_xlabel("qa (1/m)")
ax.set_xlim(0, π)
ax.set_ylim(0)
ax.grid(true)
ax.legend()
ax.set_ylabel("Frequency (THz)")
ax_top = ax.twiny()
ax_top.set_xlim(ax.get_xlim())
# 第一ブルリアンゾーンの説明：
# https://ocw.nagoya-u.jp/files/109/lecnote_08.pdf
ax_top.set_xticks([0, π])
ax_top.set_xticklabels([L"\Gamma", "X"])
# 各特徴的な線の記述
axhline((2b*(1/M_Na + 1/M_Cl))^(1/2)/(2π)*1e-12,color="black",linestyle="--")
axhline((2b*(1/M_Na))^(1/2)/(2π)*1e-12,color="black",linestyle="--")
axhline((2b*(1/M_Cl))^(1/2)/(2π)*1e-12,color="black",linestyle="--")
text(y=(2b*(1/M_Na + 1/M_Cl))^(1/2)/(2π)*1e-12, x=3.14, s=L"$\left[ 2b \left( \frac{1}{M} + \frac{1}{m} \right) \right]^{1/2}$")
text(y=(2b*(1/M_Na))^(1/2)/(2π)*1e-12, x=3.14, s=L"$\left(\frac{2b}{m}\right)^{1/2}$")
text(y=(2b*(1/M_Cl))^(1/2)/(2π)*1e-12, x=3.14, s=L"$\left(\frac{2b}{M}\right)^{1/2}$")
# display
fig.tight_layout()
display(fig)
# -> this plot is simliller to this arcticle:
# https://www.researchgate.net/figure/Phonon-band-structure-of-NaF-NaCl-NaBr-and-NaI_fig2_267048720今回は正の第一ブルリアンゾーンの半分までを計算した結果が以下のようになる。またJuliaのグラフについては過去記事も参考にしてほしい。
また、この結果について先行研究[3]があるので、こちらも確認してみる。
数値の細かいところには差分はあるがΓ-Lでは大きな差分がないことが確認できる。
詳細な計算は第一原理計算などのほうが向いていると考えるが、基礎的な理解としては十分な検討ができたと考えている。
光学フォノン、音響フォノン(光学的振動、音響的振動)
今回の計算で得た結果について分析していく。上の分枝(blanch);\(\omega_+\)を光学的振動(Optical vibration);光学モード(Optical mode);光学フォノン(Optical phonon)といい、下の分枝を;\(\omega_-\)を音響的振動(Acoustical vibration);音響モード(Optical mode);音響フォノン(Acoustic phonon)と呼ぶ。
using PyPlot
# パラメータ
a = 1.0
N = 20
x = a .* (0:N-1)
A = 0.2
ω = 1.0
t = 0
# 振動パターン
function get_displacement(mode::Symbol)
u = zeros(N)
for n in 1:N
phase = ω*t - 0.1π*n
if mode == :TO
u[n] = A * (-1)^n * cos(phase)
elseif mode == :TA
u[n] = A * cos(phase)
end
end
return u
end
model_sin(x, p) = p[1] * sin.(p[2] * x .+ p[3])
function plot_smooth_curve(x_raw, u_raw; color="black", label="")
x_dense = range(minimum(x_raw), stop=maximum(x_raw), length=200)
# 関数する：
p0 = [0.2, 0.1, 0.0] # 初期パラメータの推定値
fit = curve_fit(model_sin, x_raw, u_raw, p0)
p = fit.param
# fittingとplotting
u_dense = model_sin(x_dense, p)
plot(x_dense, u_dense, color=color, linestyle="--", label=label)
end
# 描画関数
function plot_mode(mode::Symbol, subplot_idx)
u = get_displacement(mode)
colors = [n % 2 == 0 ? "blue" : "red" for n in 1:N]
sizes = [n % 2 == 0 ? 300 : 150 for n in 1:N] # matplotlibのサイズはポイント^2
subplot(2, 1, subplot_idx)
scatter(x, u, s=sizes, c=colors)
even_idx = 2:2:N
odd_idx = 1:2:N
plot_smooth_curve(x[even_idx], u[even_idx], color="blue", label="Cl wave")
plot_smooth_curve(x[odd_idx], u[odd_idx], color="red", label="Na wave")
ylim(-0.3, 0.3)
xlim(-0.5, N-0.5)
xlabel("格子位置")
title(mode == :TO ? L"$\omega_+$:TOモード" : L"$\omega_-$:TAモード")
grid(true)
end
# 図を描画
figure(figsize=(6, 4))
plot_mode(:TO, 1)
plot_mode(:TA, 2)
tight_layout()
display(gcf())
NaCl結晶を例に\(\omega_+\)と\(\omega_-\)について横波を描いた場合を上図に示す。
このように、振動に光学的、音響的と呼ぶ理由について説明をしていく。\(q=0\)の場合を考えると、\(\omega_+\)と\(\omega_-\)は次の式に代入する
\begin{equation} -M\omega^2A = bB(1 + e^{iqa})-2bA\\ -m\omega^2B = bA(1 + e^{-iqa})-2bB \end{equation} こののち、(\(A,B\) )~(重い原子の振幅,軽い原子の振幅)について整理をすると \begin{equation} \begin{split} &\omega_+の場合： &\frac{A}{B} = -\frac{m}{M}\\ &\omega_-の場合： &\frac{A}{B} = 1 \end{split} \end{equation} この関係は重要な意味を示唆している。
\(\omega_+\)：光学的振動の場合は重い原子と軽い原子において、逆位相＆逆の重み付けで変位していることがわかる。これが意味するところは重心位置が0変位であるといったことである。また、\(m\)と\(M\)の変位が逆向きであり、正と負の電荷の変位が逆向きで振動していることなる。よって、この振動は分極を起こしており、この時の振動数と同じ電磁波(光)の吸収を起こすようになる。このため、この振動が”光学的”と呼ばれる。
一方、\(\omega_-\)：音響的振動の場合は原子は一体となって変位し、振幅は等しく変位している。これが意味していることは全体で並進あるいは停止していることを示す。この場合は正負が同方向に振動するため分極は発生しない。また、結晶中を音波が伝わるときの振動様式と同じであるため、”音響的”と呼ばれる。
参考文献
[1]青木恭弘. 光ファイバの誘導ラマン、ブリュアン散乱現象. 光学 18巻6号, 303–304 (1989).
[2]坂田亮. 理工学基礎 物性科学. (培風館, 1989).
[3]Messaoudi, I., Zaoui, A. & Ferhat, M. Band-gap and phonon distribution in alkali halides. physica status solidi (b) 252, (2014).
コメント