三角形格子のフォトニック結晶欠陥導波路の固有モード計算

これまで正方格子のフォトニック結晶導波路について固有値問題、伝達問題を計算しましたが、三角形格子ではどうなるかやってみました。

三角形格子の場合、正方格子の解析で用いた周期構造導波路の固有値問題の定式化では文献と比較してうまくいかないケースがありました。

正方格子の固有値解析では、界を

　　Φ(x, y) = φ(x, y) exp(-jβx)

とおき、φに関する固有方程式を解く方法をとりました。この方法は固有値問題をβの２次の固有値問題に帰着できるので、周波数 f を与えて固有値βを求めることができました。

ところで、ブロッホ・フロケの定理で課される境界条件は、

　　Φ(d, y) = Φ(0, y) exp(-jβd)
　　Φ' (d, y) = Φ' (0, y) exp(-jβd)
　　ここで d：周期構造周期(距離)
　　　　　 Φ' = dΦ/dx : Φのxに関する偏微分

なので、この条件を直接適用しΦに関する固有方程式を解く方法もあります。こちらの方が一般的な方法なのかもしれません。この方法は、exp(-jβd)の項が固有方程式に残るため、βを与えて固有値 f (周波数) を求めることになります。

上記２種類の定式化については下記が参考になります。

Φ(x, y) = φ(x, y) exp(-jβx)と置く方法

(1) Yasuhide Tsuji and Masanori Koshiba,

“Finite element method using port truncation by perfectly matched layer boundary conditions for optical waveguide discontinuity problems,”

Journal of Lightwave Technology, vol. 20, p. 463 - 468, March 2002

Φ(x, y) を直接解く方法

(2) 小柴正則, “光・波動のための有限要素法の基礎”, 森北出版,

　　7. 周期構造導波路問題, p.169 - p.186

(3) Naoya Kono, Yasuhide Tsuji

"A Novel Finite-Element Method For Nonreciprocal Magneto-Photonic Crystal Waveguides"

https://kitir.lib.kitami-it.ac.jp/dspace/bitstream/10213/937/1/4226.pdf

Journal of Lightwave Technology, vol. 22, no. 7, July 2004

さて、計算した三角形格子フォトニック欠陥導波路は、正方格子を45°回転させたものです。引用元はこちら。

Attila Mekis, Shanhui Fan, and J. D. Joannopoulos

"Bound states in photonic crystal waveguides and waveguide bends"

http://www.stanford.edu/group/fan/publication/Mekis_PRB_58_4809_1998.pdf

Physical Review B Volume 58, no 8, 15 august 1998-II

引用元のFig.1 (b)に45°の三角形格子のロッドを３列取り除いた導波路が取り上げられています。ロッドの半径r = 0.18a (a:格子定数)、屈折率n = 3.4です。

結論から言うとΦ(x, y) = φ(x, y) exp(-jβx)と置く方法ではうまく解けませんでした。

まず、Φ(x, y) を直接解く方法で計算した結果を示します。

Cad図面

f:id:ryujimiya:20130122110856p:plain

計算結果　(Φ(x, y) を直接解く方法)

分散特性　(Φ(x, y) を直接解く方法)

　βを与えて規格化周波数 a/λ を求めました。以下、横軸は βd /(2π) (d = √2 a: x方向周期構造周期)、縦軸は規格化周波数a/λです。

偶モードの分散特性　(Φ(x, y) を直接解く方法)

f:id:ryujimiya:20130122110915p:plain

奇モードの分散特性　(Φ(x, y) を直接解く方法)

f:id:ryujimiya:20130122110932p:plain

高次偶モードの分散特性　(Φ(x, y) を直接解く方法)

f:id:ryujimiya:20130122110948p:plain

　文献と比較してほぼ同じ結果が得られていると思います。

界分布

　カットオフ付近の界分布を示します。

偶モード　(βd / (2π) = 0.32)　(Φ(x, y) を直接解く方法)

f:id:ryujimiya:20130116071521p:plain

奇モード　(βd / (2π) = 0.0)　(Φ(x, y) を直接解く方法)

f:id:ryujimiya:20130116071640p:plain

高次偶モード　(βd / (2π) = 0.0)　(Φ(x, y) を直接解く方法)

f:id:ryujimiya:20130116071721p:plain

一方、Φ(x, y) = φ(x, y) exp(-jβx)と置いた場合の計算結果は次の通りでした。

計算結果　(Φ(x, y) = φ(x, y) exp(-jβx)と置く方法)

偶モード　(Φ(x, y) = φ(x, y) exp(-jβx)と置く方法)

f:id:ryujimiya:20130122111100p:plain

奇モード　(Φ(x, y) = φ(x, y) exp(-jβx)と置く方法)

f:id:ryujimiya:20130122111114p:plain

高次偶モード　(Φ(x, y) = φ(x, y) exp(-jβx)と置く方法)

f:id:ryujimiya:20130122111137p:plain

文献と比較すると、βd/(2π) >= 0.2以降で差異がみられます。

なぜ一致しなくなるのかよくわかりませんが、Φ(x, y) = φ(x, y) exp(-jβx)と置くことの妥当性を検証する必要がありそうです。

【追記】特に高次偶モードの差異が大きいのですが、これは周期構造方向(x方向)の分割数が足りないからだと思われます。βd/(2π)が0.5の界分布を見ると、Φ(x, y)はcos(πx/d)のような分布です（Φを直接解く方法の分散特性の図を参照）。一方、φ(x, y)はcos(2πx/d)のような分布で（Φ = φ exp( - jβz)と置く方法の分散特性の図を参照)、Φよりφの方がより振動の激しい分布になっています。したがって、φの近似にはΦの近似より細かな分割が必要となることが予想できます。

なお、正方格子誘電体ロッドの場合はどちらの方法でもほぼ同じ結果が得られました。

ではまた。

【追記２】

三角形格子の周期構造領域を変えて計算した結果を示します。

＜180°周期方向にずらした場合＞

f:id:ryujimiya:20130122200810p:plain