ビーム伝搬法を活用した光導波路デバイスの構造最
適設計
メタデータ

言語: Japanese
出版者:
公開日: 2019-06-25
キーワード (Ja):
キーワード (En):作成者: 井口, 亜希人
メールアドレス:所属:

URL

https://doi.org/10.15118/00009913

博士学位論文

ビーム伝搬法を活用した
光導波路デバイスの構造最適設計

井口

亜希人

目次
第 1 章 序論
1.1 背景 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 本論文の構成 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1
1
3

第 2 章 ビーム伝搬法の定式化
2.1 Maxwell 方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 2 次元有限差分ビーム伝搬法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.1 波動方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.2 近軸式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.3 広角式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.4 Crank-Nicolson 法に基づく定式化 . . . . . . . . . . . . . . . . . . . . . . .
2.2.5 電力の評価方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3 3 次元セミベクトル有限差分ビーム伝搬法 . . . . . . . . . . . . . . . . . . . . . .
2.3.1 波動方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3.2 Crank-Nicolson 法と交互方向陰解法に基づく定式化 . . . . . . . . . . . . .
2.3.3 電力の評価方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.4 3 次元フルベクトル有限差分ビーム伝搬法 . . . . . . . . . . . . . . . . . . . . . .
2.4.1 波動方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.4.2 Crank-Nicolson 法と交互方向陰解法に基づく定式化 . . . . . . . . . . . . .
2.4.3 電力の評価方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.5 2 次元時間領域有限差分ビーム伝搬法 . . . . . . . . . . . . . . . . . . . . . . . . .
2.5.1 波動方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.5.2 狭帯域時間領域ビーム伝搬法 . . . . . . . . . . . . . . . . . . . . . . . . .
2.5.3 Crank-Nicolson 法と交互方向陰解法に基づく定式化 . . . . . . . . . . . . .
2.5.4 電力の評価方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.6 虚軸ビーム伝搬法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5
6
7
7
8
8
9
10
12
12
15
17
19
19
21
24
25
25
25
26
27
28

第 3 章 感度解析に基づく光導波路デバイスの構造最適設計
3.1 感度解析に基づく構造適設計 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 構造表現手法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3 構造単純化のための平滑化フィルタ . . . . . . . . . . . . . . . . . . . . . . . . . .

29
29
32
34

第 4 章 2 次元ビーム伝搬法を用いた光導波路デバイスのトポロジー最適化
4.1 感度解析手法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2 曲がり導波路 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3 パワー分岐導波路 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

35
35
37
42

i

パワー 2 分岐導波路 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
パワー任意分岐導波路 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

42
45

3 次元セミベクトルビーム伝搬法を用いた光導波路デバイスのトポロジー最適化
感度解析手法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
曲がり導波路 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
パワー 2 分岐導波路 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
モード次数変換素子 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

48
48
50
56
59

第 6 章 2 次元時間領域ビーム伝搬法を用いた光導波路デバイスのトポロジー最適化
6.1 感度解析手法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.2 曲がり導波路 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.3 導波路型反射器 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

63
64
65
69

第 7 章 フルベクトルビーム伝搬法を用いたトポロジー最適化
7.1 感度解析手法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.2 偏波回転素子のトポロジー最適化 . . . . . . . . . . . . . . . . . . . . . . . . . . .

72
72
74

第8章

77

4.3.1
4.3.2
第5章
5.1
5.2
5.3
5.4

結論

参考文献

79

研究業績リスト

87

謝

辞

89

付 録 A 境界条件

90付 録 B Stern の差分式

93付 録 C 偏波結合項の差分式

96付 録 D 2 次の高精度差分式 (IFD2)

100

ii

第1章

1.1

序論

背景

近年，ブロードバンドサービスの普及，スマートフォンの台頭，データのリッチ化や高度化に伴い，インターネットの総通信量が増加傾向にある．光通信システムのさらなる高速・大容量化を達成するため，光通信システムの最も基本的な要素である光導波路デバイスの小型化，低損失化，低コスト化が要求されている．計算機の計算処理能力の向上や電磁界シミュレーション技術の発展によって，近年は光導波路デバイスの設計に数値シミュレーションが積極的に導入されており，開発コストの低減に大きく貢献している．
従来，様々な光導波路デバイスは，設計者の電磁気現象に関する深い理解や，これまでに提案されてきた動作原理の知識等に基づいて設計が行われてきた．しかし，計算機の高性能化に伴い，そのような設計者の知識や経験をほとんど要しない，数値シミュレーションを活用した数理的な手法に基づく最適設計手法の利用が広まりつつある．構造の最適化手法は図 1.1 に示すように，寸法最適化，形状最適化，トポロジー最適化の 3 つに大別可能である．寸法最適化や形状最適化では，設計者が特定の構造を指定し，その構造を基礎として構造の寸法・形状最適化を行う．形状最適化は寸法最適化より設計自由度は高いが，構造のトポロジーを変化させることはできない．ト
ポロジー最適化は指定した設計領域内の材料分布そのものを最適化できるため，場合によっては，特別な初期材料分布を与えなくても自動的に最適な構造を発現させることが可能である．先に挙げた 3 つの最適化手法の中で，トポロジー最適化は最も設計自由度が高い最適化手法であり，それゆえデバイスの大幅な性能改善，および新たな動作原理を持つデバイスの創生の可能性を秘めている．トポロジー最適化法は構造力学分野での発展ののち，その応用は熱，音響，電磁波へと様々な分野に展開されている [1][2]．光導波路デバイスに関しては，その適用例は 2004 年のフォトニック結晶導波路における 90 度曲がりの設計例 [3] をはじめとして，主にフォトニック結晶導波

図 1.1: 構造最適化の分類1

路デバイスが数多く報告されており，60 度，90 度，120 度曲がり [4][5][7][11]，Y 分岐 [6]，T 字導
波路 [8]，スローライト [12][14][17]，波長分離器 [19][21]，モード次数変換器 [20] が作製例も含めて報告されている．フォトニック結晶導波路デバイスの他にも，平面光波回路 (Planner Lightwave
Circuit, PLC) デバイスやシリコン細線導波路デバイス等にも応用されており，Y 分岐 [9]，T 字導
波路 [10]，90 度曲がり [10][13]，90 度交差導波路 [15]，パワー 3 分岐 [16]，3 波長分離器 [18] と，
多岐にわたり設計例が報告されている．また，トポロジー最適化と銘打たれてはいないが，トポロ
ジーの変化が生じる自由度の高い設計手法として Objective first 法 [23][22] が提案されており，シ
リコン細線導波路デバイスの設計・作製を通して小型で高性能なデバイスの実現可能性が示されている．
トポロジー最適化では，設計領域内部の材料分布を何らかの設計変数を用いて表現し，数値解析を繰り返しながら所望の特性を最大化あるいは最小化するような設計変数の組を探索する．そのためトポロジー最適化では，数値解析による反復計算が計算時間の大半を占める．従来の光導波路デバイスのトポロジー最適化において，文献 [5]〜[18] では数値解析手法に有限要素法 (Finite Element
Method， FEM)，文献 [19]〜[21] では時間領域有限差分法 (Finite Diﬀerence Time Domain Method，
FDTD 法)，文献 [23] では周波数領域有限差分法 (Finite Diﬀerence Frequency Domain Method, FDFD法) が利用されている．ところで，PLC デバイス等，長手方向に数 mm 程度となるような長いデバイスを解析する場合に対しては，ビーム伝搬法 (Beam Propagation Method， BPM)[24]〜 [33]が効率的な数値解析手法として幅広く利用されている．BPM は高速フーリエ変換 (Fast Fourier
Transformation， FFT) に基づく FFT-BPM[24] を初めとして，横方向を有限差分法 (Finite Diﬀerence
Method， FDM) により離散化する FD-BPM[25]〜 [31] や，有限要素法を用いる FE-BPM[32]，[33]などの各種 BPM が開発されている．BPM は光の伝搬方向に沿って逐次的に電磁界分布を算出す
る数値解析手法であり，伝搬方向に緩慢変化する電磁界を未知変数とするため，後方散乱波を考慮することはできないが FEM や FDTD 法と比較して伝搬方向の離散点数を大幅に低減できる特徴を持つ．周波数領域の解析手法である BPM を時間領域手法に拡張した，時間領域 BPM(TimeDomain BPM, TD-BPM) も既に提案されている [34]〜[38]．この手法は後方散乱波が考慮可能であり，2 次元問題において FDTD 法よりも計算時間が短縮可能なことが報告されている [38]．BPM
を解析手法として用いたトポロジー最適設計例はこれまで報告されておらず，これが実現できればトポロジー最適化の計算コストが低減可能となり，トポロジー最適化の PLC デバイスへの適用
や，FEM や FDTD 法では計算コスト上困難な設計問題への応用が期待できる．
BPM の原理に基づいた自由度の高い最適設計手法として波面整合法 [39]〜 [43] が既に提案されており，PLC デバイスの低損失化に大きく貢献している．波面整合法では，デバイスの出力界分
布および位相が所望のものとなるように，BPM の反復計算過程を経て設計領域内の屈折率分布が自動的に合成される．文献 [39] で提案された初期の波面整合法は，いわゆる「モザイク」状の複雑な屈折率分布が出現し，実験結果との不一致を生んでいた．その後手法の改善が進み，「ソリッ
ド」型の波面整合法が提案された [40]．このアプローチはモザイク状の複雑な材料分布の出現回避に成功しているが，設計自由度は形状最適化程度に落ちている．
構造最適化問題では，所望の特性を目的関数で表現し，その最大化あるいは最小化問題を解く
ことによって目的の特性を最大化する．トポロジー最適化において，設計変数空間の探索手法として MMA(Method of Moving Asymptotes)[44] 等の逐次凸計画法や勾配法が広く用いられている．
ヒューリスティクスとして遺伝的アルゴリズムや粒子群最適化等の多点探索手法があるが，トポ
ロジー最適化においては設計変数の数が通常莫大であるため，現状の計算機性能において，特に3 次元解析が必要な場合は，多点探索手法の適用は計算機コスト上現実的ではない．MMA や勾配2

法に基づいて設計空間を探索する際，目的関数の勾配，すなわち感度計算が必要になるが，この感度を計算する効率的な方法として随伴変数法 (Adjoint Variable Method，AVM) があり，トポロ
ジー最適化において広く適用されている．BPM を用いた AVM に基づく光導波路デバイスの最適化法は既に報告されている [45][46] ものの，これらの報告は寸法最適化を目指しており，トポロ
ジーを目指した報告ではない．
本論文は，光導波路デバイスにおける構造最適化の計算コスト低減を目指して，感度解析に基づくビーム伝搬法を活用した光導波路デバイスのトポロジー最適化手法に関する研究を行った成果をまとめたものである．具体的に本研究では AVM を用いた感度解析に基づき，• 2 次元 BPM を用いたトポロジー最適設計• 偏波結合項を無視する 3 次元セミベクトル (スカラ)BPM を用いたトポロジー最適設計• 3 次元フルベクトル BPM を用いたトポロジー最適設計• 2 次元時間領域 BPM を用いたトポロジー最適設計の実現を新たに行っている．2 次元 BPM を用いたトポロジー最適設計，およびセミベクトル (スカラ)BPM を用いたトポロジー最適設計の実現により，長手方向に数 mm オーダーの石英系 PLC
デバイスに対する効率的なトポロジー最適化が可能となったため，今後さらなる小型化が要求されている PLC デバイスの低損失化が期待される．また，近年，コアとクラッドの比屈折率差が非常に大きな SiN 系あるいは Si 細線導波路デバイスが盛んに報告されているが，これらの導波路はBPM の近似制約上，後方散乱波が大きい場合等，原理的に適用が困難なケースが多い．そのため，
本研究では，後方散乱波が考慮可能な時間領域 BPM を用いたトポロジー最適設計手法の実現も行い，2 次元近似設計問題を通してその有効性を示している．また，フルベクトル BPM を用いたトポロジー最適設計についても検討を行い，比屈折率差の大きな強導波路デバイスの設計に対して本論文で提案するアプローチの利用可能性を検討している．比屈折率差が大きなシリコン細線デバ
イス等のトポロジー最適設計のさらなる効率化を行うにあたり，後方散乱波が考慮可能な双方向
BPM[47]〜[56] の利用などが考えられるが，本研究はそのような今後の研究の発展において先駆的な役割を果たすと思われる．本研究では特に，横方向の離散化に有限差分法を用いる FD-BPMの利用を考え，設計領域の比誘電率分布の表現手法は密度法を適用する場合に対して感度解析手法の定式化を行っている.

1.2

本論文の構成本論文の構成は以下のとおりである．

2章
本研究で利用する FD-BPM 手法の具体的な定式化を行っている．はじめに，等方・線形な誘電体材料で構成される導波路を伝搬する光波が満たす Maxwell 方程式を示している．その後，各節においてそれぞれの問題に対応する波動方程式を導いたのちに 2 次元 BPM，交互方向陰解法 (Alternating Direction Implicit Method, ADI 法) に基づく 3 次元セミベクトル (スカラ)BPM，フルベクトル BPM，2 次元時間領域 BPM の定式化を行っている．また，導波モードの計算手法として用いた虚軸 BPM についても述べている．

3

3章
光導波路デバイスの感度解析に基づくトポロジー最適設計の実装に必要な事項を述べている．はじめに，本研究におけるトポロジー最適化の流れを処理項目ごとに概説する．本研究では，材料分布の表現方法に密度法の適用を考えるが，密度法では実際の作製上問題となる
グレイ材料の出現やチェッカーボードパターンに代表される微細で複雑な構造の出現が問題となる．本章では，本研究におけるグレイ領域の出現の抑制方法および微細構造の除去を目的として導入する構造平滑化フィルタについても述べる．

4章
BPM を用いたトポロジー最適設計の実現のための初期検討として，数値解析手法に 2 次元
BPM を用いる場合に対して，AVM に基づく感度解析手法の定式化を新たに行っている．本
感度解析手法を比屈折率差が小さな導波路デバイスのトポロジー最適設計問題へ適用し，定式化の妥当性および提案アプローチによるトポロジー最適化手法の有効性を示す．具体的には，S 字曲がり導波路のトポロジー最適設計例を通して感度解析手法の妥当性を示し，続いてパワー 2 分岐および任意分岐の設計例を通して本設計手法の汎用性を確かめている．
5章
BPM を活用したトポロジー最適化手法を 3 次元設計問題へ応用するため，数値解析手法に
ADI 法に基づく 3 次元セミベクトル (スカラ)FD-BPM を適用する場合に対して感度解析手法の定式化を新たに行っている．3 章同様，比屈折率差が小さな導波路デバイスのトポロジー
最適設計問題への適用例を通して，本感度解析手法の妥当性および設計手法の有用性を示している．
6章
数値解析手法に ADI 法に基づく 2 次元時間領域 FD-BPM(2DTD-FD-BPM) を用いた場合の感度解析手法の定式化を行っている．比屈折率差が大きな導波路デバイスのトポロジー最適設計問題への適用例を通して感度解析手法の妥当性を示している．また，TD-BPM はガウシア
ンパルスの伝搬解析によって 1 度の解析で透過スペクトルが算出可能であり，本トポロジー
最適化手法によって広帯域なデバイスが効率的に設計できる可能性を示している．

7章
比屈折率差の大きな導波路においてはデバイスの偏波依存性が強くなり，偏波の適切な制御が必要になる．偏波制御を行うデバイスの設計には，偏波の結合を考慮したフルベクトル解
析が必要である．6 章では 2 次元時間領域ビーム伝搬法を活用する強導波路デバイスのトポ
ロジー最適化手法の有効性を示したが，陽解法に基づく FDTD 法に対して，3 次元問題では時間領域ビーム伝搬法の優位性が明らかにされておらず，フルベクトル時間領域ビーム伝搬法の活用が計算効率の向上につながるかは明確ではない．そのため本章では，偏波の結合を考慮したフルベクトル解析が必要なデバイスの設計問題の効率化を目的として，3 次元フルベクトル FD-BPM を用いた場合の感度解析手法の定式化を行っている．偏波回転器のトポ
ロジー最適設計問題への適用例を通して感度解析手法の妥当性を示し，本アプローチの利用可能性を検討している．

8章
本研究で得られた成果および結論と今後の課題を示している．

4

第2章

ビーム伝搬法の定式化

本章では，初めに誘電体導波路中を伝搬する光波が満たす Maxwell 方程式を示し，その後，各節において各問題に対応する波動方程式を導出したのちに有限差分ビーム伝搬法 (Finite Diﬀerence
Beam Propagation Method, FD-BPM) の具体的な定式化を行う．BPM では光波を伝搬方向にキャリ
ア周波数程度で振動する成分と緩やかに振動する成分とに分離し，そのうち緩慢変化する成分を未知変数とする手法である．通常の BPM では非常に広角に伝搬する波や後方散乱波を考慮することはできないが，緩慢変化する電磁界を未知とするため伝搬方向に対する刻み幅を大きくとることが可能であり，その計算効率は有限要素法 (Finite Element Method, FEM) や有限差分時間領域法(Finite Diﬀerence Time Domain, FDTD 法) などと比較すると飛躍的に高い．BPM は高速フーリエ変換 (Fast Fourier Transformation, FFT) を利用する FFT-BPM をはじめとして，横方向の離散化に差分法を適用する FD-BPM や FEM を適用する FE-BPM 等，各種の BPM が開発されている．特にFD-BPM について述べると，広角化 [57][58] や高精度化手法 [59]〜[61] も種々提案されている．
本研究では特に，Crank-Nicolson 法を利用する陰解法型の FD-BPM[25] の利用を考える．このタ
イプの FD-BPM は 2 次元周波数領域問題の場合，バンド幅 3 の行列演算が必要であるが，Thomasのアルゴリズムの利用によって高速に解を得ることができる．3 次元問題や 2 次元時間領域問題ではバンド幅が比較して大きくなるが，交互方向陰解法 (Alternative Direction Implicit, ADI 法)[26][28]の利用によってバンド幅 3 の行列演算に帰着させることができるため，効率的な伝搬解析が可能である．本章では，2 次元 FD-BPM，セミベクトル (スカラ)FD-BPM，フルベクトル FD-BPM お
よび BPM を時間領域手法に拡張した TD-BPM について述べ，最後に，本研究において導波路の固有モードを算出する手法として利用する虚軸 BPM について述べる．また，2 次元 FD-BPM の定式化では Padé 展開による広角化の方法についても述べる．

5

2.1 Maxwell 方程式
導波路を構成する誘電体材料が等方線形であり，磁性をもたない場合を考える．また，本章の定式化では完全整合層 (Perfectly Matching Layer, PML)[62]〜[65] によって計算領域が終端される場合を考える．この場合，Maxwell 方程式は次式のようになる．

∂ H̃
∂t
∂
Ẽ
∇′ × H̃ = ε0 εr∂t∇′ × Ẽ = −µ0

(2.1)
(2.2)

ここで，∇′ は ex ，ey ，ez を直交座標系の単位ベクトルとして

∇′ =

1 ∂
1 ∂
1 ∂
ex +
ey +
ez
sx ∂x
sy ∂y
sz ∂z

である．sα (α = x, y, z) は PML の座標ストレッチングにより生じる係数であり，次式で与える．

(
sα = 1 − j

d

)2

TPML

tan δ

(2.3)

ここで，TPML は PML の厚さ，d は計算領域端からの距離，tan δ は PML 終端 (d = TPML ) における損失角である．その他式中で使用した記号の意味を以下に示す．

Ẽ(x, y, z, t) = E˜x (x, y, z, t)ex + E˜y (x, y, z, t)ey + Ẽz (x, y, z, t)ez : 電界ベクトル [V/m]H̃(x, y, z, t) = H̃x (x, y, z, t)ex + H̃y (x, y, z, t)ey + H̃z (x, y, z, t)ez : 磁界ベクトル [A/m]ε0 : 真空における誘電率 [F/m]µ0 : 真空における透磁率 [A/m2 ]εr : 比誘電率
電磁波が一定の角周波数 ω で振動する場合，電磁界は時間に対して定常振動する成分と振幅成分とにわけてẼ(x, y, z, t) = Re {E(x, y, z) exp(jωt)}(2.4)

H̃(x, y, z, t) = Re {H(x, y, z) exp(jωt)}(2.5)と表現できる．式 (2.4) および (2.5) の表現を式 (2.1)，(2.2) に代入すると，次の周波数領域における Maxwell 方程式を得る．

∇′ × E = −jωµ0 H

(2.6)

′

∇ × H = jωε0 εr E

(2.7)

6

2.2 2 次元有限差分ビーム伝搬法
2.2.1

波動方程式

式 (2.6)，(2.7) より，y 方向に対して電磁界および導波路構造の変化がない場合 (∂/∂y = 0) を仮定すると次式が得られる．

∂Ey= −jωµ0 Hx
∂z
1 ∂Ez
∂Ex
−
= −jωµ0 Hy
∂z
sx ∂x
1 ∂Ey
= −jωµ0 Hz
sx ∂x
∂Hy
−
= jωε0 εr Ex
∂z
∂Hx
1 ∂Hz
−
= jωε0 εr Ey
∂z
sx ∂x
1 ∂Hy
= jωε0 εr Ez
sx ∂x
−

(2.8)
(2.9)
(2.10)(2.11)(2.12)(2.13)

これらの式は 2 つの独立した組に分類でき，式 (2.8)，(2.10)，(2.12) にしたがって伝搬する波はTE(Transverse Electric) 波，式 (2.9)，(2.11)，(2.13) にしたがう波は TM(Transverse Magnetic) 波とそれぞれ呼ばれる．なお，ここでは z 方向を光波の伝搬方向としており，BPM では光波の伝搬方
向に PML を付与する必要がないため，上式において sz = 1 としている．
TE 波の波動方程式は式 (2.8)，(2.10)，(2.12) から次のように導出される．
()∂ 2 Ey1 ∂Ey
1 ∂++ k02 εr Ey = 0(2.14)
sx ∂x sx ∂x∂z 2√
ここで，k0 (= ω µ0 ε0 ) は自由空間波数である．同様に TM 波の波動方程式は式 (2.9)，(2.11)，(2.13)
より次のように求められる．
(
)
(
)
1 ∂Hy
∂
1 ∂Hy
1 ∂
+
+ k02 Hy = 0(2.15)
sx ∂x sx εr ∂x∂z εr ∂z式 (2.14)，(2.15) は統一的に次のように表記できる．
(
)
(
)
1 ∂
p ∂Φ
∂
∂Φ
+
p
+ k02 qΦ = 0sx ∂x sx ∂x∂z∂z(2.16)

ここで，Φ，p，q は表 2.1 のように定義している．
表 2.1: Φ，p，q の定義

TE wave
TM wave

Φ
Ey
Hy

7

p
1
1/εr

q
εr
1

2.2.2近軸式

電磁界 Φ を伝搬方向 (+z 方向) に対して緩慢変化する振幅項とキャリア波数程度で振動する位相項に分離して次のように表現する．

Φ(x, z) = ϕ(x, z) exp(−jk0 n0 z)(2.17)

ここで，n0 は参照屈折率である．式 (2.17) を式 (2.16) に代入して次式を得る．
(
)
(
)
(
)
∂2ϕ
∂ϕ 1 ∂p ∂ϕ
1 ∂
p ∂ϕ
2
2
−
j2k
n
+
−
jk
n
ϕ
+
+
k
ε
−
n
ϕ=0
0
0
0
0
r
0
0
∂z 2
∂z
p ∂z ∂zsx p ∂x sx ∂x(2.18)式 (2.18) 左辺第 2 項は TM 波解析の場合に現れる項である．この項を無視した場合，コアとクラッドの比屈折率差が大きな場合は BPM 解析においてデバイスの相反性および電力保存性が満たされないことが指摘されている [66]〜[69] が，本研究では特に弱導波路の設計問題を扱うため，以降∂p/∂z は微小量として無視する．さらに，伝搬する光波の z 方向位相定数が k0 n0 に十分近い場合，
ϕ は z 方向に対して緩慢に変化するため，∂2ϕ∂ϕ≪ 2k0 n02∂z∂z

(2.19)

が成り立ち，∂ 2 ϕ/∂z 2 ≈ 0 とみなせる．これを近軸近似，あるいは Fresnel 近似と呼ぶ．式 (2.18)より，近軸近似に基づく 2 次元 BPM の基本方程式が次のように導出される．

∂ϕ= (Dxx + ν) ϕ
∂z
(
)
p ∂ϕ
1 ∂
Dxx ϕ =
sx p ∂x sx ∂x()
νϕ = k02 εr − n20 ϕ

j2k0 n02.2.3(2.20)(2.21)(2.22)

広角式

伝搬方向に対する位相定数が k0 n0 から大きく離れる光波 (例えば，z 軸に対して広角伝搬する光波) を扱う場合，2 階の微分項 (∂ 2 ϕ/∂z 2 ) の影響が無視できなくなり近軸式では精度が劣化する．
近軸近似を適用する際に無視した 2 階微分項を考慮する方法として，Padé 展開を用いる方法がある．ここでは，Padé 展開に基づく広角 BPM の基礎方程式を導く．
∂ 2 ϕ/∂z 2 を無視しない場合，BPM の基本方程式は式 (2.18) より次のようになる.

j2k0 n0∂ϕ ∂ 2 ϕ− 2 = (Dxx + ν) ϕ ≡ P ϕ∂z∂z(2.23)式 (2.23) を形式的に次のように表現する．
()∂∂
j2k0 n0 −ϕ = Pϕ
∂z
∂z
∂ϕ
P
ϕ
=
∂
∂z
j2k0 n0 − ∂z

(2.24)(2.25)

8ここで，式 (2.25) 両辺の ϕ に作用している演算子が∂P=∂∂z n j2k0 n0 − ∂z

(2.26)
n−1

∂
=0
∂z −1

(2.27)のような漸化式の関係にあるとみなすと，n = 1 の場合は

∂ϕ
P
=
∂
∂z
j2k0 n0 − ∂z

ϕ
0

∂ϕ
P
=
ϕ
∂z
j2k0 n0 − j2kP n(2.28)

0 0となる. 式 (2.28) 右辺の ϕ に作用する演算子は，分母分子共に P に関する 1 次式になっている. 式(2.28) は Padé(1,1) 展開に基づく BPM の基本方程式である. Fresnel 近似に基づく基本方程式 (2.20)
は n = 0 とした場合であり，Padé(1,0) 展開に基づく BPM の基本方程式に相当する.

2.2.4

Crank-Nicolson 法に基づく定式化

近軸式に基づく BPM 基本方程式 (2.20) を Crank-Nicolson 法に基づく陰的スキーム [25] により差分化して，次式を得る．

j2k0 n0

ϕk+1 − ϕkϕk+1 + ϕk= (Dxx + ν)∆z2

(2.29)

ここで，∆z は z 方向の離散点間隔，k は z 方向に対する離散点番号であり，ϕk は ϕ(x, k∆z) を意
味している．式 (2.29) より，ϕk に関する項を右辺，ϕk+1 に関する項を左辺に移行して整理すると，
BPM の逐次更新式が次のように得られる.
{
}
{
}
∆z
∆z
k+1
1−
(Dxx + ν) ϕ
= 1+
(Dxx + ν) ϕk(2.30)
j4k0 n0
j4k0 n0Padé(1,1) 展開に基づく広角式 (2.28) の場合，逐次更新式は同様の手順により次のように導出できる．
{
(
)
}
∆z
1
1+
−
(Dxx + ν) ϕk+1
4k02 n02 j4k0 n0
{
(
)
}
1
∆z
= 1+
(Dxx + ν) ϕk+(2.31)
4k02 n20 j4k0 n0
近軸式の場合と比較して，下線部の項が追加されている．式 (2.30) および (2.31) は横方向の差分化により，ともに次のように表現できる．
k+1
k
k
k
L1 ϕk+1+ L3 ϕk+1
i−1 + L2 ϕii+1 = R1 ϕi−1 + R2 ϕi + R3 ϕi+1(2.32)

ここで，∆x は x 方向における離散点間隔であり，ϕki は ϕ(i∆x, k∆z) を意味している．本研究では，横方向の差分化には Stern による方法 [70]，[71] を用いる．Stern の方法では，離散点間の中9

心に媒質境界があると仮定し，媒質境界における電磁界の境界条件を考慮して差分式が導出される．具体的な係数 Lα ，Rα (α = 1, 2, 3) の導出方法は付録 A および付録 B に示す.
解析領域終端における境界条件の処理により，式 (2.32) は[A]k {ϕ}k+1 = [B]k {ϕ}k(2.33)の形式で表現できる．ここで，行列 [A]k ，[B]k はサイズが Nx × Nx の 3 重対角行列となる．次ス
テップの界分布 {ϕ}k+1 (= [A]−1k [B]k {ϕ}k ) は行列 [A]k が 3 重対角行列のため，Thomas のアルゴ
リズムにより高速に計算できる．

2.2.5

電力の評価方法

導波モードフィールドには直交性があり，正規直交であることを仮定すると，n 番目の固有モー
ド電界フィールド E n と m 番目の固有モード磁界フィールド H m に対して次が成り立つ.
∫∫
∫∫
(
)
∗
∗
(E n × H ∗m ) · ez dxdy =Ex,n Hy,m− Ey,n Hx,m
dxdy
Ω

Ω

= δmn

(2.34)

ここで，E t,n ，H t,n は固有モードの横方向電磁界，Ω は横方向の計算領域，δnm は Kronecker のデルタであり，上添字 ∗ は複素共役をとることを表している. また，+z 方向に進む光波が運ぶ電力の実効値 P は次式で計算できる．
{∫ ∫
}
1
∗
P = Re(E × H ) · ez dxdy
2
Ω
{∫ ∫
}
(
)
1
∗
∗
= Re
(2.35)Ex Hy − Ey Hx dxdy
2Ω

2 次元問題の場合は電磁界の 1 成分のみで電力評価が可能であり，以下に 2 次元問題における TE
波と TM 波のそれぞれについて，伝搬している電磁界フィールドから任意の導波モードの電力を算出する方法を述べる.
TE 波の場合
光導波路デバイスを伝搬する光波のフィールドは固有モードフィールドの線形結合により，次のように表せる．
∑
Ey =cm Ey,m(2.36)mHx =∑cm Hx,m(2.37)m

ここで，下添字 m は m 番目のモードフィールドであることを意味しており，cm は m 番目のモー∗ を乗じて計算領域の全範囲で積分す
ドフィールドの複素振幅係数である．式 (2.36) の両辺に Hx,n
ると
∫
∑ ∫
∗
∗
Ey,m Hx,ndx
Ey Hx,ndx =cm(2.38)Ω

mΩ

10

となる．TE 波の場合 Ex = 0 であることを考慮すると，モードフィールドの直交性より式 (2.38)
右辺の積分値は m ̸= n の場合零となるため∫∫∗∗
Ey Hx,n dx = cn
(2.39)
Ey,n Hx,ndxΩΩとなる．よって，n 番目のモードの複素振幅係数 cn は次式のように表現できる．
∫∗ dx
Ey Hx,ncn = ∫ Ω∗Ω Ey,n Hx,n dx(2.40)

ここで，位相定数 βn の TE 波が +z 方向に定常伝搬している場合，式 (2.8) より

Hx,n = −

βn
Ey,n
ωµ0

(2.41)

成り立つため，式 (2.41) を式 (2.40) に代入すると，複素振幅係数は Ey 成分のみを用いて次式のように表せる.
∫∗ dx
Ey Ey,ncn = ∫Ω(2.42)
2Ω |Ey,n | dx
式 (2.35) と式 (2.42) より，導波路を伝搬する光波のうち，n 番目の固有モードが +z 方向に運ぶ y
方向に 1m あたりの電力 PTE,n は次式で計算できる．
{∫
}
1
∗ ∗
PTE,n = − Recn Ey,n cn Hx,n dx2Ω∫βn|cn |2|Ey,n |2 dx(2.43)=2ωµ0ΩTM 波の場合
TE 波の場合と同様の議論により，n 番目のモードの複素振幅係数は次式で表現できることが導
ける．
∫ 1
∗ dx
Hy Hy,ncn = ∫Ω ε1r(2.44)
2Ω εr |Hy,n | dx
よって，2 次元問題の場合，TM 波の n 番目のモードの電力は次式で計算できる．
∫
βn
1
2
PTM,n =|cn ||Hy,n |2 dx
2ωε0
Ω εr

11

(2.45)

2.3 3 次元セミベクトル有限差分ビーム伝搬法
2.3.1

波動方程式

電界 (E) を未知とする場合
式 (2.6) の両辺の回転をとると，ベクトル公式∇ × (∇ × A) = ∇ (∇ · A) − ∇2 A

(2.46)

より次を得る．

()()∇′ × ∇′ × E = ∇′ ∇′ · E − ∇′2 E = −jωµ0 ∇′ × H(2.47)式 (2.7) より，式 (2.47) は E のみを用いて()∇′ ∇′ · E − ∇′2 E − ω 2 µ0 ε0 εr E = 0(2.48)と表現できる．さらに，式 (2.7) の発散をとり，ベクトル公式∇ · (∇ × A) = 0(2.49)を用いることによって，∇′ · E = −1 ′∇ εr · E
εr(2.50)の関係が得られるため，式 (2.48) は
(
)
1 ′
′2
′
∇ E+∇
∇ εr · E + k02 εr E = 0
εr(2.51)と書き換えられる．式 (2.51) は電界の 3 成分 (Ex ,Ey ,Ez ) 全てを考慮するフルベクトル波動方程式である．
式 (2.51) を成分ごとに書くと，
(
)
(
)
1 ∂
1 ∂Ex
1 ∂
1 ∂Ey
∂ 2 Ez++sx ∂x sx ∂x
sy ∂y sy ∂y
∂z 2
{ (
)}
1 1 ∂εr
1 ∂
1 ∂εr
∂εr
+
+ k02 εr Ex = 0
Ex +
Ey +
Ez
(2.52)
sx ∂x εr sx ∂x
sy ∂y
∂z
(
)
(
)
∂ 2 Ey
1 ∂Ey
1 ∂
1 ∂Ey
1 ∂
+
+
sx ∂x sx ∂x
sy ∂y sy ∂y
∂z 2
{ (
)}
1 1 ∂εr
1 ∂
1 ∂εr
∂εr
+ k02 εr Ey = 0(2.53)
+
Ex +
Ey +
Ez
sy ∂y εr sx ∂x
sy ∂y
∂z
(
)
(
)
1 ∂Ez
1 ∂
1 ∂Ez
∂ 2 Ez1 ∂++sx ∂x sx ∂x
sy ∂y sy ∂y
∂z 2
{ (
)}
1 1 ∂εr
∂
1 ∂εr
∂εr
+
+ k02 εr Ez = 0
Ex +
Ey +
Ez
(2.54)∂z εr sx ∂x
sy ∂y
∂z

12

となる．ここで，2 次元問題の場合と同様に sz = 1 としている．導波路構造が z 方向に緩やかに変化している (∂εr /∂z ≈ 0) と仮定し，さらに，偏波の結合が生じず片方の偏波成分が無視できる
場合，すなわち，式 (2.52) は Ey ，式 (2.53) は Ex を無視できる場合，式 (2.52) と式 (2.53) はそれ
ぞれ
(
)
(
)
(
)
1 ∂εr
1 ∂
1 ∂Ex
1 ∂
1 ∂Ex
∂ 2 Ex1 ∂+Ex + k02 εr Ex = 0 (2.55)++
sx ∂x sx ∂x
sy ∂y sy ∂y∂z 2
sx ∂x sx εr ∂x
(
)
(
)
(
)
1 ∂εr
∂ 2 Ey
1 ∂Ey
1 ∂
1 ∂Ey
1 ∂
1 ∂
+
+
Ey + k02 εr Ey = 0 (2.56)+sx ∂x sx ∂x
sy ∂y sy ∂y∂z 2
sy ∂y sy εr ∂yと一つの偏波成分のみの波動方程式に変形できる．これらはそれぞれ，Ex あるいは Ey に関する
セミベクトル (スカラ) 波動方程式である．
以降の表記上の便宜のため，次の関係
{
}
(
)
(
)
1 ∂ (εr Eα )
1 ∂
1 ∂Eα
1 ∂
1 ∂εr
1 ∂
=
+
Eα
[α = x, y] (2.57)
sα ∂α sα εr∂α
sα ∂α sα ∂α
sα ∂α sα εr ∂αを用いて，式 (2.55)，(2.56) を次のように表現する．
{
}
(
)
1 ∂
1 ∂ (εr Ex )
1 ∂
1 ∂Ex
∂ 2 Ex+++ k02 εr Ex = 0sx ∂x sx εr∂x
sy ∂y sy ∂y
∂z 2
(
)
{
}
∂ 2 Ey1 ∂Ey
1 ∂1 ∂ (εr Ey )1 ∂+++ k02 εr Ey = 0sx ∂x sx ∂x
sy ∂y sy εr∂y∂z 2(2.58)(2.59)

磁界 (H) を未知とする場合
式 (2.7) の両辺の回転をとり，ベクトル公式 (2.46) を用いて次のように式を変形する．
()()∇′ × ∇′ × H = ∇′ ∇′ · H − ∇′2 H = jωε0 ∇′ × (εr E)(2.60)式 (2.6) の発散をとって，ベクトル公式 (2.49) を適用することにより∇·H =0(2.61)の関係が得られるため，式 (2.60) は次のように書き換えられる．

∇′2 H + jωε0 ∇′ × (εr E) = 0(2.62)

さらに，式 (2.62) の左辺第 2 項に次のベクトル公式∇ × (ϕA) = ∇ϕ × A + ϕ∇ × A(2.63)

を適用して，次を得る．

∇′2 H + jωε0 εr ∇′ × E + jωε0 ∇′ εr × E = 0(2.64)式 (2.64) の左辺第 2 項に式 (2.6)，第 3 項に式 (2.7) を適用すると，次の磁界に関するフルベクトル
波動方程式が得られる．

∇′2 H +()
1 ′∇ εr × ∇′ × H + k02 εr H = 0εr
13(2.65)

式 (2.65) を成分ごとに書くと
(
)
(
)
1 ∂
1 ∂Hx
1 ∂
1 ∂Hx
∂ 2 Hx++sx ∂x sx ∂x
sy ∂y sy ∂y
∂z 2
{
(
)
(
)}
1 1 ∂εr 1 ∂Hy
1 ∂Hx∂εr ∂Hx
1 ∂Hz
+
−
−
−
+ k02 εr Hx = 0
εr sy ∂y sx ∂x
sy ∂y
∂z
∂z
sx ∂x
(2.66)
(
)
(
)
2
∂ Hy
1 ∂
1 ∂Hy
1 ∂
1 ∂Hy
+
+
sx ∂x sx ∂x
sy ∂y sy ∂y
∂z 2
{
(
)
(
)}
1 ∂εr 1 ∂Hz∂Hy
1 ∂εr 1 ∂Hy
1 ∂Hx
+
−
−
−
+ k02 εr Hy = 0εr ∂z sy ∂y∂z
sx ∂x sx ∂x
sy ∂y(2.67)
(
)
(
)
1 ∂
1 ∂Hz
1 ∂
1 ∂Hz
∂ 2 Hz++sx ∂x sx ∂x
sy ∂y sy ∂y
∂z 2
{
(
)
(
)}
1 1 ∂εr ∂Hx
∂Hy
1 ∂Hz
1 ∂εr 1 ∂Hz+−−−+ k02 εr Hz = 0
εr sx ∂x∂z
sx ∂x
sy ∂y sy ∂y∂z(2.68)となる．導波路構造が z 方向に緩やかに変化している (∂εr /∂z ≈ 0) と仮定し，片方の偏波成分が無視できる場合，すなわち，式 (2.66) は Hy ，式 (2.67) は Hx を無視できる場合，式 (2.66) と式 (2.67)
はそれぞれ
(
)
(
)
1 ∂εr ∂Hx
1 ∂
1 ∂Hx
1 ∂
1 ∂Hx
∂ 2 Hx(2.69)++− 2+ k02 εr Hx = 02
sx ∂x sx ∂x
sy ∂y sy ∂y∂z
sy εr ∂y ∂y()()
1 ∂εr ∂Hy
∂ 2 Hy
1 ∂
1 ∂Hy
1 ∂
1 ∂Hy
+
+
−
+ k02 εr Hy = 0(2.70)
sx ∂x sx ∂x
sy ∂y sy ∂y∂z 2
s2x εr ∂x ∂xと一つの偏波成分のみのセミベクトル (スカラ) 波動方程式に変形できる．
表記の便宜上，次の関係
(
)
(
)
εr ∂
1 ∂Hβ
1 ∂
1 ∂Hβ
1 ∂εr ∂Hβ=− 2[(α, β) = (x, y), (y, x)]sα ∂α sα εr ∂α
sα ∂α sα ∂α
sα εr ∂α ∂αを用いて式 (2.69)，(2.70) を次のように変形しておく．
(
)
(
)
1 ∂Hx
εr ∂
1 ∂Hx
∂ 2 Hx1 ∂+++ k02 εr Hx = 0sx ∂x sx ∂x
sy ∂y sy εr ∂y
∂z 2
(
)
(
)
∂ 2 Hy
1 ∂Hy
1 ∂
1 ∂Hy
εr ∂
+
+
+ k02 εr Hy = 0sx ∂x sx εr ∂xsy ∂y sy ∂y∂z 2(2.71)(2.72)(2.73)

統一表現
上記で導出したセミベクトル波動方程式 (2.58)，(2.59)，(2.72)，(2.73) を統一的に次のように表
現する．

∂2Φ
+ k02 εr Φ = 0∂z 2
Φ，Dxx ，Dyy の定義を表 2.2 に示す．
Dxx Φ + Dyy Φ +(2.74)

14

表 2.2: Φ，Dxx ，Dyy の定義

Φ
Ex

2.3.2

1 ∂
sx ∂x

{

Dxx Φ
1 ∂(εr Ex )sx ε r
∂x

Ey

1 ∂
sx ∂x

Hx

1 ∂
sx ∂x

Hy

εr ∂
sx ∂x

(

(

(

1 ∂Ey
sx ∂x
1 ∂Hx
sx ∂x

)

1 ∂Ex
sy ∂y

1 ∂
sy ∂y
1 ∂
sy ∂y

)

1 ∂Hy
sx εr ∂x

Dyy Φ
(

}

)

{

εr ∂
sy ∂y

1 ∂(εr Ey )
sy εr
∂y

(

1 ∂
sy ∂y

)

1 ∂Hx
sy εr ∂y

(

1 ∂Hy
sy ∂y

}

)

)

Crank-Nicolson 法と交互方向陰解法に基づく定式化

電磁界 Φ を伝搬方向 (+z 方向) に対して緩慢変化する振幅項とキャリア波数程度で振動する位相項に分離して次のように表現する．

Φ(x, y, z) = ϕ(x, y, z) exp(−jk0 n0 z)(2.75)式 (2.75) を式 (2.74) に代入すると

Dxx ϕ + Dyy ϕ +

(
)
∂ϕ
∂2ϕ
− j2k0 n0+ k02 εr − n20 ϕ = 02∂z∂z(2.76)となる．近軸近似 (∂ 2 ϕ/∂z 2 ≈ 0) が成り立つ場合，次の近軸近似に基づくセミベクトル BPM の基
本方程式が得られる．

j2k0 n0∂ϕ= (Dxx + Dyy + ν) ϕ∂z()
νϕ = k02 εr − n20 ϕ(2.77)(2.78)

Crank-Nicolson 法を用いて伝搬方向 (z 方向) を差分化すると，式 (2.77) は
j2k0 n0

ϕk+1 − ϕk∆z= (Dxx + Dyy + ν)

ϕk+1 + ϕk
2(2.79)となる．ここで，ϕk は ϕ(x, y, k∆z) を意味している．式 (2.79) を整理すると，次のセミベクトル
BPM の逐次更新式を得る．
{
}
{
}
∆z
∆z
1−
(Dxx + Dyy + ν) ϕk+1 = 1 +(Dxx + Dyy + ν) ϕk(2.80)
j4k0 n0
j4k0 n0
式 (2.80) を用いて BPM 解析を行う場合，横方向の計算領域分割数を Nx = Ny = N とすると，伝搬ステップ毎にバンド幅が 2N + 1 の連立一次方程式を解く必要がある．交互方向陰的差分法(Alternating Direction Implicit Method, ADIM) を適用すると，バンド幅 3 の連立一次方程式を 2N回解くことによって次ステップの界分布が求められる．バンド幅 3 の連立一次方程式は Thomas のアルゴリズムで高速に解くことができるため，離散点数が多いほど ADIM に基づく BPM の方が計算コストの観点から有利になる．

15

ADIM に基づく BPM の逐次更新式を導出する．式 (2.80) に対して，
}
(
{
)2
(
∆z
∆z
ν)(
ν ) k+1
1−(Dxx + Dyy + ν) ϕk+1 +
Dxx +
Dyy +
ϕ
j4k0 n0
j4k0 n0
2
2
{
}
(
)2
(
∆z
∆z
ν)(
ν) k
k
= 1+
(Dxx + Dyy + ν) ϕ +
Dxx +
Dyy +
ϕ
j4k0 n0
j4k0 n022(2.81)のように，下線で示した (∆z)2 オーダーの誤差項を付加する．誤差項を付加したことによって，逐
次更新式は
{
}{
}
∆z (
∆z (
ν)
ν)
1−
Dxx +
Dyy +
1−
ϕk+1
j4k0 n02
j4k0 n0
2
}{
}
{
∆z (
∆z (
ν)
ν)
= 1+
Dxx +
1+
Dyy +
ϕk
(2.82)
j4k0 n02
j4k0 n021のように，演算子を分割した形で表現できる．中間界 ϕk+ 2 を導入して，式 (2.82) を次の 2 ステップの演算に分割する．

1st step{1−
2nd step
{
1−

(

∆z
j4k0 n0∆z
j4k0 n0

}
}
{
(
)
∆z
1
ν)
ν
Dxx +
Dyy +
ϕk+ 2 = 1 +ϕk2
j4k0 n02((2.83)

}
{
}
∆z (
1
ν)
ν)
k+1
Dyy +
ϕ
= 1+
Dxx +
ϕk+ 2
2
j4k0 n02(2.84)

横方向を Stern の方法で離散化すると，式 (2.83)，(2.84) はそれぞれ次のように表現できる．

1st step
k+ 1

k+ 1

k+ 1

1st
1st
1st k
1st k
1st k
2
2
2
L1st
1 ϕi−1,j + L2 ϕi,j + L3 ϕi+1,j = R1 ϕi,j−1 + R2 ϕi,j + R3 ϕi,j+1
nd2(2.85)

step
k+ 1

k+ 1

k+ 1

k+1
2nd k+1
2nd k+1
2nd
2nd
2nd
2
2
2
L2nd
1 ϕi,j−1 + L2 ϕi,j + L3 ϕi,j+1 = R1 ϕi−1,j + R2 ϕi,j + R3 ϕi+1,j(2.86)

ここで，ϕki,j は ∆x，∆y をそれぞれ x，y 方向における離散点間隔として，ϕ(i∆x, j∆y, k∆z) を意味している．式 (2.85)，(2.86) を行列・ベクトル形式で次のように表現する.

1st step[Γ2 ]k {ϕ}k+ 1 = [Γ1 ]k {ϕ}k(2.87)

2

2nd step[Γ4 ]k {ϕ}k+1 = [Γ3 ]k {ϕ}k+ 1(2.88)2[Γα ](α = 1, 2, 3, 4) は計算領域端の境界条件の処理によって Nx Ny × Nx Ny の正方行列となるが，ADIM の適用により式 (2.87) では Ny 回，式 (2.88) では Nx 回の 3 重対角行列の処理により高速に
次ステップの界が求められる.
16

2.3.3電力の評価方法

電磁界の横方向成分 E t = {Ex , Ey }T ，H t = {Hx , Hy }T は固有モードフィールドの線形結合∑
Et =cm E t,m(2.89)mHt =∑cm H t,m(2.90)

mで表現できる. ここで，Φt,m (Φt = E t , H t ) は m 番目のモードフィールド，cm はその複素振幅係()∗ − E H∗数である. Ex Hy,ny x,n を横方向の計算領域 Ω にわたって積分すると，モードフィールドの直交性より次が成り立つ.
∫∫∑ ∫∫ (
(
)
)
∗
∗
∗
∗
Ex Hy,n− Ey Hx,ndxdy =
cm
Ex,m Hy,n− Ey,m Hx,n
dxdy
Ω

m

∫∫

= cn

Ω

(

)
∗
∗
Ex,n Hy,n− Ey,n Hx,ndxdy

(2.91)Ω

したがって，n 番目のモードの複素振幅係数 cn は次式で計算できる.
)∫∫ (∗ − E H∗Ex Hy,ny x,n dxdyΩ)cn = ∫∫ (∗∗Ω Ex,n Hy,n − Ey,n Hx,n dxdy()
∗ − H E∗
また， Hy Ex,n
x y,n を考えると，同様に次が導ける．
)∫∫ (∗ − H E∗
dxdy
H
E
x
y
y,n
x,n
)
cn = ∫∫ Ω(∗∗Ω Hy,n Ex,n − Hx,n Ey,n dxdy

(2.92)(2.93)

電界を未知とする場合は式 (2.92)，磁界の場合は (2.93) を用いる．セミベクトル解析の場合，ある近似の下では電磁界の 1 成分のみで電力の評価ができる．以下に TE-like(Ex ，Hy ) 波およびTM-like(Ey ，Hx ) 波の場合の電力の評価方法を示す.

TE-like 波の場合
電界のみで電力を評価する方法を考える. 式 (2.6) より，+z 方向に位相定数 βn で伝搬するモードに対して

Hy,n =

1 ∂Ez,n
βn
Ex,n +
ωµ0
jωµ0 ∂x(2.94)

が成り立つ. Ex モードでは Ey = 0 であることを考慮して，式 (2.94) を式 (2.92) に代入すると∗ )∫∫ ( βn∗ − 1 E ∂Ez,n dxdy
E
E
x
x
x,n
jωµ0
∂x
Ω ωµ0
)
cn,Ex = ∫∫ ((2.95)
∗
∂E
2
βn
z,n
1
|E
|
−
E
dxdy
x,n
x,n
jωµ0
∂x
Ω ωµ0
を得る. ここで，|∂Ez,n /∂x| は 2 次の微小量であり，無視できる場合を考えると，複素振幅は電界のみを用いて次のように計算できる.
∫∫ ∞∗ dxdyEx Ex,ncn,Ex = ∫∫−∞
(2.96)
∞
2
|E
|
dxdy
x,n
−∞

17

電力は式 (2.35) より，次式で計算できる．
∫∫
βn
2|Ex,n |2 dxdy
Pn,Ex =
|cn,Ex |2ωµ0Ω(2.97)

次に，磁界のみで電力を評価する場合を考える. 式 (2.7) より，+z 方向に位相定数 βn で伝搬するモードに対して

Ex,n =

βn1 ∂Hz,nHy,n +
ωε0 εrjωε0 εr ∂y(2.98)

が成り立つ. Hy モードでは Hx = 0 であることを考慮して，式 (2.98) を式 (2.93) に代入すると∗ )∫∫ ( βn
∂Hz,n
1
∗ −
H
H
H
dxdy
y
y
y,n
jωε0 εr∂yΩ ωε0 εr)cn,Hy = ∫∫ ((2.99)
∗
∂Hz,n
2
βn
1
|H
|
−
H
dxdy
y,n
y,n
jωµ0 εr∂yΩ ωε0 εrを得る. ここで，|∂Hz,n /∂y| は 2 次の微小量であり，無視できる場合を考えると，複素振幅は磁界のみを用いて次のように計算できる.
∫∫ 1∗ dxdyHy Hy,ncn,Hy = ∫∫Ω ε1r(2.100)
2
|H
|
dxdy
y,n
Ω εr
電力は式 (2.35) より，次式で計算できる．
∫∫
βn
1
2
Pn,Hy =
cn,Hy|Hy,n |2 dxdy
2ωε0
ε
r
Ω

(2.101)TM-like 波の場合
TE-like 波の場合と同様の手順により，Ey モードに対して次式が導ける．
∫∫∗ dxdyEy Ey,ncn,Ey = ∫∫Ω2Ω |Ey,n | dxdy

Pn,Ey =

βn
2
cn,Ey
2ωµ0

∫∫
|Ey,n |2 dxdy

(2.103)Ω

同様に Hx モードに対して，∫∫ 1∗ dxdyHx Hx,ncn,Hx = ∫∫Ω ε1r
2Ω εr |Hx,n | dxdy

βnPn,Hx =
|cn,Hx |22ωε0(2.102)∫∫(2.104)

1|Hx,n |2 dxdy
Ω εr

が導ける.

18

(2.105)

2.4 3 次元フルベクトル有限差分ビーム伝搬法
2.4.1

波動方程式

電界 (E) を未知とする場合
電界に関するフルベクトル波動方程式 (2.51) より，∂εr /∂z ≈ 0 を仮定すると次が成り立つ.

)
1 ∂Ex∂ 2 Ex
+
sy ∂y
∂z 2
{ (
)}
1 1 ∂εr
1 ∂εr
1 ∂
+
Ex +
Ey
+ k02 εr Ex = 0sx ∂x εr sx ∂x
sy ∂y
(
)
(
)
∂ 2 Ey
1 ∂
1 ∂Ey
1 ∂
1 ∂Ey
+
+
sx ∂x sx ∂x
sy ∂y sy ∂y
∂z 2
{ (
)}
1 1 ∂εr
1 ∂εr
1 ∂
Ex +
Ey
+ k02 εr Ey = 0+sy ∂y εr sx ∂x
sy ∂y

1 ∂
sx ∂x

(

1 ∂Ex
sx ∂x

)

1 ∂
+
sy ∂y

(

式 (2.106)，(2.107) は式 (2.57) を用いると次のように表現できる.
{
}
(
)
1 ∂
1 ∂ (εr Ex )
1 ∂
1 ∂Ex
∂ 2 Ex++sx ∂x sx εr∂x
sy ∂y sy ∂y
∂z 2
(
)
1 ∂εr
1 ∂
Ey + k02 εr Ex = 0+
sx sy ∂x εr ∂y(){}∂ 2 Ey1 ∂Ey
1 ∂1 ∂ (εr Ey )1 ∂++sx ∂x sx ∂x
sy ∂y sy εr
∂y
∂z 2
(
)
1 ∂εr
1 ∂
Ex + k02 εr Ey = 0+
sx sy ∂y εr ∂xさらに，式 (2.108)，(2.109) に次の関係
{
}
(
)
∂
1 ∂ (εr Eβ ) ∂Eβ
∂
1 ∂εr
−
=
Eβ
∂α εr
∂β
∂β
∂α εr ∂β(2.106)(2.107)(2.108)(2.109)[(α, β) ∈ {(x, y), (y, x)}](2.110)を適用して，最終的に次の表現を得る.
{
}
(
)
1 ∂
1 ∂ (εr Ex )
1 ∂
1 ∂Ex
∂ 2 Ex++sx ∂x sx εr∂x
sy ∂y sy ∂y
∂z 2
{
}
1 ∂
1 ∂ (εr Ey ) ∂Ey++ k02 εr Ex = 0−
sx sy ∂x εr
∂y
∂y
(
)
{
}
∂ 2 Ey
1 ∂
1 ∂Ey
1 ∂
1 ∂ (εr Ey )++sx ∂x sx ∂x
sy ∂y sy εr
∂y
∂z 2
{
}
1 ∂ (εr Ex ) ∂Ex
1 ∂+ k02 εr Ey = 0+−
sx sy ∂y εr∂x∂x(2.111)19(2.112)

磁界 (H) を未知とする場合
磁界に関するフルベクトル波動方程式 (2.65) より，∂εr /∂z ≈ 0 を仮定すると次が成り立つ.

(

)

(

)

∂ 2 Hx∂z 2()
1 ∂εr 1 ∂Hy
1 ∂Hx
+
−
+ k02 εr Hx = 0 (2.113)
sy εr ∂y sx ∂x
sy ∂y
(
)
(
)
∂ 2 Hy
1 ∂
1 ∂Hy
1 ∂
1 ∂Hy
+
+
sx ∂x sx ∂x
sy ∂y sy ∂y
∂z 2
(
)
1 ∂Hx
1 ∂εr 1 ∂Hy−+ k02 εr Hy = 0 (2.114)−
sx εr ∂x sx ∂x
sy ∂y

1 ∂
sx ∂x

1 ∂Hx
sx ∂x

+

1 ∂
sy ∂y

1 ∂Hx
sy ∂y

+

式 (2.113)，(2.114) は式 (2.71) を用いて次のように表現できる．
(
)
(
)
1 ∂Hx
εr ∂
1 ∂Hx
∂ 2 Hx1 ∂1 ∂εr ∂Hy+++ k02 εr Hx = 0+2
sx ∂x sx ∂x
sy ∂y sy εr ∂y∂zsx sy εr ∂y ∂x()()∂ 2 Hy
1 ∂Hy
1 ∂Hy
εr ∂
1 ∂
1 ∂εr ∂Hx+++
+ k02 εr Hy = 0sx ∂x sx εr ∂xsy ∂y sy ∂y∂z 2
sx sy εr ∂x ∂yさらに，式 (2.115)，(2.116) に次の関係
(
)
∂Hα
∂
1 ∂Hα
1 ∂εr ∂Hα
− εr
=
∂α∂β
∂α εr ∂βεr ∂α ∂β[(α, β) ∈ {(x, y), (y, x)}]を適用して，最終的に次の表現を得る.
(
)
(
)
1 ∂
1 ∂Hx
εr ∂
1 ∂Hx
∂ 2 Hx++sx ∂x sx ∂x
sy ∂y sy εr ∂y
∂z 2
{
(
)}
∂Hy
∂
1 ∂Hy
1
− εr
+ k02 εr Hx = 0+
sx sy ∂y∂x
∂y εr ∂x()()
∂ 2 Hy
εr ∂
1 ∂Hy
1 ∂
1 ∂Hy
+
+
sx ∂x sx εr ∂xsy ∂y sy ∂y
∂z 2
{
(
)}
1
∂Hx
∂
1 ∂Hx
+
− εr
+ k02 εr Hy = 0
sx sy ∂x∂y∂x εr ∂y(2.115)(2.116)(2.117)(2.118)(2.119)統一表現
式 (2.111)，(2.112)，(2.118)，(2.119) を統一的に以下のように表現する.

∂ 2 Φx+ BΦy + k02 εr Φx = 0∂z 2∂ 2 ΦyÃxx Φy + Ãyy Φy ++ B̃Φx + k02 εr Φy = 0∂z 2

Axx Φx + Ayy Φx +式中の演算子記号と未知変数との対応は表 2.3 に示す.

20

(2.120)(2.121)

表 2.3: 式 (2.120)，(2.121) に使用する演算子の定義

Φ
1 ∂
sx ∂x

Axx ΦxE
1 ∂(εr Ex )
sx εr
∂x

Ayy Φx

1 ∂sy ∂yÃxx Φy

1 ∂
sx ∂x
1 ∂
sy ∂y

Ãyy Φy

2.4.2

{

BΦy

1 ∂
sx sy ∂xB̃Φx

1 ∂sx sy ∂y

{

{

{

(

(

1 ∂Ex
sy ∂y
1 ∂Ey
sx ∂x

}

)

(

)

1 ∂
sx ∂x

)

εr ∂
sy ∂y

)

1 ∂(εr Ey )
sy εr
∂y

H
(

εr ∂
sx ∂x

}

∂E
1 ∂(εr Ey )
− ∂yy
εr
∂y
1 ∂(εr Ex )
x
− ∂E
εr
∂x
∂x

}
}

1
sx sy
1
sx sy

{
{

(

1 ∂
sy ∂y

1 ∂Hx
sx ∂x

1 ∂Hx
sy εr ∂y1 ∂Hy
sx εr ∂x

(

1 ∂Hy
sy ∂y

∂Hy
∂
∂y∂x − εr ∂y∂Hx∂∂x∂y − εr ∂x

(

(

)

)

1 ∂Hy
εr ∂x
1 ∂Hx
εr ∂y

)}
)}

Crank-Nicolson 法と交互方向陰解法に基づく定式化

横方向電磁界 Φt (x, y, z) を +z 方向に緩慢変化する項とキャリア波数程度で振動する項にわけて[] []Φx (x, y, z)
ϕx (x, y, z)
Φt (x, y, z) ==exp(−jk0 n0 z) ≡ ϕt (x, y, z) exp(−jk0 n0 z)
Φy (x, y, z)
ϕy (x, y, z)(2.122)と表現する. 式 (2.122) を式 (2.120)，(2.121) に代入し，近軸近似 (∂ 2 ϕt /∂z 2 ≈ 0) を適用すると次のフルベクトル BPM の基本方程式が得られる.
[]∂ϕt
Axx + Ayy + νBϕt(2.123)
j2k0 n0=∂z
B̃Ãxx + Ãyy + ν()
νϕt = k02 εr − n20 ϕt(2.124)式 (2.123) 右辺の行列を

j2k0 n0∂ϕt= ([X] + [Y ]) ϕt∂z][
Axx + ν20[X] =B̃Ãxx + ν2[]
Ayy + ν2
B
[Y ] =
0Ãyy + ν2(2.125)(2.126)(2.127)と分割する. 式 (2.125) に Crank-Nicolson 法を適用して z 方向を離散化すると

j2k0 n0

ϕk+1− ϕkt
ϕk+1 + ϕktt= ([X] + [Y ]) t∆z
2となり，さらに整理すると
{
}
{
}
∆z
∆z
k+1
[I] −
([X] + [Y ]) ϕt= [I] +([X] + [Y ]) ϕt k
j4k0 n0
j4k0 n021(2.128)(2.129)

となる. ここで，[I] は単位行列，ϕkt は ϕt (x, y, k∆z) を意味している.
式 (2.129) の両辺に (∆z)2 オーダーの誤差項を加えて演算子を分割すると
{
}
(
)2
∆z
∆z
[I] −
([X] + [Y ]) +[X] [Y ] ϕt m+1
j4k0 n0
j4k0 n0
{
}
(
)2
∆z
∆z
= [I] +([X] + [Y ]) +[X] [Y ] ϕt m
j4k0 n0
j4k0 n0

(
[I] −

∆z
[X]
j4k0 n0

)(

)
∆z
[I] −
[Y ] ϕt m+1
j4k0 n0
(
)(
)
∆z
∆z
[X]
[I] +
[Y ] ϕt m= [I] +
j4k0 n0
j4k0 n0k+ 12となる. したがって，中間界 ϕtが計算できる.

1st step
(

∆z
[I] −
j4k0 n0[(2.130)(2.131)を導入して次の 2 ステップで次ステップの横方向電磁界 ϕk+1t

Axx + ν20B̃Ãxx + ν2

])

(
ϕt

k+ 12

=

∆z
[I] +
j4k0 n0[

Ayy + ν2B
0Ãyy + ν2])
ϕt k(2.132)

2nd step
(

∆z
[I] −
j4k0 n0[

B
Ayy + ν20Ãyy + ν2])(
ϕt k+1 =∆z[I] +j4k0 n0[0
Axx + ν2
B̃Ãxx + ν2])
1

ϕt k+ 2(2.133)式 (2.132)，(2.133) は，具体的に次の手順で計算できる.

1st step
{
1−
{
1−

}
{
}
∆z (
ν)
∆z (
ν)
∆z
k+ 12
Axx +
ϕx = 1 +
Ayy +
ϕkx +
Bϕky
j4k0 n02
j4k0 n02
j4k0 n0(}ν)∆zÃxx +
j4k0 n0
2

k+ 12

ϕy

=

∆z
B̃ϕx
j4k0 n0

k+ 21

{
+ 1+

(

}
ν)

∆z
Ãyy +
j4k0 n02(2.134)ϕky(2.135)

2nd step
{
1−

}
{
}
∆z (
ν)
ν)
∆z
∆z (
k+ 21
k+ 1
Ãyy +
Ã
+
ϕk+1
=
B̃ϕ
+
1
+
ϕy 2
x
xx
y
j4k0 n02
j4k0 n0
j4k0 n02(2.136)
{
}
{
}
(
)
(
)
∆z
∆z
ν
ν
∆z
k+ 1
1−
ϕk+1
= 1+
ϕx 2 +
Bϕk+1
Ayy +
Axx +
x
y
j4k0 n02
j4k0 n02
j4k0 n0(2.137)22

Axx ϕ，Ãxx ϕ，Ayy ϕ，Ãyy ϕ の離散化に Stern の方法，偏波の結合に寄与する項 Bϕ，B̃ϕ の離散化に文献 [72][73] の方法を用いると，式 (2.134)〜(2.137) は次のように表現できる.
1st step
k+ 1

k+ 1

k+ 1

2
2
L1st,1ϕx,i−1,j+ L1st,1
ϕx,i,j2 + L1st,1
ϕx,i+1,j
=
1
2
3
{ }
{
}
T
R11st,1 ϕkx,i,j−1 + R21st,1 ϕkx,i,j + R31st,1 ϕkx,i,j+1 + T 1st,1
ϕky

k+ 1

k+ 1

(2.138)k+ 1

2
2+ L1st,2=
L1st,2
ϕy,i,j2 + L31st,2 ϕy,i+1,jϕy,i−1,j2
1}{{ 1st,2 }T
k+ 12
1st,2 k1st,2 k1st,2 kR1 ϕy,i,j−1 + R2 ϕy,i,j + R3 ϕy,i,j+1 + T
ϕx(2.139)2nd step
2nd,1 k+1
L2nd,1ϕk+1ϕy,i,j + L2nd,1
ϕk+1
1
3y,i,j−1 + L2
y,i,j+1 =
{
}
{
}
1
1
1
T
k+ 12
2nd,1 k+ 2
2nd,1 k+ 2
2nd,1 k+ 2
2nd,1R1
ϕy,i−1,j + R2
ϕy,i,j + R3
ϕy,i+1,j + T
ϕx(2.140)

2nd,2 k+1
ϕk+1
ϕx,i,j + L2nd,2
ϕk+1
L2nd,23
1
x,i,j+1 =x,i,j−1 + L2
}
}
{
1
1
1
T {
k+ 2
k+
k+ 2
R12nd,2 ϕx,i−1,j
ϕk+1+ R22nd,2 ϕx,i,j2 + R32nd,2 ϕx,i+1,j+ T 2nd,2y(2.141)

ここで，式中のベクトル {T }，{ϕ} は 6 つの成分からなるベクトルであり，その具体的な導出方法
は付録 C に示す.
式 (2.138)〜(2.141) を行列・ベクトル表現で次のように表現する.

1st step1[Γ2 ] {ϕt }k+ 2 = [Γ1 ] {ϕt }k(2.142)

2nd step1[Γ4 ] {ϕt }k+1 = [Γ3 ] {ϕt }k+ 2(2.143)

ここで，{ϕt } および 2Nx Ny × 2Nx Ny の正方行列 [Γα ] (α = 1, 2, 3, 4) は次のように表現できる．
[
]
{ϕx }
{ϕt } =
(2.144)
{ϕy }
[
]
(
)
δz [Ayy ] + 1 + δz2ν [I]δz [B]()[Γ1 ] =
(2.145)0δz [Ãyy ] + 1 + δz2ν [I]
[
]
)
(
0
−δz [Axx ] + 1 − δz2ν [I]()(2.146)[Γ2 ] =
−δz [B̃]−δz [Ãxx ] + 1 − δz2ν [I][]()δz [Axx ] + 1 + δz2ν [I]0()[Γ3 ] =
(2.147)
δz [B̃]δz [Ãxx ] + 1 + δz2ν [I][]()−δz [Ayy ] + 1 − δz2ν [I]−δz [B]()(2.148)[Γ4 ] =0−δz [Ãyy ] + 1 − δz2ν [I]
ここで，δz = j4k∆zと定義している．
0 n0

23

2.4.3電力の評価方法

コアとクラッドの比屈折率差が大きな強導波路デバイスを解析する場合，電力を電界あるいは磁界のみで評価する方法では電力が保存されないことが指摘されている [74]〜[76]．電界と磁界の両方を用いて電力評価を行うことにより，より正確な電力評価が可能となることが考えられる．
モードの直交性より，n 番目のモードの複素振幅係数 cn は次のように表現できる.
)∫∫ (∗ − H E∗Hy Ex,n
x y,n dxdyΩ)
(2.149)cn = ∫∫ (∗∗Ω Hy,n Ex,n − Hx,n Ey,n dxdy
式 (2.149) の複素振幅係数は電界と磁界の固有モードフィールドをあらかじめ求めておき，未知変数を磁界として BPM 解析を実行することにより算出する．磁界を未知変数としているのは，電界は媒質境界で不連続となるためであり，磁界を未知変数として BPM 解析を実行した方が精度の向上が見込めるためである．よって，式 (2.35) より，n 番目のモードの実効電力は次式で計算できる.
{∫ ∫
}
(
)
1
∗
∗
Pn = |cn |2 Re
Ex,n Hy,n− Ey,n Hx,ndxdy
(2.150)2Ω

24

2.5 2 次元時間領域有限差分ビーム伝搬法
2.5.1

波動方程式

∂/∂y = 0 とする 2 次元問題の場合，Maxwell の方程式 (2.1)，(2.2) を成分ごとに書くと次のよう
になる．
1 ∂Ey
∂Hx
= −µ0
sz ∂z
∂t
∂Hy
1 ∂Ey
1 ∂Ez
−
= −µ0
sz ∂z
sx ∂x
∂t
1 ∂Ey
∂Hz
= −µ0
sx ∂x
∂t
∂H
1
∂Ex
y
−
= ε0 ε r
sz ∂z
∂t
∂Ey
1 ∂Hz
1 ∂Hy
−
= ε0 εr
sz ∂z
sx ∂x
∂t
1 ∂Hy
∂Ez
= ε0 εrsx ∂x∂t−(2.151)(2.152)(2.153)(2.154)(2.155)(2.156)

このとき，TE 波 (主成分は Ey ) の波動方程式が式 (2.151)，(2.153)，(2.155) から，TM 波 (主成
分は Hy ) の波動方程式が式 (2.152)，(2.154)，(2.156) からそれぞれ次のように求められる．
()()
εr ∂ 2 Ey
1 ∂
1 ∂Ey
1 ∂
1 ∂Ey
=
+
(2.157)c2 ∂t2sz ∂z sz ∂z
sx ∂x sx ∂x()()
1 ∂ 2 Hy
1 ∂
1 ∂Hy
1 ∂
1 ∂Hy
=
+
(2.158)c2 ∂t2sz ∂z sz εr ∂zsx ∂x sx εr ∂x
式 (2.157)，(2.158) は統一的に次のように書ける．

εr ∂ 2 Φ= (Dzz + Dxx ) Φc2 ∂t2(2.159)

ここで，演算子 Dαα (α = x, z) は次のように定義している．

(
)
1 ∂Φ
 1 ∂
s ∂α s ∂α ) (Φ = Ey )Dαα Φ = εα ( α
1 ∂Φ
 r ∂
sα ∂α sα εr ∂α (Φ = Hy )2.5.2

(2.160)狭帯域時間領域ビーム伝搬法

電磁界 Φ の次のような表現を考える．

Φ(x, z, t) = Re {ϕ(x, z, t) exp (jω0 t)}(2.161)ω0 は解析する光波の中心周波数であり，光波に含まれる周波数成分が ω0 に非常に近い場合，ϕ は時間方向に緩やかに振動する関数となる．式 (2.161) を式 (2.159) に代入し，SVEA
∂2ϕ
∂ϕ
≪ 2ω0
2
∂t
∂t

(2.162)25

を適用して ∂ 2 ϕ/∂t2 ≈ 0 とすると，次の狭帯域 (Narrow Band, NB)TD-BPM の基本方程式が導ける．

∂ϕ= {ζ (Dzz + Dxx ) + ξ} ϕ∂t

(2.163)

ここで，ζ ，ξ は次のように定義している．

c2
j2ω0 εr
ω0
ξ=
j2

ζ=

2.5.3

Crank-Nicolson 法と交互方向陰解法に基づく定式化

式 (2.163) に Crank-Nicolson 法を適用して整理すると次のようになる．

{1 − δt (Dzz + Dxx + ν)} ϕk+1 = {1 + δt (Dzz + Dxx + ν)} ϕk(2.164)

ここで，時間軸の刻み幅を ∆t として，ϕk は ϕ(x, z, k∆t) を表し，δt ，ν は次のように定義している．

c2 ∆t
c∆t
=
j4ω0 εr
j4k0 εr
ω02 εrν = 2 = k02 εrcδt =式 (2.164) に ADIM を適用して{((
ν )} k+1ν )} {
1 − δt Dxx +ϕ
1 − δt Dzz +
2
{2
(
)}
{
(
ν
ν )} k= 1 + δt Dzz +1 + δt Dxx +ϕ2
2(2.165)1を得る．式 (2.165) は中間界 ϕk+ 2 を導入して次のように計算できる．

1st step{((ν )} k+ 1 {ν )} k1 − δt Dzz +ϕ 2 = 1 + δt Dxx +ϕ2
2
2nd step{((ν )} k+1 {ν )} k+ 1
1 − δt Dxx +ϕ= 1 + δt Dzz +ϕ 2
22(2.166)(2.167)式 (2.166)，(2.167) は，横方向の差分化により次のように表現できる．

1st step
k+ 1

k+ 1

k+ 1

1st
1st
1st k
1st k
1st k
2
2
2
L1st
1 ϕi−1,j + L2 ϕi,j + L3 ϕi+1,j = R1 ϕi,j−1 + R2 ϕi,j + R3 ϕi,j+1(2.168)

2nd step
k+ 1

k+ 1

k+ 1

k+1
2nd
2nd
2nd k+1
2nd k+1
2nd
2
2
2
L2nd
1 ϕi,j−1 + L2 ϕi,j + L3 ϕi,j+1 = R1 ϕi−1,j + R2 ϕi,j + R3 ϕi+1,j(2.169)

ここで，ϕki,j は ϕ(i∆x, j∆z, k∆t) を意味している．本研究では，TD-BPM 解析において横方向の
1st
2nd
2nd
差分化には 2 次精度の高精度差分式 [59] を用いている．係数 L1st
α ，Rα ，Lα ，Rα (α = 1, 2, 3)の具体的な導出方法は付録 D に示す.

26

式 (2.168)，(2.169) は行列・ベクトル形式で次のように表現できる.

1st step[Γ2 ]k {ϕ}k+ 1 = [Γ1 ]k {ϕ}k(2.170)2

2nd step[Γ4 ]k {ϕ}k+1 = [Γ3 ]k {ϕ}k+ 1(2.171)2[Γα ](α = 1, 2, 3, 4) は計算領域端の境界条件の処理によって Nz Nx × Nz Nx の正方行列となる．セミベクトル BPM の場合と同様，3 重対角行列の処理のみで次ステップの界が求められる.

2.5.4

電力の評価方法

TD-BPM ではガウシアンパルスの伝搬解析によって，1 度の解析により透過電力スペクトルが算出可能である．導波路を伝搬するパルスの観測面を入射側と出力側の 2 箇所に設置し，観測パ
ルスをそれぞれ ϕin (α, t)，ϕout (α, t) [α ∈ {x, z}] とする．観測されるパルスのうち，任意のポートの固有モードフィールド成分はモードフィールドの直交性より，次式で計算できる．
∫
ϕX,n (t) =ψn∗ (α)ϕX (α, t)dα(2.172)Ω

ここで，X ∈ {in, out}，ψn は n 番目のポートにおける固有モードフィールドである．
このとき，port1 を入射ポートとすると，正規化透過電力スペクトルは次式で計算できる．

∫T
P (λ) =0∫T

ϕout,n (t) exp(−j 2πcλ t)dt2(2.173)2πc0 ϕin,1 (t) exp(−j λ t)dt

ここで，λ は波長，T は TD-BPM における最大観測時間である．式 (2.173) の分子分母はそれぞれ
FFT により効率的に計算できる．

27

2.6

虚軸ビーム伝搬法

固有モード電磁界を数値計算により算出する手法として，虚軸 BPM(Imaginaly Distance BPM，
ID-BPM)[77]〜[79] がある. ID-BPM は BPM 解析ソルバとの適合性に優れているため，本研究では固有モードの算出に ID-BPM を用いる. 以下に，ID-BPM の概要を述べる.
+z 方向に伝搬する電磁界の横方向成分 Φt は，固有モードフィールドの重ね合わせとして次のように表現できる.

Φt (x, y, z) =∑cn Φt,n (x, y) exp(−jβn z)(2.174)n

ここで，cn ，Φt,n ，βn はそれぞれ n 番目の固有モードの複素振幅，フィールドパターン，位相定数である. 伝搬軸を z 方向から虚軸方向に変更し，z → jτ とすると

Φt (x, y, z) =∑cn Φt,n (x, y) exp(βn τ )(2.175)nとなる. 入射電磁界を虚軸方向に十分長い距離伝搬させると，位相定数が最も大きな基本モードの振幅が支配的となるため，その固有モードフィールドが抽出できる. 実際のプログラムでは，∆z → j∆τ の変更を行うのみで ID-BPM による固有モード解析が実行できる.
ID-BPM では，ステップ毎に伝搬フィールドから平均実効屈折率 ne (= βn /k0 ) を見積り，参照
屈折率 n0 を ne に動的更新することによって，精度良く固有モードフィールドおよびその位相定数が算出できる. フルベクトル解析の場合を例にとると，+z に伝搬する光波の位相定数が βn である場合，波動方程式 (2.120)，(2.121) より次が導ける．
[]
Axx + Ayy + k0 εrB2βn Φt =
Φt ≡ [P ]Φt(2.176)B̃
Ãxx + Ãyy + k0 εr
両辺に (Φ∗t ·) を作用させ，横方向計算領域 Ω にわたって積分して式を整理すると次が得られる．
√ ∫∫∗Ω Φt · [P ]Φt dxdy∫∫
βn =(2.177)∗Ω Φt · Φt dxdy
式 (2.177) を用いて伝搬フィールドからその平均的な位相定数 βn を見積り，これを参照屈折率に反映させながら解析を行うと，βn は目的の固有モードの位相定数に漸近していく．
高次の固有モードを求める場合，固有モードの偶奇対称性を利用するか，あるいは Gram-Schmidtの直交化法を利用する [80]，[81]. n 番目の高次モードを求める場合，それよりも低次のモードを全て先に求めておき，低次モード成分を除去しながら ID-BPM 解析を行う．すなわち，伝搬ステッ
プ毎にΦt ← Φt −n−1∑

cm Φt,m(2.178)m=0の操作をすることによって，n 番目の高次モードフィールド Φt,n が求められる.

28

第3章

感度解析に基づく光導波路デバイスの構造
最適設計

本章では，光導波路デバイスの感度解析に基づく構造最適化手法の流れおよび最適化の実装に必要な種々の技法について述べる．トポロジー最適化では設計領域内部の材料分布を何らかの設計
変数により表現する必要があり，初歩的な構造表現手法として密度法がある．密度法は設計領域をいくつかのドメインに分割し，それぞれのドメインに密度パラメータと呼ばれる設計変数を割り当て，その密度パラメータを媒介変数として材料分布を表現する方法である．また，構造最適化問題
では，所望のデバイス出力特性を目的関数で表現し，それを最小化 (あるいは最大化) する設計変数の組を探索する．感度解析に基づいて最小化問題を解く方法として，線形計画法，MMA(Methodof Moving Asymptotes)[44] に代表される逐次凸計画法，勾配法が利用できる．本研究では，構造表
現手法および設計変数空間の探索手法はそれぞれ初歩的な手法の利用を考え，構造表現手法として密度法，探索手法として勾配法の中で最も初歩的な最急降下法を適用した場合を検討する．密度法を適用した場合，コア材料とクラッド材料の中間的な屈折率を持つ，いわゆるグレイ領域の出現やチェッカーボードパターンに代表される複雑な構造の出現が問題となる．本章では，0 から1 への遷移領域を設けた，改良した Heaviside 関数を利用したグレイ領域の除去方法および画像処理の分野で利用されている平滑化フィルタを利用した複雑な構造の除去方法についても述べる．

3.1

感度解析に基づく構造適設計

図 3.1 に本研究における構造最適化のフローチャートを示す. 本研究におけるトポロジー最適化
では，初期構造を決定したのち，BPM による光波伝搬解析，目的関数の収束判定，設計変数に対する目的関数の感度計算，設計変数の更新を繰り返すことによって最適化構造を得る. それぞれのプロセスの概要を以下に説明する.

初期構造の決定トポロジー最適化では，所望の特性を最大限満たすような構造を自動的に発現させることができるため，問題によっては一様媒質を初期構造として与えても特性の良い構造が得られる場合がある. しかし，目的関数が多峰性を有する等，目的関数の性質によっては初期構造依存性が強くなる場合があり，この場合は特別な初期構造を考える必要がある．本研究では，• 一様媒質• 既に提案されている原理に基づくデバイス構造のいずれかを初期構造として与える.

29

startInitialize design parameters

Wave propagation analysis
with BPM

Is the objective fuction
converged?

Yes

No
Sensitivity analysis

Update design parametersend
図 3.1: 勾配法に基づくトポロジー最適化のフローチャート．

30

BPM による光波伝搬解析
このプロセスでは，BPM 解析によってデバイスを伝搬する光波の振る舞いを解析し，目的関数の値を計算する.
本研究では特に，デバイスの特性が出力パワーで表現される場合を考える．出力パワーは S パ
ラメータ (正規化振幅透過係数) により表現できる. BPM 解析によってデバイスの出力ポートにおける出力電磁界を計算し，特定の固有モードに対する S パラメータを求め，その透過パワーを計算する. 固有モードの正規直交性を仮定すると，入射ポート (port 1) から n 番目のポートへの S パ
ラメータは次式で計算できる.
∫∫
)
( ∗
∗
(3.1)
Sn1 =
ψy,n ϕx,Nz −1 − ψx,nϕy,Nz −1 dxdyΩ

ここで，P1 は入射電磁界の電力であり，ϕ，ψ は (ϕ, ψ) ∈ {(E, H), (H, E)} である. ψn は n 番目のポートの固有モード電磁界フィールドであることを表しており，Nz は z 方向における分割数，
ϕNz −1 は BPM により計算された Nz − 1 ステップ目の電磁界 (出力ポートにおける電磁界) である.
横方向を有限差分法で離散化する場合，式 (3.1) の積分は区分求積法を適用して次のように表現で
きる.
(
)
(3.2)
Sn1 = ∆x∆y {ψy,n }† {ϕx }Nz −1 − {ψx,n }† {ϕy }Nz −1 ≡ {gn }† {ϕt }Nz −1[]{ψy,n }
{gn } = ∆x∆y
(3.3)
− {ψx,n }[]{ϕx }Nz −1(3.4){ϕt }Nz −1 ={ϕy }Nz −1
ここで，† はエルミート転置 (共役転置) をとることを表している．n 番目のポートにおける固有
モードの正規化出力パワーは |Sn1 |2 により計算できる．

目的関数の収束判定
収束判定は，反復回数が α 回目の目的関数の値を Cα として，

δ > |Cα − Cα−1 |(3.5)

により判定する場合が多い. ここで，δ は任意の微小量である. ただし，本研究では直線探索を行
わない最急降下法により設計空間の探索を行う設計例があり，その例では規定の回数最適化を反復し，その中で最も特性の良い構造を最適構造としている場合もある.

設計パラメータの感度解析
感度解析では，目的関数 C の設計変数 ρ に対する勾配{∂C ∂C ∂C∂C
∂ρ1 ∂ρ2 ∂ρ3 · · · ∂ρN}T≡ ∇ρ Cを計算する. ここで，N は設計パラメータの数である.

31

トポロジー最適化法では通常設計パラメータの数が非常に多く，目的関数の勾配を中心差分近似する場合，設計パラメータの数を N とすると，追加で 2N 回の BPM 伝搬解析が必要となる. 随伴変数法 (Adjoint Variable Method, AVM) は構造最適化法において広く用いられている感度解析手法であり，1 度の追加の光波伝搬解析により目的関数の勾配が計算できる. BPM を伝搬解析手法とする場合の随伴変数法に基づく感度解析の方法は以降の章で述べる．

設計パラメータの更新
本研究では，勾配法の中で最も単純な最急降下法を用いる. 最急降下法では設計変数を次式により更新する.

ρα+1 = ρα − K∇ρ C(3.6)

ここで，ρα は反復回数 α 回目における設計変数ベクトルである. また，K は任意の比例係数であり，目的関数が確実に最小化されるように直線探索によって決定することが多い. しかし，直線探索は追加の BPM 伝搬解析を必要とし，いくつかの設計例では最適化計算の効率化のため，K を定数とする場合がある.

3.2

構造表現手法

密度法では，設計領域内を正方セル等のドメインで分割し，各ドメインに正規化密度パラメー
タ ρ を割り当てて比誘電率分布を表現する. 導波路構造が 2 媒質 (比誘電率:εr,1 ，εr,2 ) で構成され，かつコアの高さが一様である場合，離散セルで分割された設計領域内部における比誘電率分布は
ρ を用いて次のように表現できる.

εr,i,k = εr,2 + (εr,1 − εr,2 ) ρi,k(3.7)

ここで，εr,i,k = εr (i∆x, k∆z)，ρi,k = ρ(i∆x, k∆z) である．
しかし，式 (3.7) で材料分布を表現する場合，ρ が 0(ホワイト) と 1(ブラック) の中間的な値 (グレイな値) をとりうる. 実際のデバイスではコアかクラッド材料に 2 値化する必要があり，このグ
レイ領域をしきい値を基準に 2 値化したとき，特性が劣化する場合がある.
本研究では，グレイ領域の出現を最終的に除去するため，次のように比誘電率を表現する.

εr,i,k = εr,2 + (εr,1 − εr,2 ) H (ρi,k )(3.8)

ここで，H(ρ) は改良した Heaviside 階段関数であり，次のように定義する.
1 (2ρ)r(0 ≤ ρ ≤ 0.5)
H(ρ) = 21 1 − (2 − 2ρ)r (0.5 < ρ ≤ 1)
2

(3.9)

ここで，r はグレイ領域の出現を抑制するペナルティ係数であり，関数 H(ρ) の値は r によって図3.2 のように変化する. 図より，r が大きくなるほど H(ρ) が階段関数に近づき，グレイ領域の出現範囲が小さくなることがわかる.

32最適化の初期段階から r が大きい値の場合，ρ の分布があまり大きく変動せず，最適化が進まな
いおそれがある. そのため，本研究では，最適化の反復回数を重ねるにつれて r の値を大きくする
方法をとる. 具体的には，(1) 最適化の総反復回数および r の最小値と最大値を指定し，反復毎に r を線形に増加させる
方法(2) r の最小値と最大値のみを指定し，目的関数の変化量が十分小さくなった場合に r を増加させる方法のいずれかの方法を採用している．すなわち (1) の方法では，r の最大値，最小値をそれぞれ rmax ，rmin ，最適化の総反復回数を N として，反復 i 回目におけるペナルティ係数を ri を

ri = rmin + (rmax − rmin )iN(3.10)

により決定する．(2) の方法では，r0 = rmin ，任意の定数を α，微小量を δ として{
αri (|Ci − Ci−1 | < δ)
ri+1 =
ri(else)(3.11)

により r を決定する．

1.0

H(ρ)

0.8
0.6
r=1

0.4

r = 100.2

r = 64
r=∞

0.0
0.0

0.2

0.4

0.6

ρ
図 3.2: 遷移領域を設けた Heaviside 階段関数．

33

0.8

1.0

3.3

構造単純化のための平滑化フィルタ

光導波路デバイスのトポロジー最適化において，実際の作成が困難である複雑な構造がしばしば出現する. 本研究では，反復毎あるいは任意の反復回数を経る毎に，正規化密度パラメータに対して 3 × 3 単純移動平均フィルタを適用し，構造の単純化を目指す.

ρfiltered
=
i,k

1
1
∑
∑
ρi+n,k+mn=−1 m=−1(3.12)

9

設計例によっては，次のように感度空間に対して 3 × 3 単純移動平均フィルタを適用する場合も
ある．

(

∂C
∂ρi,k

)filtered
=

(
)
1
1
∑
∑
1
∂C
9 ∂ρi+n,k+m

n=−1 m=−1

ここで，C は目的関数を意味している．

34(3.13)

第4章

2 次元ビーム伝搬法を用いた光導波路デバ
イスのトポロジー最適化

本章では，光波伝搬解析手法として 2 章で述べた 2 次元ビーム伝搬法，構造表現手法として 3 章で述べた密度法，感度解析手法として随伴変数法を適用する場合に対して，具体的な感度解析の方法を述べる．定式化した感度解析手法を曲がり導波路のトポロジー設計問題へ適用し，感度解析手法の妥当性を確認する．さらに，本手法を入力信号の電力を分岐するパワー分岐導波路のト
ポロジー最適設計問題へ適用し，本手法の汎用性を確認する．

4.1

感度解析手法

正規化透過電力の正規化密度パラメータに対する感度は，∂ |Sn1 |2∂S ∗∗ ∂Sn1= Sn1+ Sn1 n1
∂ρ
∂ρ
∂ρ
(
)∗
∗ ∂Sn1∗ ∂Sn1
+ Sn1
= Sn1
∂ρ
∂ρ
{
}
∗ ∂Sn1= 2Re Sn1
∂ρ

(4.1)

となり，正規化透過電力の感度は S パラメータの感度を求めることによって計算できる. 1 番目のポートから n 番目の出力ポートへの S パラメータは式 (3.2) より，

Sn1 = {gn }† {ϕ}Nz −1(4.2)となる. ここで，ベクトル {gn } は n 番目の出力ポートの固有モードフィールド {ψn } に関するベ
クトルであり，{gn } = ∆x {ψn }(4.3)

である. 離散点ごとに正規化密度パラメータが割り当てられ，座標 (i∆x, k∆z) における正規化密度パラメータを ρi,k とする．式 (4.2) の両辺を ρi,k で微分して次を得る．

∂ {ϕ}Nz −1∂Sn1= {gn }†
∂ρi,k
∂ρi,k

(4.4)

ただし，ここでは出力ポートの導波路形状が設計空間に含まれない場合を仮定している.

35

ところで，2DFD-BPM の逐次更新式 (2.33) より，出力ポートにおける電磁界フィールド {ϕ}Nz −1は次のように表現できる.

()(){ϕ}Nz −1 = [A]−1 [B] Nz −2 [A]−1 [B] Nz −3 · · ·
(
)
(
) (
)
· · · [A]−1 [B] k+1 [A]−1 [B] k [A]−1 [B] k−1 · · ·() ()· · · [A]−1 [B] 2 [A]−1 [B] 1 {ϕ}0

(4.5)

ここでは，正規化密度パラメータは離散点ごとに割り当てられると仮定しているため，ρi,k に依存する行列は [A]k ，[B]k である. そのため，式 (4.5) を ρi,k で微分すると次のようになる．
()( −1 )∂ {ϕ}Nz −1 ( −1 )∂[A]−1−1 ∂[B]kk= [A] [B] Nz −2 · · · [A] [B] k+1[B]k + [A]k
{ϕ}k
∂ρi,k
∂ρi,k
∂ρi,k
(
)
( −1 )( −1 )−1 ∂[A]k−1−1 ∂[B]k= [A] [B] Nz −2 · · · [A] [B] k+1 −[A]k[A]k [B]k + [A]k
{ϕ}k
∂ρi,k
∂ρi,k
)
(
(
)
(
)
∂[B]k
∂[A]k
{ϕ}
+
{ϕ}
= [A]−1 [B] Nz −2 · · · [A]−1 [B] k+1 [A]−1
−
(4.6)
k+1
k
k
∂ρi,k
∂ρi,k
式 (4.6) を式 (4.4) に代入して，次の式を得る.
(
)
∂Sn1
∂[A]k
∂[B]k
T
= {λn,k }−{ϕ}k+1 +
{ϕ}k
∂ρi,k
∂ρi,k
∂ρi,k

(4.7)

ここで，T は転置をとることを表しており，随伴ベクトル {λn,k }T は次のように定義している．

(
)
(
)
(
)
{λn,k }T = {gn }† [A]−1 [B] Nz −2 [A]−1 [B] Nz −3 · · · [A]−1 [B] k+1 [A]−1kこの随伴ベクトルは次の随伴方程式を解くことにより算出できる．
)
)
(
)
(
(
T
T
T
[B]T [A]−1· · · [B]T [A]−1[A]Tk {λn,k } = [B]T [A]−1
Nz −3

k+1

Nz −2

(4.8)

{gn }∗(4.9)

式 (4.7) は()∂ [A]−1∂Sn1Tk [B]k= {Ψ}k+1
{ϕ}k
∂ρi,k
∂ρi,k

(4.10)のように表現を書き換えることができる．すなわち，z = k∆z の断面に属する正規化密度パラメー
タに対する感度は，図 4.1 の感度解析の概念図に示すように，z = k∆z まで +z 方向に伝搬してきた電磁界 ϕk ，その伝搬ステップにおける伝搬演算子の微係数，および出力ポート (port n) における固有モードフィールドを z = (k + 1)∆z の断面まで「逆伝搬」させたフィールド Ψk+1 をそれぞれ求め，これらの積をとることで計算できる．Ψk+1 は，正確には逆伝搬させて得られる界分布ではなく，式 (4.9) に示すように，出力ポートにおける固有モードフィールドに対して伝搬演算子を逆方向から作用させていくことによって計算する．単純に S パラメータを中心差分により計算することも可能だが，この場合，設計変数の数の 2 倍だけ追加の BPM 伝搬解析が必要になる．上記の随伴変数法に基づく感度計算方法は，追加の伝搬解析は「逆伝搬」解析の 1 度のみで良く，効率よく S パラメータの感度が計算できる．

36

z

y
x

Ψk+1

φk
φ0

ψn∗
[A]k , [B]k

図 4.1: BPM を用いた場合の随伴変数法に基づく感度解析

4.2

曲がり導波路S 字曲がり導波路の設計モデルを図 4.2 に示す. 設計モデルにおける構造パラメータは表 4.1 にまとめて示す. ここで，入出力ポートの導波路幅 w は単一モード動作となるような導波路幅に設定する.
目的関数は，波長 1.55 µm の TE 波が port 1 に入射されたとき，port 2 への TE 波の出力が最大化，すなわち挿入損失が最小化されるように，次のように設定する.
()2(4.11)
Minimize C = 1 − |S21 |2
伝搬解析には近軸近似に基づく 2DFD-BPM を用いる．参照屈折率 n0 は入射モードの実効屈折率に固定し，離散点間隔は ∆x = 0.2 µm，∆z = 1 µm と設定する．解析領域端の境界条件は透明境界条件 (Transparent Boundary Condition, TBC)[82][83] を適用する．最速降下法における設計パ∂C
ラメータの移動量は K = 1/ maxi,j ∂ρ
に固定する．また，密度法におけるペナルティ係数 r はi,j
反復 200 回の間に 2 から 64 まで線形に変化させて設計を行う．
初期構造を ρ = 0.3 の一様媒質としてトポロジー最適設計を行い，構造平滑化フィルタの適用の有無による違いを確認する．構造平滑化フィルタは 3 × 3 移動平均フィルタを利用し，最適化の反復ごとに適用する．図 4.3 に反復回数に対する目的の出力特性のグラフを示す．図より，構造平
滑化フィルタを用いない場合は反復回数を経るにつれて特性が改善するように最適化が進み，反復 50 回程度で特性が収束している．構造平滑化フィルタを用いた場合は，フィルタを適用した影響によりわずかな出力パワーの振動がみられるが，おおよそ特性が改善するように最適化が進んでいる．2 値化した最適化構造および伝搬波形を構造平滑化フィルタを用いなかった場合を図 4.4に，用いた場合を図 4.5 に示す．図 4.4(a) と図 4.5(a) を比較すると，どちらの最適化構造も入出力
ポート付近にモード不整合を低減するような導波路オフセットが付加されており，さらにクラッドへの放射を抑制するグレーティングのような構造が現れている．構造平滑化フィルタを用いなかった図 4.4(a) の最適化構造はモザイク状の微細な構造の出現が確認できるが，構造平滑化フィルタを用いた図 4.5(a) の構造はそのような微細構造が除去され，より単純な構造が得られている．
最適化構造における正規化出力パワーは，図 4.4(a) では 0.941，図 4.5(a) では 0.955 であった．

37

図 4.6 は正規化出力パワーの波長依存特性である．構造平滑化フィルタを用いなかった場合は1.45 µm〜1.65 µm の帯域でおよそ 0.028 のパワーの変動に対し，構造平滑化フィルタを用いた場
合はおよそ 0.021 の変動がある．構造平滑化フィルタを用いた場合の方が波長依存性も小さく，出力特性が良い構造が得られていることがわかる．リブ型導波路における交差導波路の 2 次元近似によるトポロジー最適化例 [15] では，複数の動作波長におけるパワー出力を考慮した目的関数を設定することにより，単一波長のみを考慮した場合と比較して構造が単純化することが報告されている．平滑化フィルタを適用した設計により得られた図 4.5(a) の構造は，微細な構造が除去されたことによって比較的広帯域動作が実現されたのではないかと考えられる．

38

y

z
d

x port 1
w
n1

Design regionDport 2

w
d
n2
l

L

l

図 4.2: S 字曲がり導波路の 2 次元近似設計モデル

表 4.1: 2 次元 S 字曲がり導波路の設計モデルにおける構造パラメータ

n1
1.45

n2
1.445

w
5 µm

L
800 µm

39

d
30 µm

l
100 µmD
40 µm

Normalized output power

1.0
0.8
0.6

without smoothing filter
with smoothing filter

0.4
0.2
0.0

0

50

100

150

200

Number of iterations
図 4.3: 反復回数に対する S 字曲がり導波路の正規化出力パワー(b) 伝搬波形 (|Ey |)(a) 最適化構造

図 4.4: 構造平滑化フィルタを適用しなかった場合の S 字曲がり導波路の最適化構造と伝搬波形(b) 伝搬波形 (|Ey |)(a) 最適化構造

図 4.5: 構造平滑化フィルタを適用した場合の S 字曲がり導波路の最適化構造と伝搬波形

40Normalized output power

1.00

without smoothing filter
with smoothing filter

0.98
0.96
0.94
0.92
0.90
0.88
0.86
0.84
1.45

1.5

1.55

1.6

1.65

Wavelength [µm]図 4.6: S 字曲がり導波路の最適化構造における正規化出力パワーの波長依存特性

41

4.3

パワー分岐導波路

4.3.1

パワー 2 分岐導波路

パワー 2 分岐素子の設計モデルを図 4.7 に示す. 設計モデルにおける構造パラメータは表 4.2 にまとめて示す. ここで，入出力ポートの導波路幅 w は単一モード動作となるような導波路幅に設定している.
目的関数は，波長 1.55 µm の TE 波が port 1 に入射されたとき，port 2，port 3 へパワーが 1:1 で等分岐するように，次のように設定する.

(
Minimize C =1− |S21 |2
2

(

)2
+

1
− |S31 |2
2)2(4.12)

伝搬解析には近軸近似に基づく 2DFD-BPM を用いる．参照屈折率 n0 は入射モードの実効屈折率に固定し，離散点間隔は ∆x = 0.2 µm，∆z = 1 µm と設定する. 解析領域端の境界条件は TBCを適用する. また，密度法におけるペナルティ係数 r は反復 200 回の間に 2 から 64 まで線形に変∂Cに固定化させて設計を行う. 最速降下法における設計パラメータの移動量は K = 1/ maxi,j ∂ρi,jする．さらに，port 2，3 にパワーが均等に分岐するように，x = 50 µm を対称軸とする構造の対称条件を課して設計を行う. 加えて，より単純な構造を得るため，構造平滑化フィルタ (3 × 3 移動
平均フィルタ) を最適化の反復ごとに適用して設計を行う.
初期構造を ρ = 0.3 の一様媒質としてトポロジー最適設計を行った場合の，反復回数に対する目
的の出力特性のグラフを図 4.8 に示す. 図より，構造平滑化フィルタの影響で特性がわずかに振動しているが，反復回数を経るにつれて出力特性が目標値の 0.5 に向かって改善するように最適化が進んでいることがわかる.
2 値化した最適化構造と伝搬波形を図 4.9 に示す. 図 4.9(a) より，最適化構造では，出力ポート付近にモード不整合を抑制するオフセットが付加され，放射を抑制するグレーティングのような構造も現れている. また，導波路の分岐部分では，直進方向への放射を抑制する構造が現れており，このような構造は文献 [84] などで提案されている. 最適化構造における port 2，3 での出力特性はともに 0.490 であった.

42

y

z

x

d2

d1

port 2w

port 1w

Design region D

port 3

w

n1
d2

n2
l

L

l

図 4.7: パワー 2 分岐素子の 2 次元近似設計モデル

表 4.2: パワー 2 分岐素子の設計モデルにおける構造パラメータ

n1
1.45

n2
1.445

w
5 µm

L
500 µm

d1
50 µm

43

d2
30 µm

l
50 µm

D
40 µm

Normalized output power

0.5
0.4
0.3

port2
port3

0.2
0.1
0.0

0

50

100

150

Number of iterations
図 4.8: 反復回数に対するパワー 2 分岐導波路の正規化出力パワー(b) 伝搬波形 (|Ey |)(a) 最適化構造

図 4.9: パワー 2 分岐導波路の最適化構造と伝搬波形

44

200

4.3.2

パワー任意分岐導波路(X : Y ) の任意パワー分岐導波路の最適化では，分岐比と過剰損失の項を分けて目的関数を次のように設定する．
Minimize C = W1,α C1 + W2,α C2()2
C1 = X |S31 |2 − Y |S21 |2()2
C2 = 1 − |S21 |2 − |S31 |2(4.13)(4.14)(4.15)

ここで，W1,α ，W2,α は重みであり，次のように選ぶ．

W1,α =
W2,α =

C1,α(4.16)

C1,α + C2,αC2,α(4.17)

C1,α + C2,α

ここで，Cn,α は反復 α 回目における目的関数の値である．この重みの導入により，C1 と C2 のうち，特性のより悪い方の目的関数が優先的に最小化されると考えられる.
解析条件および最適化の設計条件はパワー等分岐素子の時と同様とし，初期構造を図 4.9(a) に示すパワー 2 分岐導波路の最適化構造として，(3:2)，(2:1)，(3:1) パワー分岐素子の最適設計をそれぞれ行った. 図 4.10 にそれぞれのパワー分岐素子の 2 値化した最適化構造と伝搬波形を示す. 表
4.3 は最適化波長 (1.55 µm) におけるそれぞれのパワー分岐導波路の過剰損失と分岐比である. な
お，過剰損失は次の式によって計算している．

()
Excess loss = −10 log10 |S21 |2 + |S31 |2(4.18)

表より，最適化構造における各パワー分岐導波路の分岐比をみると，3:2 および 2:1 分岐導波路では目標値の 1.5，2.0 を達成しており，3:1 分岐導波路においても 2.92 と，目的関数に設定した特性がほぼ実現されていることがわかる. 過剰損失を見ると，最も過剰損失が大きい 3:1 分岐導波路においてもその損失が 0.286 dB であり，低損失な分岐導波路が実現できている．図 4.11 は最適化構
造における挿入損失と分岐比の波長依存特性である. 図より，最適化構造は挿入損失，分岐比ともに顕著な波長依存性は見られず波長平坦な動作が実現されている.
表 4.3: 最適化したパワー分岐導波路の波長 1.55 µm における過剰損失と分岐比
過剰損失 [dB]

分岐比

0.163
0.108
0.286

1.50
2.00
2.92

(3:2) power divider(2:1) power divider(3:1) power divider

45

(ii) 伝搬波形 (|Ey |)(i) 最適化構造(a) (3:2) 分岐素子(ii) 伝搬波形 (|Ey |)(i) 最適化構造(b) (2:1) 分岐素子(ii) 伝搬波形 (|Ey |)(i) 最適化構造(c) (3:1) 分岐素子

図 4.10: パワー任意分岐導波路の最適化構造と伝搬波形

46

0.8

1:1
3:2
2:1
3:1

Excess loss [dB]

0.7
0.6
0.5
0.4
0.3
0.2
0.1
0.0
1.45

1.50

1.55

1.60

1.65

Wavelength [µm](a) 過剰損失

4.5

1:1
3:2
2:1
3:1

Splitting ratio

4.0
3.5
3.0
2.5
2.0
1.5
1.0
0.5
1.45

1.50

1.55

1.60

1.65

Wavelength [µm](b) 分岐比図 4.11: パワー 2 分岐導波路の最適化構造における過剰損失と分岐比の波長依存特性

47

第5章

3 次元セミベクトルビーム伝搬法を用いた
光導波路デバイスのトポロジー最適化

前章では 2 次元問題に対する BPM を用いた光導波路のトポロジー最適化法を提案したが，上部
クラッドや基板への光波の放射を考慮したより実際的な設計を行うには 3 次元解析が必要である．
本章では偏波の結合が無視できる場合に対して有効な，3 次元セミベクトル BPM を用いた感度解析に基づくトポロジー最適化法について検討を行う．はじめに，構造表現手法に密度法，光波の伝搬解析手法に ADI 法を適用したセミベクトル BPM を用いた場合に対して，AVM に基づく感度解析手法の定式化を行う．提案アプローチを曲がり導波路の設計問題に適用し，全ての未知変数(Ex ，Ey ，Hx ，Hy ) においてトポロジー最適設計が妥当に行われることを示す．また，パワー分
岐導波路，モード次数変換素子への設計問題へと適用し，提案手法の有効性および汎用性を確認
する．

5.1

感度解析手法

セミベクトル問題の場合，1 番目のポートから n 番目の出力ポートへの S パラメータは式 (3.2)
より，

Sn1 = {gn }† {ϕ}Nz −1(5.1)

で算出でき，固有モードフィールドに関するベクトル {gn } は次のように表現される．

{gn } = ∆x∆y {ψn }(5.2)

ここでは，設計領域における屈折率分布は x，z 方向のみに構造の変化があり，コアの高さが一様な平面光波回路素子の設計を考える．このとき，正規化密度パラメータ分布は 2 次元問題同様，設計領域内部において x と z に依存する関数となる．離散点ごとに正規化密度パラメータを割り当てる場合を考え，座標 (i∆x, k∆z) における正規化密度パラメータを ρi,k とする．式 (5.1) の両辺
を ρi,k で微分して次を得る．

∂ {ϕ}Nz −1∂Sn1= {gn }†
∂ρi,k
∂ρi,k

(5.3)

ここで，出力ポートの導波路構造は設計領域に含まれないことを仮定している.
ADI 法に基づくセミベクトル FD-BPM の逐次更新式 (2.170)，(2.171) より，出力ポートにおけ48

る電磁界フィールド {ϕ}Nz −1 は次のように表現できる.

()(){ϕ}Nz −1 = [Γ4 ]−1 [Γ3 ] [Γ2 ]−1 [Γ1 ][Γ4 ]−1 [Γ3 ] [Γ2 ]−1 [Γ1 ]
···
Nz −3
( Nz −2) ()()· · · [Γ4 ]−1 [Γ3 ] [Γ2 ]−1 [Γ1 ][Γ4 ]−1 [Γ3 ] [Γ2 ]−1 [Γ1 ][Γ4 ]−1 [Γ3 ] [Γ2 ]−1 [Γ1 ]
···
k+1
k
k−1
(
) (
)
· · · [Γ4 ]−1 [Γ3 ] [Γ2 ]−1 [Γ1 ][Γ4 ]−1 [Γ3 ] [Γ2 ]−1 [Γ1 ] {ϕ}0
2

1

(5.4)
離散点ごとに正規化密度パラメータを割り当てる場合を考えているため，ρi,k に依存する行列は[Γα ]k (α = 1, 2, 3, 4) である. そのため，式 (5.4) を ρi,k で微分して次を得る．

)()∂ {ϕ}Nz −1 (= [Γ4 ]−1 [Γ3 ] [Γ2 ]−1 [Γ1 ]· · · [Γ4 ]−1 [Γ3 ] [Γ2 ]−1 [Γ1 ]
∂ρi,k
Nz −2
k+1
(
∂ [Γ4 ]k∂ [Γ3 ]k− [Γ4 ]−1{ϕ}k+1 + [Γ4 ]−1{ϕ}k+ 1
k
k
2
∂ρi,k
∂ρi,k
)
−1 ∂ [Γ2 ]k−1 ∂ [Γ1 ]k− [Γ4 ]−1{ϕ}k+ 1 + [Γ4 ]−1{ϕ}kk [Γ3 ]k [Γ2 ]kk [Γ3 ]k [Γ2 ]k
2
∂ρi,k
∂ρi,k

(5.5)

式 (5.5) を式 (5.3) に代入して，次の式を得る.

{{}T ∂ [Γ4 ]k}T ∂ [Γ3 ]k∂Sn1= − λn,k+1/2{ϕ}k+1 + λn,k+1/2{ϕ}k+ 1
2
∂ρi,k
∂ρi,k
∂ρi,k
1
∂ [Γ1 ]k∂ [Γ2 ]k{ϕ}k+ 2 + {λn,k }T{ϕ}k
− {λn,k }T
∂ρi,k
∂ρi,k
{
}
ここで，随伴変数ベクトル λn,k+1/2 ，{λn,k } は次のように定義している．
(){}T
λn,k+1/2= {gn }† [Γ4 ]−1 [Γ3 ] [Γ2 ]−1 [Γ1 ]{{λn,k }T = λn,k+1/2

}T

Nz −2

(
)
· · · [Γ4 ]−1 [Γ3 ] [Γ2 ]−1 [Γ1 ]

(5.6)

k+1

[Γ4 ]−1k(5.7)[Γ3 ]k [Γ2 ]−1
k

(5.8)

すなわち，随伴変数ベクトルは次の随伴変数方程式を解くことにより算出できる.
)
{
} (
T
T
···
[Γ4 ]Tk λn,k+1/2 = [Γ1 ]T [Γ2 ]−1 [Γ3 ]T [Γ4 ]−1
k+1
(
)
T
T
{gn }∗· · · [Γ1 ]T [Γ2 ]−1 [Γ3 ]T [Γ4 ]−1
Nz −2
{
}
[Γ2 ]Tk {λn,k } = [Γ3 ]Tk λn,k+1/2(5.9)(5.10)

2 次元問題の場合と同様，光波伝搬の順方向と「逆方向」の 2 度の BPM 解析により S パラメータに対する設計変数の感度を効率よく計算できる．

49

5.2

曲がり導波路曲がり導波路の設計問題を通して，前節で定式化した感度解析手法の妥当性を確認する．S 字曲
がり導波路の設計モデルを図 5.1 に示す．設計モデルにおける構造パラメータは表 5.1 にまとめて示す．目的関数は，port 1 に波長 1.55 µm の基本モードが入射されたとき，port 2 におけるその出力が最大化されるように，次のように設定する．
()2
Minimize C = 1 − |S21 |2(5.11)
伝搬解析には近軸近似と ADI 法に基づくセミベクトル FD-BPM を用いる．参照屈折率 n0 は入射モードの実効屈折率に固定し，離散点間隔は ∆x = ∆y = 0.2 µm，∆z = 1 µm と設定する．解析
領域端には厚さ 10 µm の PML を装荷し，PML 終端は電気壁あるいは磁気壁とする．
はじめに，全ての未知変数 (Ex ，Ey ，Hx ，Hy ) に対して感度解析が妥当に実行されるか確認する．
感度解析の妥当性の確認のため，最速降下法におけるステップ幅は十分小さく K = 0.01/ maxi,j

∂C
∂ρi,j

とし，平滑化フィルタを適用せずに最適設計を行う．初期構造は ρ = 0.3 の一様媒質を与え，密度法におけるペナルティ係数 r は反復 200 回の間に 8 から 1024 まで線形に変化させて設計を行う．
図 5.2 にそれぞれの未知変数に対する，最適化の反復回数に対する正規化出力パワーのグラフを示す．図 5.2 およびその挿入図より，全ての未知変数において反復回数が進むにつれて特性が改善していることが確認できる．したがって，本トポロジー最適化手法がどの未知変数においても妥当に実行できることが確認できた．図 5.3 に未知変数の違いによる最適化構造とその伝搬波形を示す．図 5.3(a)(c)(e)(f) に示す最適化構造を見ると，未知変数によって最適化構造に大きな違いは無く，入出力ポート部においてモード不整合を低減する導波路オフセットが付加されていることが確認できる．また，2 次元近似設計における曲がり導波路の最適化構造とは異なり，放射を抑制するようなグレーティング構造は現れておらず，代わりに導波路幅を広くしてコアへの閉じ込めを強くし，深さ方向への放射を抑制するような導波路構造となっているように見える．これは，あえてクラッド部へ放射させるような構造だと，実際には放射損失が大きくなってしまうためと考えられ，3 次元最適設計によって実際的な最適化構造が得られることがわかる．各最適化構造において，正規化出力パワーはそれぞれ 0.929(Ex ), 0.930(Hy ), 0.924(Ey ), 0.929(Hx ) であり，挿入損失が 0.343 dB 以下の曲がり導波路が実現されている．

50

y

z

y

z

PML

h

x x port 1wdDesign region
n1 or n2

n1
H

H

D

w
d

n2
l

l

L

図 5.1: S 字曲がり導波路の最適設計モデル

n1
1.45

表 5.1: S 字曲がり導波路の設計モデルにおける構造パラメータ
n2
w
h
l
d
D
H
L

1.445

8 µm

4 µm

port 2100 µm

51

30 µm

20 µm

15 µm

900 µm

Ex
Hy
Ey
Hx

Normalized power

Normalized power

1.0
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0.0

1.0
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0.0

Ex
Hy
Ey
Hx
0

2

4

6

8

10

Number of iterations

ρ = 0.3

0

50

100

150

200

Number of iterations
図 5.2: 各未知変数の違いによる反復回数に対する S 字曲がり導波路の正規化出力パワー

52

(i)
(ii)

(i)
(iii)

(ii)

(iii)

(b) 伝搬波形 (ϕ = Ex )(a) 最適化構造 (ϕ = Ex )
(i)
(ii)

(i)
(iii)

(ii)

(iii)

(d) 伝搬波形 (ϕ = Ey )(c) 最適化構造 (ϕ = Ey )
(i)
(ii)

(i)
(iii)

(ii)

(iii)

(f) 伝搬波形 (ϕ = Hx )(e) 最適化構造 (ϕ = Hx )
(i)
(ii)

(i)
(iii)

(ii)

(iii)

(h) 伝搬波形 (ϕ = Hy )(g) 最適化構造 (ϕ = Hy )図 5.3: 未知変数の違いによる S 字曲がり導波路の最適化構造とその伝搬波形

53続いて，平滑化フィルタ (3 × 3 移動平均フィルタ) を最適化の反復毎に適用して最適設計を行った結果を示す．また，ここでは，初期構造による最適化構造の違いについても議論する．最速降∂C
下法におけるステップ幅は K = 1/ maxi,j ∂ρ
に固定し，未知変数は Hy とする．初期構造はi,j

ρ = 0.3 の一様媒質あるいは入出力ポートをつなぐ傾斜直線導波路 (コア部:ρ = 0.8，クラッド部
ρ = 0.3) を与え，密度法におけるペナルティ係数 r は反復 200 回の間に 8 から 1024 まで線形に変化させて設計を行う．図 5.4 に最適化の反復回数に対する正規化出力パワーのグラフを示す．図より，平滑化フィルタを適用しているため出力特性がわずかに振動しているが，おおむね反復毎に特性が改善していることがわかる．また，特性の改善の様子は異なるが，初期構造にかかわらず
どちらもほぼ同じ出力値に収束している．図 5.5 に最適化構造とその伝搬波形を示す．図 5.5(a)(c)の挿入図を見ると，平滑化フィルタを適用しない場合の最適化構造と比較して微細な構造が除去されており，より実際的な構造が得られている．最適化構造はどちらの構造も導波路オフセットが付加されており，かつ幅の太い導波路が形成されている．放射損失を低減する原理としてはどちらも共通であると考えられるため，この問題においては，最適化構造の顕著な初期構造依存性は見られない．正規化出力パワーは図 5.5(a) において 0.947，図 5.5(c) において 0.941 であった．挿入損失が 0.254 dB 以下であり，平滑化フィルタを適用しない場合より損失が改善している．これは，微細な構造による散乱損失の影響が構造平滑化フィルタの影響により除去されたことが一つの理由と考えられる．

54

Normalized power

1.0
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0.0

Uniform medium
Tilted waveguideρ = 0.8
ρ = 0.3

ρ = 0.3

Uniform medium
050

Tilted waveguide
100

150

200

Number of iterations
図 5.4: 構造平滑化フィルタを適用した場合の反復回数に対する S 字曲がり導波路の正規化出力
パワー

(i)
(ii)

(i)
(iii)

(ii)

(iii)

(b) (a) における伝搬波形 (|Hy |)(a) 最適化構造 (初期構造:一様媒質)
(i)
(ii)

(i)
(iii)

(ii)

(iii)

(d) (c) における伝搬波形 (|Hy |)(c) 最適化構造 (初期構造:直線傾斜導波路)

図 5.5: 構造平滑化フィルタを適用した場合の S 字曲がり導波路の最適化構造と伝搬波形

55

5.3

パワー 2 分岐導波路

入力パワーを 2 等分岐する導波路の設計問題を通して，本手法のさらなる汎用性を確認する．パワー 2 分岐導波路の 3 次元設計モデルを図 5.6 に示す．設計モデルにおける構造パラメータは表5.2 にまとめて示す．目的関数 C は，port 1 に波長 1.55 µm の x 偏波基本モードが入射されたとき，正規化出力パワーが port 2，port 3 へ等分岐されるように，次のように設定する．
(
)2 (
)2
1
1
− |S21 |2 +− |S31 |2
Minimize C =(5.12)
22
伝搬解析には近軸近似と ADIM に基づくセミベクトル FD-BPM を用いる．参照屈折率 n0 は入射モードの実効屈折率に固定し，離散点間隔は ∆x = ∆y = 0.2 µm，∆z = 1 µm と設定する．解析領域端には厚さ 10 µm の PML を装荷し，PML 終端は電気壁あるいは磁気壁とする．最速降下∂C法におけるステップ幅は K = 1/ maxi,j ∂ρ
に固定し，未知変数は Hy とする．密度法における
i,jペナルティ係数 r は反復 200 回の間に 8 から 1024 まで線形に変化させて設計を行う．さらに，出力パワーのばらつきが生じないようにするため，x = 50 µm を対称面とする構造対称条件を課して設計を行う．加えて，微細な構造の出現を回避するため，構造平滑化フィルタ (3 × 3 移動平均フィルタ) を最適化の反復ごとに適用して設計を行う.
本設計においても，異なる 2 つの初期構造を与えて，最適化構造の初期構造依存性を確かめる．
初期構造として設計領域において ρ = 0.3 の一様媒質および直線 Y 分岐導波路 (コア領域:ρ = 0.8，
クラッド領域 ρ = 0.3) を与える．図 5.7 に最適化の反復回数に対する正規化出力パワーのグラフを示す．port 2 と port 3 の出力は対称条件を課しているため同じであり，ゆえに図では片方の portの正規化出力パワーのみを示している．図より，平滑化フィルタの影響による特性の振動が見られるが，どちらの場合も最適化の反復毎に目標値である 0.5 に向かって特性が改善していることが確認できる．図 5.8 に最適化構造および伝搬波形を示す．図を見ると，どちらの最適化構造も出力ポート付近にモード不整合を抑制するオフセットが付加されている．また，図 5.8(a)(c) の拡大図に示すように，導波路の分岐部では類似した構造が現れている．この構造は文献 [84] 等で提案されているような，分岐部分における正面方向への放射損失を低減する構造になっていると考えられる．そのため，動作原理の観点から見ると，最適化構造の顕著な初期構造依存性は見られない．
最適化構造における port 2(port 3) における正規化出力パワーは，図 5.8(a) では 0.471，図 5.8(c) では 0.478 であり，過剰損失にしてそれぞれ 0.259 dB，0.195 dB である．初期構造として直線 Y 分
岐導波路を用いたが，この構造をコア・クラッドに 2 値化した構造において port 2(port 3) の出力
パワーは 0.413，過剰損失にして 0.830 dB であり，この構造と比較して大幅に特性が改善されて
いる．

56

y

z

y

z

PML

h
d1

x x

d2

w

port 1w

Design region
n1 or n2n1H

port 2HD

port 3

w
d2

n2
l

l

L

図 5.6: パワー 2 分岐導波路の 3 次元最適設計モデル

表 5.2: パワー 2 分岐導波路の 3 次元設計モデルにおける構造パラメータ
n1
n2
w
h
l
d1
d2
H
L
D

1.45

1.445

8 µm

4 µm

50 µm

40 µm

57

20 µm

15 µm

900 µm

40 µm

Normalized power

0.50
0.45
0.40
0.35
0.30
0.25
0.20
0.15
0.10
0.05
0.00

Uniform medium
Straight Y−branchρ = 0.3
ρ = 0.3
ρ = 0.8

Uniform medium
050

Straight Y−branch
100

150

200

Number of iterations
図 5.7: 反復回数に対するパワー 2 分岐導波路の正規化出力パワー(a) 最適化構造 (初期構造:一様媒質)(b) (a) における伝搬波形 (|Hy |)(c) 最適化構造 (初期構造:直線 Y 分岐導波路)(d) (c) における伝搬波形 (|Hy |)図 5.8: パワー 2 分岐素子の最適化構造と伝搬波形

58

5.4

モード次数変換素子

既存の波長分割多重 (Wavelength Division Multiplexing, WDM) システムの伝送容量の限界を越えるため，モード分割多重 (Mode Division Multiplexing, MDM) システムが次世代システムの候補として考案されている．MDM システムにおいてモード多重された信号の生成・分離を行う集積光回路ベースのモードマルチプレクサにおいて，モード次数変換器は必須の基本要素である．本章のこれまでの設計例では，単一モード動作の導波路素子を考えてきたが，ここでは x 偏波第 1 高次
モード (TE1 モード) を基本モード (TE0 ) に変換するモード次数変換素子の設計を行い，本設計手法が複数のモードを扱う素子に対しても有効であることを示す．
図 5.9 にモード次数変換素子の 3 次元設計モデルを示す．図中の構造パラメータは表 5.3 にまとめて示す．目的関数 C は，port 1 に波長 1.55 µm の TE1 モードが入射されたとき，port 2 における基本モードの正規化出力パワーが最大化されるように，次のように設定する．
()
2 2
Minimize C = 1 − S21(TE1→TE0),(5.13)
2

ここで， S21(TE1→TE0) は port 1 に TE1 モードが入射されたときの port 2 における TE0 モードの正規化出力パワーである．伝搬解析には近軸近似と ADIM に基づくセミベクトル FD-BPM を用いる．参照屈折率 n0 は入射モードの実効屈折率に固定し，離散点間隔は ∆x = ∆y = 0.2 µm，∆z = 2 µm と設定する．解析領域端には厚さ 10 µm の PML を装荷し，PML 終端は電気壁あるい∂Cは磁気壁とする．最速降下法におけるステップ幅は K = 1/ maxi,j ∂ρ
に固定し，未知変数は
i,j

Hy とする．密度法におけるペナルティ係数 r は反復 1500 回の間に 1 から 1024 まで線形に変化させて設計を行う．また，構造平滑化フィルタ (3 × 3 移動平均フィルタ) を最適化の反復ごとに適用して設計を行う.
初期構造に ρ = 0.3 の一様媒質と入出力ポートを結ぶ直線テーパ (コア領域:ρ = 0.8，クラッド
領域 ρ = 0.3) を与え，最適化結果の初期構造依存性を比較する．図 5.10 に最適化の反復回数に対する正規化出力パワーを示す．図より，初期構造によって出力パワーの収束値に多きな差が現れており，一様媒質を初期構造として与えた場合の方が良い特性が得られている．本章におけるこれまでの曲がり導波路やパワー分岐導波路の設計では，最適化構造の初期構造依存性に顕著な差は見られなかったが，このように問題よっては初期構造の選択は一つの重要な要素になることがわかる．図 5.11 に最適化構造と伝搬波形を示す．TE0 モードの正規化出力パワーはそれぞれ，図5.11(a) では 0.911，図 5.11(c) では 0.800 である．図 5.11(a) の構造は [85] で提案されているような
非対称 Y 分岐型のような構造に近く，モード間の断熱変化を引き起こすような構造になっていると考えられる．図 5.11(b) は導波路幅の違いを利用した位相調節，および MMI の動作によってモード変換が実現されていると考えられる．断熱変化をモード変換の原理として利用する場合，通常構造が伝搬方向に緩やかに変化するため，導波路構造の変化に対する光波の散乱損失が比較的低く抑えられる．図 5.11(a) の構造が比較的良い特性を持つのは，そのような断熱構造を持つのが理由と考えられる．
最後に，図 5.11(a) の最適化構造が非対称直線 Y 分岐構造よりもコンパクトかつ特性が良いことを示す．非対称直線 Y 分岐として図 5.12 のようなモデルを考え，w11 ，w12 ，および d12 の寸法最適値と結果を比較する．その他図の構造パラメータは表 5.3 と一致させる．図 5.13 に寸法最適化
構造と伝搬波形を示す．寸法の最適値はそれぞれ w11 = 15 µm，w12 = 8.6 µm，d12 = 12 µm であった．寸法最適化構造における TE0 モードの正規化出力パワーは 0.846 であり，トポロジー最適化構造の方が確かに寸法最適化構造よりも良い特性が得られている．
59

z

y

y

z

PML

d

d

port 1

w1

h
H

x x
n1

D

port 2

Design region
n1 or n2

w2

n2

H

L

l

l

図 5.9: モード次数変換素子の 3 次元最適設計モデル

表 5.3: モード次数変換素子の 3 次元設計モデルにおける構造パラメータ
n1
n2
w1
w2
h
d
L
l
H
D

1.45

1.445

18 µm

6 µm

4 µm

50 µm

60

40 µm

15 µm

900 µm

80 µm

Normalized power

1.0
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0.0

Uniform medium
Linear taperρ = 0.8
ρ = 0.3
ρ = 0.3
Uniform medium

0

300

600

Linear taper

900

1200

1500

Number of iterations
図 5.10: 最適化の反復回数に対する TE0 モードの正規化出力パワー(a) 最適化構造 (初期構造:一様媒質)(b) (a) における伝搬波形 (|Hy |)(c) 最適化構造 (初期構造:線形テーパ導波路)(d) (c) における伝搬波形 (|Hy |)図 5.11: モード次数変換素子の最適化構造とその伝搬波形

61

z

y

y

z

PML

x x

d

port 1

w1

h
H

n1
D

H

d12

w11

d

w12
port 2
w2

n2

L

l

l

図 5.12: 非対称直線 Y 分岐構造に基づくモード次数変換素子の寸法最適設計モデル(b) 伝搬波形 (|Hy |)(a) 寸法最適化構造

図 5.13: 非対称直線 Y 分岐構造に基づくモード次数変換素子の寸法最適化構造とその伝搬波形

62

第6章

2 次元時間領域ビーム伝搬法を用いた光導波路デバイスのトポロジー最適化

本章では，伝搬解析手法に 2 次元 TD-BPM を用いたトポロジー最適設計について，感度解析手法およびその有効性を検討する．これまでの章では，コアとクラッドの比屈折率差が非常に小さな導波路デバイスのトポロジー最適化について検討を行ってきた．しかし，通常の BPM は後方反射波を考慮できないため，これまでのアプローチを比屈折率差が大きな導波路デバイスの設計に適用することは難しい．BPM を時間領域手法へと拡張した TD-BPM は時間方向に緩慢変化する電磁界を未知量とし，広角伝搬ないし後方散乱する光波が考慮可能である．また，TD-BPM は通常陰解法であり，同じく時間領域手法である陽解法の FDTD 法と比較して時間方向の刻み幅を大きくとることが可能である．さらに，TD-BPM はガウシアンパルスの伝搬解析によって，一度の伝搬解析で透過スペクトルが算出できる特徴を持つ．それゆえ，TD-BPM を用いたトポロジー最適化の実現によって，広帯域動作を目的とした比屈折率差の大きな導波路デバイスの効率的なトポ
ロジー最適化が期待できる．
本章では，はじめに TD-BPM を用いた場合の感度解析手法を示し，その後，曲がり導波路のト
ポロジー最適設計を通して感度解析手法の妥当性を確認する．さらに，導波路型反射器の設計を通して，本アプローチが広帯域動作などの透過・反射スペクトル自体を最適化したい場合に便宜な手法であることを示す．

63

6.1

感度解析手法導波路デバイスの性能を表す量として，次の重なり積分を考える．
∫∫
ηn1 (t) =ψn∗ (x, z)ϕ(x, z, t)dxdz,

(6.1)

Ω

ここで，Ω は計算領域，ψn は port n における所望の出力固有モードフィールドである．|ηn1 (t)| のピークが時刻 t = τ ∆t に現れるとし，そのピーク値をτ
max {|ηn1 (t)|} = |ηn1 (τ ∆t)| ≡ |ηn1
|

(6.2)

t

τ は次の形式で表現できる．
と定義する．式 (6.1) を区分求積法により算出する場合，ηn1
τηn1= {gn }† {ϕ}τ ,(6.3)

τ を ρ(i∆x, j∆z)
ここで，{gn } = ∆x∆z {ψn } である．入出力ポートが設計領域に含まれない場合，ηn1(≡ ρij ) で微分して次を得る．
τ∂ {ϕ}τ
∂ηn1= {gn }†
∂ρij
∂ρij

(6.4)

TD-BPM の逐次更新式 (2.170)，(2.171) より，式 (6.4) は次のように書き換えられる．
τ −1
τ
)
∂ηn1
∂ ∏(
= {gn }†
[Γ4 ]−1 [Γ3 ][Γ2 ]−1 [Γ1 ] k {ϕ}0
∂ρij
∂ρij

(6.5)

k=0

総乗の微分を書き換えて，最終的に以下の式を得る．
τ [
τ
∑
{
}T ∂[Γ4 ]k{}T ∂[Γ3 ]k∂ηn1=− λn,k−1/2{ϕ}k + λn,k−1/2{ϕ}k−1+ 1
2
∂ρij
∂ρij
∂ρij
k=1
]
T ∂[Γ2 ]kT ∂[Γ1 ]k− {λn,k−1 }{ϕ}k−1+ 1 + {λn,k−1 }{ϕ}k−1
2
∂ρij
∂ρij

(6.6)

随伴行列ベクトル {λn,k } は，次の随伴方程式を解くことにより算出できる．

{
}
(
)T
(
)T
[Γ4 ]Tk λn,k+1/2 = [Γ1 ]Tk+1 · · · [Γ2 ]−1[Γ3 ]Tτ−1 [Γ4 ]−1{gn }∗
τ −1
τ −1
{
}
[Γ2 ]Tk {λn,k } = [Γ3 ]Tk λn,k+1/2

(6.7)
(6.8)

ADI 法に基づくセミベクトル BPM の場合の感度計算式 (5.6) と比較すると，右辺に総和記号が追τ の感度は順方向と「逆加されている点が異なる．しかし，これまでの感度解析手法と同様に，ηn1方向」の 2 回の伝搬解析のみで効率よく計算可能である．

64

6.2

曲がり導波路コアとクラッドの比屈折率差が大きな強導波路における曲がり導波路の 2 次元近似設計を通して，提案する感度解析手法の妥当性を確認する．曲がり導波路の 2 次元近似設計モデルを図 6.1 に示す．図中の構造パラメータは表 6.1 にまとめて示す．時刻 t = 0 における初期電磁界として，port
1 に次のガウシアンパルスを与える．

ϕ(x, z, t = 0) = ϕ1 (x)g(z − li ) exp [−jβ1 (z − li )][]g(z) = exp − (z/a)2(6.9)(6.10)

ここで，a = 2 µm とし，ϕn および βn は波長 1.55 µm における port n の固有モードフィールドお
よび位相定数である．目的関数は，port 1 に波長 1.55 µm の基本モードが入射されたとき，port 2
におけるその出力が最大化されるように，次のように設定する．
)2
(
τ |2
|η21
Minimize C = 1 −(6.11)
ideal 2
η21
ideal は port 2 における理想の出力パルス ψ を用いて次のように定義する．
ここで，η21
2
∫∫
ideal
η21
=
|ψ2 (x, z)|2 dxdz

(6.12)Ωψ2 (x, z) = ϕ2 (x)g(z − lr2 ) exp [−jβ2 (z − lr2 )](6.13)

伝搬解析には近軸近似と ADI 法に基づく 2 次元 TDFD-BPM を用いる．中心波長は 1.55 µm とし，離散点間隔は ∆x = ∆z = 0.035 µm，∆t = 1 fs と設定する．解析領域端には厚さ 1.05 µm の PMLを装荷し，PML 終端は電気壁とする．最速降下法におけるステップ幅は二分法を用いた直線探索により決定し，最適化の反復毎に感度平滑化フィルタ (3 × 3 移動平均フィルタ) を適用して最適設計を行う．密度法におけるペナルティ係数 r は初期値を 8 とし，目的関数の値が収束条件 |Cα−1 − Cα |< 10−4 (Cα : 反復 α 回目の目的関数の値) を満たした場合に r を 1.5 倍の値へ更新し，最大 1024 まで変化させて設計を行う．
図 6.2 に最適化の反復回数に対する目的関数の値のグラフを示す．図中の挿入図は初期構造であり，単純な S 字曲がり導波路を与えている．図より，ペナルティ係数 r の更新に伴う特性の劣化が生じている箇所が見られるが，それらの箇所を除けば，目的関数の値が最小化するように最適化が進んでおり，感度解析手法の妥当性が確認できる．最適化の反復が終了した反復 103 回目における最適化結果を図 6.3 に示す．図 6.3(a) は port 2 における正規化透過パワースペクトルである．図では，比較のため，FEM を用いた密度法に基づくトポロジー最適化構造および入出力部に導波路
オフセットを付加した S 字曲がり導波路の特性も併せて示している．オフセット付き曲がり導波路の特性と比較すると，本設計アプローチにより得られた最適化構造は損失の大幅な改善が達成されている．本設計アプローチと FEM を用いた場合の最適化構造における特性を比較すると，ほぼ同程度の特性となっているが，TD-BPM を用いた設計の方が比較的波長平坦な特性が得られている．図 6.3(b) に示す本設計アプローチにより得られた最適化構造を見ると，曲がり導波路に固定の周期を持たないグレーティング構造が現れている．グレーティング導波路において，グレー
ティング周期を変調することによって波長平坦な反射特性が得られ，サイドローブレベルも低減
することが知られている．図 6.3(b) の構造は中心波長から離れた波長帯域において反射スペクト
ルが平坦になるようにグレーティング周期が変調され，結果，比較的波長平坦な特性特性が得られたのではないかと考えられる．

65

y

z

Reference plane 1

Reference plane 2

PMLx

Input pulsedw
port 1

Design region

n1

n1 or n2LxDw port 2

n2
d

Lz

ld

li
lr1

lr2
Wz
図 6.1: 強導波路における曲がり導波路の 2 次元近似設計モデル

表 6.1: 強導波路における曲がり導波路の 2 次元近似設計モデルにおける構造パラメータ
n1
n2
w
Lz
Lx
d

2.2

1.445

0.5 µm5 µm

5 µm

1.5 µm

D
4 µm

Wz
45 µm

li
5 µm

lr1
10 µm

ld
10 µm

lr2
35.525 µm

66

0.7
n1

Objective function

0.6
0.5
0.4
0.3
Initial structure

n2

0.2
0.1
0.0

0

20

40
60
80
Number of iterations

100

120

図 6.2: 最適化の反復回数に対する強導波路における曲がり導波路の出力パワー

67

Normalized output power

1.00
0.95
0.90

Optimized structure (TD−BPM)
Optimized structure (FEM)
Bending waveguide with offsets

0.85
0.80
1.53 1.54 1.55 1.56 1.57 1.58 1.59 1.6 1.61 1.62 1.63
Wavelength [ µ m](a) 正規化出力パワースペクトル(b) 最適化構造 (TD-BPM を用いた場合)(c) 最適化構造 (FEM を用いた場合)(d) オフセット付き曲がり導波路
図 6.3: 強導波路における曲がり導波路の 2 次元近似設計におけるトポロジー最適設計結果

68

6.3

導波路型反射器

本設計アプローチのさらなる有用性を確認するため，図 6.4 に示す設計モデルを用いて導波路型
反射器の設計を行う．設計モデルにおける構造パラメータは表 6.2 にまとめて示す．中心波長，刻み幅，PML 厚さ，最速降下法におけるステップ幅の決定方法および目的関数は前節の曲がり導波路の設定と同様とし，port 1 から入射される TE 基本モードの光波が port 2 へ反射されるような反射器の実現を目的とする．同じく，入射パルスには式 (6.9) のガウシアンパルスを与える．
初期構造に ρ = 0.498 の一様媒質を与えてトポロジー最適設計を行った結果を図 6.5 に示す．図6.5(a) は port 2 における正規化反射パワースペクトルである．図では比較のため，FEM を用いた密度法に基づく従来設計手法による最適化構造の特性も併せて示している．図より，FEM を用いた場合の最適化構造では中心波長のみを考慮した設計であるため，中心波長をピークとする波長依存性の強い反射スペクトルが得られている．本アプローチはガウシアンパルスの伝搬解析により比較的広い波長帯域の透過特性が考慮されているため，本アプローチを用いて得られた最適化
構造では比較して非常に波長平坦な反射スペクトルが得られている．そのため，本設計アプロー
チのように，時間領域手法を用いたトポロジー最適設計は波長平坦な特性など，透過・反射スペ
クトル自体を最適化したい場合に比較的便宜な手法であると考えられる．

69

yReference plane 2
Reference plane 1

z
PML

x

d

Input pulsew
port 1

Design region

n1

n1 or n2w
port 2Lxdn2li(= lr2)
lr1

D

Lz

ld
Wz

図 6.4: 強導波路における導波路型反射器の 2 次元近似設計モデル

表 6.2: 強導波路における導波路型反射器の 2 次元近似設計モデルにおける構造パラメータ
n1
n2
w
Lz
Lx
d

2.2

1.445

0.5 µm5 µm

5 µm

2.5 µmD
2 µm

Wz
25.9 µm

li
5 µm

lr1
10 µm

ld
10 µm

lr2
5 µm

70

Normalized output power

1.00
0.95
0.90
0.85
0.80
0.75

with TD−BPM
with FEM

0.70
1.53 1.54 1.55 1.56 1.57 1.58 1.59 1.6 1.61 1.62 1.63
Wavelength [ µ m](a) 正規化反射パワースペクトル(b) 最適化構造 (TD-BPM を用いた場合)(c) 最適化構造 (FEM を用いた場合)

図 6.5: 強導波路における導波路型反射器の 2 次元近似設計におけるトポロジー最適化結果

71

第7章

フルベクトルビーム伝搬法を用いたトポロ
ジー最適化

本章では，フルベクトルビーム伝搬法を活用した，コアとクラッドの屈折率差が比較的大きい
強導波路デバイスのトポロジー最適設計について検討を行う．
シリコンフォトニクスにおける偏波無依存化やデジタルコヒーレント通信において偏波回転・分離は必須の要素であり，これまでに方向性結合や断熱変換などを原理とする種々の偏波回転素子や偏波分離素子が提案されている [86]〜[97]．偏波回転素子に代表される偏波の結合が生じる素子の数値解析にはフルベクトル解析が必要であるが，特に断熱変換を原理とする素子は比較的長い素子長が必要になるため，現在の計算機性能においても FDTD 法や有限要素法計算等を利用した場
合は計算コストが問題となる．ビーム伝搬法は計算効率が非常に高い計算手法であるものの，通常のビーム伝搬法は後方散乱波を考慮できないため，汎用的に強導波路デバイスの解析に適用す
ることは困難である．しかし，伝搬方向に十分構造が緩やかに変化する場合は強導波路においても適用可能であり，実際に偏波回転素子の設計に適用した例も報告されている [76]．方向性結合や導波モードの断熱変換を原理とする素子は光波の伝搬方向に対して導波路構造が緩やかに変化するため，強導波路の最適設計においてもビーム伝搬法の活用が可能であると考えられる．本章
では，はじめにフルベクトルビーム伝搬法を用いた構造最適化について感度計算方法を示し，[89]で提案された 2 層構造を有する偏波回転素子をコンセプトとした設計モデルを用いて本手法の利用可能性を検討する．

7.1

感度解析手法

ここでは，設計領域における屈折率分布は x，z 方向のみに構造の変位がある場合を考え，離散点ごとに割り当てられる座標 (i∆x, k∆z) における正規化密度パラメータを ρi,k とする．フルベク
トル問題の場合，1 番目のポートから n 番目の出力ポートへの S パラメータは式 (3.2) で表現され，以下に再掲する．
()
Sn1 = ∆x∆y {ψy,n }† {ϕx }Nz −1 − {ψx,n }† {ϕy }Nz −1 ≡ {gn }† {ϕt }Nz −1[]{ψy,n }
{gn } = ∆x∆y
− {ψx,n }][{ϕx }Nz −1{ϕt }Nz −1 ={ϕy }Nz −1

72

ADI 法に基づくセミベクトル BPM の場合と同様の導出手順で，正規化密度パラメータに対する感度が次のように求められる.
{}T ∂ [Γ4 ]k{}T ∂ [Γ3 ]k∂Sn1= − λn,k+1/2{ϕt }k+1 + λn,k+1/2{ϕt }k+ 1
2
∂ρi,k
∂ρi,k
∂ρi,k
∂ [Γ2 ]k∂ [Γ1 ]k
− {λn,k }T{ϕt }k+ 1 + {λn,k }T{ϕt }k
2
∂ρi,k
∂ρi,k
{
}
ここで，随伴変数ベクトル λn,k+1/2 ，{λn,k } は次のように定義している．
(){}T= {gn }† [Γ4 ]−1 [Γ3 ] [Γ2 ]−1 [Γ1 ]
λn,k+1/2{{λn,k }T = λn,k+1/2

}T

Nz −2

(
)
· · · [Γ4 ]−1 [Γ3 ] [Γ2 ]−1 [Γ1 ]

(7.1)

k+1

[Γ4 ]−1k(7.2)[Γ3 ]k [Γ2 ]−1
k

(7.3)

随伴変数ベクトルは次の随伴変数方程式を解くことにより算出できる.
)
{
} (
T
T
[Γ4 ]Tk λn,k+1/2 = [Γ1 ]T [Γ2 ]−1 [Γ3 ]T [Γ4 ]−1
···
k+1
(
)
T
T
· · · [Γ1 ]T [Γ2 ]−1 [Γ3 ]T [Γ4 ]−1{gn }∗
Nz −2{}[Γ2 ]Tk {λn,k } = [Γ3 ]Tk λn,k+1/2

(7.4)
(7.5)

上記の式 (7.1)〜(7.5) は，ADI 法に基づくセミベクトル BPM を用いた場合の感度計算式と同じ形式であり，行列 [Γα ](α ∈ {1, 2, 3, 4}) と電磁界振幅ベクトルが式 (2.144)〜(2.148) に示した，フルベクトル解析に対応する形に入れ替わっている．なお，本研究では順方向および感度解析のための逆方向伝搬解析では，伝搬解析の精度を考慮して，ともに磁界を未知変数として伝搬解析を行う．

73

7.2

偏波回転素子のトポロジー最適化

図 7.1 に示す設計モデルを用いて，前節で示した感度計算方法の妥当性を確かめる．設計モデ
ルに示す材料定数および寸法を表 7.1 に示す．密度法に基づくトポロジー最適化により，偏波回転素子の上部層を設計領域として偏波回転素子の最適化を行う．目的関数は，port 1 に波長 1.55 µmの TM-like 基本モードが入射されたとき，port 2 の TE-like 基本モード出力が最大化されるように，
次のように設定する．
()2
Minimize C = 1 − |S21,TM→TE |2(7.6)
ここで，|S21,TM→TE |2 は port 1 に入射した TM-like モードが port 2 から TE-like モードに変換されて出力される正規化出力パワーを意味している．伝搬解析には近軸近似と ADI 法に基づくフルベクトル FD-BPM を用いる．参照屈折率 n0 は入射モードの実効屈折率に固定し，離散点間隔は∆x = ∆y = 0.01 µm，∆z = 0.02 µm と設定する．解析領域端には厚さ 0.5 µm の PML を装荷し，PML 終端は電気壁あるいは磁気壁とする．構造単純化のため，最適化の反復毎に 3x3 単純移動平均フィルタによる感度平滑化フィルタを適用し，さらに z = (L + l)/2 を対称面とする対称条件を∂Cとし，二分課して設計を行う．最速降下法におけるステップ幅は最大値を K = 0.1/ maxi,j ∂ρi,j法による直線探索により移動量を決定する．初期構造は長手方向の長さが 5 µm のテーパが付加された図 7.2(a) に示す構造 (コア部:ρ = 0.8，クラッド部 ρ = 0.4) を与え，密度法におけるペナル
ティ係数 r は 8 から 1024 まで変化させて設計を行う．初期構造の屈折率分布をコア，クラッド材
料に 2 値化した場合，その正規化出力パワーは TM-like 基本モードが 0.387，TE-like 基本モード
が 0.624 であり，初期構造の段階では十分な偏波回転機能が実現されていない．
図 7.3 に最適化の反復回数に対する正規化出力パワーのグラフを示す．図より，反復毎に port 2からの TE-like モード出力パワーが増加しており，感度解析手法の妥当性が確認できる．図 7.4 に最適化構造とその伝搬波形を示す．伝搬界分布を見ると，TM-like 基本モードが TE-like 基本モードに変換されている様子が確認できる．
ここで用いた偏波回転素子の設計モデルは前述のとおり [89] で提案されている動作原理を想定している．コア断面形状が x 軸と y 軸の構造対称性を持たない場合，x 偏波成分と y 偏波を併せ持つ 2 つのハイブリッドモードが励起し，その 2 つのモードの相互作用により偏波回転が実現される．2 つのハイブリッドモードの実効屈折率を高い方からそれぞれ ne,1 ，ne,2 とすると，偏波変換に必要な素子長 Lπ は次式で求められる．

Lπ =λ2(ne,1 − ne2 )(7.7)

式 (7.7) に示すように，2 つのハイブリッドモードの実効屈折率差が大きいほど短い素子長で偏波回転素子が実現できる．図 7.4(a) に示す最適化構造は，コア断面形状において励起される 2 つのハイブリッドモードの実効屈折率差が大きくなるように上部層の構造が最適化されてたのではな
いかと考えられる．最適化構造における正規化出力パワーは TM-like 基本モードが 0.064，TE-like
基本モードが 0.953 であり，正規化パワーの和が 1.017 と 1 を越えており，パワー保存則が満たされていない．これは最適化構造がビーム伝搬法の近似制約を越えた構造であるためと考えられる．
十分な精度の最適設計のためには，ビーム伝搬法で取り扱える構造のみが出現するような構造制
約条件の考案が課題である．

74

l

z
x

w

port 2

L
l

y

h
D

port 1h2

n1hDesign Regionh1(n1 or n2 )n2w

図 7.1: 偏波回転素子の最適設計モデル

n1
2.2

表 7.1: 偏波回転素子の設計モデルにおける構造パラメータ
n2
w
h
l
h1
h2
D

1.445

0.8 µm

0.8 µm1 µm

0.4 µm

0.4 µm

3 µm

L
20 µm

(a) 初期構造(b) 2 値化した初期構造(c) 伝搬波形 (|Hx |)(d) 伝搬波形 (|Hy |)図 7.2: 偏波回転素子の最適設計における初期構造と 2 値化した初期構造における伝搬波形

75Normalized power in port 2

1.00
0.95
0.90
0.85
0.80
0.75
0.70
0.65

0

5

10
15
20
25
Number of iterations30

図 7.3: 最適化の反復回数に対する TE-like モードの正規化出力パワー(a) 最適化構造(c) 伝搬波形 (|Hy |)(b) 伝搬波形 (|Hx |)図 7.4: 偏波回転素子のトポロジー最適化構造とその伝搬波形

76

35

第8章

結論

光導波路の解析手法として広く利用されている計算効率の高いビーム伝搬法を活用し，感度解析手法として随伴変数法を用いる，光導波路素子のための密度法に基づくトポロジー最適設計ア
プローチを新たに提案した．
はじめに 4 章では，2 次元ビーム伝搬法を用いる 2 次元近似設計問題の場合に対して，密度法における設計変数である正規化密度パラメータに対する正規化出力パワーの感度解析方法の定式化を行った．コアとクラッドの屈折率が小さい弱導波路における曲がり導波路の設計に本設計アプ
ローチを適用し，最適化の反復回数を経る毎に特性のほぼ単調な改善を確認できたため，感度解析手法の妥当性が確かめられた．さらに各種パワー分岐素子 (分岐比 1:1，3:2，2:1，3:1) の設計問
題を通して，本設計手法の汎用性を確かめた．次に 5 章にて，平面光波回路における面外方向の放射を考慮する，より実際的な 3 次元設計問題のために，交互方向陰解法に基づくセミベクトル
ビーム伝搬法を用いた場合の感度計算方法を定式化した．本設計手法を弱導波路における曲がり
導波路の 3 次元設計問題に適用し，未知電磁界に依らず最適化が妥当に実行されることを確認した．さらにパワー等分岐，モード次数変換素子の設計問題を通して，本手法のさらなる汎用性を示した．これらの研究成果によりトポロジー最適設計の効率的な利用が可能になったため，現行
あるいは今後もさらなる発展が予想される石英系光導波路デバイスに対して，各種光回路の基本
要素の小型化・低損失化が期待される．
通常のビーム伝搬法は後方散乱波を考慮しない本質的に近似的な解析手法であるため，コアとクラッドの屈折率差が非常に大きな，シリコンフォトニクスデバイスを扱う場合は上記設計アプ
ローチの適用が困難であった．そこで 6 章では，後方散乱を考慮できる時間領域ビーム伝搬法を用いた場合の感度計算方法を新たに提案し，導波路曲がりや反射器の設計問題を通して本設計手法の妥当性・有効性を確かめた．また，時間領域ビーム伝搬法を活用したトポロジー最適化は広帯域動作を目的としたデバイス等，スペクトル自体の最適化を行う場合に比較的便宜である可能性を示した．同じ時間領域手法である，現在広く利用されている陽解法に基づく FDTD 法と比較すると，時間領域ビーム伝搬法は陰的解法に基づくため，時間方向に対する刻み幅を比較して大きくとることができ，2 次元問題では比較的計算コスト面で優位であることが示されている [38]．
しかし，3 次元問題ではその優位性が明らかではなく，そのためセミベクトルあるいはフルベクト
ル時間領域ビーム伝搬法を活用して強導波路デバイスのトポロジー最適設計問題が効率化されるかは明確ではない．そのため 7 章では，3 次元問題における偏波回転素子などの偏波結合を考慮する必要がある設計問題のために，通常の後方散乱を無視する交互方向陰解法に基づくフルベク
トルビーム伝搬法を用いた場合の感度計算方法を示し，本設計アプローチの利用可能性を検討した．偏波回転素子の設計問題に本設計アプローチを適用し，偏波変換効率が最大化するように最適化が進むことを確認し，感度計算の妥当性を示した．方向性結合器を原理とするデバイスや導波モードの断熱変換を原理とするデバイスに対しては，後方散乱を無視する通常のビーム伝搬法が活用可能と考えられるが，トポロジー最適化では設計領域に任意の構造が出現する問題がある．
以上の結果を受けて，本研究における今後の課題として次が挙げられる・77• 導波モードが断熱的に変化するような導波路構造のみが出現するような構造制約条件の定式化• 双方向ビーム伝搬法 [47]〜[56] の活用
導波モードの断熱変換を原理とする光導波路デバイスは総じて比較的高い製造トレランスを有するが，素子長が長く，小型化が課題となっている．導波モードが断熱的に変化するような導波路構造のみが出現する構造制約条件が定式化できれば，通常のビーム伝搬法を活用したトポロジー
最適化の活用が可能であると考えられ，断熱変化を原理とする光導波路デバイスの小型化・低損失化の貢献が期待できる．また，さらに広汎の問題に対する計算効率の高いトポロジー最適化を実現する手段として，双方向ビーム伝搬法の活用が考えられる．双方向ビーム伝搬法は不連続面における反射・透過波を算出するプロセスと伝搬方向に構造が一様な領域の伝搬プロセスを交互に行う，後方散乱が考慮可能な伝搬解析手法であり，強導波路デバイスの設計問題に広く適用可能である．双方向ビーム伝搬法を利用したトポロジー最適化手法の有用性について，手法の計算効率の検討を含め今後の課題とする．

78

参考文献
[1] M. P. Bendsøe and O. Sigmund, “Topology Optimization: Theory, Methods, and Applications,”
Springer-Verlag Berlin Heidelberg New York, 2003.
[2] 西脇眞二, 泉井一浩, 菊池昇, “トポロジー最適化,” 丸善出版, 2013.
[3] J. S. Jensen and O. Sigmund, “Systematic design of photonic crystal structures using topology
optimization: Low-loss waveguide bends,” Appl. Phys. Lett., Vol. 84, No. 12, pp. 2022-2024, Mar.
2004.
[4] P. I. Borel, A. Harpøth, L. H. Frandsen, M. Kristensen, P. Shi, J. S. Jensen, and O. Sigmund,“Topology optimization and fabrication of photonic crystal structures,” Opt. Express, Vol. 12, No.
9, pp. 1996-2001, May 2004.
[5] L. H. Frandsen, A. Harpøth, P. I. Borel, M. Kristensen, J. S. Jensen, and O. Sigmund, “Broadband
photonic crystal waveguide 60◦ bend obtained utilizing topology optimization,” Opt. Express, Vol.
12, No. 24, pp. 5916-5921, Nov. 2004.
[6] P.I. Borel, L.H. Frandsen, A. Harpøth, M. Kristensen, J. S. Jensen, and O. Sigmund, “Topology
optimised broadband photonic crystal Y-splitter,” Electron. Lett., Vol. 41, No. 2, pp. 69-71, Jan.
2005.
[7] J. S. Jensen, O. Sigmund, L. H. Frandsen, P. I. Borel, A. Harpøth, and M. Kristensen, “Topology
design and fabrication of an eﬃcient double 90◦ photonic crystal waveguide bend,” IEEE Photon.
Technol. Lett., Vol. 17, No. 6, pp. 1202-1204, June 2005.
[8] J. S. Jensen and O. Sigmund, “Topology optimization of photonic crystal structures: a highbandwidth low-loss T-junction waveguide,” J. Opt. Soc. Am. B, Vol. 22, No. 6, pp. 1191-1198,
June 2005.
[9] L. H. Frandsen, P. I. Borel, J. S. Jensen, and O. Sigmund, “Topology optimized photonic wire
splitters,” In CLEO/QELS 2006 Technical Digest CD-Rom., pp. 1-2, 2006.
[10] Y. Tsuji, K. Hirayama, T. Nomura, K. Sato, and S. Nishiwaki, “Design of optical circuit devices
based on topology optimization,” IEEE Photon. Technol. Lett., Vol. 18, No. 7, pp. 850-852, Apr.
2006.
[11] Y. Watanabe, N. Ikeda, Y. Sugimoto, Y. Takata, Y. Kitagawa, A. Mizutani, N. Ozaki, and K.
Asakawa, “Topology optimization of waveguide bends with wide, flat bandwidth in air-bridgetype photonic crystal slabs,”J. Appl. Physics, Vol. 101, No. 11, pp. 850-852, June 2007.

79

[12] R. Stainko and O. Sigmund, “Tailoring dispersion properties of photonic crystal waveguides by
topology optimization,” Waves in Random and Complex Media, Vol. 17, No. 4, pp. 477-489, Nov.
2007.
[13] Y. Tsuji and K. Hirayama, “Design of optical circuit devices using topology optimization methodwith function-expansion-based refractive index distribution,” IEEE Photon. Technol. Lett., Vol. 20,
No. 12, pp. 982-984, June 2008.
[14] R. Matzen, J. S. Jensen, and O. Sigmund, “Systematic design of slow-light photonic waveguides,”J. Opt. Soc. Amer. B, Vol. 28, No. 10, pp. 2374-2382, Oct. 2011.
[15] K. Fujimoto, Y. Tsuji, K. Hirayama, T. Yasui, S. Sato, and R. Kijima, “A study on topology
optimization of optical circuits consisting of multi-materials,” J. Lightw. Technol., Vol. 30, No. 13,
pp. 2210-2215, July 2012.
[16] T. Yasui, Y. Tsuji, J. Sugisaka, and K. Hirayama, “Design of three-dimensional optical circuit
devices by using topology optimization method with function-expansion-based refractive index
distribution,” J. Lightw. Technol., Vol. 31, No. 23, pp. 3765-3770, Dec. 2013.
[17] H. Goto, Y. Tsuji, T. Yasui, and K. Hirayama, “A study on optimization of waveguide dispersion
property using function expansion based topology optimization method,” IEICE Trans. Electron.,Vol. E97-C, No. 7, pp. 670-676, July 2014.
[18] Z. Zhang, Y. Tsuji, T. Yasui, and K. Hirayama, “Design of ultra-compact triplexer with functionexpansion based topology optimization,” Opt. Express, Vol. 23, No. 4, pp. 3937-3950, Feb. 2015.
[19] L.H. Frandsen, Y. Elesin, O. Sigmund, J.S. Jensen, K. Yvind, “Wavelength selective 3D topology
optimized photonic crystal devices,” in Tech. Dig. CLEO, 2013, CTh4L.6.
[20] L. H. Frandsen, Y. Elesin, L. F. Frellsen, M. Mitrovic, Y. Ding, O. Sigmund, and K. Yvind,“Topology optimized mode conversion in a photonic crystal waveguide fabricated in silicon-oninsulator material,” Opt. Express, Vol. 22, No. 7, pp. 8525-8532, Apr. 2014.
[21] L. H. Frandsen, Y. Elesin, X. Guan, O. Sigmund, and K. Yvind, “Flat-top drop filter based on asingle topology optimized photonic crystal cavity,” in Tech. Dig. CLEO, 2015, SF2H.7.
[22] A. Y. Piggott, J. Lu, K. G. Lagoudakis, J. Petykiewicz, T. M. Babinec, and J. Vuc̆ković, “Inverse
design and demonstration of a compact and broadband on-chip wavelength demultiplexer,” Nature
Photon., Vol. 9, pp. 374-377, May 2015.
[23] J. Lu and J. Vuc̆ković, “Nanophotonic computational design,” Opt. Express., Vol. 21, No. 11, pp.
13351-13367, June 2013.
[24] M. D. Feit and J. A. Fleck, “Light propagation in graded-index optical fibers,” Appl. Opt., Vol. 17,
No. 24, pp. 3990-3998, Dec. 1978.
[25] Y. Chung and N. Dagli, “An assessment of finite diﬀerence beam propagation method,” IEEE J.
Quantum. Electron., Vol. 26, No. 8, pp. 1335-1339, Aug. 1990.
80

[26] J. Yamauchi, T. Ando, H. Nakano, “Beam-propagation analysis of optical fibers by alternating
direction implicit method,” Electron. Lett., Vol. 27, No. 18, pp. 1663-1665, Aug. 1991.
[27] P. L. Liu and B. J. Li, “Study of form birefringence in waveguide devices using the semivectorial
beam propagation method,” IEEE Photon. Technol. Lett., Vol. 3, No. 10, pp. 913-915, Oct. 1991.
[28] P. L. Liu and B. J. Li, “Semivectorial beam-propagation method for analyzing polarized modes of
rib waveguides,” IEEE J. Quantum. Electron., Vol. 28, No. 4, pp. 778-782, Apr. 1992.
[29] W. P. Huang and C. L. Xu, “Simulation of three-dimensional optical waveguides by a full-vector
beam propagation method,” IEEE J. Quantum. Electron., Vol. 29, No. 10, pp. 2639-2649, Oct.
1993.
[30] J. Yamauchi, G. Takahashi, and H. Nakano, “Full-vectorial beam-propagation method based on
the McKee-Mitchell scheme with improved finite-diﬀerence formulas,” J. Ligthw. Technol., Vol.
16, No. 12, pp. 2458-2464, Dec. 1998.
[31] Y. Hsueh, M. Yang, and H. Chang, “Three-dimensional noniterative full-vectorial beam propagation method based on the alternating direction implicit method,” J. Ligthw. Technol., Vol. 17, No.
11, pp. 2389-2397, Nov. 1999.
[32] Y. Tsuji and M. Koshiba, “A finite element beam propagation method for strongly guiding andlongitudinally varying optical waviguides,” J. Ligthw. Technol., Vol. 14, No. 2, pp. 217-222, Feb.
1996.
[33] Y. Tsuji, M. Koshiba, and T. Shiraishi, “Finite element beam propagation method for threedimensional optical waveguide structures,” J. Lightw. Technol., Vol. 15, No. 9, pp. 1728-1734,
Sept. 1997.
[34] P-L. Liu, Q. Zhao, F-S. Choa, “Slow-wave finite-diﬀerence beam propagation method,” IEEE
Photon. Technol. Lett., Vol. 7, No. 8, pp. 890-892, Aug. 1995.
[35] G. H. Jin, J. Harari, J.P. Vilcot, D. Decoster, “An improved time-domain beam propagation method
for integrated optics components,” IEEE Photon. Technol. Lett., Vol. 9, No. 3, pp. 348-350, Mar.
1997.
[36] M. Koshiba, Y. Tsuji, M. Hikari, “Time-domain beam propagation method and its application tophotonic crystal circuits,” J. Lightw. Technol., Vol. 18, No. 1, pp. 102-110, Jan. 2000.
[37] J. Shibayama, A. Yamahira, T. Mugita, J. Yamauchi, H. Nakano, “A finite-diﬀerence time-domain
beam-propagation method for TE- and TM-wave analyses,” J. Lightw. Technol., Vol. 21, No. 7, pp.
1709-1715, Aug. 2003.
[38] J. Shibayama, M. Muraki, J. Yamauchi, H. Nakano, “Comparative study of several time-domain
methods for optical waveguide analyses,” J. Lightw. Technol., Vol. 23, No. 7, pp. 2285-2293, July
2005.

81

[39] T. Hashimoto, T. Saida, I. Ogawa, M. Kohtoku, T. Shibata, and H. Takahashi, “Optical circuitdesign based on a wavefront-matching method,” Opt. Lett., Vol. 30, No. 19, pp. 2620-2622, Oct.
2005.
[40] Y. Sakamaki, T. Saida, T. Hashimoto, and H. Takahashi, “New optical waveguide design based on
wavefront matching method,” J. Lightw. Technol., Vol. 25, No. 11, pp. 3511-3518, Nov. 2007.
[41] Y. Sakamaki, T. Saida, T. Shibata, Y. Hida, T. Hashimoto, M. Tamura, and H. Takahashi, “Y-branchwaveguides with stabilized splitting ratio designed by wavefront matching method,” IEEE Photon.
Technol. Lett., Vol. 18, No. 7, pp. 817-819, Apr. 2006.
[42] Y. Sakamaki, T. Saida, T. Hashimoto, S. Kamei, and H. Takahashi, “Loss reduction of waveguidecrossings by wavefront matching method and their application to integrated optical circuits,” J.
Lightw. Technol., Vol. 27, No. 13, pp. 2257-2263, July 2009.
[43] Y. Sakamaki, S. Kamei, T. Hashimoto, T. Kitoh, and H. Takahashi, “Loss uniformity improvementof arrayed-waveguide grating with mode-field converters designed by wavefront matching method,”
J. Lightw. Technol., Vol. 27, No. 24, pp. 5170-5715, Dec. 2009.
[44] K. Svanberg, “The method of moving asymptotes–a new method for structural optimization,” Int.
J. Numer. Meth. Eng., Vol. 24, No. 2, pp. 359-373, Feb. 1987.
[45] M. A. Swillam, M. H. Bakr, and X. Li, “Eﬃcient adjoint sensitivity analysis exploiting the FDBPM,” J. Ligthw. Technol., Vol. 25, No. 7, pp. 1861-1869, July 2007.
[46] M. A. Swillam, M. H. Bakr, and X. Li, “Full vectorial 3-D sensitivity analysis and design optimization using BPM,” J. Lightw. Technol., Vol. 26, No. 5, pp. 528-536, Mar. 2008.
[47] P. Kaczmarski and P.E. Lagasse, “Bidirectional beam propagation method,” Electron. Lett., Vol.
24, No. 11, pp. 675-676, May 1988.
[48] H. Rao, R. Scarmozzino, R. M. Osgood, “A bidirectional beam propagation method for multiple
dielectric interfaces,” IEEE Photon. Technol. Lett., Vol. 11, No. 7, pp. 830-832, July 1999.
[49] H. El-Refaei, D. Yevick, I. Betty, “Stable and noniterative bidirectional beam propagation method,”IEEE Photon. Technol. Lett., Vol. 12, No. 4, pp. 389-391, Apr. 2000.
[50] P. L. Ho and Y. Y. Lu, “A bidirectional beam propagation method for periodic waveguides,” IEEE
Photon. Technol. Lett., Vol. 14, No. 3, pp. 325-327, Mar. 2002.
[51] Y. Y. Lu, S. H. Wei, “A new iterative bidirectional beam propagation method,” IEEE Photon.
Technol. Lett., Vol. 14, No. 11, pp. 1533-1535, Dec. 2002.
[52] H. Zhang, J. Mu, W.-P. Huang, “Improved bidirectional beam-propagation method by a fourthorder finite-diﬀerence scheme,” J. Lightw. Technol., Vol. 25, No. 9, pp. 2807-2813, Sept. 2007.
[53] H. Zhang, J. Mu, W.-P. Huang, “Assessment of rational approximations for square root operator
in bidirectional beam propagation method,” J. Lightw. Technol., Vol. 26, No. 5, pp. 600-607, Mar.
2008.
82

[54] H. A. Jamid and M. Z. M. Khan, “A numerical approach for full-vectorial analysis of 3-D guided
wave structures with multiple and strong longitudinal discontinuities,” IEEE J. Quantum. Electron.,Vol. 45, No. 2, pp. 117-124, Feb. 2009.
[55] J. Xiao, S. Wu, X. Sun, “A stable and accurate preconditioner for bidirectional beam propagation
method,” IEEE J. Quantum. Electron., Vol. 51, No. 4, 6100107 , Apr. 2015.
[56] S. Wu, J. Xiao, “An eﬃcient semivectorial bidirectional beam propagation method for 3-D opticalwaveguide structures,” J. Lightw. Technol., Vol. 34, No. 4, pp. 1313-1321, Feb. 2016.
[57] G. R. Hadley, “Wide-angle beam propagation using Padé approximant operators,” Opt. Lett.,
Vol.17, No. 20, pp. 1426-1428, Oct. 1992.
[58] J. Shibayama, T. Takahashi, J. Yamauchi, H. Nakano, “A three-dimensional horizontally wideangle noniterative beam-propagation method based on the alternating-direction implicit scheme,”
J. Lightw. Technol., Vol. 18, No. 5, pp. 661-663, Mar. 2006.
[59] J. Yamauchi, M. Sekiguchi, O. Uchiyama, J. Shibayama, and H. Nakano, “Modified finitediﬀerence formula for the analysis of semivectorial modes in step-index optical waveguides,”IEEE Photon. Technol. Lett., Vol. 9, No. 7, pp. 961-963, July 1997.
[60] J. Yamauchi, G. Takahashi, and H. Nakano, “Modified finite-diﬀerence formula for semivectorial
H-field solutions of optical waveguides,” IEEE Photon. Technol. Lett., Vol. 10, No. 8, pp. 11271129, Aug. 1998.
[61] J. Yamauchi, M. Toshiyuki, and H. Nakano, “Semivectorial H-field analysis of rib waveguides bya modified beam-propagation method based on the generalized Douglas scheme,” Opt. Lett., Vol.
25, No. 24, pp. 1771-1773, Dec. 2000.
[62] J. P. Bérenger, “A perfectly matched layer for the absorption of electromagnetic waves,” J. Comput.
Phys., Vol. 114, No. 2, pp. 185-200, Oct. 1994.
[63] W. C. Chew and W. H. Weedon, “A 3-D perfectly matched medium from modified Maxwell’sequations with stretched coordinates,” Micro. Opt. Technol. Lett., Vol. 7, No. 13, pp. 599-604,
Sept. 1994.
[64] W. P. Huang, C. L. Xu, W. Lui, and K. Yokoyama, “The Perfectly Matched Layer (PML) Boundary
Condition for the Beam Propagation Method,” IEEE Photon. Technol. Lett., Vol. 8, No. 5, pp.
649-651, May 1996.
[65] C. Vassallo and F. Collino, “Highly eﬃcient absorbing boundary conditions for the beam proDarzation method,” J. Lightw. Technol., Vol. 14, No. 6, pp. 1570-1577, June 1996.
[66] C. Vassallo, “Diﬃculty with vectorial BPM,” Electron. Lett., Vol. 33, No. 1, pp. 61-62, Jan. 1997.
[67] H. J. W. M. Hoekstra, “On beam propagation methods for modelling in integrated optics,” Opt.
and Quantum Electoron., Vol. 29, pp. 157-171, 1997.

83

[68] L. Poladian and F. Ladouceur, “Unification of TE and TM beam propagation algorithms,” IEEE
Photon. Technol. Lett., Vol. 10, No. 1, pp. 105-107, Jan. 1998.
[69] J. Yamauchi, K. Matsubara, T. Tsuda, and H. Nakano, “Norm-conserving finite-diﬀerence beampropagation method for TM wave analysis in step-index optical waveguides,” J. Lightw. Technol.,Vol. 18, No. 5, pp. 721-728, May 2000.
[70] M. S. Stern, “Semivectorial polarised finite diﬀerence method for optical waveguides with arbitraryindex prifiles,” IEE Proc. J, Vol. 135, No. 1, pp. 56-63, Feb. 1988.
[71] M. S. Stern, “Semivectorial polarised H field solutions for dielectric waveguides with arbitraryindex profiles,” IEE Proc. J, Vol. 135, No. 5, pp. 333-338, Oct. 1988.
[72] Y. He and F.G. Shi, “Improved full-vectorial beam propagation method with high accuracy forarbitrary optical waveguides,” IEEE Photon. Technol. Lett., Vol. 15, No. 10, pp. 1381-1383, Sept.
2003.
[73] J. Xiao and X. Sun, “A Modified full-vectorial finite-diﬀerence beam propagation method basedon H-fields for optical waveguides with step-index profiles,” Opt. Commun., Vol. 266, No. 2, pp.
505-511, Oct. 2006.
[74] J. Yamauchi, Y. Nito, and H. Nakano, “A modified semivectorial BPM retaining the eﬀects of thelongitudinal field component and its application to the design of a spot-size converter,” J. Lightw.
Technol., Vol. 27, No. 13, pp. 2470-2476, July 2009.
[75] Y. Nito, T. Oda, S. Takase, J. Yamauchi, and H. Nakano, “A beam-propagation method using bothelectric and magnetic fields,” IEEE Photon. Technol. Lett., Vol. 23, No. 7, pp. 429-431, Apr. 2011.
[76] Y. Nito, J. Shibayama, J. Yamauchi, and H. Nakano, “Full-vectorial beam-propagation methods
based on a fundamental scheme–design of a short polarization converter,” J. Lightw. Technol., Vol.
32, No. 21, pp. 4111-4118, Nov. 2014.
[77] D. Yevick and W. Bardyszewski, “Correspondence of variational finite-diﬀerence (relaxation)and imaginary-distance propagation methods for modal analysis,” Opt. Lett., Vol.17, No. 5, pp.
329-330, Mar. 1992.
[78] C. L. Xu, W. P. Huang, and S. K. Chaudhuri, “Eﬃcient and accurate vector mode calculations bybeam propagation method,” J. Lightw. Technol., Vol. 11, No. 7, pp. 1209-1215, July 1993.
[79] S. Jüngling and J. C. Chen, “A study and optimization of eigenmode calculations using theimaginary-distance beam-propagation method,” IEEE J. Quantum. Electron., Vol. 30, No. 9, pp.
2098-2105, Sept. 1994.
[80] J. C. Chen and S. Jüngling, “Computation of higher-order waveguide modes by imaginary-distancebeam propagation method,” Opt. and Quantum Electoron., Vol. 26, pp. S199-A205, 1994.
[81] F. Wijnands, H. J. W. M. Hoekstra, G. J. M. Krijnen, and R. M. de Ridder, “Modal fields calculation
using the finite diﬀerence beam propagation method,” J. Lightw. Technol., Vol. 12, No. 12, pp.
2066-2072, Dec. 1994.
84

[82] G. R. Hadley, “Transparent boundary condition for beam propagation,” Opt. Lett., Vol.16, No. 9,
pp. 624-626, May. 1991.
[83] G. R. Hadley, “Transparent boundary condition for the beam propagation method,” IEEE J.
Quantum. Electron., Vol. 28, No. 1, pp. 363-370, Jan. 1992.
[84] Z. Weissman, E. Marom, and A. Hardy, “Very low-loss Y-junction power divider,” Opt. Lett.,
Vol.14, No. 5, pp. 293-295, Mar. 1989.
[85] W. Chen, P. Wang, and J. Yang, “Optical mode interleaver based on the asymmetric multimode Y
junction,” IEEE Photon. Technol. Lett., Vol. 26, No. 20, pp. 2043-2046, Oct. 2014.
[86] K. Mertens, B. Opitz, R. Hövel, K. Heime, and H. J. Schmitt, “First realized polarization converter
based on hybrid supermodes,” IEEE Photon. Technol. Lett., Vol. 10, No. 3, pp. 388-390, Mar.
1998.
[87] J. Z. Huang, R. Scarmozzino, G. Nagy, M. J. Steel, and R. M. Osgood, Jr., “Realization of acompact and single-mode optical passive polarization converter,” IEEE Photon. Technol. Lett.,Vol. 12, No. 3, pp. 317-319, Mar. 2000.
[88] H. Fukuda, K. Yamada, T. Tsuchizawa, T. Watanabe, H. Shinojima, and S. Itabashi, “Polarization
rotator based on silicon wire waveguides,” Opt. Express, Vol. 16, No. 4, pp. 2628-2635, Feb. 2008.
[89] Z. Wang and D. Dai, “Ultrasmall Si-nanowire-based polarization rotator,” J. Opt. Soc. Amer. B,
Vol. 25, No. 5, pp. 747-753, May. 2008.
[90] D. Vermeulen, S. Selvaraja, P. Verheyen, P. Absil, W. Bogaerts, D. V. Thourhout, and G. Roelkens,“Silicon-on-insulator polarization rotator based on a symmetry breaking silicon overlay,” IEEE
Photon. Technol. Lett., Vol. 24, No. 6, pp. 482-484, Mar. 2012.
[91] M. Aamer, A. M. Gutierrez, A. Brimont, D. Vermeulen, G. Roelkens, J.-M. Fedeli, A. Håkansson,
and P. Sanchis, “CMOS compatible silicon-on-insulator polarization rotator based on symmetry
breaking of the waveguide cross section,” IEEE Photon. Technol. Lett., Vol. 24, No. 22, pp.
2031-2034, Nov. 2012.
[92] A. Xie, L. Zhou, J. Chen, and X. Li, “Eﬃcient silicon polarization rotator based on modehybridization in a double-stair waveguide,” Opt. Express, Vol. 23, No. 4, pp. 3960-3970, Feb.
2015.
[93] J. Feng, R. Akimoto, and H. Zeng, “Asymmetric silicon slot-waveguide-assisted polarizing beam
splitter,” IEEE Photon. Technol. Lett., Vol. 28, No. 12, pp. 1294-1297, June 2016.
[94] X. Sun, M. Z. Alam, J. S. Aitchison, and M. Mojahedi, “Polarization rotator based on augmented
low-index-guiding waveguide on silicon nitride/silicon-on-insulator platform,” Opt. Lett., Vol. 41,
No. 14, pp. 3229-3232, July 2016.
[95] C. Sun, Y. Yu, G. Chen, and X. Zhang, “A low crosstalk and broadband polarization rotator andsplitter based on adiabatic couplers,” IEEE Photon. Technol. Lett., Vol. 28, No. 20, pp. 2253-2256,
Oct. 2016.
85

[96] W. Jiang, “Fabrication-tolerant polarization splitter and rotator based on slanted silicon waveguides,” IEEE Photon. Technol. Lett., Vol. 30, No. 7, pp. 614-617, Apr. 2018.
[97] H.-C. Chung and S.-Y. Tseng, “Ultrashort and broadband silicon polarization splitter-rotator using
fast quasiadiabatic dynamics,” Opt. Express, Vol. 26, No. 8, pp. 9655-9665, Apr. 2018.

86

研究業績リスト
1. 論文(1) A. Iguchi, Y. Tsuji, T. Yasui, and K. Hirayama, “Topology optimization of optical waveguidedevices based on beam propagation method with sensitivity analysis,” J. Lightw. Technol., Vol. 34
No. 18, pp. 4214-4220, Sept. 2016.
(2) A. Iguchi, Y. Tsuji, T. Yasui, and K. Hirayama, “Eﬃcient topology optimization of optical waveguide devices utilizing semi-vectorial finite-diﬀerence beam propagation method,” Opt. Express,
Vol. 25, No. 3, pp. 28210-28222, Nov. 2017.
(3) S. Kawai, A. Iguchi, and Y. Tsuji, “Study on high precision and stable finite element beam propagation method based on incomplete third order hybrid edge/nodal element,” J. Lightw. Technol.,Vol. 36, No. 11, pp. 2278-2285, June 2018.
(4) A. Iguchi, Y. Tsuji, T. Yasui, and K. Hirayama, “Topology optimal design for optical waveguides
using time domain beam propagation method,” Electron. Express, Vol. 15, No. 11, #20180417,
May 2018.
(5) K. Mori, K. Morimoto, T. Tanaka, A. Iguchi, and Y. Tsuji, “Topology optimization of nonlinear
optical waveguide devices considering output signal phase,” Opt. Commun., Vol. 439, pp. 290-294,
May 2019.

2. 国際会議(1) A. Iguchi, Y. Tsuji, T. Yasui, and K. Hirayama, “Topology optimal design of optical waveguides
in consideration of polarization dependence using BPM and AVM,” Progress In Electromagnetics
Research Symposium (PIERS 2016) in Shanghai, Aug. 2016.(中国, 上海市)
(2) A. Iguchi, Y. Tsuji, T. Yasui, and K. Hirayama, “Topology optimization using time domain beam
propagation method for design of optical waveguide devices,” Progress In Electromagnetics Research Symposium (PIERS 2017) in Singapore, Nov. 2017.(シンガポール)(3) A. Iguchi, Y. Tsuji, T. Yasui, and K. Hirayama, “Topology optimization using beam propagation
method for fabrication tolerant optical waveguide devices,” Asia-Pacific Microwave Conference(APMC2018) in Kyoto, Nov. 2018.(京都府, 京都市)

3. 研究会(1) 辻 寧英, 井口亜希人, 安井 崇, 平山浩一, “関数展開法に基づく 3 次元光導波路デバイスのトポロジー最適設計と構造単純化に関する検討,” 電子情報通信学会技術研究報告, Vol. 114,
No. 143, EST2014-38, pp. 175-179, 2014 年 7 月.(北海道, 室蘭市)87

(2) 井口亜希人, 辻 寧英, 安井 崇, 平山浩一, “関数展開法に基づくモード変換素子の 3 次元ト
ポロジー最適設計,” 電子情報通信学会技術研究報告, Vol. 115, No. 268, EST2015-3, pp. 13-18,
2015 年 5 月.(東京都, 小金井市)(3) 辻 寧英, 井口亜希人, 安井 崇, 平山浩一, “ビーム伝搬解析を用いた光導波路デバイスのトポ
ロジー最適設計,” 電子情報通信学会技術研究報告, Vol. 115, No. 213, EST2015-64, pp. 57-61,
2015 年 9 月.(沖縄県, 石垣市)
(4) 井口亜希人, 辻 寧英, 安井 崇, 平山浩一, “FD-BPM と随伴変数法に基づく感度解析を用いた 3 次元光導波路デバイスのトポロジー最適設計,” 電子情報通信学会技術研究報告, Vol. 115,
No. 434, EST2015-113, pp. 141-146, 2016 年 1 月.(神奈川県, 神戸市)
(5) 井口亜希人, 辻 寧英, 安井 崇, 平山浩一, “有限差分ビーム伝搬法を用いた強導波路デバイ
スの感度解析に基づくトポロジー最適設計,” 電子情報通信学会技術研究報告, Vol. 116, No.
56, EST2016-5, pp. 21-26, 2016 年 5 月.(東京都, 港区)

4. 全国大会(1) 井口亜希人, 辻 寧英, 安井 崇, 平山浩一, “関数展開法に基づく偏波回転素子の 3 次元トポロジー最適設計,” 電子情報通信学会総合大会, C-15-10, 2015 年 3 月.(滋賀県, 草津市)
(2) 井口亜希人, 辻 寧英, 安井 崇, 平山 浩一, “フルベクトル有限差分ビーム伝搬法を利用した
光導波路素子のトポロジー最適設計,” 電子情報通信学会総合大会, C-15-23, 2017 年 3 月.(愛知県, 名古屋市)(3) 井口亜希人, 辻 寧英, 安井 崇, 平山浩一, “ビーム伝搬法を活用した関数展開法に基づく光導波路素子のトポロジー最適設計” 電子情報通信学会総合大会, CS-2-5, 2018 年 3 月.(東京都,足立区)

5. 北海道支部大会・インターネットシンポジウム(1) 井口亜希人, 辻 寧英, 安井 崇, 平山浩一, “関数展開法に基づくトポロジー最適化を用いた
3 次元光導波路デバイス設計の効率化に関する検討,” 平成 26 年度電気・情報関係学会北海道支部連合大会, 102, 2014 年 10 月．
(2) 井口亜希人, 辻 寧英, 安井 崇, 平山浩一, “多層構造を有する光回路デバイスの関数展開法に基づく 3 次元トポロジー最適設計,” IEICE 北海道支部学生会インターネットシンポジウム,
2015 年 2 月.
(3) 井口亜希人, 辻 寧英, 安井 崇, 平山浩一, “ビーム伝搬法と随伴変数法による感度解析を用いたモードインターリーバのトポロジー最適設計,” 平成 27 年度電気・情報関係学会北海道支部連合大会, 145, 2015 年 11 月.
(4) 井口亜希人, 辻 寧英, 安井 崇, 平山浩一, “セミベクトル 3 次元ビーム伝搬法と随伴変数法
による感度解析を用いたモード分離・変換素子のトポロジー最適設計,” 平成 28 年度電気・情報関係学会北海道支部連合大会, 160, 2016 年 11 月.

88

謝

辞

本論文をまとめるにあたり，特に指導教員である室蘭工業大学波動エレクトロニクス研究室の辻寧英教授には，研究課題のきっかけを賜り，また終始研究に関する御指導ならびに御討論，論文についての御助言等を頂き，深く感謝致します．また，本論文の査読を御承諾いただき，学部から博士後期課程を通して本研究について有益な御助言・御指摘を頂いた長谷川弘治教授，渡邊浩太准教授に深く感謝致します．
北見工業大学工学部の平山浩一教授および安井崇准教授には，自身の学会発表内容や出版論文の質を向上させる数多くのコメントを頂きました．博士後期課程においてインターンの受け入れを御快諾頂いた古河電工株式会社研究開発本部の長谷川淳一様およびフォトニックデバイス開発部の皆様に深く感謝致します．普段の研究生活では得がたい光回路製品の研究・開発の現場および最前線事情等，自身の研究にとって大変有益な体験をさせて頂きました．また，法政大学理工学部の柴山純教授には，エレクトロニクスシミュレーション研究会を通して，有限差分ビーム伝搬法に関するいくつかの御助言を頂きました．神奈川大学工学部電気電子情報工学科の特別助教である張沢君先生や鐘正氏には，私の学部時代から博士前期課程に至るまで，国際会議の発表原
稿や論文を執筆するに当たり，多くの御援助ならびに御助言を頂きました．
波動エレクトロニクス研究室の卒業生の方々および博士前期課程の方々には日頃の研究の他，様々なことに関して共に議論をして頂き，感謝致します．最後になりますが，自身が研究に邁進できたのは両親の理解があってのことでした．
ここに改めて，以上の方々に深く感謝致します．

89

付 録A

境界条件

光波が位相定数 β で +z 方向に伝搬している場合，PML 領域を考慮した周波数領域の Maxwell
方程式 (2.1)，(2.2) およびガウスの法則を成分毎に書くと，次のようになる．

1 ∂Ez+ jβEy = −jωµ0 Hx
sy ∂y
1 ∂Ez
−jβEx −= −jωµ0 Hy
sx ∂x
1 ∂Ey
1 ∂Ex
−
= −jωµ0 Hz
sx ∂x
sy ∂y
1 ∂Hz
+ jβHy = jωε0 εr Ex
sy ∂y
1 ∂Hz
−jβHx −= jωε0 εr Ey
sx ∂x
1 ∂Hx
1 ∂Hy
−
= jωε0 εr Ez
sx ∂x
sy ∂y

(A.1)
(A.2)
(A.3)
(A.4)
(A.5)
(A.6)

1 ∂ (εr Ex )1 ∂ (εr Ey )+− jβεr Ez = 0
sx
∂x
sy
∂y
1 ∂Hy
1 ∂Hx
+
− jβHz = 0
sx ∂x
sy ∂y

(A.7)
(A.8)

媒質間の電磁界の境界条件は，Stern の差分式 [70],[71] および 2 次以降の高精度差分式 [59]〜[61]の導出に必要である．以降に，セミベクトル解析の場合における誘電体媒質間の電磁界の境界条件を求める．ここで，α，β を (α, β) ∈ {(x, y), (y, x)} と定義する．
図 A.1 に示す，α 軸に対して垂直な媒質境界を考える．このとき，ni と ni+1 の媒質境界において以下の境界条件がそれぞれ成り立つ．

(for Eα mode)
εr,i+1 Eα,R = εr,i Eα,L

(A.9)

Ez,R = Ez,L
1 ∂Eα,R1 ∂Eβ,R
1 ∂Eα,L1 ∂Eβ,L+=+sα,i+1 ∂αsβ,i+1 ∂βsα,i ∂αsβ,i ∂β1 ∂Eα,R
1 ∂Eα,L=sα,i+1 ∂αsα,i ∂α90(A.10)

ni−1
φi−1

ni
∆α
2

∆α
2

φi

ni+1
∆α
2

∆α
2

φi+1

α

φL φR

i−1

i

i+1

図 A.1: 差分グリッドと媒質境界(for Eβ mode)
Eβ,R = Eβ,L

(A.11)

Hz,R = Hz,L
1 ∂Eβ,R
1 ∂Eα,R1 ∂Eβ,L
1 ∂Eα,L−=−sα,i+1 ∂αsβ,i+1 ∂βsα,i ∂αsβ,i ∂β
1 ∂Eβ,R1 ∂Eβ,L=sα,i+1 ∂αsα,i ∂α(A.12)(for Hα mode)
Hα,R = Hα,L

(A.13)

Hz,R = Hz,L
1 ∂Hα,R1 ∂Hβ,R
1 ∂Hα,L1 ∂Hβ,L+=+sα,i+1 ∂αsβ,i+1 ∂βsα,i ∂αsβ,i ∂β1 ∂Hα,R
1 ∂Hα,L=sα,i+1 ∂αsα,i ∂α91(A.14)

表 A.1: θE ，θH の定義

ϕ=E

ϕ=H

T

N

T

N

θE

1

εr,i
εr,i+111θHsα,i+1sα,isα,i+1sα,isα,i+1 εr,i+1
sα,i εr,isα,i+1sα,i媒質境界に対して T: 接線成分，N: 法線成分

(for Hβ mode)
Hβ,R = Hβ,L

(A.15)

Ez,R = Ez,L
∂Hβ,R
∂Hα,R
1
1
1 ∂Hβ,L
1 ∂Hα,L−=−sα,i+1 εr,i+1 ∂αsβ,i+1 εr,i+1 ∂βsα,i εr,i ∂αsβ,i εr,i ∂β∂Hβ,R1
1 ∂Hβ,L=sα,i+1 εr,i+1 ∂αsα,i εr,i ∂α(A.16)

ここで，下添字 L および R はそれぞれ媒質境界の左極限および右極限の電磁界であることを意味している．これらの境界条件の導出には，式 (A.3)，式 (A.6)，式 (A.7)，および式 (A.8) を用いている．
式 (A.9)〜(A.16) の境界条件は統一的に，次のように表現できる．

ϕR = θ E ϕL
∂ϕR
∂ϕL
= θH
∂α
∂α

(A.17)(A.18)

ここで，θE ，θH は未知電磁界成分によって異なり，表 A.1 のように定義される．
高精度差分式の導出のためには，追加の境界条件が必要になる．媒質境界極限にて，次の 1 次
元波動方程式
)
( 2
∂
2
β2ϕ =
+
k
ε
(A.19)
0 r ϕ∂α2が満足されることを用いて，高次の境界条件を導く．式 (A.17) に式 (A.19) を用いると，次が導ける．
( 2)∂ 2 ϕR∂=θ+η
ϕL(A.20)E∂α2∂α2
ここで，η = k02 (εr,i − εr,i+1 ) である．式 (A.18) に式 (A.19) を用いると，次が導ける．
( 3)∂∂ 3 ϕR
∂
ϕL
= θH
+η
∂α3
∂α3
∂α
式 (A.20)〜(A.21) の境界条件は，2 次精度差分式 (IFD2) の導出に利用する．

92(A.21)

付 録B

Stern の差分式

図 A.1 において，隣り合う離散点における屈折率が一様 (ni−1 = ni = ni+1 ) である場合を考える．このとき，α = i∆α を展開中心とした場合の ϕ の 2 階導関数の差分式は次のようになる．

ϕ̃i+1 − 2ϕi + ϕ̃i−1
∂2ϕ
(B.1)
=
+ O(∆α2 )2∂α∆α2
ここで，表記の便宜上 (∆α)j を ∆αj と表現している．また，O(∆αj ) は打ち切り誤差が j 次のオーダーであることを表している．隣接点において屈折率が異なる場合，式 (B.1) の ϕ̃i+1 および ϕ̃i−1は実際の界振幅とは異なる値をとる．Stern の方法では，この見掛け上の振幅を媒質間の境界条件を考慮することによって，実際の振幅に置き換える方法をとる．
ϕi+1 と ϕi について，次のような Taylor 展開を考える．
()∞∑
1 ∆α n ∂ n ϕL
(B.2)
ϕ̃i+1 =
n!
2
∂αn
n=0
(
)
∞
∑
1
∆α n ∂ n ϕL
ϕi =
−
(B.3)
n!
2
∂αn
n=0
(
)
∞
∑
1 ∆α n ∂ n ϕR
ϕi+1 =
(B.4)
n!
2
∂αn
n=0
(
)
∞
∑
∆α n ∂ n ϕR
1
−
(B.5)
ϕ̃i =
n!
2
∂αn
n=0

式 (B.2) は α = (i + 12 )∆α における左極限を展開中心とした ϕi+1 の Taylor 級数表現であるが，ϕi+1が含まれる媒質の外を展開中心としているため，見掛け上の振幅となっている．同様の理由により，式 (B.5) の ϕ̃i も見掛け上の振幅である．式 (B.2) と式 (B.3) の和と差，および式 (B.3) と式 (B.4)の和と差を計算すると，それぞれ以下のようになる．
)
1(
ϕ̃i+1 + ϕi = ϕL + O(∆α2 )
(B.6)
2
) ∂ϕ
1 (
L
+ O(∆α2 )(B.7)
ϕ̃i+1 − ϕi =
∆α
∂α
(
)
1
(B.8)
ϕi+1 + ϕ̃i = ϕR + O(∆α2 )
2
) ∂ϕ
1 (
R
(B.9)
ϕi+1 − ϕ̃i =+ O(∆α2 )∆α∂αO(∆α2 ) の項を落とし，式 (B.6) と式 (B.8) に境界条件の式 (A.17) を，式 (B.7) と式 (B.9) に式 (A.18)
をそれぞれ適用して以下の関係を得る．
()
ϕi+1 + ϕ̃i = θE ϕ̃i+1 + ϕi
(B.10)()
ϕi+1 − ϕ̃i = θH ϕ̃i+1 − ϕi(B.11)
93

式 (B.10) と式 (B.11) の和を計算し，見掛け上の電磁界 ϕ̃i+1 に関して次の関係を得る．

ϕ̃i+1 = −θE − θH2
ϕi +ϕi+1 ≡ Fi ϕi + Fi+1 ϕi+1
θE + θH
θE + θH(B.12)

見掛け上の電磁界 ϕ̃i−1 に関しても，同様の手順で次の関係が得られる．

ϕ̃i−1 =′ − θ′
θE
2
H
ϕ
−
i−1
′ + θ′′ + θ ′ ϕi ≡ Gi−1 ϕi−1 + Gi ϕi
θE
θ
H
E
H

(B.13)′ ，θ ′ は表 A.1 において i + 1 → i − 1 の入れ替えを行った量である．式 (B.12)，(B.13)
ここで，θEH
を式 (B.1) に代入すると，最終的に次の差分式を得る．

∂2ϕ
Fi+1
Fi + Gi − 2
Gi−1≈
ϕi+1 +
ϕi +
ϕi−1
2
2
2
∂α
∆α
∆α
∆α2
表 2.2 に示した導関数の Stern の方法による差分表現を以下に列挙する．

1 ∂ 2 Eα
≈ AEα,i+1 + BEα,i + CEα,i−1s2α ∂α22εr,i+1
1A=2sα,i ∆α sα,i εr,i + sα,i+1 εr,i+1
(
)
2εr,i
2εr,i
1
B=−
+
sα,i ∆α2 sα,i εr,i + sα,i+1 εr,i+1 sα,i εr,i + sα,i−1 εr,i−12εr,i−11C=2sα,i ∆α sα,i εr,i + sα,i−1 εr,i−1

Dαα Eα =1 ∂ 2 Eβ
≈ AEβ,i+1 + BEβ,i + CEβ,i−1s2α ∂α212A=sα,i ∆α2 sα,i + sα,i+1
(
)
1
2
2
B=−
+
= −A − Csα,i ∆α2 sα,i + sα,i+1 sα,i + sα,i−1
12C=sα,i ∆α2 sα,i + sα,i−1

Dαα Eβ =1 ∂ 2 Hα
≈ AHα,i+1 + BHα,i + CHα,i−1s2α ∂α212A=2sα,i ∆α sα,i + sα,i+1
(
)
1
2
2
B=−
+
= −A − Csα,i ∆α2 sα,i + sα,i+1 sα,i + sα,i−112C=2sα,i ∆α sα,i + sα,i−1Dαα Hα =94(B.14)

1 ∂ 2 Hβ≈ AHβ,i+1 + BHβ,i + CHβ,i−1s2α ∂α2
2εr,i
1
A=
2
sα,i ∆α sα,i εr,i + sα,i+1 εr,i+1
(
)
2εr,i
2εr,i
1
B=−
+
= −A − Csα,i ∆α2 sα,i εr,i + sα,i+1 εr,i+1 sα,i εr,i + sα,i−1 εr,i−1
2εr,i
1
C=
2
sα,i ∆α sα,i εr,i + sα,i−1 εr,i−1

Dαα Hβ =

95

付 録C

偏波結合項の差分式

フルベクトル解析の場合に現れる偏導関数
{
}
1 ∂ϕ
∂
∂α f (α) ∂βの座標 (α, β) = (i∆α, j∆β) を差分中心としたときの差分式を考える. 本節で導出する差分式は論文 [72][73] の差分式と一致する．以降，座標 (i∆α, j∆β) を (i, j)，ϕ(i∆α, j∆β) を ϕi,j と略記する．
初めに，
(
)
1 ∂ϕ
f ∂β i± 1 ,j
2について，(i, j) を展開中心として Taylor 展開すると，次のようになる．
(
)
(
)
(
)
∞
∑
1 ∆α n ∂ n
1 ∂ϕ
1 ∂ϕ
=
f ∂β i+ 1 ,jn!
2
∂αn f ∂β i,j
n=0
2
)
(
)
(
)
(
∞
∑
1
∆α n ∂ n
1 ∂ϕ
1 ∂ϕ
=
−
f ∂β i− 1 ,jn!
2
∂αn f ∂β i,j
2

(C.1)
(C.2)

n=0

式 (C.1) と式 (C.2) の差をとり，次を得る．
{(
}
(
)
)
)
(
(
)
1 ∂ϕ
1
1 ∂ϕ
1 ∂ϕ
∂
=
−
+ O ∆α2∂α f ∂β i,j∆αf ∂β i+ 1 ,jf ∂β i− 1 ,j
2

(C.3)

2

次に，ϕi± 1 ,j±1 について，(i ± 12 , j) を展開中心として Taylor 展開すると，次のようになる．
2

ϕi± 1 ,j+1 =2

ϕi± 1 ,j−1 =
2

∞
∑
1
n=0
∞
∑
n=0

n!

(∆β)n

∂ n ϕi± 1 ,j
2

(C.4)

∂β n
n

∂ ϕi± 1 ,j12(−∆β)n
n!
∂β n

(C.5)

式 (C.4) と式 (C.5) の差をとり，次を得る．
( )
)
(
)
∂ϕ
1 (
=
ϕi± 1 ,j+1 − ϕi± 1 ,j−1 + O ∆β 222
∂β i± 1 ,j
2∆β

(C.6)

2

続いて，ϕi,j±1 ，ϕi+1,j±1 について，(i + 12 , j ± 1) を展開中心として Taylor 展開すると，次のよう
になる．
(
) n
∞
∑
∆α n ∂ ϕi+ 12 ,j±1
1
−
(C.7)
ϕi,j±1 =
n!
2
∂αn
n=0
(
) n
∞
∑
1 ∆α n ∂ ϕi+ 12 ,j±1
ϕi+1,j±1 =
(C.8)
n!
2
∂αn
n=0

96

式 (C.7) と式 (C.8) の和をとり，次を得る．

ϕi+ 1 ,j±1 =2()
1(ϕi+1,j±1 + ϕi,j±1 ) + O ∆α2
2

(C.9)

同様に，ϕi−1,j±1 ，ϕi,j±1 について，(i − 12 , j ± 1) を展開中心として Taylor 展開すると，次のよう
になる．

(
) n
∞
∑
∆α n ∂ ϕi− 12 ,j±11−ϕi−1,j±1 =
n!
2
∂αn
n=0
(
) n
∞
∑
1 ∆α n ∂ ϕi− 12 ,j±1ϕi,j±1 =n!
2∂αn(C.10)(C.11)

n=0

式 (C.10) と式 (C.11) の和をとり，次を得る．

ϕi− 1 ,j±1 =2()
1(ϕi−1,j±1 + ϕi,j±1 ) + O ∆α22(C.12)

次に，fi,j ，fi+1,j について，(i + 12 , j) を展開中心として Taylor 展開すると，次のようになる．

(
) n
∞
∑
1
∆α n ∂ fi+ 12 ,j
−
n!
2
∂αn
n=0
(
) n
∞
∑
1 ∆α n ∂ fi+ 12 ,jfi+1,j =n!
2∂αnfi,j =(C.13)(C.14)

n=0

式 (C.13) と式 (C.14) の和をとり，次を得る．

fi+ 1 ,j =2()
1(fi,j + fi+1,j ) + O ∆α22(C.15)

さらに，fi−1,j ，fi,j について，(i − 12 , j) を展開中心として Taylor 展開すると，次のようになる．

(
) n
∞
∑
1
∆α n ∂ fi− 12 ,jfi−1,j =
−
n!
2
∂αn
n=0
(
) n
∞
∑
1 ∆α n ∂ fi− 12 ,jfi,j =n!
2∂αn(C.16)(C.17)

n=0

式 (C.16) と式 (C.17) の和をとり，次を得る．

fi− 1 ,j =2()
1
(fi−1,j + fi,j ) + O ∆α22(C.18)()()
これまでに求めた 5 つの式 (C.6)，(C.9)，(C.12)，(C.15)，(C.18) において，O ∆α2 および O ∆β 2の項を落とした上でこれらを式 (C.3) に代入すると，最終的に次の差分式が得られる．
(
)
1 ∂ϕ
∂
≈ A1 ϕi−1,j−1 + A2 ϕi,j−1 + A3 ϕi+1,j−1∂α f ∂β i,j+A4 ϕi−1,j+1 + A5 ϕi,j+1 + A6 ϕi+1,j+1

97(C.19)

ここで，係数 A1 〜A6 はそれぞれ次のようになる．

112∆α∆β fi−1,j + fi,j
)
(
1
1
1
−
A2 =
2∆α∆β fi−1,j + fi,jfi,j + fi+1,j1
1
A3 = −
2∆α∆β fi,j + fi+1,j1
1
A4 = −
2∆α∆β fi−1,j + fi,j
(
)
1
1
1
A5 = −−
2∆α∆β fi−1,j + fi,jfi,j + fi+1,j1
1
A6 =2∆α∆β fi,j + fi+1,j

A1 =(C.20)(C.21)(C.22)(C.23)(C.24)(C.25)

式 (C.19) を応用すると，電界および磁界のそれぞれに対する偏波結合項の差分式は次のようになる．
{}
1 ∂ (εr Eβ ) ∂Eβ1 ∂−≈ T1 Eβ,i−1,j−1 + T2 Eβ,i,j−1 + T3 Eβ,i+1,j−1
sα sβ ∂α εr∂β∂β+T4 Eβ,i−1,j+1 + T5 Eβ,i,j+1 + T6 Eβ,i+1,j+1()
2εr,i−1,j−1
1
1
T1 =
−1
sα,i,j sβ,i,j 4∆α∆β εr,i−1,j + εr,i,j()
2εr,i,j−12εr,i,j−1
1
1
T2 =
−
sα,i,j sβ,i,j 4∆α∆β εr,i−1,j + εr,i,j
εr,i,j + εr,i+1,j()
2εr,i+1,j−1
1
1
T3 = −−1sα,i,j sβ,i,j 4∆α∆β εr,i,j + εr,i+1,j()
2εr,i−1,j+1
1
1
T4 = −−1sα,i,j sβ,i,j 4∆α∆β εr,i−1,j + εr,i,j()
2εr,i,j+12εr,i,j+1
1
1
T5 = −−
sα,i,j sβ,i,j 4∆α∆β εr,i−1,j + εr,i,j
εr,i,j + εr,i+1,j()
2εr,i+1,j+1
1
1
−1
T6 =
sα,i,j sβ,i,j 4∆α∆β εr,i,j + εr,i+1,j98(C.26)(C.27)(C.28)(C.29)(C.30)(C.31)(C.32)

1
sα sβ

{

∂Hα
∂
− εr
∂α∂β
∂α

(

1 ∂Hα
εr ∂β

)}
≈ T1 Hα,i−1,j−1 + T2 Hα,i,j−1 + T3 Hα,i+1,j−1+T4 Hα,i−1,j+1 + T5 Hα,i,j+1 + T6 Hα,i+1,j+1}{
2εr,i,j1
1T1 =1−sα,i,j sβ,i,j 4∆α∆βεr,i−1,j + εr,i,j{}2εr,i,j2εr,i,j1
1−
T2 = −
sα,i,j sβ,i,j 4∆α∆β εr,i−1,j + εr,i,j
εr,i,j + εr,i+1,j{}2εr,i,j1
1
T3 = −
1−sα,i,j sβ,i,j 4∆α∆βεr,i,j + εr,i+1,j{}2εr,i,j1
1
T4 = −
1−sα,i,j sβ,i,j 4∆α∆βεr,i−1,j + εr,i,j{}2εr,i,j2εr,i,j1
1−T5 =sα,i,j sβ,i,j 4∆α∆β εr,i−1,j + εr,i,j
εr,i,j + εr,i+1,j{}2εr,i,j1
1
T6 =1−sα,i,j sβ,i,j 4∆α∆βεr,i,j + εr,i+1,j99(C.33)(C.34)(C.35)(C.36)(C.37)(C.38)(C.39)

付 録D

2 次の高精度差分式 (IFD2)

α = (i + 21 )∆α の媒質境界における右極限を展開中心として ϕi+1 を，α = i∆α を展開中心として左極限電磁界 ϕL をそれぞれ Taylor 展開する．
(
)
(
)
( 4)
∆α ∂ϕR1 ∆α 2 ∂ 2 ϕR
1 ∆α 3 ∂ 3 ϕR
ϕi+1 = ϕR +
+
+
+
O
∆α
(D.1)
2 ∂α
2!
2
∂α2
3!
2
∂α3
(
)
(
)
( 4)
∆α ∂ϕi1 ∆α 2 ∂ 2 ϕi
1 ∆α 3 ∂ 3 ϕi
(D.2)
ϕL = ϕi +
+
+
+
O
∆α
2 ∂α
2!
2
∂α2
3!
2
∂α3
式 (D.1)，(D.2) に式 (A.17)〜(A.21) の境界条件を適用すると，次が導ける．

ϕi+1 = f0 ϕi + f1)(
∂ϕi
∂ 2 ϕi∂ 3 ϕi+ f2+
f+ O ∆α4
3
2
3
∂α
∂α
∂α

}
{
η
f0 = θE 1 + ∆α2
8
{
}
{
}
1
η
1
η
3
3
f1 = θE∆α + ∆α + θH
∆α + ∆α
2
16
2
48
1
1
f2 = θE ∆α2 + θH ∆α2441
1
f3 = θE ∆α3 + θH ∆α3
12
12

(D.3)

(D.4)
(D.5)
(D.6)
(D.7)

式 (D.3) は，媒質間の境界条件を考慮して，ϕi+1 を α = i∆α を展開中心として Taylor 展開した式とみなせる．これは一般化 Taylor 展開と呼ばれている．
同様に ϕi−1 を α = i∆α を展開中心として一般化 Taylor 展開すると，次が導ける．

ϕi−1 = e0 ϕi + e1()
∂ 2 ϕi∂ 3 ϕi∂ϕi+ e2+e+ O ∆α4
3
2
3
∂α
∂α
∂α

(D.8)

ここで，係数 e0 〜e3 は式 (D.4)〜(D.7) において ∆α → −∆α，εr,i+1 → εr,i−1 ，si+1 → si−1 の入れ替えを行うことで得られる．
FD-BPM では，次の α = i∆α における ϕ の 2 階導関数が波動方程式に現れる．

1 ∂ 2 ϕi2sα,i ∂α2(D.9)

この 2 階導関数の 2 次精度差分式を求める．式 (D.3) の両辺に e1 を，式 (D.8) の両辺に f1 をそれ
ぞれ乗じる．

()∂ 3 ϕi∂ 2 ϕi∂ϕi+fe+ O ∆α5+ f2 e1
3
1
2
3
∂α
∂α
∂α
(
)
∂ϕi
∂ 2 ϕi∂ 3 ϕi
f1 ϕi−1 = f1 e0 ϕi + f1 e1+ f1 e2+ f1 e3+ O ∆α5
2
3
∂α
∂α
∂α
e1 ϕi+1 = f0 e1 ϕi + f1 e1

100(D.10)(D.11)

式 (D.10)，(D.11) より 1 階導関数の項を消去して整理すると，次が得られる．

( 2)
f1 ϕi−1 + (f0 e1 − f1 e0 ) ϕi − e1 ϕi+1
∂ 2 ϕi f1 e3 − f3 e1 ∂ 3 ϕi
=
+
+
O
∆α
f1 e2 − f2 e1∂α2
f1 e2 − f2 e1 ∂α3(D.12)

式 (D.12) 右辺第 2 項は 3 次の打ち切り誤差を持つ．以上より，次の差分式が得られる．

( 2)
1 ∂ 2 ϕi
1 f1 ϕi−1 + (f0 e1 − f1 e0 ) ϕi − e1 ϕi+1+O∆α=f1 e2 − f2 e1
s2α,i ∂α2
s2α,i

101

(D.13)