(57)【要約】【課題】本発明は、電波、光波、Ｘ線などを含む電磁波工学の分野に属する。フォトニック結晶偏光分離素子に代表される誘電率周期構造における遮断域の減衰定数を精度良く計算する手法を提供する。【解決手段】  全方位周期境界型ＦＤＴＤ法などの既存の数値解法によって求めた伝搬域・遮断域境界での周期構造内の電磁界分布の試験関数を、固有周波数の変分表現に代入して数値積分を評価することにより、遮断域での減衰定数と固有周波数との関係を従来手法に比べ少ない計算量で精度良く計算することができる。【選択図】        図１
2007064911.tif 000002
【特許請求の範囲】
【請求項１】
誘電率が３次元空間内で周期的に変化している構造体が、ある特定の周波数を持つ電磁波に対して呈する遮断特性において、その電磁波の周波数と減衰量との関係を計算する手法であって、電磁波の固有周波数の変分表現を利用することを特徴とする計算手法
【請求項２】
請求項１に記載の計算手法であって、電磁波の固有周波数の変分表現を計算する際に用いる試験関数として、空間微分の有限差分近似を行うことを特徴とする数値解法で求められたものを用いることを特徴とする計算手法
【請求項３】
請求項２に記載の計算手法であって、電磁波の固有周波数の変分表現を計算する際に用いる試験関数として、有限差分時間領域法を含む手法で求められた関数を用いることを特徴とする計算手法
【請求項４】
請求項１から請求項３に記載の計算手法を電子計算機上で実施するためのプログラム、及びプログラムの記録された記憶媒体
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、電波、光波、Ｘ線などを含む電磁波工学の分野に属する。特に特定の周波数および偏波面および進行方向を持つ電磁波を、透明もしくはほぼ透明な誘電体材料からなる周期構造体を用いて効率よく反射させる反射型素子において、その反射量を数値計算により求める際の手法に関する。
【背景技術】
【０００２】
誘電率が空間中で周期的に変化する誘電率周期構造の中では、電磁波が位相の揃った多重反射を起こす。特に周期構造の周期をd、構造中での電磁波の平均波長をλaとしたとき、λa＝ｄ／２付近のある有限な幅の波長帯または周波数帯では、全ての反射波が同位相で強めあういわゆるブラッグ反射が生じる。ブラッグ反射の生じる周波数帯では、周期構造を構成する個々の誘電体がその周波数の電磁波に対して無損失であっても、電磁波の電力は周期構造によって完全に反射される。
【０００３】
この現象はＸ線物理の分野においては結晶のＸ線回折分析、光エレクトロニクスの分野においては誘電体多層膜反射鏡、そして電波工学の分野では八木・宇田アンテナなど、電磁波工学の広い分野にわたって応用されてきた。
特に光エレクトロニクスの分野においては、誘電率が空間の一方向にのみ周期変調された、伝統的な交互多層膜型反射鏡の他に、近年、誘電率が空間中で２つ以上の方向に周期変調された構造体である「フォトニック結晶」の応用が活発に検討されている。
【０００４】
フォトニック結晶の示す立体的なブラッグ反射の特性を利用することで、従来にはない光制御素子を実現することができると期待されている。これらの素子の一例としては、薄膜型偏光分離素子（特許文献１、非特許文献１）、超小型光集積回路、無閾値レーザなどがある。
【０００５】
誘電率周期構造を応用する際、多くの場合は周期をどこかで打ち切り、その端面を電磁波の入射面とする。このような半無限周期構造にその端面からブラッグ反射の周波数帯にある周波数の波を入射させると、波の電界及び磁界の振幅は周期構造体の端面から内部に侵入するに従い、指数関数的に減少する。
【０００６】
１周期進むごとに振幅が指数関数的に、ｅｘｐ（−α）倍になるとすると、この実数αは波の周波数の関数であり、一般にブラッグ反射帯域の中央付近の周波数において最も大きく、そこから周波数が離れるに従って徐々に小さくなり、やがてブラッグ反射の帯域から外れるとゼロになる。現実の応用においてはこのような半無限の周期構造が用いられることは少なく、素子の生産性、及び素子寸法の小型化を重視して、入射側端面から数えてある周期数のところで周期構造を再び打ち切る場合がほとんどである。
【０００７】
このような有限周期の周期構造においては、原理的には周期を打ち切った地点から先に入射波の電力の一部が流れていく。すなわち入射波の一部が周期構造を透過することになり、入射側から見て電力の反射率が１００％とならない。ブラッグ反射の一周期あたりの減衰定数αが大きい場合、周期を打ち切った地点では電磁界の振幅が十分に減衰しているので、反射率は１００％に十分に近い値となる。一方、減衰定数αが小さいと透過率が増し、反射率が低下する。
【０００８】
このように動作周波数において元の周期構造の減衰定数αがいかなる値を持っているかは、所望の反射率および透過率を得るのに必要な周期数を決定するのに極めて重要な情報であり、これを高精度かつ簡便に計算することが上記の電磁波工学の各応用分野で必要となる。
光エレクトロニクス分野で広く用いられている交互多層膜のような一次元の周期構造については、そのブラッグ反射の減衰定数を数値丸め誤差の範囲内程度で厳密に求める手法が確立されており、減衰定数と周波数の関係がどのような曲線を描くかも広く知られている。
【０００９】
また一部の特殊な二次元誘電率周期構造についても、一次元構造での手法を拡張することで厳密解を得ることができることが知られている（非特許文献２）。その一方で、一般の二次元以上の多次元周期構造では、その減衰定数を解析的ないし厳密に求めることは困難とされており、もっぱらマックスウェルの方程式もしくは波動方程式を数値的に解くことで周期構造内の電磁界分布を求め、そこから減衰定数を計算することが行われている。代表的な数値解法と、それらの特質は次の通りである。
１）有限要素法（非特許文献３）。マクスウェルの方程式と変分原理を組み合わせた手法で、計算量が同程度なら他の手法よりも精度は高い傾向にある。他方、３次元計算に必要な立体空間メッシュの自動生成方法が難しく、計算プログラムの開発は簡単ではない。
２）有限差分時間領域法（ＦＤＴＤ法）と周期境界条件の組み合わせ（非特許文献４）。       
【００１０】
周期構造の一周期分を解析空間とし、解析空間の端同士を周期境界条件で終端することで、伝搬域の電磁波の分散関係、すなわち波数ベクトルと周波数の関係を計算することができる。この型の計算方法を、本文章中では便宜的に「全方位周期境界型のＦＤＴＤ法」と呼ぶことにする。遮断域を取り扱うことは原理的にできない。      
３）同じくＦＤＴＤ法と吸収境界条件の組み合わせ（非特許文献５）。十分に周期数の多い半無限周期構造を解析空間中に構築し、そこに平面波が入射した場合の反射・透過波の振る舞いを計算する。
【００１１】
反射スペクトルから遮断域の各周波数における減衰定数を知ることができる。３次元空間のxyz座標においてｚ≧ｚ0（ただしｚ0＞０）の空間に周期構造が広がっており、z=0の位置からz軸の正の方向に向けて平面波を入射させる場合、ｚ＝ｚ1（ただし  ｚ1＜０）の位置とｚ＝ｚ2（ただしｚ2＞０）の位置にそれぞれ吸収境界を設置する。なおｚ2＞＞ｚ1である。ｚ＝ｚ1の吸収境界は周期構造の端面に当たって反射してきた波を吸収するためのものである（非特許文献６）。
【００１２】
一方ｚ＝ｚ2の吸収境界は周期構造の奥深くに侵入した波を徐々に吸収させ、ｚ＝０の地点から見るとあたかも周期構造の向こう側の端が無限の彼方にあるかのようにみなせるようにするためのものである。一方ｘ方向には、ｘ＝０からｘ＝ａxまでを、またｙ方向にはｙ＝０からｙ＝ａyまでを解析空間とする。ここでａx，ａyはそれぞれ周期構造のｘ方向及びｙ方向の１周期の長さである。ｘ＝０とｘ＝ａxの２つの面同士を周期境界条件で終端する。
【００１３】
ｙ＝０とｙ＝ａyの面同士も周期境界条件で終端する。入射した波のパワーが徐々に吸収境界によって減少することから、この型の計算方法を本文章中では便宜的に「吸収境界型のＦＤＴＤ法」と呼ぶことにする。ｚ＝ｚ1からｚ＝ｚ2までの間には周期構造の単位構造を非常に多数置く必要があるので、大きな計算空間と長い計算時間が必要とされる。
４）平面波展開法（非特許文献７）。これも２）の方法と同じく、取り扱えるのは伝搬域に限られる。
上記の各方法のうち、２）と３）のＦＤＴＤ法はアルゴリズムが分かり易くプログラムも比較的容易に作成できるため、フォトニック結晶の研究開発分野では最も一般的に利用されている。
【００１４】
しかし上述３）のように、遮断域中での光波の振る舞いを求めるのには非常に多くの計算量が必要となるというな問題点を抱えていた。
ところで従来型の誘電体光導波路においては伝搬定数の変分表現が古くから知られており、導波モードの分散関係の計算など、導波路理論において解析的な処理の難しい種々の問題を取り扱うための便利なツールとして用いられてきた（非特許文献８）。
【００１５】
一方フォトニック結晶のような周期構造体においても固有周波数の変分表現が近年見出され、そのパスバンド、すなわち伝搬域における分散関係の計算のためのツールとして用いる手法が提案されている（非特許文献９、非特許文献１０）。しかし遮断域の減衰定数をこの変分表現を利用して求める方法はこれまで知られていなかった。
【特許文献１】特許第 ３２８８９７６号
【非特許文献１】Ｙ．Ｏｈｔｅｒａ，  Ｔ．  Ｓａｔｏ，  Ｔ．  Ｋａｗａｓｈｉｍａ，  Ｔ．  Ｔａｍａｍｕｒａ，  ａｎｄ  Ｓ．  Ｋａｗａｋａｍｉ，“Ｐｈｏｔｏｎｉｃ  Ｃｒｙｓｔａｌ  Ｐｏｌａｒｉｚａｔｉｏｎ  Ｓｐｌｉｔｔｅｒｓ，”  Ｅｌｅｃｔｒｏｎｉｃｓ  Ｌｅｔｔｅｒｓ，  ｖｏｌ．３５，ｎｏ．１５，  ｐｐ．１２７１−１２７２，  １９９９．
【非特許文献２】Ｓ．  Ｋａｗａｋａｍｉ，  “Ａｎａｌｙｔｉｃａｌｌｙ  Ｓｏｌｕｂｌｅ  ｍｏｄｅｌ  ｏｆ  ｐｈｏｔｏｎｉｃ  ｃｒｙｓｔａｌ  ｓｔｒｕｃｔｕｒｅｓ  ａｎｄ  ｎｏｖｅｌ  ｂｅｈａｖｉｏｒ  ｏｆ  ｒｅｌａｔｅｄ  ｒｅｓｏｎａｔｏｒｓ  ａｎｄ  ｗａｖｅｇｕｉｄｅｓ，”  Ｊｏｕｒｎａｌ  ｏｆ  Ｌｉｇｈｔｗａｖｅ  Ｔｅｃｈｎｏｌｏｇｙ，  ｖｏｌ．２０，ｎｏ．８，ｐｐ．１６４４−１６５０，  ２００２．
【非特許文献３】橋本尚人、小柴正則、「無限長周期構造導波路の有限要素法解析の一手法」、平成３年度電気関係学会北海道支部連合大会、講演番号５５
【非特許文献４】Ｃ．Ｔ．Ｃｈａｎ，  Ｑ．Ｌ．Ｙｕ，  ａｎｄ  Ｋ．Ｍ．Ｈｏ，  “Ｏｒｄｅｒ−Ｎ  ｓｐｅｃｔｒａｌ  ｍｅｔｈｏｄ  ｆｏｒ  ｅｌｅｃｔｒｏｍａｇｎｅｔｉｃ  ｗａｖｅｓ，”    Ｐｈｙｓｉｃａｌ  Ｒｅｖｉｅｗ  Ｂ，  ｖｏｌ．５１，ｎｏ．２３，ｐｐ．１６６３５−１６６４２，１９９５．
【非特許文献５】Ｓ．Ｆａｎ，  Ｐ．Ｒ．Ｖｉｌｌｅｎｅｕｖｅ，  ａｎｄ  Ｊ．Ｄ．Ｊｏａｎｎｏｐｏｕｌｏｓ，  “Ｌａｒｇｅ  ｏｍｉｎｉｄｉｒｅｃｔｉｏｎａｌ  ｂａｎｄ  ｇａｐｓ  ｉｎ  ｍｅｔａｌｌｏｄｉｅｌｅｃｔｒｉｃ  ｐｈｏｔｏｎｉｃ  ｃｒｙｓｔａｌｓ，”  Ｐｈｙｓｉｃａｌ  Ｒｅｖｉｅｗ  Ｂ，  ｖｏｌ．５４，ｎｏ．１６，ｐｐ．１１２４５−１１２５１，１９９６．
【非特許文献６】Ｍ．Ｋｏｓｈｉｂａ，  Ｙ．Ｔｓｕｊｉ，  ａｎｄ  Ｓ．Ｓａｓａｋｉ，“Ｈｉｇｈ−ｐｅｒｆｏｒｍａｎｃｅ  ａｂｓｏｒｂｉｎｇ  ｂｏｕｎｄａｒｙ  ｃｏｎｄｉｔｉｏｎｓ  ｆｏｒ  ｐｈｏｔｏｎｉｃ  ｃｒｙｓｔａｌ  ｗａｖｅｇｕｉｄｅ  ｓｉｍｕｌａｔｉｏｎｓ，”  ＩＥＥＥ  Ｍｉｃｒｏｗａｖｅ  ａｎｄ  Ｗｉｒｅｌｅｓｓ  Ｃｏｍｐｏｎｅｎｔｓ  Ｌｅｔｔｅｒｓ，ｖｏｌ．１１，ｎｏ．４，ｐｐ．１５２−１５４，２００１．
【非特許文献７】Ｋ．Ｍ．Ｌｅｕｎｇ  ａｎｄ  Ｙ．Ｆ．Ｌｉｕ，“Ｆｕｌｌ  ｖｅｃｔｏｒ  ｗａｖｅ  ｃａｌｃｕｌａｔｉｏｎ  ｏｆ  ｐｈｏｔｏｎｉｃ  ｂａｎｄ  ｓｔｒｕｃｔｕｒｅｓ  ｉｎ  ｆａｃｅ−ｃｅｎｔｅｒｅｄ−ｃｕｂｉｃ  ｄｉｅｌｅｃｔｒｉｃ  ｍｅｄｉａ，”Ｐｈｙｓｉｃａｌ  Ｒｅｖｉｅｗ  Ｌｅｔｔｅｒｓ，ｖｏｌ．６５，ｎｏ．２１，ｐｐ．２６４６−２６４９，１９９０．
【非特許文献８】Ａ．Ｗ．Ｓｎｙｄｅｒ  ａｎｄ  Ｊ．Ｄ．Ｌｏｖｅ，“Ｏｐｔｉｃａｌ  Ｗａｖｅｇｕｉｄｅ  Ｔｈｅｏｒｙ，”  Ｋｌｕｗｅｒ  Ａｃａｄｅｍｉｃ  Ｐｕｂｌｉｓｈｅｒｓ，Ｌｏｎｄｏｎ，１９８３,  ｐ．６４３．
【非特許文献９】Ｙ．Ｏｈｔｅｒａ，  Ｋ．Ｋｕｒｏｋａｗａ，  ａｎｄ  Ｓ．Ｋａｗａｋａｍｉ，“Ｖａｒｉａｔｉｏｎａｌ  ｅｘｐｒｅｓｓｉｏｎ  ｆｏｒ  ｔｈｅ  ａｎａｌｙｓｉｓ  ｏｆ  ｐｈｏｔｏｎｉｃ  ｃｒｙｓｔａｌ  ｄｅｖｉｃｅｓ，”  ＩＥＥＥ  Ｊｏｕｒｎａｌ  ｏｆ  Ｑｕａｎｔｕｍ  Ｅｌｅｃｔｒｏｎｉｃｓ，ｖｏｌ．３８，ｎｏ．７，ｐｐ．９１９−９２６，２００２．
【非特許文献１０】Ｙ．Ｏｈｔｅｒａ，  Ｙ．Ｓａｓａｋｉ，  ａｎｄ  Ｓ．Ｋａｗａｋａｍｉ，“Ｐｏｓｔｐｒｏｃｅｓｓｉｎｇ  ｏｆ  ＦＤＴＤ  ｓｏｌｕｔｉｏｎｓ  ｆｏｒ  ｐｒｅｃｉｓｅ  ｃａｌｃｕｌａｔｉｏｎｓ  ｏｆ  ｅｉｇｅｎｆｒｅｑｕｅｎｃｉｅｓ  ｏｆ  ｐｈｏｔｏｎｉｃ  ｐｅｒｉｏｄｉｃ  ｓｔｒｕｃｔｕｒｅｓ  ｕｔｉｌｉｚｉｎｇ  ｔｈｅ  ｖａｒｉａｔｉｏｎａｌ  ｅｘｐｒｅｓｓｉｏｎ，”  Ｊｏｕｒｎａｌ  ｏｆ  Ｌｉｇｈｔｗａｖｅ  Ｔｅｃｈｎｏｌｏｇｙ，ｖｏｌ．２２，ｎｏ．６，ｐｐ．１６２８−１６３６，２００４．
【発明の開示】
【発明が解決しようとする課題】
【００１６】
本発明は上述した従来の数値解析手法が有していた、誘電率周期構造の遮断周波数域における電磁波の減衰定数を求めるのに多大な計算量が必要とされる、あるいは求めること自体が困難であった、という課題の解決を目的としている。
【課題を解決するための手段】
【００１７】
本発明によれば、誘電率が３次元空間内で周期的に変化している構造体が、ある特定の周波数を持つ電磁波に対して呈する遮断特性において、その電磁波の周波数と減衰量との関係を計算する手法であって、電磁波の固有周波数の変分表現を利用することを特徴とする計算手法が得られる。
【００１８】
また本発明は、上記の計算手法であって、電磁波の固有周波数の変分表現を計算する際に用いる試験関数として、空間微分の有限差分近似を行うことを特徴とする数値解法で求められたものを用いることを特徴とする計算手法を提供する。
また本発明は、上述の計算手法であって、電磁波の固有周波数の変分表現を計算する際に用いる試験関数として、ＦＤＴＤ法を含む手法で求められた関数を用いることを特徴とする計算手法を提供する。
また本発明は、上述の計算手法を電子計算機上で実行するのに必要なプログラムのアルゴリズムを手供する。
【発明の効果】
【００１９】
本発明によれば、伝搬域の分散関係の計算に多用されている、全方位周期境界型のＦＤＴＤ法と変分表現とを連携させることで、全方位周期境界型のＦＤＴＤ法単体ではこれまで原理的に求めることが不可能であった遮断域中での分散関係、すなわち減衰定数と周波数の関係を計算することができるようになる。
またこれまで吸収境界型のＦＤＴＤ法を用いて遮断特性が求められることも行なわれてきたが、これを本発明の方法に置き換えることによって、計算時間を大幅に短縮する効果が得られる。また吸収境界型のＦＤＴＤ法では伝搬域と遮断域の境界付近の周波数において、減衰特性を計算することが難しかったが、本発明は遮断域中の全ての周波数において、減衰定数を精度よく計算することができるという効果を持つ。
【発明を実施するための最良の形態】
【００２０】
以下、本発明の実施の形態について図面及び式を参照しながら説明する。
本発明は上記の問題を解決するために、「誘電率周期構造における電磁波の固有周波数の変分表現」を既存の数値解法と連携させて用いることを特徴とする。
以下では簡単のため、周期構造のブラベ格子型が斜方晶系、すなわち３つの基本併進ベクトルがいずれも直交するような場合について説明する。その他の格子型でも議論の根幹は同じである。
誘電率及び透磁率が３次元空間内でx, y, zの各方向にそれぞれ周期ａx，ａy，ａzをもって周期的に変化するような、誘電率周期構造中における電磁波の固有周波数の変分表現は以下の式で与えられることが知られている（非特許文献９）。
【００２１】
【数１】
2007064911.tif 000003
ただし、
【００２２】
【数２】
2007064911.tif 000004
である。εrは比誘電率、μr は比透磁率のテンソルであり、Rを任意の格子ベクトルとして、
【００２３】
【数３】
2007064911.tif 000005
【００２４】
【数４】
2007064911.tif 000006
で定義される周期性を持つ。
【００２５】
【数５】
2007064911.tif 000007
は位置ベクトルである。式（１）でwは電磁波の固有角周波数、c0は真空中の光速であり、Ｅsは電界の試験関数、
【００２６】
【数６】
2007064911.tif 000008
はＥsと共役な試験関数、yは伝搬定数ベクトルを表す。以後の文章ではyのことをＢｌｏｃｈ波数ベクトルと呼ぶこともある。式（２５）はyの複素共役y*に対応する試験関数と考えることができるので、実用上は単にＥsの複素共役とすればよい。式（１）の分子、分母の積分は周期構造の各方向１周期分の体積中で計算するものとし、Vuはその単位格子領域を表す。式（２）の
【００２７】
【数７】
2007064911.tif 000009
は領域Vuの表面における表面積分項であり、nは各面において外側を向く法線ベクトルである。試験関数Ｅs、式６（２５）として以下の式（５）、式（６）に示すＢｌｏｃｈ関数型のものを採用した場合は式７（２６）はゼロとなることが容易にわかる。
【００２８】
【数８】
2007064911.tif 000010
【００２９】
【数９】
2007064911.tif 000011
R'も任意の格子ベクトルである。伝搬定数ベクトルyは伝搬域中では純虚数であるが、遮断域中では一般には複素ベクトルとなる。すなわち、
【００３０】
【数１０】
2007064911.tif 000012
βは通常の意味での、すなわち波の空間的な平均位相変化を表すところの伝搬定数ベクトルであり、αは減衰定数ベクトルである。βもαもともに実数の成分を持つ。
【００３１】
式（１）の変分表現の重要な特徴は、「試験関数の１次の変分からは、固有周波数ω2の２次以上の誤差しか生じない」という点である。注目する固有周波数における電磁界の固有関数の形が完全にわかっている場合、その固有関数を上式の試験関数Ｅs,式６（２５）に代入して右辺を計算すると、原理的に固有周波数が求められる。ところが固有関数の形が厳密には不明である場合でも、その概略を知っていればそれを試験関数に用いることにより、高い精度で固有周波数の近似値を得ることができる。
【００３２】
さて、１次元周期構造でよく知られている知見と同様に、２次元、３次元の周期構造においても注目する遮断域の中での電磁界の固有関数を、遮断域の下端及び上端における固有関数の適当な線形結合でよく近似することができる。この性質を利用して、式（５）、式（６）の周期関数uに相当するこれらの固有関数をそれぞれｕL，ｕHとし、a ,bを複素数として以下の式で試験関数を定義する。
【００３３】
【数１１】
2007064911.tif 000013
ここで、
【００３４】
【数１２】
2007064911.tif 000014
【００３５】
【数１３】
2007064911.tif 000015
である。遮断域の端、すなわち逆格子空間におけるブリルアン・ゾーンの端ではβは既知である。一方αの方は適当に値を設定する。遮断域と伝搬域の境界付近を計算したい場合は、減衰が弱いはずなのでαとして小さな値を設定する。他方、遮断域の中心周波数付近では減衰が大きいので、適度に大きな値をαに設定する。なおこの時点では設定したαに対し意味のある固有周波数が計算できるかどうかは分からない。実際には起こり得ないほど大きな値のαに対応する固有値は存在しないからである。
【００３６】
【数１４】
2007064911.tif 000016
(*は複素共役を表す)の関係を使うと、式（８）と式（１）より遮断域中で使用する変分表現として下式を得る。
【００３７】
【数１５】
2007064911.tif 000017
この式を次のように整理する。
【００３８】
【数１６】
2007064911.tif 000018
【００３９】
【数１７】
2007064911.tif 000019
【００４０】
【数１８】
2007064911.tif 000020
p,θを実数として、
【００４１】
【数１９】
2007064911.tif 000021
のように係数を変換すると、式（１２）は更に以下に変形できる。
【００４２】
【数２０】
2007064911.tif 000022
与えられた複素数のyに対してｕL，ｕHを使って固有周波数の自乗、ω2を求めるのが本手法の主眼である。真の固有周波数においては、式（１５）の値はｕL，ｕHの係数変数p, qの変分に対し停留値を取る。すなわち次式が成り立つ。
【００４３】
【数２１】
2007064911.tif 000023
【００４４】
【数２２】
2007064911.tif 000024
途中の式変形の詳細は省略するが、式（１６）、（１７）を計算していくと次式に到達する。
【００４５】
【数２３】
2007064911.tif 000025
【００４６】
【数２４】
2007064911.tif 000026
【００４７】
【数２５】
2007064911.tif 000027
式（１９）には複号があるので、一般には与えられた複素伝搬定数ベクトルyに対し周波数の固有値が２つ、存在する。根号の中が0になると重根となり、それは遮断域中で減衰定数が最大値ａmaxを取る周波数を表すことになる。なおａmaxより大きな減衰定数から出発して形式上得られるwは物理的に意味のない解である。
以上では電界Eを用いた定式化について説明したが、磁界Hを用いた定式化も可能である。変分表現の式（１）、積分項の式（１３）、式（１４）は以下の形となる。その他の式については省略する。
【００４８】
【数２６】
2007064911.tif 000028
【００４９】
【数２７】
2007064911.tif 000029
【００５０】
【数２８】
2007064911.tif 000030
式（１５）から式（２０）はそのまま使うことができる。
以上を踏まえ、電界を使用する場合の実際の計算手順を次のように整理することができる。
手順１．まず全方位周期境界型のＦＤＴＤ法など既存の数値解法を用いて、注目する遮断域の下端及び上端の周波数を求める。ＦＤＴＤ法における計算精度は解析空間の分割数と時間ステップ数に依存するが、この段階では大まかな周波数がわかれば十分であるので空間分割数及び時間ステップ数は上端、下端二つの周波数を決定するのに必要な最小限の値でよい。
【００５１】
手順２．手順１で求めた２つの周波数における固有関数を、やはり既存の適当な数値解法で計算する。全方位周期境界型のＦＤＴＤ法でこのような固有関数を計算する手順は公知である。これらを式（５）の形式で表した際の周期関数部分をそれぞれｕL（ｒ），ｕH（ｒ）とする。
【００５２】
手順３．ｕL（ｒ），ｕH（ｒ）は空間格子上の離散的な点でしか値が与えられていない。この後の手順にある数値積分を精度よく計算するために、これらを実空間上で補間し、十分に細かい空間格子上での関数値を生成しておく。なお減衰曲線の大まかな形のみを知りたい場合など、計算精度をさほど要求しない用途においては、この手順を省略してもよい。
【００５３】
手順４．減衰定数aを適当に決める。手順１を計算する時点で伝搬定数ベクトルの虚部βは既に決めてあるので、これらと式（９）及び式（１０）とから遮断域下端と上端における試験関数ＥL（ｒ），ＥH（ｒ）が定まる。
手順５．式（１３）、式（１４）の８つの体積積分を数値積分で計算する。
手順６．式（２０）から係数Ｂ12，Ｂ23，Ｂ14，Ｂ34を計算する。
手順７．式（１８）、式（１９）からｕLとｕHの係数変数p, qを求める。
手順８．式（１５）から固有周波数の推定値wを得る。
手順９．減衰定数aを変えながら手順４から８を繰り返すことで、aとwの関係を示す曲線が得られる。
【００５４】
以上に示した計算の各手順は、電子計算機上にて適当なプログラムを記述して実行するか、もしくは既存の数値計算ソフトウェアのマクロ機能などを組み合わせて実行するのが現実的である。その場合の計算手順のフローチャートを図１に示した。また本発明の計算手法は、以上に列挙した数式の他に、計算機用のプログラムを記録した記憶媒体の形態で提供することができる。
以下では実際の計算例を示す。
【実施例１】
【００５５】
（２次元正方格子、air hole型周期構造）
簡単のため比誘電率、比透磁率をいずれもスカラー関数とし、比透磁率を至るところ1とする。格子定数がx,zのいずれの方向にもaであり、比誘電率の平方根であるところの屈折率分布が図２で表される構造中をxz面内で伝搬する波に注目する。すなわちこの構造は、屈折率が２．１の背景媒質２０１の中に、屈折率が１．５の円柱２０２がｘｚ面内で正方格子状に配列しているもので、円柱はｙ方向には無限に長く伸びている。屈折率がｘとｚの２方向にのみ周期的に変化しており、ｚ方向には一様であるため、このような構造は一般に２次元周期構造と呼ばれる。
【００５６】
Ｂｌｏｃｈ波数ベクトルyのx成分及びy成分を0とし、z成分の実部αzも0に固定する。z成分の虚部βzのみを変化させてＴＥ波、すなわち電界ベクトルが紙面に垂直なy成分のみからなる波の伝搬域での分散関係を全方位周期境界型のＦＤＴＤ法で計算した。結果を図３に示す。この計算における空間分割数はx,z各方向共に２０、時間ステップ数は６５５３６、縦軸である正規化された周波数の分解能は２．８８×１０-4である。また横軸の波数の分割数は１００とした。
【００５７】
なおここでは分散関係の形状を説明するために、時間ステップ数、波数の刻みをともに十分に細かくして計算した結果を示している。計算手順の詳細は例えば（非特許文献４）に開示されている。図３において、横軸が0及び0.5の位置を慣習に従ってそれぞれG点、Z点と呼ぶことにする。Z点では縦軸が0.28付近に第一のバンドギャップ、すなわち遮断域の存在が見て取れる。縦軸0.8付近の位置も遮断域である。一方、縦軸0.58から0.68の間で途切れているモード同士も遮断域を形成している。
【００５８】
他方、G点ではまず縦軸0.8付近に遮断域がある。縦軸0.54近傍では３つのモードが結合した遮断域が存在する。
以上を踏まえ、Z点第一遮断域の減衰定数を求める。図３に３０１及び３０２で示す点における電界のy成分Ｅyの固有関数を、同じく全方位周期境界型のＦＤＴＤ法によって計算した。結果を図４に示す。同図（ａ）は実部、（ｂ）は虚部である。周期構造の単位セルの分割数はx, z両方向ともに２０である。また時間ステップ数は６５５３６とした。時間軸上の刻みΔｔは下式で決まる上限値Δｔmaxの９９．９％の値とした。ＦＤＴＤ法では計算が安定に行なわれるためにΔｔ≦Δｔmaxであることが必要とされ、これをＣＦＬ条件という。
【００５９】
【数２９】
2007064911.tif 000031
ここでｖP, maxは解析空間中での電磁波の位相速度の最大値をあらわす。またΔｘ＝ａ／Ｎx，Δｚ＝ａ／Ｎzである。ＣＦＬ条件はわかりやすく言えば、「時間刻みΔｔは電磁波が空間グリッドの隣の点まで到達するのにかかる時間よりも短く設定しなければならない」という制約である。Ｅyから抽出された周期関数成分ｕyにx, z平面上で３次自然スプライン補間を施し、等価的に分割数１００×１００の詳細なグリッド上のデータを生成した。結果を図５に図５（ａ）、（ｂ）は補間前のｕyの実部及び虚部を、また同図（ｃ）、（ｄ）はそれらをｘとｚの各方向に元の５分の１の刻みにまで補間した後の分布を示す。次にこの空間補間後のｕyを用いて各数値積分を評価した。具体的には例えば式（１４）のＩd1及びＩd4の評価は次式に従って行った。
【００６０】
【数３０】
2007064911.tif 000032
【００６１】
【数３１】
2007064911.tif 000033
ただし（２６）において
【００６２】
【数３２】
2007064911.tif 000034
ｕL，ｕHはそれぞれ低周波側の端３０１, 高周波側の端３０２におけるｕyを表す。
【００６３】
次にαzを0から少しずつ増加させて式（１５）〜式（２０）から固有周波数を計算した。結果を図６に示す。固有値がほぼ重根となったところでαzの変化を打ち切っている。
【００６４】
１次元交互多層膜中を膜法線に沿って進む波の示す減衰曲線が、このαz−ω平面内で楕円を描くことは良く知られている。特にαzの最大値は遮断域の幅に比例する。この計算の構造は２次元周期を持つが、今着目しているZ点上の第一遮断域は一次元交互多層膜とほぼ同じく、主にz方向の前進波と後進波の結合によって生じているとみなせるので、図６のような上下ほぼ対称の楕円形の減衰曲線は合理的な結果である。
【００６５】
続いて偶モードの第二遮断域と奇モードの第一遮断域について同様の手順で計算を行った。これは図３に３０３で示した領域である。結果を図７に示す。この計算における空間分割数もx,z各方向共に２０、時間ステップ数は６５５３６とした。偶モードの試験関数の構築に使用したのは図中に７０１と７０２で示した点の固有関数である。同様に奇モードについては７０３と７０４における固有関数を用いた。
【００６６】
同様にして残りの各地点での減衰特性も計算した。結果を図８にまとめて示す。
【００６７】
また同じ構造中でのＴＭ波、すなわち磁界Hがy成分のみからなる波についても同様の手順で計算できる。この場合は磁界を使った表示式（２１）、（２２）、（２３）を用いるのが便利である。結果を図９に示す。ＴＥ波の場合と同様、各遮断域における減衰曲線が計算できているのがわかる。
【００６８】
以上の例では周期構造の格子定数を特に具体的な値に限定せず、aという変数のまま表現しているが、具体的なフォトニック結晶構造の遮断特性はそれぞれの結晶の寸法に応じて、例えばa=0.2mmなどとすれば直ちに求めることができる。またその他の波長域についても、例えばミリ波領域用の周期構造体ならa=1mm、周波数２．４ＧＨｚの超短波帯用ならa=6cmなどと適宜aに具体的な値を代入して図８、図９の縦横の軸を実次元の量に読み替えればただちに具体的な波長域に対する減衰特性のデータが得られる。
【実施例２】
【００６９】
ここでは本手法による計算に必要な計算時間を、遮断域での減衰定数を見積もるのに良く用いられる吸収境界終端型ＦＤＴＤ法の計算時間を比較することを念頭に置いた実施例を示す。吸収境界終端型ＦＤＴＤ法における解析空間の設定の概念図を図１０に示した。周期構造体の単位格子の寸法及び誘電率分布は実施例１に示したものと同じである。また構造はx方向には無限に周期が繰り返されている一方、遮断特性を求めようとするz方向にはある周期数の地点で打ち切られており、そこで一様媒質と接している。
【００７０】
この計算における空間分割数はx方向に２０、ｚ方向には１２２０である。ｚ方向の手前側の完全整合層の層数は１１０、向こう側の完全整合層の層数は４００である。時間ステップ数は３２７６８とした。吸収境界終端型ＦＤＴＤ法で実施例１に示したＴＥ波のZ点遮断域の減衰特性を計算した結果を図１１及び図１２に示す。それぞれは偶モードの第一遮断域、奇モードの第一遮断域についての結果である。
【００７１】
図において、本発明の手法で求めた減衰曲線を太い実線で、また吸収境界終端型ＦＤＴＤ法で計算した減衰曲線を細い点線で示している。吸収境界条件型のＦＤＴＤ法では減衰の弱い、伝搬域の端近くで曲線が大きく振動しているが、これは入射端面と出射端面の間での多重反射による共振の効果である。遮断域では両者はよく一致しており、特に図１１ではほとんど重なっている。
【００７２】
この吸収境界条件型のＦＤＴＤ法によれば、原理的には誘電率の周期構造がどの様な形であっても遮断域の減衰特性を求めることができる。ＦＤＴＤ法の計算量は（単位セルあたりのグリッド分割数）×（周期数）×（時間ステップ数）で決まる。上の例では単位セルあたりのグリッド点数は標準値（一辺１２〜２０）でよいが、長手方向には２０ないし５０程度の周期数が必要である。また励振波のエネルギーが反射波もしくは透過波として解析空間から出て行くまでに長い時間ステップ数が必要である。
【００７３】
これに対し本発明の方法で用いる全方位周期境界型のＦＤＴＤ法では、計算空間は縦横一周期分の単位セルでよい。また時間ステップ数は、時間波形にＦＦＴを施した際に、注目するモードの周波数を隣のモードのそれと分離できる程度であればよいので、前者の数分の１程度でよい。一方、なるべく正確に試験関数を合成するために、単位セルあたりのグリッド点数は前者より多くしなければならないことが多い。
【００７４】
具体的には単位格子の各辺の方向あたり４０点程度とする。以上を整理して、特にＺ点の第一遮断域の減衰特性に対する計算時間を比較したのが図１３である。それぞれの手法の詳細は次の通りである。なお計算にはクロック周波数2.8GHzのPentium(登録商標)-4チップを搭載した計算機を用い、プログラムはＣ＋＋言語で記述し、カーネルバージョン２．６．９のＬｉｎｕｘ上でバージョン３．３．２のｇｃｃ（ＧＮＵ  Ｃ  Ｃｏｍｐｉｌｅｒ）でコンパイルしたものを使用した。
【００７５】
<本手法>
１）第１手順（全方位周期境界条件型FDTD法による固有周波数の概算）
(空間グリッド数 40x40)x (時間ステップ数 8192)．
平均計算時間＝  ２．５９秒
２）第２手順（全方位周期境界条件型FDTD法による固有関数の計算）
(空間グリッド数 40x40)x (時間ステップ数 8192) ．
平均計算時間＝  ２．６３秒
３）第３手順（固有関数の補間と変分表現の計算）
（空間グリッド数  120x120）x(減衰定数の刻み 84)．
平均計算時間＝  １．８８秒
１）〜３）合計      約  ７．１秒
<吸収境界条件型FDTD法>
１）(横方向グリッド数 20)x(縦方向グリッド数 非PML 710 + (PML 510)x2)x(時間          ステップ数 32768) ．
平均計算時間＝  １７４．２秒
すなわち本手法によって、この例計算時間を吸収境界型のＦＤＴＤ法の２４分の１以下に短縮できた。
【実施例３】
【００７６】
ここでは本発明の計算手法を、実際の周期構造型光デバイスの一つであるフォトニック結晶偏光分離素子に適用する実施例を示す。
【００７７】
解析の対象とする周期構造の概念図を図１４に示す。この構造体は周期的な凹凸の列の刻まれた石英の基板１４０１上に、高屈折率媒質である屈折率２．１の五酸化タンタル(Ｔａ2Ｏ5)の層１４０２と、低屈折率媒質である屈折率１．５の石英の層１４０３とを交互に、ｒｆバイアス・スパッタ法によって成膜したものである。ｙ方向には屈折率の変化はない。屈折率がｘｚ面内で周期的に変化しているため、２次元周期構造である。このような構造体の積層表面から、ｚの正の方向に向かって光を入射する。
【００７８】
このとき光の電界がｙ成分のみを持つ波（ＴＥ波）と磁界がｙ成分のみを持つ波（ＴＭ波）の間に異方性が生じる。そしてある特定の光波長帯においては、ＴＥ波とＴＭ波のどちらか一方のみが遮断域となることがある。このとき遮断域となる偏波は周期構造体中に侵入できずに反射し、一方伝搬域となる偏波は構造体を透過する。このようにして光波長を適切に選ぶことで、本素子は反射型の偏光分離素子として機能する。この素子の製造方法及び動作原理の詳細は特許文献１及び非特許文献１に開示されている。
【００７９】
具体的な寸法などを規定して減衰量を計算した例を以下に示す。まず周期構造の単位格子１４０４の寸法はｘ方向、ｚ方向ともに２５０ｎｍとし、高屈折率媒質と低屈折率媒質の厚みはいずれも１２５ｎｍとした。また三角波状の周期的な凹凸の斜面の角度は３５°とした。この寸法でのＴＥ波のΓ−Ｚ方向の分散関係を計算した結果を図１５に示した。伝搬域は全方位周期境界型のＦＤＴＤ法で計算し、遮断域は本発明の手法で計算した。また実用上はあまり重要ではない、奇対称モードの分散曲線は表示していない。
【００８０】
ＴＥ波を反射させ、ＴＭ波を透過させるタイプの偏光分離素子では、ＴＥ波の遮断域の減衰量を知ることが必要である。そこで図１５に１５０１で示したＺ点の遮断域、及び１５０２で示したΓ点の遮断域をより詳細に計算した。それぞれの結果を図１６及び図１７に示す。ここでは減衰量の単位を、フォトニック結晶中を１周期、つまり五酸化タンタル層と石英層を１層ずつ通過する際のパワーの減衰量に変換して示している。
【００８１】
Ｚ点遮断域の中心周波数ではおおよそ１．８ｄＢ／周期、Γ点遮断域のうち上側の遮断域の中心周波数ではおおよそ４．２ｄＢ／周期の減衰が得られることが分かる。前者を利用する場合、例えばＴＥ波の透過率を−４０ｄＢまで減衰させるには、少なくとも２３周期程度の積層周期が必要であることがこの結果からわかる。このようにして本発明の手法を、素子の設計に利用することができる。
なお本発明の手法が取り扱うことのできる構造の寸法や屈折率は、当然ここで示した値に限定されるものではない。
【産業上の利用可能性】
【００８２】
本発明にかかる遮断域中の減衰定数の計算方法は、電波、光波、Ｘ線などを含む電磁波工学の分野において、透明もしくはほぼ透明な誘電体材料からなる周期構造体を用いた反射型素子の数値設計に幅広く適用することができる。特に既存の電磁界数値解析ソフトウェアの拡張機能として提供するなどの形態で広く産業的に利用することができる。
【図面の簡単な説明】
【００８３】
【図１】本発明の計算手法の流れをフローチャートとして示した図である。
【図２】本発明の計算手法が取り扱うことのできる誘電率周期構造の一例を示す概念図である。
【図３】図２の構造中を伝搬する電磁波のうち、ＴＥ波の分散曲線を示す図である。
【図４】図３の分散曲線において、点３０１と点３０２での電界の固有関数を示す図である。同図（ａ）は実部、（ｂ）は虚部を表す。
【図５】図３の分散曲線において、点３０１と点３０２での電界の固有関数に含まれる空間的周期関数成分を示す図である。同図（ａ）は実部、（ｂ）は虚部を表す。また同図（ｃ）、（ｄ）は実部と虚部それぞれのデータを空間的に５分の１の刻みで補間したものを示す。
【図６】図３の分散曲線において、点３０１と点３０２の間の減衰曲線を本発明の手法によって計算した結果を示す図である。
【図７】図３の分散曲線において、３０３で示すＺ点上の遮断域の減衰曲線を本発明の手法によって計算した結果を示す図である。
【図８】図２の構造中を伝搬するＴＥモードの電磁波の分散関係を、伝搬域中では全方位周期境界型ＦＤＴＤ法単独で、遮断域中では本発明の手法でそれぞれ計算した結果をまとめて示した図である。
【図９】図２の構造中を伝搬するＴＭモードの電磁波の分散関係を、伝搬域中では全方位周期境界型ＦＤＴＤ法単独で、遮断域中では本発明の手法でそれぞれ計算した結果をまとめて示した図である。
【図１０】遮断域の減衰特性を吸収境界型ＦＤＴＤ法で計算する場合の解析空間の設定方法の一例を示す図である。
【図１１】図３の分散曲線の点３０１と点３０２の間の減衰曲線を、本発明の手法と吸収境界型ＦＤＴＤ法のそれぞれで計算した結果を合わせて示した図である。
【図１２】図３の分散曲線に３０３で示す遮断域のうち、奇対称モードの減衰曲線を本発明の手法と吸収境界型ＦＤＴＤ法のそれぞれで計算した結果を示した図である。
【図１３】吸収境界型ＦＤＴＤ法と本発明の手法それぞれの計算時間の一例を示す図である。
【図１４】自己クローニング型フォトニック結晶による反射型偏光分離素子の構造の一例を示す図である。
【図１５】図１４の偏光分離素子において、ＴＥ波のΓ―Ｚ方向の分散曲線を示す図である。
【図１６】図１４の偏光分離素子において、１５０１に示すＺ点遮断域での減衰特性を示す図である。
【図１７】図１４の偏光分離素子において、１５０２に示すΓ点遮断域での減衰特性を示す図である。
【符号の説明】
【００８４】
２０１    屈折率２．１の背景媒質
２０２    屈折率１．５の円柱状媒質
３０１    Ｚ点上の第１モード
３０２    Ｚ点上の第２モード
３０３    Ｚ点上の遮断域
７０１    Ｚ点上の偶対称モードのひとつ
７０２    Ｚ点上の偶対称モードのひとつ
７０３    Ｚ点上の奇対称モードのひとつ
７０４    Ｚ点上の奇対称モードのひとつ
１４０１    周期的凹凸溝付き石英基板
１４０２    高屈折率層
１４０３    低屈折率層
１４０４    周期構造の単位格子
【図１】
2007064911.tif 000035
【図２】
2007064911.tif 000036
【図３】
2007064911.tif 000037
【図４】
2007064911.tif 000038
【図５】
2007064911.tif 000039
【図６】
2007064911.tif 000040
【図７】
2007064911.tif 000041
【図８】
2007064911.tif 000042
【図９】
2007064911.tif 000043
【図１０】
2007064911.tif 000044
【図１１】
2007064911.tif 000045
【図１２】
2007064911.tif 000046
【図１３】
2007064911.tif 000047
【図１４】
2007064911.tif 000048
【図１５】
2007064911.tif 000049
【図１６】
2007064911.tif 000050
【図１７】
2007064911.tif 000051