RCWA法のメモ書き-備考-

はじめに

こちらのサイトはRCWAの数値解析説明の補足になります。
全体の流れは元サイトを参考にしてください。

補足1

フロケの定理とは

本論文中、最初にRCWAが提案された論文中[1]においても\(k_{xi}\)の導出は”from the Floquet conditions(theorism)”よりの表記しかなく、他にも具体的な解説があまり多く存在しなかったので追記します。
周期的な構造をもった媒質に対する電磁場のふるまいを表すのによく使われるのがfloquetの定理といい、floquetの定理自体は一般的な数学の定理の一つです。
光電磁場に適用するように特にわかりやすい記事[2][3]があったのでこちらを参考にします。
この定理は行列の常微分方程式の知識が前提として必要であり、その復習も行うと
\begin{equation}
\dot{x}(t) = A(t)x(t)
\end{equation}
を考えます。ここで\(A(t),x(t)\)はそれぞれn×nの正方行列、n次のベクトルを表しています。
このときのn個の基本解をもち、この基本解を\( x_1(t),x_2(t)…x_n(t)\)として並べたn次の正方行列
\begin{equation}
G(t)=( x_1(t),x_2(t)…x_n(t))
\end{equation}
をこの微分方程式の基本行列解と呼びます。
\(x_n(t)\)は解なので
\begin{equation}
\dot{G(t)}=A(t)G(t)
\end{equation}
を満たします。また、当たり前ですが基本解行列は一意に決まらず無数に存在し、その中の\(G_1(t),G_2(t)\)の関係を見てみると任意の\(t_0\)に対して
\begin{equation}
G_2(t) = G_1(t)G_1^{-1}(t_0)G_2(t_0)
\end{equation}
が成り立つかを確認します。
ここで\(G_1^{-1}(t_0)G_2(t_0)\)(\(=C\)としたときのこの\(C\)をモノドロミー行列という。)の時間依存性を確認します。
\begin{equation}
\frac{d(G_1^{-1}(t)G_2(t))}{dt}=0
\end{equation}
を確認できれば時間依存しないことがわかり、\(G_2(t) = G_1(t)G_1^{-1}(t_0)G_2(t_0)\)の関係式が成り立つことを証明できｍうぇう。
逆行列の微分公式
\begin{equation}
\frac{dG^{-1}}{dt}=-G_1^{-1}\dot{G_1}G_1^{-1}
\end{equation}
と定義から任意の\(a\)に対して\(\dot{G_a}=AG_a\)を利用すると
\begin{align}
\frac{d(G_1^{-1}(t)G_2(t))}{dt}
&=\dot{G_1^{-1}}G_2+G_1^{-1}\dot{G_2}\\
&=-G_1^{-1}\dot{G_1}G_1^{-1}G_2+G_1^{-1}\dot{G_2}\\
&=-G_1^{-1}AG_1G_1^{-1}G_2+G_1^{-1}AG_2\\
&=0
\end{align}
よって時間依存しないことが示されます。

これより実際にフロケの定理を確認します。ここからは係数行列に以下のように周期性がある場合を考えと。
\begin{equation}
A(t+T)=A(t)
\end{equation}

proof.
\begin{equation}
\dot{G}(t+T)=A(t+T)G(t+T)=A(t)G(t+T)
\end{equation}が基本解行列の性質より成り立つ、さきほどのモノドロミー行列の確認における式\(G_2(t) = G_1(t)G_1^{-1}(t_0)G_2(t_0)\)に対して\(G_1\to G(t),G_2\to G(t+T),t_0\to 0\)とすれば定理1が成り立つことがわかる。

proof.
\(P(t)=G(t)e^{-tB}\)、ある行列\(B\)について\(G^{-1}(0)G(T)=e^{TB}\)と定義して定理2より
\begin{align}
P(t+T)&=G(t+T)e^{-(t+T)B}\\
&=G(t)[G^{-1}(0)G(T)]e^{-(t+T)B}\\
&=G(t)[e^{TB}]e^{-(t+T)B}\\
&=G(t)e^{-tB}\\
&=P(t)
\end{align}となり、周期性が示される。

電磁場におけるフロケの定理

こちらについてはかなり苦労して”なんとなく”納得はできたので残しておこうと思います。
言い訳となってしまいますが、自分は数学を専門としていないため、高校-大学初年度レベルの数学しか知らず不足があれば教えてください。
まずは\(k_{xi}\)の導出及びこれを含んだ電場表記に対して、文献によっては”floquetの定理”、”floquet-blochの定理”、”floquet-fourieの定理”などと呼ばれています。
この背景には、これらの表記はfloquetの定理の多様な応用性にあります。実際に調べていく中でどの文献の表現も間違っていないと思いました。
さらに、常微分方程式では一般的な理解なのかは知りませんが導出までおこなっている文献が非常に少なかったです。
このように表記ゆれや導出している文献の少なさ、古さ(そもそもfloquetの定理がドイツ語で書かれている)が大きな混乱のもとになっていましたので、自分なりに整理しました。

\(k_{xi}\)を含む数式がとある文献[4](数10以上の文献を読んでまともにrefrenceがあった文献がこれのみ)で導出が行われていました。
この文献内ではTE波を考えていて、回折格子内の\(y\)電場\(E_y\)を\(E_y = X(x)Z(z)\)に変数分離後\(x\)成分(\(X(x)\)について次のように与えられているとしています。
\begin{equation}
X(x) = \sum_{l=\infty}^{\infty}\beta_l \exp[j2\pi(\xi+lb)x]
\end{equation}
ここで変数は\(\xi,l,b\)はそれぞれ\(\sin\theta /\lambda, i,1/\Lambda\)と置き換えることができ、
\begin{equation}
X(x) = \sum_{i=\infty}^{\infty}\beta_i \exp[j2\pi/\lambda_0(n_1 \sin\theta+i(\lambda_0/\Lambda))x]
\end{equation}
と求めたい\(k_{xi}\)の形が見えています。
ここからはこの数式の導出を考えていきます。これを導くのにはmathieu方程式を解くと求まるらしいので教科書[5]をもとにこれを解きます。
mathieu方程式は2階の常微分方程式で以下のように表せる。
\begin{equation}
\frac{\partial^2 V}{\partial x^2}+\frac{\partial^2 V}{\partial y^2} = \frac{1}{c^2}\frac{\partial^2 V}{\partial t^2}
\end{equation}
ここで\(V=u(x,y)\cos(pt+\epsilon)\)と置くことで以下のように書き換えができます。
\begin{equation}
\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2} + \frac{p^2}{c^2}u = 0
\end{equation}
続いて、\(x=h\cosh \xi\cos\eta, y=h\sinh \xi \sin \eta\)と置くことで
\begin{equation}
\frac{\partial^2 u}{\partial \xi^2} + \frac{ \partial^2 u}{ \partial \eta^2 } + \frac{h^2 p^2}{c^2} ( \cosh^2 \xi-\cosh^2 \eta ) = 0
\end{equation}
このように書くことで、
\begin{equation}
u=F(\xi)G(\eta)
\end{equation}
と変数分離できます。\(V\)は電場の進行方向が\(y\)軸に平行だとした(TE波)ときの電場\(E\)とと書き換えることができて、maxwell方程式も上記のように変数分離できることが示される。Mathieu方程式とmaxwell方程式の関係がわかったところで実際にこれを解いていきます。
ここで以降の議論がわかりやすいように、
\begin{equation}
F(\xi)→X(x)
\end{equation}
と書き換える。(ここでの\(\xi\)と\(x\)は直接関係あるわけではなく、筆者が書きやすいようにこうしました。わかりづらくてごめんなさい。)
この数式は天文学の問題として1883年にfloquetによって解析的な調査が行われ、この解析の一環として上述したfloquetの定理が導かれた。\(g(x)，h(x)\)をmathieu方程式(係数が周期\(2\pi\)を持つ任意の一次方程式)の解の基本系とすると
\begin{equation}
X(x)=Ag(x)+Bh(x)
\end{equation}
であり、\(A,B\)は定数です。\(g(x+2\pi),h(x+2\pi)\)もこの方程式の解となるので、以下のような関係となります。
\begin{equation}
g(x+2\pi)=\alpha_1 g(x)+\alpha_2 h(x), h(x+2\pi)=\beta_1 g(x)+\beta_2 h(x)
\end{equation}
\(\alpha_i,\beta_j,(i,j=1,2)\)は定数で、このとき
\begin{equation}
X(x+2\pi)=(A\alpha_1+B\beta_1)g(x)+(A\alpha_2+B\beta_2)h(x)
\end{equation}
この結果\(X(x+2\pi)=kX(x)\)(\(k\)は定数)となり、
\begin{equation}
A\alpha_1+B\beta_1 = kA, A\alpha_2+B\beta_2=kB
\end{equation}
ここで\(A=B=0\)以外の根を持つとき、以下の式が成り立ちます。
\begin{equation}
\left|
\begin{array}{c c}
\alpha_1-k & \beta_1 \\
\alpha_2 & \beta_2-k
\end{array}
\right|=0
\end{equation}
\(k\)をこの方程式のいずれかの根とすると、関数\(X(x)\)は次のような微分方程式の解となるように構成することができます。
\begin{equation}
X(x+2\pi)=kX(x)
\end{equation}
ここで\(k=e^{2\pi\mu}\)となるように\(\mu\)をおいて、\(\phi(x)=e^{-\mu z}X(x)\)と定義すると
\begin{equation}
\phi(x+2\pi)=e^{-\mu(x+2\pi)}X(x+2\pi ) = \phi(x)
\end{equation}
したがって、微分方程式は
\begin{equation}
u = e^{\mu z}\phi(x)
\end{equation}
となるような解をもちます。(\(\phi(x)\)は\(2\pi\)の周期をもつ関数である)
これはfloquetの定理3を表しているのがわかります。
上記の解の形は一般的にはHillの方程式の解の形であり、Hillの方程式についてより考察していきます。
Hillの方程式は以下のような微分方程式の形をとります。
\begin{equation}
\frac{d^2 u}{dz^2} + J(z) u = 0
\end{equation}
最初に置いたように\(V=u(x,y)\cos(pt+\epsilon)\)であり、maxwell方程式を想定する場合、\(V\)は電場\(E\)を表しています。(つまり、\(u\)は電場の空間成分である)
ここで\(J(z)\)は周期\(\pi\)を持つとして、教科書中[5]では以下のような条件を考えていました。
(1)\(z\)は複素数
(2)\(J(z)\)は実軸を含む複素平面上にあり、実軸に平行な領域

ちなみに、上記条件について我々が解こうとするmaxwell方程式に置き換えてみると
(1)maxwell方程式において\(z = x\)であり、\(x\)は空間座標を表し、実数となる。よって、(1)の条件は満たす。
(2)maxwell方程式において\(J(z)\)は元論文[4]で\(\nabla^2 E_y + k^2\epsilon_{r0}(1+a \cos 2\pi b x )E_y = 0\)とされており、\(J(x) = k^2\epsilon_{r0}(1+a \cos 2\pi b x) \)なので\(J(z)\)は周期\(\pi\)を持つ実数である。よって、(2)を満たしている。

このとき、\(J(z)\)は離散フーリエ展開すると\(\theta_0 + 2\sum_{n=1}^{\infty} \theta_n \cos n z \)であり、\(\theta_n = \theta_{-n}\)であるとすると、
\begin{equation}
u = e^{ \mu z} \sum_{n=- \infty }^{\infty } b_n e^{2niz}
\end{equation}
と微分方程式の解の形が\(e^{\mu z}\phi(x)\)となることからも、以上のような解が推測されます。(元論文[4]ではこの数式を指して、引用としていた)
ここまでが元論文[4]、教科書[5]の引用をもとにいろいろ計算をしてきたましたが、最終的には解の形が推測できる/みなされる(=”assume”)とのことで明確に示されたわけではないです。
ただ、自分としては数学的にも上記のような解の形になることが不自然ではないことが分かったうえで、実験事実としても説明できていることから納得ができました。

[1]M. G. Moharam et.al., “Rigorous coupled-wave analysis of planar-grating diffraction,” J. Opt. Soc. Am. 71, 811-818 (1981)

[2]池田 達彦,”フロケ理論と光学応答：２準位量子系を例にして“,物性若手夏の学校テキスト,2758-2159(2024)

[3]解析学B 講義ノート6

[4]C. B. Burckhardt, “Diffraction of a Plane Wave at a Sinusoidally Stratified Dielectric Grating,” J. Opt. Soc. Am. 56, 1502-1508 (1966)

[5]E.T. whittaker et.al., “A cource of modern analysis” 4th ed. (1940)

備考2

\begin{equation}
\frac{\partial U_{xi}(z)}{\partial z} = -k_0 \sum_{p}\varepsilon_{i-p}S_{yp}+\left( \frac{k_{xi}^2}{k_0} \right) S_{yi}
\end{equation}
こちらの式の導出を考えていく。
\begin{equation}
E_{gr,y} = \sum_{i} S_{yi} (z) \exp(-jk_{xi}x)\\
H_{gr,x} = -j\left(\frac{ \epsilon_0}{\mu_0}\right)^{1/2} \sum_i U_{xi} (z) \exp(-jk_{xi}x)
\end{equation}
上記2式を
\begin{align}
\frac{\partial H_{gr,x}}{\partial z}
= -j\left(\frac{\epsilon_0}{\mu_0}\right)^{1/2}\left(-k_0 \varepsilon (x)E_{gr,y}-\frac{1}{k}\frac{\partial^2 E_{gr,y}}{\partial x^2}\right)
\end{align}
こちらに適用することを考える。
回折格子中の誘電率は以下のように与えられる。ここで\(h\)はフーリエ次数を表す。
\begin{equation}
\varepsilon (x) = \sum_h \varepsilon_h \exp \left( j \frac{2\pi h}{\Lambda} \right)
\end{equation}

ここでバイナリの回折格子を考えていて、回折格子の材質部分(凸部分)の屈折率\(n_{rd}\)と非材質部分(凹部分)の屈折率\(n_{gr}\)から
\begin{equation}
\begin{cases}
\varepsilon_h = (n_{rd}^2 – n_{gr}^2)\frac{\sin (\pi h f)}{\pi h} & \rm{if} h\neq 0, \\
\varepsilon_0 = n_{rd}^2 f+n_{gr}^2 (1-f) & \rm{if} h= 0,
\end{cases}
\end{equation}
ここで具体的に\(\varepsilon (x)\)を展開していく
\begin{align}
\varepsilon (x) E_{gr,y}
&= \sum_h \varepsilon_h \exp \left( j \frac{2\pi h}{\Lambda} \right) \sum_{p} S_{yp} (z) \exp(-jk_{xp}x)\\
&= \sum_p \sum_h \varepsilon_h S_{yp} (z) \exp(-j(k_{xp}-\frac{2\pi h}{\Lambda})x)\\
&= \sum_p \sum_h \varepsilon_h S_{yp} (z) \exp(-j(k_{xp}-\frac{2\pi h}{\Lambda})x)\\
&= \sum_p \sum_h \varepsilon_h S_{yp} (z) \exp(-jk_{xp+h}x)\\ & \uparrow k_{xp+h}=k_{xp}-\frac{2\pi h}{\Lambda}\\
&= \sum_p \sum_i \varepsilon_{i-p} S_{yp} (z) \exp(-jk_{xi}x)\\
&\uparrow i=p+h\\
&= \sum_i \exp(-jk_{xi}x) \sum_p \varepsilon_{i-p} S_{yp} (z)
\end{align}
よってこれを代入することで、今回の導出ができる。

[1]M. G. Moharam and T. K. Gaylord, “Diffraction analysis of dielectric surface-relief gratings,” J. Opt. Soc. Am. 72, 1385-1392 (1982)

備考3

個人的にこのあたりの行列が出てきたあたりでわからなくなったので補足しよう思う。
\begin{equation}
\mathbb{Q}^2 \mathbb{S}_y =\mathbb{A} \mathbb{S}_y
\end{equation}
まず、この式のお気持ちは左が求めたい数値で、右辺が既知の行列となっている。そこで\(\mathbb{A}\)を深堀していくと
\begin{align}
\mathbb{A} &= \mathbb{K}_x^2-\mathbb{E}\\
&= \left(
\begin{array}{c c c c c}
(k_{x(-n)}/k_0)^2 & 0 & … & 0 & 0 \\
0 & (k_{x(-n+1)}/k_0)^2 & … & 0 & 0 \\
& & … & &\\
0 & 0 & … & (k_{x(n-1)}/k_0)^2 & 0 \\
0 & 0 & … & 0 & (k_{x(n)}/k_0)^2
\end{array}
\right)\\
&- \left(
\begin{array}{c c c c c}
\varepsilon_{-n-n} & \varepsilon_{-n-(n-1)} & … & \varepsilon_{-n+(n-1)} & \varepsilon_{-n+n} \\
\varepsilon_{(-n+1)-n} & \varepsilon_{(-n+1)-(n-1)} & … & \varepsilon_{(-n+1)+(n-1)} & \varepsilon_{(-n+1)+n} \\
& & … & &\\
\varepsilon_{(n-1)-n} & \varepsilon_{(n-1)-(n-1)} & … & \varepsilon_{(n-1)+(n-1)} & \varepsilon_{(n-1)+n} \\
\varepsilon_{n-n} & \varepsilon_{n-(n-1)} & … & \varepsilon_{n+(n-1)} & \varepsilon_{n+n} \\
\end{array}
\right)
\end{align}
となる。こうやってみるとわりと意味が分かるようになり、\(\mathbb{K}_x^2, \mathbb{E}\)ともに実数のみで構成された行列であることがわかる。
このことからともにエルミネート行列であり、エルミート行列の和である\(\mathbb{A}\)の固有値が実数であることがわかる。