RCWA(厳密結合波解析)法のメモ書き

はじめに
RCWA法について
論文の内容について
1. 前提条件
2. 数学的理解
  1. maxwell方程式による立式
  2. 式の解析
参考
1. 関連

はじめに

RCWA(Rigorous Coupled-Wave Analysis)法は回折格子の光学特性の解析法でMoharamらによって提案された。

現在、光学素子の技術のうち光メタサーフェイス技術というのが注目されている。波長オーダー未満の形状を利用することで光の波面を制御が可能となる。これを利用した光メタサーフェイス技術は従来の幾何光学素子では達成できない光学素子の光学特性を得ることができる。しかし、この光学特性のシミュレーションにおいては利用していた幾何光学的解析は不可能なため、波動光学的観点からのシミュレーションが必要となる。そのシミュレーションの一つとしてRCWA法が注目されている。

その中で特に日本語の解説記事[1][2][3]などが多くある論文[4]について、自分なりに勉強したのでそのメモ書きを兼ねて書き残しておこうと思う。
特に参考文献[1]がわかりやすかったので元論文と合わせて一通り読んでほしい。

RCWA法について

RCWA法とは光電磁場解析の一つであり、マクスウェル方程式に対して空間離散フーリエ変換した方程式を解いて光電磁場を計算する。メタサーフェイスの代表的なシミュレーション方法にあるFDTD法(有限差分時間領域法)は時間ドメインなのに対して、周波数ドメインなことが特徴の一つである。また、maxwell方程式はMathieu方程式の一種であり、このMathieu方程式は解析的に解くのが困難である。しかし、周期的な条件(フロケの条件)を満たすとき(Hillの方程式と呼ばれる)に解析的に解くことができるようになり、RCWA法ではこの解法を利用している。そのため、いわゆる全領域において有限要素法にて計算を行っているFDTDは計算が重いのに対して、RCWA法では比較的軽い計算リソースで計算が可能である。一方、前述したとおり周期構造がない場合はRCWA法を適用できないので柔軟な計算はできない。

論文の内容について

前提条件

今回、上記のようなバイナリ型のグレーティングを考える。また、TE偏光の場合を今回考える。

数学的理解

maxwell方程式による立式

光電場は次のように表せる。

$$ E = \exp \{ -j(k・r-\omega t) \} $$

TE波とは、入射電場( $ E_{inc} $:添え字はincidentから)のy軸に平行な成分のみを持つ波であり、空間部分に注目すると

$$ E_{inc,y} = \exp \{ -j(k_0 n_1 (\sin \theta x + \cos \theta z) \} $$

また、波数$ k_0 $は$ k_0 =\frac{2 \pi}{\lambda_0} $ が成り立つ。ここで$ \lambda_0 $は真空中の波長。
入射前の領域$ \rm{I} $と回折格子の媒質内での領域$ \rm{I}\hspace{-1pt}\rm{I} $とすると電場は

\begin{align}
E_{I,y} = E_{inc,y} + \sum_{i} R_i \exp \{ -j(k_{x_i}x – k_{I,z_i}z ) \} \\
E_{II,y} = \sum_{i} T_i \exp \{ -j(k_{x_i}x – k_{II,z_i} (z-d) ) \}
\end{align}

と表せる。ここで領域$ \rm{I} $の左項は入射電場、右項は反射成分。領域$ \rm{II} $は透過成分を表している。
また、$ R,T $は規格化された反射率、透過率を表しており、添え字の$ i $は$ i $次の反射/透過成分を表している。
フロケの定理*(補足1)から
\begin{align}
k_{xi} = k_0 [n_1 \sin \theta – i (\lambda_0 / \Lambda)]
\end{align}
そして、エネルギー保存の法則から波数は保存されるので
\begin{align}
(n_l k_0)^2 &= k_{xi}^2+k_{zi}^2\\
k_{zi}^2 &= k_0^2 [n_l^2-k_{xi}^2/k_0^2]
\end{align}
との関係が成り立つ。よって
\begin{equation}
k_{\mathfrak{L},zi}
\begin{cases}
+k_0 [n_l^2 – (k_{xi}/k_0)^2 ]^{1/2} & \text{if $k_0 n_l > k_{xi}$,} \\
-jk_0 [(k_{xi}/k_0)^2 – n_l ^2]^{1/2} & \text{if $k_0 n_l < k_{xi}$,}
\end{cases}
\end{equation}

ここで、$ \mathfrak{L},l $は領域を表しており$ k_{\mathfrak{L},zi} $は各領域における波数を表している。

$ k_{\mathfrak{L},zi} $の上式の場合は波数が実数となり、電場の位相部分に影響を与えているのがわかる。そのためこの項が回折格子の回折作用を与える。
また、下式は虚数となるので、電場の振幅成分に影響を与える。指数関数的に電場振幅が減少していき、一般的にはこれはエバネッセント成分となる。
$i$の次数が大きければ大きいほど$ k_{\mathfrak{L},zi} $が虚数を示すことがわかり、一般的な回折格子が高次の回折光がなかなか見えないことがわかる。

また、$k_{xi}$について、
$ k_{\mathfrak{L},zi} $と$ k_{xi} $はエネルギー保存の関係から以上の図のような関係となっており、$k_{xi} = k_0 [n_1 \sin \theta – i (\lambda_0 / \Lambda)]$の式は凹凸一周期につき、何個の波が入るかを表している。
($i$が1のときは$k_0(\lambda_0 / \Lambda)=2\pi/\Lambda$となり、上図のように波が一周期入っているような形になる。)

続いて回折格子領域$gr$の電場を考える。まずは、マクスウェル方程式から磁場は電場を使って
\begin{equation}
H= \left(\frac{j}{\omega \mu} \right) \nabla \times E
\end{equation}
となり、$ \mu $は真空の透磁率、$ \omega $は角周波数である。
回折格子領域では、tangential電場(y平面)と磁場(x平面)について、正規化された$ S,U$の振幅を用いて空間フーリエ展開すると
\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}
ここでは添え字の$i$はさきほどとは異なり$ i $次の空間離散フーリエ成分を表している。この式を見てみると$ S,U$ は$z$依存の項で深度方向の電場分布を表すことがわかり、回折格子の透過率、反射率については$ S,U$の項が重要である。また、回折格子の水平方向は$\exp(-jk_{xi}x)$の項で表されている。
この式をより解析するために電場と磁場はマクスウェル方程式を思い出すと
\begin{equation}
\nabla \times E = – \frac{\partial B}{\partial t} = -j\omega H\\
\nabla \times H = i + \frac{\partial D}{\partial t} = j\omega \varepsilon (x) \epsilon_0 E
\end{equation}
となり、$\epsilon_0, \varepsilon_0$は真空中での透磁率、誘電率を表す。$ i = 0$と仮定して(回折格子内で電流が通らないとして)回折格子中での電場、磁場は
\begin{equation}
\frac{\partial E_{gr,y}}{\partial z} = j\omega \mu_0 H_{gr,x}\\
\frac{\partial H_{gr,x}}{\partial z} = j\omega \mu_0 H_{gr,x} \varepsilon (x) \epsilon_0 E_{gr,y} + \frac{\partial H_{gr,z}}{\partial x}
\end{equation}
上式について$S_{yi},U_{xi}$を用いて表すと
\begin{equation}
\frac{\partial S_{yi}(z)}{\partial z} = k_0 U_{xi}
\end{equation}
が得られる。ここで$k_0 = \omega/c = \omega\sqrt{\epsilon_0\mu_0} $を表す。下式について
\begin{align}
\frac{\partial H_{gr,x}}{\partial z}
&= j\omega \mu_0 H_{gr,x} \varepsilon (x) \epsilon_0 E_{gr,y} + \frac{\partial H_{gr,z}}{\partial x} \\
&= j\omega \mu_0 H_{gr,x} \varepsilon (x) \epsilon_0 E_{gr,y} + \frac{\partial}{\partial x}\left(\frac{j}{\omega \mu_0}\frac{\partial E_{gr,y}}{\partial x}\right)\\
&= j\omega \mu_0 H_{gr,x} \varepsilon (x) \epsilon_0 E_{gr,y} + \left(\frac{j}{\omega \mu_0}\frac{\partial^2 E_{gr,y}}{\partial x^2}\right) \\
&= -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}
ここまでの式変形から$S_{yi},U_{xi}$を用いて表す*(補足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}
$z′=k_0z$と定義して、ここまでの$S_{yi},U_{xi}$について$n$次の空間離散フーリエ展開を行い、これを($n×n$)の行列式でまとめると
\begin{equation}
\left(
\begin{array}{c}
\frac{\partial \mathbb{S}_{y}}{\partial z’} \\
\frac{\partial \mathbb{U}_{x}}{\partial z’}
\end{array}
\right)
=
\left(
\begin{array}{c c}
0 & \mathbb{I} \\
\mathbb{A} &0
\end{array}
\right)
\left(
\begin{array}{c}
\mathbb{S}_y \\
\mathbb{U}_x
\end{array}
\right)
\end{equation}
ここで$\mathbb{I}$は単位行列、$\mathbb{A}$は$\mathbb{A}=\mathbb{K}_x^2-\mathbb{E}$で定義される。$\mathbb{K}_x$は対角成分を$k_{xi}/k_0$とする対角行列で$\mathbb{E}$は$\varepsilon_{i-p}$を$ip$成分とする行列である。上の行列式を計算すると
\begin{equation}
\frac{\partial \mathbb{S}_y(z’)}{\partial z’} = \mathbb{I}\mathbb{U}_x, \frac{\partial \mathbb{U}_x(z’)}{\partial z’} = \mathbb{A}\mathbb{S}_y\\
\end{equation}
左式に$\frac{\partial}{\partial z’}$を作用させると
\begin{equation}
\frac{\partial \mathbb{S}_y^2(z’)}{\partial z’^2} = \mathbb{I}\frac{\partial \mathbb{U}_x(z’)}{\partial z’} = \mathbb{A}\mathbb{S}_y\\
\end{equation}
となる。ここで$\mathbb{S}_y$の$z$依存性を$\exp (q_m k_0z)$($m$はフーリエ次数)の形で仮定すると
\begin{equation}
\mathbb{Q}^2 \mathbb{S}_y =\mathbb{A} \mathbb{S}_y
\end{equation}
ここで$\mathbb{Q}$は成分$q_m$を対角成分にもつ対角行列である。上記の関係から行列$\mathbb{A}$の$m$番目の固有値が$g_m$のとき対応関係から$q_m = \pm \sqrt{g_m}$ になる。*(補足3)すなわち、固有値$g_m$に対応する固有ベクトルを$\mathbb{w}_m$とすると、$m$次のフーリエ次数のグレーティング領域を表す電磁場の固有モードは
\begin{equation}
E_m =\sum_i w_{i,m} \exp(\pm q_m k_0 z) \exp(−jk_{xi}x)
\end{equation}
のように得られる。ここで$w_{i,m}$は固有ベクトル$\mathbb{w}_m$の$i$成分を表す。固有モード$E_m$はx軸方向の波数が異なる$w_{i,m}$を重ね合わせた($=\sum_i w_{i,m}$)電場である。また、z軸方向に関しては$\exp(\pm q_m k_0 z)$を見ればわかるように$q_m$が実数ならエバネッセント波、純虚数なら伝搬する波になっている。また、このときの磁場はマクスウェル方程式から
\begin{equation}
H_m =-j \left( \frac{\epsilon_0}{\mu_0}\right)^{1/2} \sum_i w_{i,m} \pm q_m \exp(\pm q_m k_0 z) \exp(−jk_{xi}x)
\end{equation}
と計算できる。$S_{y,i},U_{x,i}$について特定の固有値$q_m$に対してz軸の正負両方向に進む波($\exp[-k_0q_mz],\exp[-k_0q_m(z-d)]$)があること、波数が異なる$w_{i,m}$を重ね合わせる($=\sum_i w_{i,m}$)ことを考えると次のように展開できる。
\begin{equation}
S_{y,i}(z) = \sum_{m=1}^n w_{i,m} (c_m^+\exp[-k_0q_mz]+c_m^-\exp[k_0q_m(z-d)])\\
U_{x,i}(z) = \sum_{m=1}^n v_{i,m} (-c_m^+\exp[-k_0q_mz]+c_m^-\exp[k_0q_m(z-d)])
\end{equation}
今回置いた$c_m^+,c_m^-$はのちに決定する係数で、$v_{i,m}$は$v_{i,m}=w_{i,m}q_m$で定義され、$v_{i,m}$を係数とする行列$\mathbb{V}$を用いて行列形式にすると$\mathbb{V}=\mathbb{WQ}$となる。

ここで再度各電場や磁場についてまとめてみる。
Region I(入射領域)
\begin{align}
E_{I,y} &= E_{inc,y} + \sum_{i} R_i \exp \{ -j(k_{x_i}x – k_{I,z_i}z ) \}\\
E_{inc,y} &= \exp \{ -j(k_0 n_1 (\sin \theta x + \cos \theta z) \}\\
H_{I,x} &= H_{inc,x} – j \sum_{i} R_i \frac{jk_{I,z_i}}{k_0} \exp \{ -j(k_{x_i}x – k_{I,z_i}z ) \}\\
H_{inc,x} &= -j\sqrt{\frac{\epsilon_0}{\mu_0}}(-jn_I\cos\theta)\exp \{ -j(k_0 n_1 (\sin \theta x + \cos \theta z) \}
\end{align}

Region grating(回折格子領域)
\begin{align}
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)\\
S_{y,i}(z) &= \sum_{m=1}^n w_{i,m} (c_m^+\exp[-k_0q_mz]+c_m^-\exp[k_0q_m(z-d)])\\
U_{x,i}(z) &= \sum_{m=1}^n v_{i,m} (-c_m^+\exp[-k_0q_mz]+c_m^-\exp[k_0q_m(z-d)])
\end{align}

Region II(回折格子の媒質内)
\begin{align}
E_{II,y} &= \sum_{i} T_i \exp \{ -j(k_{x_i}x – k_{II,z_i} (z-d) )\}\\
H_{II,x} &= -j \sqrt{\frac{\epsilon_0}{\mu_0}}\sum_{i} T_i \frac{jk_{I,z_i}}{k_0}\exp \{ -j(k_{x_i}x – k_{II,z_i} (z-d) )\}
\end{align}

式の解析

ここで境界条件を考えていく。$\delta_{i0}\equiv\exp \{ -j(k_0 n_1 (\sin \theta x + \cos \theta z) \}/\exp(-jk_{xi}x)$と定義して、region Iとregion gratingの境界部($z=0$)において連続であることを用いると、
\begin{align}
\delta_{i0}+R_i &= \sum_{m=1}^n w_{i,m} (c_m^++c_m^-\exp[-k_0q_md])\\
j\left(n_I\cos\theta \delta_{i0}+\frac{k_{I,zi}}{k_0}R_i\right) &= \sum_{m=1}^n v_{i,m} (c_m^+-c_m^-\exp[-k_0q_md])
\end{align}
ここで対角行列$\mathbb{X},\mathbb{Y}_I,\mathbb{Y}_{II}$(対角成分をそれぞれ$\exp(-k_0q_md),k_{I,zi}/k_0,k_{II,zi}/k_0$;$m$がフーリエ次数、$i$が回折次数)を定義して、行列成分に書き直すと
\begin{equation}
\left(
\begin{array}{c}
\delta_{i0} \\
jn_I\cos\theta \delta_{i0}
\end{array}
\right)
+
\left(
\begin{array}{c}
\mathbb{I} \\
-j\mathbb{Y}_I
\end{array}
\right)\mathbb{R}
=
\left(
\begin{array}{c c}
\mathbb{W} & \mathbb{WX}\\
\mathbb{V} & -\mathbb{VX}
\end{array}
\right)
\left(
\begin{array}{c}
\mathbb{c}^+ \\
\mathbb{c}^-
\end{array}
\right)
\end{equation}
region IIとregion gratingの境界部($z=d$)において連続であることを用いると、
\begin{align}
T_i &= \sum_{m=1}^n w_{i,m} (c_m^+\exp[-k_0q_md]+c_m^-)\\
\frac{jk_{II,zi}}{k_0}T_i &= \sum_{m=1}^n v_{i,m} (c_m^+\exp[-k_0q_md]-c_m^-)
\end{align}
同様に行列形式に直して
\begin{equation}
\left(
\begin{array}{c}
\mathbb{I} \\
j\mathbb{Y}_{II}
\end{array}
\right)\mathbb{T}
=
\left(
\begin{array}{c c}
\mathbb{WX} & \mathbb{W}\\
\mathbb{VX} & \mathbb{-V}
\end{array}
\right)
\left(
\begin{array}{c}
\mathbb{c}^+ \\
\mathbb{c}^-
\end{array}
\right)
\end{equation}
ここで改めて$i=0のとき\delta_{i0}=1,i\neq 0のとき\delta_{i0}=1$である。これは入射波が平面波を考えていることからも自明である。2つの行列式は展開の次数を$n$(最大フーリエ次数)と$4n$個の未知数$\mathbb{R},\mathbb{T},\mathbb{c}^+,\mathbb{c}^-$に関する連立一次方程式となる。ここで$\mathbb{X}$が$\exp(-k_0q_md)$の対角行列であることを思い出すと一つ問題がでてくる。それはフーリエ次数$m$が大きいときは一般的に$q_m$が大きくなること、回折格子領域の厚さ$d$は波長のオーダーと比較して大きいことである。これは$\exp(-k_0q_md)$が0に近づくことを意味していて、一般的な計算機に使われているガウスの消去法を用いて計算しても解の不安定性につながる。そこでは先に$\mathbb{R},\mathbb{T}$を解くことで解決できるとしている。元論文には個々の計算はなかったが[1]の文献には途中式があったので、そちらを参考にすると反射率について$z=0$のときの電場の境界条件を使って
\begin{equation}
\mathbb{R} = \mathbb{Wc}^+ + \mathbb{Wc}^- – \delta_{i0}
\end{equation}
$z=0$のときの磁場の境界条件の式に代入して
\begin{align}
n_I\cos\theta \delta_{i0}-j\mathbb{Y_IR} &= \mathbb{Vc}^+ – \mathbb{Vc}^-\\
n_I\cos\theta \delta_{i0}-j\mathbb{Y_I}(\mathbb{Wc}^+ + \mathbb{Wc}^- – \delta_{i0}) &= \mathbb{Vc}^+ – \mathbb{Vc}^-\\
(\mathbb{V}+j\mathbb{Y_I W})\mathbb{c}^+ + (-\mathbb{VX}+j\mathbb{Y_I WX})\mathbb{c}^- &= n_I\cos\theta \delta_{i0} + j\mathbb{Y_I}\delta_{i0}
\end{align}
透過率においても同様に$z=d$のときの電場と磁場の境界条件を使って
\begin{equation}
(\mathbb{VX}-j\mathbb{Y_{II} W X})\mathbb{c}^+ + (-\mathbb{V}-j\mathbb{Y_{II} W})\mathbb{c}^- = 0
\end{equation}
これらをまとめて行列にすると
\begin{equation}
\left(
\begin{array}{c c}
\mathbb{V}+j\mathbb{Y_I W} & \mathbb{VX}-j\mathbb{Y_{II} W X}\\
-\mathbb{VX}+j\mathbb{Y_I WX} & -\mathbb{V}-j\mathbb{Y_{II} W}
\end{array}
\right)
\left(
\begin{array}{c}
\mathbb{c}^+ \\
\mathbb{c}^-
\end{array}
\right)
=
\left(
\begin{array}{c}
n_I\cos\theta \delta_{i0} + j\mathbb{Y_I}\delta_{i0} \\
0
\end{array}
\right)
\end{equation}
この方程式を見てみると確かに、対角成分に$\mathbb{X}\sim 0$を含まないので安定的に解けることがわかる。よって、この方程式を解くことによって$\mathbb{X},\mathbb{Y}$を得ることができる。最後に絶対回折効率($反射側DE_{ri},透過側DE_{ti}$ )はポインティングベクトルの計算から次のように導かれる。
\begin{align}
DE_{ri} &= R_iR_i^* \Re \left( \frac{k_{I,zi}}{k_0 n_I \cos \theta}\right) \\
DE_{ti} &= T_iT_i^* \Re \left( \frac{k_{II,zi}}{k_0 n_I \cos \theta}\right)
\end{align}