はじめに
以前よりメタサーフェイスのシミュレーションに関する記事をあげてきた。この記事内では2次元フォトニクスであるメタサーフェイスのシミュレーション手法の１つ、RCWAの計算を行っていた。これらの計算はモジュールを用いたものであり、自分の手でプログラムを実装したことはなかった。
この理由の一つにRCWAの計算をプログラムに落とし込むのに労力がかかるうえに、すでに充実したモジュールがあるからである。
しかし、実際に自分でプログラムを実装することでより計算の理解が深まると考え、今回は実装を行った。
今回は、フォトニクス結晶入門[迫田2004]といった教科書に沿って解説を行い、一部作者プログラムを参考にしながら(Fortranのためちゃんとは読めなかった…)を実装した。
1次元フォトニクス結晶(誘電体多層膜)について
1次元フォトニクス(誘電体多層膜)について考えてみる。誘電体多層膜は様々な光学素子に用いられ、カメラのレンズや様々な光学素子の反射対策/高反射率達成のために利用されている。
図1. 1次元フォトニクス結晶の例
図1に１次元フォトニクス結晶の模式図を示す。各層の厚みが\(l_a\)と\(l_b\)となるような同じ構造が繰り返しているものとする。ここでその比誘電率を\(\epsilon_a\)と\(\epsilon_b\)とする。このときの透過率は以下のような周波数スペクトルとなる。
図2. 1次元フォトニクス結晶の透過率スペクトル
図2からわかることは光の干渉に起因する周期的な透過率の変動に加えて、透過率が0になる箇所が複数存在する。このような周波数帯をフォトニクスバンドギャップという。フォトニクスバンドギャップを利用した応用は多岐にわたり、誘電体多層膜ミラーやバンドパスフィルタなどに利用される。フォトニクスバンドギャップにおいて光は反射された光の位相はほとんどそろっており、互いに強め合うように反射する。
本記事では、この透過率をJuliaの計算によって行っていく。
透過率の計算法
誘電体多層膜の一層の立式
図3. 誘電体多層膜の透過と反射
図3が示すような場合を考える。入射波の変更方向に\(y\)軸をとることにして、各領域の電場の\(y\)成分は次のように表すことができる。
\begin{align*}
E_{Ay} &= E_i e^{i(k_A x-\omega t)} + E_r e^{-i(k_A x+\omega t)} \\
E_{By} &= E_1 e^{i(k_B x-\omega t)} + E_2 e^{-i(k_A x+\omega t)} \\
E_{Cy} &= E_t e^{i(k_C x-\omega t)}
\end{align*}
ここでexpの乗数が正の場合は進行波、負の場合は逆行波となっている。電場同様に磁場の\(z\)成分は
\begin{align*}
H_{Az} &= \frac{k_A}{\mu_0 \omega} \left( E_i e^{ik_Ax} – E_r e^{-ik_Ax} \right) e^{-i\omega t}\\
H_{Bz} &= \frac{k_B}{\mu_0 \omega} \left( E_1 e^{ik_Ax} – E_2 e^{-ik_Ax} \right) e^{-i\omega t}\\
H_{Cz} &= \frac{k_C}{\mu_0 \omega} E_t e^{i(k_C x-i\omega t)}
\end{align*}
である。また、境界条件(電場及び磁場の接戦成分が各境界で連続)から
\begin{align*}
E_i e^{ik_A x_0} + E_r e^{-ik_A x_0} &= E_1 e^{ik_B x_0} + E_2 e^{-ik_B x_0}\\
E_1 e^{ik_B x_1} + E_2 e^{-ik_B x_1} &= E_t e^{i k_C x_1} \\
k_A (E_i e^{ik_A x_0} – E_r e^{-ik_A x_0}) &= k_B (E_1 e^{ik_B x_0} – E_2 e^{-ik_B x_0})\\
k_B (E_1 e^{ik_B x_1} – E_2 e^{-ik_B x_1}) &= k_C E_t e^{i k_C x_1}
\end{align*}
を得る。ここで第一式と第三式からAとBの境界での電場は
\begin{align*}
\left(
\begin{array}{c}
E_i \\
E_r
\end{array}
\right)
=
\left(
\begin{array}{cc}
\frac{k_A+k_B}{2k_A}e^{-i(k_A-k_B)x_0} & \frac{k_A-k_B}{2k_A}e^{-i(k_A+k_B)x_0} \\
\frac{k_A-k_B}{2k_A}e^{-i(k_A+k_B)x_0} & \frac{k_A+k_B}{2k_A}e^{-i(k_A-k_B)x_0}
\end{array}
\right)
\left(
\begin{array}{c}
E_1 \\
E_2
\end{array}
\right)
\end{align*}
を得る。上式はAとBの境界条件から得られたが、これをそのままBとCの境界条件として読み替えると
\begin{align*}
\left(
\begin{array}{c}
E_1 \\
E_2
\end{array}
\right)
=
\left(
\begin{array}{cc}
\frac{k_B+k_C}{2k_B}e^{-i(k_B-k_C)x_1} & \frac{k_B-k_C}{2k_B}e^{-i(k_B+k_C)x_1} \\
\frac{k_B-k_C}{2k_B}e^{-i(k_B+k_C)x_1} & \frac{k_B+k_C}{2k_B}e^{-i(k_B-k_C)x_1}
\end{array}
\right)
\left(
\begin{array}{c}
E_t \\
0
\end{array}
\right)
\end{align*}
これらの数式から入射電場\(E_i\)、反射電場\(E_r\)、出射電場\(E_t\)の関係性がわかる。
誘電体多層膜の立式
一層の立式を\(N\)層の多層膜に拡張する。任意の層：\(j\)層目から\(j+1\)層目での境界を跨ぐ行列は次のように表せる。
\begin{align*}
M_j
=
\left(
\begin{array}{cc}
\frac{k_j+k_{j+1}}{2k_j}e^{-i(k_j-k_{j+1})x_j} & \frac{k_j-k_{j+1}}{2k_j}e^{-i(k_j+k_{j+1})x_j} \\
\frac{k_j-k_{j+1}}{2k_j}e^{i(k_j+k_{j+1})x_j} & \frac{k_j+k_{j+1}}{2k_j}e^{i(k_j-k_{j+1})x_j}
\end{array}
\right)
\end{align*}
\( 0<j<N \)において
\begin{align*}
\left(
\begin{array}{c}
A_j \\
B_j
\end{array}
\right)
=
M_j
\left(
\begin{array}{c}
A_{j+1} \\
B_{j+1}
\end{array}
\right)
\end{align*}
また、\(j=0\)と\(j=N\)において
\begin{align*}
\left(
\begin{array}{c}
E_i \\
E_r
\end{array}
\right)
=
M_0
\left(
\begin{array}{c}
A_{1} \\
B_{1}
\end{array}
\right)
\end{align*}
\begin{align*}
\left(
\begin{array}{c}
A_N \\
B_N
\end{array}
\right)
=
M_N
\left(
\begin{array}{c}
E_t \\
0
\end{array}
\right)
\end{align*}
である。そこで行列\(M\)を
\begin{align*}
M = M_0 M_1 \cdots M_N
\end{align*}
で定義すると
\begin{align*}
\left(
\begin{array}{c}
E_i \\
E_r
\end{array}
\right)
=
M
\left(
\begin{array}{c}
E_t \\
0
\end{array}
\right)
\end{align*}
が得られる。ここまでが教科書にある説明だが、実際の計算をする上では層内での位相変化を計算する必要がある。実際の計算では層厚\(l_j\)を考慮して
\begin{align*}
M’_j
=
\left(
\begin{array}{cc}
\frac{k_j+k_{j+1}}{2k_j}e^{-i(k_j-k_{j+1})x_j}e^{ik_jl_j} & \frac{k_j-k_{j+1}}{2k_j}e^{-i(k_j+k_{j+1})x_j}e^{ik_jl_j} \\
\frac{k_j-k_{j+1}}{2k_j}e^{i(k_j+k_{j+1})x_j}e^{-ik_jl_j} & \frac{k_j+k_{j+1}}{2k_j}e^{i(k_j-k_{j+1})x_j} e^{-ik_jl_j}
\end{array}
\right)
\end{align*}
で計算を行う。
実際の透過率計算
パラメータ名 シンボル 値 説明 層の総数 N28 シミュレーション対象の層の総数 材料1の厚さ m1_length0.5 周期構造を構成する第一の材料の厚さ 材料2の厚さ m2_length0.5 周期構造を構成する第二の材料の厚さ 欠陥層の厚さ def_length0.5 フォトニック結晶内の欠陥層の厚さ 材料1の誘電率 m1_eps2.1025 第一の材料の比誘電率 材料2の誘電率 m2_eps5.29 第二の材料の比誘電率 欠陥層の誘電率 def_eps2.1025 欠陥層の比誘電率 開始周波数 ω10.05 シミュレーションの開始周波数 終了周波数 ω21.00 シミュレーションの終了周波数 周波数ステップ数 Δω9000 周波数範囲を分割するステップ数
上記のパラメータを用いて計算を行った。第一材料の層数15層と第二材料の層数14層で計算を行っている。また、のちの拡張のために真ん中には欠陥層を入れているが全層数を28層+第一材料の層とすることで計算を行っている。以下にJuliaコードを示す。
using LinearAlgebra
# 伝送係数と反射係数を計算する関数
function calculate_trans_ref_coeff(Q::Matrix{ComplexF64})
QR = -Q[2, 1] / Q[2, 2] # 反射係数
QT = Q[1, 1] + Q[1, 2] * QR # 伝送係数
return QT, QR
end
# ========== 物理パラメータの設定 ==========
const N = 28
const m1_length = 0.5
const m2_length = 0.5
const def_length = 0.5
const m1_eps = 2.1025 + 0.0im
const m2_eps = 5.29 + 0.0im
const def_eps = 2.1025 + 0.0im
const m1_eta = sqrt(m1_eps)
const m2_eta = sqrt(m2_eps)
const def_eta = sqrt(def_eps)
const ω1 = 0.05
const ω2 = 1.00
const Δω = 9000
# ========== 行列の設定 ==========
# 入射境界の転送行列(真空→材料n)
function inc_matrix(ω, eta_n)
QK0 = 2.0 * π * ω
QK1 = QK0 * eta_n
result = [(QK1 + QK0)/(2.0 * QK1) (QK1 - QK0)/(2.0 * QK1)
(QK1 - QK0)/(2.0 * QK1) (QK1 + QK0)/(2.0 * QK1)]
return Matrix{ComplexF64}(result)
end
# 材料nでの位相変化行列
function phase_matrix(ω, eta_n, length_n)
QK1 = 2.0 * π * ω * eta_n
result = [exp(1im * QK1 * length_n) 0.0 + 0.0im
0.0 + 0.0im exp(-1im * QK1 * length_n)]
return Matrix{ComplexF64}(result)
end
# 材料n→材料mの転送行列
function change_material_matrix(ω, eta_n, eta_m, length_n)
QK1 = 2.0 * π * ω * eta_n
QK2 = 2.0 * π * ω * eta_m
Q1 = (QK1 + QK2)/(2.0 * QK1)
Q2 = (QK1 - QK2)/(2.0 * QK1)
Q3 = exp(1im * QK1 * length_n)
Q4 = exp(-1im * QK1 * length_n)
result = [Q1 * Q3 Q2 * Q3
Q2 * Q4 Q1 * Q4]
return Matrix{ComplexF64}(result)
end
# 出射境界の転送行列(材料n→真空)
function out_matrix(ω, eta_n, length_n)
QK0 = 2.0 * π * ω
QK1 = QK0 * eta_n
Q1 = (QK0 + QK1) / (2.0 * QK0)
Q2 = (QK0 - QK1) / (2.0 * QK0)
Q3 = 1im * (QK0 - QK1) * length_n
Q4 = 1im * (QK0 + QK1) * length_n
result = [Q1 * exp(-Q3) Q2 * exp(-Q4)
Q2 * exp(Q4) Q1 * exp(Q3)]
return Matrix{ComplexF64}(result)
end
# 転送行列を計算する関数（メイン計算部分）
function calculate_transfer_matrix(ω)
# 初期化
transfer_matrix = zeros(ComplexF64, 2, 2)
transfer2_matrix = zeros(ComplexF64, 2, 2)
exit_position = N/2.0 + def_length
# ========各行列の準備========
# 入射光
incident_matrix = inc_matrix(ω, m1_eta)
# 材料1での位相変化行列
phase_m1_inv = phase_matrix(ω, m1_eta, m1_length)
# 材料1→材料2境界での転送行列
m1_to_m2_matrix = change_material_matrix(ω, m1_eta, m2_eta, m1_length)
# 材料2→材料1境界での転送行列
m2_to_m1_matrix = change_material_matrix(ω, m2_eta, m1_eta, m2_length)
# 周期構造の単位セル転送行列
unit_cell_matrix = m2_to_m1_matrix * m1_to_m2_matrix
# 材料1での位相変化（出射側）
phase_exit = conj.(phase_matrix(ω, m1_eta, exit_position))
# 出射境界での転送行列
exit_matrix = out_matrix(ω, m1_eta, exit_position)
# 材料2→欠陥層の境界
m2_to_defect_matrix = change_material_matrix(ω, def_eta, m2_eta, def_length)
# 欠陥層→材料2の境界
defect_to_m2_matrix = change_material_matrix(ω, m2_eta, def_eta, m2_length)
# ========転送行列の段階的構築========
# 入射境界から開始
transfer2_matrix = phase_m1_inv * incident_matrix # 入射境界 + 最初の材料1層
transfer_matrix = m2_to_m1_matrix * transfer2_matrix # + 最初の材料2層
# ========== 前半の周期構造を構築 ==========
num_periods_before_defect = div(N-4, 8) # 欠陥層前の周期数
for i in 1:num_periods_before_defect
transfer2_matrix = unit_cell_matrix * transfer_matrix # 周期構造を追加
transfer_matrix = unit_cell_matrix * transfer2_matrix # さらに1周期追加
end
# 欠陥層の処理
transfer2_matrix = m2_to_defect_matrix * transfer_matrix # 欠陥層を追加
transfer_matrix = defect_to_m2_matrix * transfer2_matrix # 欠陥層後の材料2層を追加
# ========== 後半の周期構造 ==========
for i in 1:num_periods_before_defect
transfer2_matrix = unit_cell_matrix * transfer_matrix
transfer_matrix = unit_cell_matrix * transfer2_matrix
end
# 最終的な転送行列の完成
transfer2_matrix = m1_to_m2_matrix * transfer_matrix # 最後の材料1層を追加
transfer_matrix = phase_exit * transfer2_matrix # 位相調整
transfer2_matrix = exit_matrix * transfer_matrix # 出射境界を追加
return transfer2_matrix
end
# メイン計算関数
function main()
println("一次元フォトニック結晶の伝送特性計算")
println("パラメータ:")
println("N = $N, m1_length = $m1_length, m2_length = $m2_length, def_length = $def_length")
println("m1_eps = $m1_eps, m2_eps = $m2_eps, def_eps = $def_eps")
println("ω1 = $ω1, ω2 = $ω2, Δω = $Δω")
println()
# 結果を格納する配列
omega_array = Float64[]
transmission = ComplexF64[]
reflection = ComplexF64[]
# 周波数ループ
for i in 0:Δω
ω = ω1 + i * (ω2 - ω1) / Δω
# 転送行列の計算
Q = calculate_transfer_matrix(ω)
# 伝送係数と反射係数の計算
QT, QR = calculate_trans_ref_coeff(Q)
# 結果の記録
push!(omega_array, ω)
push!(transmission, QT)
push!(reflection, QR)
# 透過率と反射率の計算
T = abs(QT)^2
R = abs(QR)^2
end
return omega_array, transmission, reflection
end
# 実行
println("計算を開始します...")
omega_array, transmission, reflection = main()
# 透過率と反射率の計算
transmittance = abs.(transmission).^2
reflectance = abs.(reflection).^2
println("\n計算完了！")
println("透過率の最大値: $(maximum(transmittance))")
println("透過率の最小値: $(minimum(transmittance))")
println("反射率の最大値: $(maximum(reflectance))")
println("反射率の最小値: $(minimum(reflectance))")
# エネルギー保存の確認
energy_conservation = transmittance + reflectance
println("エネルギー保存 (T+R) の範囲: $(minimum(energy_conservation)) ~ $(maximum(energy_conservation))")
using PyPlot
plt = PyPlot
plt.figure(figsize=(10, 5))
plt.plot(omega_array, transmittance, label="Transmittance", color="blue")
#plt.plot(omega_array, reflectance, label="Reflectance", color="red")
plt.xlabel("Frequency (ω)")
plt.ylabel("Transmittance")
plt.legend()
plt.grid()
plt.ylim(0,1)
plt.xlim(0.1,1.0)
display(plt.gcf())この結果が上図になる。
コメント