(57)【特許請求の範囲】
【請求項１】
金属部材を含む構造物に電磁波が入射した場合の電磁界を、ＦＤＴＤ法に従って解析するためのシミュレーション方法であって、
前記構造物の解析モデルが配置されるＦＤＴＤ領域の空間離散化条件を、前記金属部材の誘電率を基にコンピュータにより決定する離散化条件決定工程と、
前記離散化条件決定工程で決定された前記空間離散化条件に基づいて、前記解析モデルに前記電磁波が入射した場合の電磁界を、帰納的畳み込み法を用いた前記ＦＤＴＤ法に従って前記コンピュータにより算出する解析工程と、
を備え、
前記離散化条件決定工程では、前記帰納的畳み込み法において前記誘電率をドルーデの分散式で表すことで得られる帰納的関係式内の２つの変数ｃ１，Ｆ，ｃ２，Ｆであって、前記ＦＤＴＤ領域において前記電磁波が有する所定波長を表すための格子点の数及び前記ＦＤＴＤ領域における前記電磁波の解析用速度並びに前記金属部材の前記誘電率を含むプラズマ角周波数ωｐ、Ｆ及び衝突周波数νｃ，Ｆを用いて下記式（１），（２）で表される前記変数ｃ１，Ｆ，ｃ２，Ｆ：
【数１】
0005113572.tif 000059
【数２】
0005113572.tif 000060
において、
前記プラズマ角周波数ωｐ，Ｆ及び前記衝突周波数νｃ，Ｆに含まれる前記誘電率を、前記所定波長の前記電磁波に対する実測値又は実測値に基づいて決まる値とした場合に、前記変数ｃ１，Ｆ，ｃ２，Ｆが、前記ＦＤＴＤ法における安定化条件に含まれる下記式（３）及び式（４）を満たすと共に、
【数３】
0005113572.tif 000061
【数４】
0005113572.tif 000062
前記解析用速度が０．７より小さくなるように、前記格子点の数及び前記解析用速度を、前記コンピュータによって決定することを特徴とするシミュレーション方法。
【請求項２】
実空間における前記所定波長をλ、前記電磁波の速度をｖ及び前記所定波長に対応する角周波数をωλとし、前記ＦＤＴＤ領域における前記格子点の数をλＦ、前記解析用速度をｖＦとし、前記角周波数ωλに対する前記誘電率の実部及び虚部の値をεｒ及びεｉとし、前記角周波数ωλに対応する前記ＦＤＴＤ領域における解析用角周波数ωＦを、
【数５】
0005113572.tif 000063
としたとき、前記衝突周波数νｃ，Ｆ及び前記プラズマ角周波数ωｐ、Ｆは、
【数６】
0005113572.tif 000064
【数７】
0005113572.tif 000065
として与えられることを特徴とする請求項１に記載のシミュレーション方法。
【請求項３】
前記解析工程では、
前記電磁波におけるＴＥモードに対しては、マクスウェル方程式から導かれる波動方程式に前記帰納的畳み込み法を適用して得られる第１の電磁界解析用の式を、前記帰納的関係式を利用して前記コンピュータにより解くことを特徴とする請求項１又は２に記載のシミュレーション方法。
【請求項４】
前記解析工程では、
前記電磁波におけるＴＭモードに対しては、マクスウェル方程式に前記帰納的畳み込み法を適用して得られる第２の電磁界解析用の式を、前記帰納的関係式を利用して前記コンピュータにより解くことを特徴とする請求項１〜３の何れか一項に記載のシミュレーション方法。
【請求項５】
前記解析工程では、
前記電磁波がＴＭモード及びＴＥモードを共に含む場合には、前記ＴＥモードに対しては、マクスウェル方程式から導かれる波動方程式に前記帰納的畳み込み法を適用して得られる第１の電磁界解析用の式を、前記帰納的関係式を利用して解き、前記ＴＭモードに対しては、マクスウェル方程式に前記帰納的畳み込み法を適用して得られる第２の電磁界解析用の式を、前記帰納的関係式を利用して解き、前記ＴＭモード及び前記ＴＥモードに対して得られた電磁界に基づいて前記ＦＤＴＤ領域における電磁界を算出する、ことを特徴とする請求項１又は２に記載のシミュレーション方法。
【請求項６】
金属部材を含む構造物に電磁波が入射した場合の電磁界を、ＦＤＴＤ法に従って解析するためのシミュレーションをコンピュータに実行させるためのシミュレーションプログラムであって、
前記構造物の解析モデルが配置されるＦＤＴＤ領域の空間離散化条件を、前記金属部材の誘電率を基にコンピュータにより決定する離散化条件決定工程と、
前記離散化条件決定工程で決定された前記空間離散化条件に基づいて、前記解析モデルに前記電磁波が入射した場合の電磁界を、帰納的畳み込み法を用いた前記ＦＤＴＤ法に従って前記コンピュータにより算出する解析工程と、
を前記コンピュータに実行させ、
前記離散化条件決定工程では、前記帰納的畳み込み法において前記誘電率をドルーデの分散式で表すことで得られる帰納的関係式内の２つの変数ｃ１，Ｆ，ｃ２，Ｆであって、前記ＦＤＴＤ領域において前記電磁波の所定波長を表すための格子点の数及び前記ＦＤＴＤ領域における前記電磁波の解析用速度並びに前記金属部材の前記誘電率を含むプラズマ角周波数ωｐ及び衝突周波数νｃ, Fを用いて下記式（８），（９）で表される前記変数ｃ１，Ｆ，ｃ２，Ｆ：
【数８】
0005113572.tif 000066
【数９】
0005113572.tif 000067
において、
前記プラズマ角周波数ωｐ，Ｆ及び前記衝突周波数νｃ，Ｆに含まれる前記誘電率を、前記所定波長の前記電磁波に対する実測値又は実測値に基づいて決まる値とした場合に、前記変数ｃ１，Ｆ，ｃ２，Ｆが、前記ＦＤＴＤ法における安定化条件に含まれる下記式（１０）及び式（１１）を満たすと共に、
【数１０】
0005113572.tif 000068
【数１１】
0005113572.tif 000069
前記解析用速度が０．７より小さくなるように、前記格子点の数及び前記解析用速度を、前記コンピュータに決定せしめることを特徴とするシミュレーションプログラム。
【請求項７】
実空間における前記電磁波の前記所定波長をλ、前記電磁波の速度をｖ及び前記所定波長に対応する角周波数をωλとし、前記ＦＤＴＤ領域における前記格子点の数をλＦ、前記解析用速度をｖＦとし、前記角周波数ωλに対する前記誘電率の実部及び虚部の値をεｒ及びεｉとし、前記角周波数ωλに対応する前記ＦＤＴＤ領域における解析用角周波数ωＦを、
【数１２】
0005113572.tif 000070
としたとき、前記衝突周波数νｃ，Ｆ及び前記プラズマ角周波数ωｐ，Ｆは、
【数１３】
0005113572.tif 000071
【数１４】
0005113572.tif 000072
として与えられることを特徴とする請求項６に記載のシミュレーションプログラム。
【請求項８】
前記解析工程では、
前記電磁波におけるＴＥモードに対しては、マクスウェル方程式から導かれる波動方程式に前記帰納的畳み込み法を適用して得られる第１の電磁界解析用の式を、前記帰納的関係式を利用して前記コンピュータに解かさしめることを特徴とする請求項６又は７に記載のシミュレーションプログラム。
【請求項９】
前記解析工程では、
前記電磁波におけるＴＭモードに対しては、マクスウェル方程式に前記帰納的畳み込み法を適用して得られる第２の電磁界解析用の式を、前記帰納的関係式を利用して前記コンピュータに解かさしめることを特徴とする請求項６〜８の何れか一項に記載のシミュレーションプログラム。
【請求項１０】
前記解析工程では、
前記電磁波がＴＭモード及びＴＥモードを共に含む場合には、前記ＴＥモードに対しては、マクスウェル方程式から導かれる波動方程式に前記帰納的畳み込み法を適用して得られる第１の電磁界解析用の式を、前記帰納的関係式を利用して前記コンピュータに解かさしめて、前記ＴＭモードに対しては、マクスウェル方程式に前記帰納的畳み込み法を適用して得られる第２の電磁界解析用の式を、前記帰納的関係式を利用して前記コンピュータに解かさしめて、前記ＴＭモード及び前記ＴＥモードに対して得られた電磁界に基づいて前記ＦＤＴＤ領域における電磁界を算出する、ことを特徴とする請求項６又は７に記載のシミュレーションプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、ＦＤＴＤ法を利用したシミュレーション方法及びシミュレーションプログラムに関するものである。
【背景技術】
【０００２】
金属部材を含む構造物の一例としては、複数の金属ストリップ（金属細線）がほぼ平行に配列されて構成される金属回折格子がある。このような金属回折格子は、ワイヤグリッド偏光子、またはディスプレイ関連光学装置における偏光コントローラとしてしばしば利用される。各金属ストリップの幅または周期が、金属回折格子に入射する光の波長に比べて小さい場合、回折格子特性は、金属ストリップの幅および周期、入射ビームの偏光および入射角に依存する。そのため、所望の回折格子特性が得られるように、金属ストリップの幅及び周期等を入射ビームの特性などに応じて設計する必要がある。
【０００３】
金属部材を含む構造物へ電磁波が入射した場合の電磁界の解析には、解析的な手法を利用したシミュレーションが適用される。例えば、厳密結合波解析（RCWA :rigorous coupled wave analysis）などの解析的方法が無限の大きさを仮定した回折格子に対して使用されることがある。しかしながら、ＲＣＷＡに例示される解析的方法では、（１）有限の大きさを仮定すべき構造又は有限の幅を有する入射ビームの場合には使用することができない、（２）不規則な境界面を扱うことができない、（３）媒質の誘電率を複素数としているが、その実部は入射ビームの角周波数に依存していないことを仮定している、等の問題点がある。よって、例えば、上述した配列方向の長さが有限である回折格子、すなわち、実際の回折格子により近い形態を正確に解析できず、また、金属部材が、銀、アルミニウム、金などの貴金属からなる場合には、可視領域においてシミュレーションの結果で得られる電磁界分布と、実空間での電磁界分布との差が大きくなる。
【０００４】
一方、上記のような問題点を解決可能なシミュレーション手法として、不規則な境界面を有しており、有限の大きさの仮定が有効な構造物又は無限の仮定が可能な構造物のどちらにも適用可能なＦＤＴＤ法（Finite Differential Time Domain method）が知られている（例えば、特許文献１，２参照）。ＦＤＴＤ法は、ＦＤＴＤ領域を複数のセルに分割して、マクスウェル方程式内の空間および時間微分を有限差分によって近似することでマクスウェル方程式を直接解くものであり、様々な金属構造を解析するために使用されてきている。そして、ＦＤＴＤ法の場合、例えば、デバイ（Debye）、ドルーデ(Drude)、またはローレンツ（Lorentz）など、誘電率についての様々な解析的モデルが、周波数領域畳み込みを使用することによって、ＦＤＴＤに組み込むことができる。このような畳込みを利用した方法は、帰納的畳み込み（Recursive Convolution）法（以下、単にＲＣ法ともいう）として知られている。
【特許文献１】特開２０００−１０５２５９号公報
【特許文献２】特開２０００−２２７４５０号公報
【発明の開示】
【発明が解決しようとする課題】
【０００５】
しかし、従来のＦＤＴＤ法を利用したシミュレーションでは、空間的離散化の値等の条件設定によっては、シミュレーションが不安定になる傾向にあった。
【０００６】
そこで、本発明は、ＦＤＴＤ法によるシミュレーションをより安定して実施可能なシミュレーション方法及びシミュレーションプログラムを提供すること目的とする。
【課題を解決するための手段】
【０００７】
本発明のシミュレーション方法は、金属部材を含む構造物に電磁波が入射した場合の電磁界を、ＦＤＴＤ法に従って解析するためのシミュレーション方法であって、構造物の解析モデルが配置されるＦＤＴＤ領域の空間離散化条件を、金属部材の誘電率を基にコンピュータにより決定する離散化条件決定工程と、離散化条件決定工程で決定された空間離散化条件に基づいて、解析モデルに電磁波が入射した場合の電磁界を、帰納的畳み込み法を用いたＦＤＴＤ法に従ってコンピュータにより算出する解析工程と、を備え、離散化条件決定工程では、帰納的畳み込み法において上記誘電率をドルーデの分散式で表すことで得られる帰納的関係式内の２つの変数ｃ１，ｃ２であって、ＦＤＴＤ領域において所定波長を表すための格子点の数及びＦＤＴＤ領域における電磁波の解析用速度並びに金属部材の上記誘電率を含むプラズマ角周波数ωｐ，Ｆ及び衝突周波数νｃ，Ｆを用いて下記式（１），（２）で表される変数ｃ１，Ｆ，ｃ２，Ｆ：
【数１】
0005113572.tif 000002
【数２】
0005113572.tif 000003
において、プラズマ角周波数ωｐ，Ｆ及び衝突周波数νｃ，Ｆに含まれる誘電率を、所定波長の電磁波に対する実測値又は実測値に基づいて決まる値とした場合に、変数ｃ１，Ｆ，ｃ２，Ｆが、ＦＤＴＤ法における計算の安定化条件としての下記式（３）及び式（４）を満たすと共に、
【数３】
0005113572.tif 000004
【数４】
0005113572.tif 000005
解析用速度が０．７より小さくなるように、格子点の数及び解析用速度をコンピュータによって決定することを特徴とする。
【０００８】
また、本発明に係るシミュレーションプログラムは、金属部材を含む構造物に電磁波が入射した場合の電磁界を、ＦＤＴＤ法に従って解析するためのシミュレーションをコンピュータに実行させるためのシミュレーションプログラムであって、構造物の解析モデルが配置されるＦＤＴＤ領域の空間離散化条件を、金属部材の誘電率を基にコンピュータにより決定する離散化条件決定工程と、離散化条件決定工程で決定された空間離散化条件に基づいて、解析モデルに電磁波が入射した場合の電磁界を、帰納的畳み込み法を用いたＦＤＴＤ法に従ってコンピュータにより算出する解析工程と、をコンピュータに実行させ、離散化条件決定工程では、帰納的畳み込み法において上記誘電率をドルーデの分散式で表すことで得られる帰納的関係式内の２つの変数ｃ１，Ｆ，ｃ２，Ｆであって、ＦＤＴＤ領域において所定波長を表すための格子点の数及びＦＤＴＤ領域における電磁波の解析用速度並びに金属部材の上記誘電率を含むプラズマ角周波数ωｐ，Ｆ及び衝突周波数νｃ，Ｆを用いて下記式（５），（６）で表される前記変数ｃ１，Ｆ，ｃ２，Ｆ：
【数５】
0005113572.tif 000006
【数６】
0005113572.tif 000007
において、プラズマ角周波数ωｐ，Ｆ及び衝突周波数νｃ，Ｆに含まれる上記誘電率を、所定波長の電磁波に対する実測値又は実測値に基づいて決まる値とした場合に、変数ｃ１，Ｆ，ｃ２，Ｆが、ＦＤＴＤ法における計算の安定化条件としての下記式（７）及び式（８）を満たすと共に、
【数７】
0005113572.tif 000008
【数８】
0005113572.tif 000009
解析用速度が０．７より小さくなるように、格子点の数及び解析用速度をコンピュータに決定せしめることを特徴とする。
【０００９】
上記本発明に係るシミュレーション方法及びシミュレーションプログラムでは、所定波長λを有する電磁波に対する金属部材の誘電率の実測値及び実測値に基づいて決まる値を使用して、ＦＤＴＤ法における計算の安定化条件を満たすようにＦＤＴＤ領域の空間離散化を図っている。そのため、より正確なシミュレーションを安定して実施することが可能である。
【００１０】
また、本発明に係るシミュレーション方法では、実空間における電磁波の所定波長をλ、電磁波の速度をｖ及び所定波長に対応する角周波数をωλとし、ＦＤＴＤ領域における格子点の数をλＦ、解析用速度をｖＦとし、角周波数ωλに対する誘電率の実部及び虚部の値をεｒ及びεｉとし、角周波数ωλに対応するＦＤＴＤ領域における解析用角周波数ωＦを、
【数９】
0005113572.tif 000010
としたとき、衝突周波数νｃ，Ｆ及びプラズマ角周波数ωｐ，Ｆは、
【数１０】
0005113572.tif 000011
【数１１】
0005113572.tif 000012
として与えられることが好ましい。
【００１１】
同様に、本発明に係るシミュレーションプログラムでは、実空間における電磁波の所定波長をλ、電磁波の速度をｖ及び所定波長に対応する角周波数をωλとし、ＦＤＴＤ領域における格子点の数をλＦ、解析用速度をｖＦとし、角周波数ωλに対する誘電率の実部及び虚部の値をεｒ及びεｉとし、角周波数ωλに対応するＦＤＴＤ領域における解析用角周波数ωＦを、
【数１２】
0005113572.tif 000013
としたとき、衝突周波数νｃ，Ｆ及びプラズマ角周波数ωｐ，Ｆは、
【数１３】
0005113572.tif 000014
【数１４】
0005113572.tif 000015
として与えられることが好ましい。
【００１２】
この場合、νｃ、Ｆ及びωＰ，Ｆは、λＦ及びｖＦの関数となる。そのため、変数ｃ１，Ｆ及びｃ２，Ｆが満たす条件と共に、ｖＦが０．７より小さい値をとるように、λＦ及びｖＦを決定することができる。
【００１３】
本発明に係るシミュレーション方法において、上記解析工程では、電磁波におけるＴＥモードに対しては、マクスウェル方程式から導かれる波動方程式に上記帰納的畳み込み法を適用して得られる第１の電磁界解析用の式を、帰納的関係式を利用してコンピュータにより解くことが好適である。同様に、本発明に係るシミュレーションプログラムにおいて、上記解析工程では、電磁波におけるＴＥモードに対しては、マクスウェル方程式から導かれる波動方程式に上記帰納的畳み込み法を適用して得られる第１の電磁界解析用の式を、帰納的関係式を利用してコンピュータに解かさしめることが好適である。
【００１４】
この場合、波動方程式を利用しているため計算量を低減することが可能である。
【００１５】
また、本発明に係るシミュレーション方法における解析工程では、電磁波におけるＴＭモードに対しては、マクスウェル方程式に帰納的畳み込み法を適用して得られる第２の電磁界解析用の式を、帰納的関係式を利用して解く態様とすることもできる。同様に、本発明に係るシミュレーションプログラムにおける解析工程では、電磁波におけるＴＭモードに対しては、マクスウェル方程式に帰納的畳み込み法を適用して得られる第２の電磁界解析用の式を、帰納的関係式を利用してコンピュータに解かさしめる態様とすることもできる。
【００１６】
更に、本発明に係るシミュレーション方法における上記解析工程では、電磁波がＴＭモード及びＴＥモードを共に含む場合には、ＴＥモードに対しては、マクスウェル方程式から導かれる波動方程式に上記帰納的畳み込み法を適用して得られる第１の電磁界解析用の式を、帰納的関係式を利用して解き、ＴＭモードに対しては、マクスウェル方程式に帰納的畳み込み法を適用して得られる第２の電磁界解析用の式を、帰納的関係式を利用して解き、ＴＭモード及びＴＥモードに対して得られた電磁界に基づいてＦＤＴＤ領域における電磁界を算出する態様とすることもできる。同様に、本発明に係るシミュレーションプログラムにおける上記解析工程では、電磁波がＴＭモード及びＴＥモードを共に含む場合には、ＴＥモードに対しては、マクスウェル方程式から導かれる波動方程式に上記帰納的畳み込み法を適用して得られる第１の電磁界解析用の式を、帰納的関係式を利用してコンピュータに解かさしめて、ＴＭモードに対しては、マクスウェル方程式に帰納的畳み込み法を適用して得られる第２の電磁界解析用の式を、帰納的関係式を利用してコンピュータに解かさしめて、ＴＭモード及びＴＥモードに対して得られた電磁界に基づいてＦＤＴＤ領域における電磁界を算出する態様とすることができる。
【発明の効果】
【００１７】
本発明によれば、ＦＤＴＤ法に従ってシミュレーションを安定して実施することができる。
【発明を実施するための最良の形態】
【００１８】
以下、図面を参照して本発明に係るシミュレーション方法及びシミュレーションプログラムの実施形態について説明する。なお、図面の説明においては同一要素には同一符号を付し、重複する説明を省略する。また、図面の寸法比率等は、説明のものと必ずしも一致していない。
【００１９】
本発明に係るシミュレーション方法及びシミュレーションプログラムは、金属部材を含む構造物に電磁波が入射した場合における電磁界を、ＦＤＴＤ法（Finite Difference Time Domain Method）に従って解析するためのものであり、上記構造物の設計などに好適に利用される。本発明に係る方法及びプログラムを利用した解析対象としての構造物の一例としては金属回折格子が挙げられる。本明細書では、特に断らない限り、光学素子は金属回折格子として説明する。
【００２０】
図１は、金属回折格子の模式図である。金属回折格子１０は、一方向（図１では、紙面に垂直な方向）に延在している複数の金属部材としての金属細線（金属ストリップ）１１が、金属細線１１の延在方向に対して略直交する方向に周期Λで配置されてなる回折格子である。説明の簡略化のために、図１に示すように、金属細線１１の配列方向をｘ軸方向とし、金属細線１１の延在方向をｚ軸方向とし、ｘ軸方向及びｚ軸方向に直交する方向をｙ軸方向とする。金属細線１１の材料としては、銀、金、アルミニウムなどが例示される。
【００２１】
金属回折格子１０の周期Λは、電磁波としての入射ビーム２０が有する所定波長λに対して十分小さく、例えば、５λ以下である。金属回折格子１０に入射ビーム２０が入射すると、入射ビーム２０は金属回折格子１０により回折される。周期Λが所定波長λに対して十分小さい場合、金属回折格子１０に入射ビーム２０が入射すると、ほぼ０次の回折光が生じることになる。以下の説明では、図１に示したように、金属回折格子１０に対して入射ビーム１０の入射側と反対側に回折される光を透過ビーム２１とも称し、入射ビーム２０の入射側と同じ側に回折される光を反射ビーム２２とも称す。なお、所定波長λは、電磁スペクトルにおける可視領域の波長であるとすることができる。
【００２２】
金属回折格子１０の特性は、周期Λ並びに金属細線１１の幅（ｘ軸方向の長さ）ｗ及び厚さ（ｙ軸方向の長さ）ｈや、回折格子法線Ｎと入射ビーム２０とのなす角である入射角θｉｎや、入射ビーム２０の偏光状態などに依存する。そのため、所望の回折格子特性を実現する金属回折格子１０を設計するために本発明に係るシミュレーション方法及びプログラムが有効である。特に、金属回折格子１０のｘ軸方向の幅Ｗが３０λより小さい場合に本発明が有効である。これは、本発明に係るシミュレーション方法及びシミュレーションプログラムが、その手法として、有限な大きさを持つ物体を扱うときに有効なＦＤＴＤ法を採用していることによる。
【００２３】
図２は、図１に示した金属回折格子のＦＤＴＤ領域における解析モデルの模式図である。ＦＤＴＤ法では、図２に示すように、ＦＤＴＤ領域１００を複数のセルに分割することによって、マクスウェル方程式を空間的及び時間的に差分化して直接数値的に解く手法である。本実施形態では、図２に示すように、解析モデルとしての金属回折格子に対しても図１に示した金属回折格子１０と同様の符号を付すこととする。また、図２中のハッチングは、セルで分割されたＦＤＴＤ領域１００内において金属細線１１を示すための便宜的なものである。
【００２４】
通常、図１及び図２に示した金属回折格子１０のように一方向（図１，２中ではｚ軸方向）に延在しており、その延在方向に変化のないものを解析対象とする場合、２次元問題として扱うことができる。そこで、２次元のＦＤＴＤ法のアルゴリズムについて説明する。
【００２５】
ＦＤＴＤ法では、以下のマクスウェル方程式を利用する。
【数１５】
0005113572.tif 000016
【数１６】
0005113572.tif 000017
μは透磁率である。また、Ｅは電界、Ｈは磁界である。更に、「∇」は、ｘｕ及びｙｕをｘ軸方向及びｙ軸方向の単位ベクトルとしたとき、ｘｕ(∂／ｘ)＋ｙｕ(∂／ｙ)で定義される微分演算子である。
【００２６】
式（１６）中のＤは電束密度を表しており、次式で与えられる。
【数１７】
0005113572.tif 000018
式（１７）中のε（ｒ，ω）は金属細線１１の誘電率である。
【００２７】
時間領域における電束密度Ｄは、式（１７）の両辺でフーリエ変換をとり、畳み込み定理を適用することによって次式で表される。
【数１８】
0005113572.tif 000019
式（１８）において、ε（ｔ）＝Ｆ｛ε（ω）｝であり、Ｆはフーリエ変換を表す。
【００２８】
誘電率ε（ω）が角周波数ωのある特殊な関数で表されるとき、式（１８）の畳み込み積分は帰納的に計算可能である。これはＲＣ法として知られている。
【００２９】
本実施形態では、金属細線１１の誘電率ε（ω）を表す関数として、ドルーデの一次近似式（ドルーデの分散式）を採用する。この場合、誘電率ε（ω）は式（１９）のように表される。
【数１９】
0005113572.tif 000020
式（１９）において、ωｐはプラズマ周波数、νｃは衝突周波数である。ε∞は、周波数が無限大のときの誘電率であり１に等しい。χ（ω）は電気比感受率である。
【００３０】
χ（ω）をフーリエ変換したもの、すなわち、Ｆ｛χ（ω）｝をχ（ｔ）と表した場合、χ（ｔ）は次式で表される。
【数２０】
0005113572.tif 000021
Ｕ（ｔ）は、時間の単位ステップ関数であって、ｔ＝０の場合、Ｕ（ｔ）＝０であり、ｔ＞０の場合、Ｕ（ｔ）＝１である。
【００３１】
Δｔを時間ステップ間隔としたとき、時刻ｔ＝ｎΔｔ（ｎは０以上の整数）において、式（１８）の畳み込み積分を評価した場合、
【数２１】
0005113572.tif 000022
が得られる。より具体的には、式（２１）は、区間［ｍΔｔ，（ｍ＋１）Δｔ］にわたってＥを定数としてとることによって、畳み込み積分から得られる。なお、ｍは０以上の整数である。
【００３２】
式（２１）は、以下のように書き表すことができる。
【数２２】
0005113572.tif 000023
ここで、Ψは、累積変数であり下記式（２３）が成り立つ。
【数２３】
0005113572.tif 000024
【００３３】
式（２３）は、
【数２４】
0005113572.tif 000025
と表すことができる。この式（２４）において、
【数２５】
0005113572.tif 000026
である。式（２５）より、
【数２６】
0005113572.tif 000027
が成立する。式（２６）において、
【数２７】
0005113572.tif 000028
【数２８】
0005113572.tif 000029
である。また、Ψ０＝０、Ψ１＝Ｅ１χ０、Ψ２＝Ｅ２χ０＋ｃ１Ｅ１＋ｃ２Ψ１とする。
【００３４】
式（２４）〜式（２８）より、累積変数Ψに関する以下の帰納的関係式が得られる。
【数２９】
0005113572.tif 000030
【数３０】
0005113572.tif 000031
式（２９）及び式（３０）では、ε∞＝１、Δｔ＝１である。また、式（２９）及び式（３０）中、χ０は式（２５）より、
【数３１】
0005113572.tif 000032
で与えられる。更に、ｎ≦０であるすべてのｎについてΨの値は０とする。
【００３５】
なお、後述するように、ＦＤＴＤ法を実施する際には、式（２９）及び式（３０）に含まれておりωｐ及びνｃを用いて表されるｃ１，ｃ２は、ωｐ及びνｃに対応するＦＤＴＤ領域１００におけるプラズマ角周波数ωｐ，Ｆ及び衝突周波数νｃ，Ｆの関数としてのｃ１，Ｆ及びｃ２，Ｆを採用する。
【００３６】
通常、２次元においては、任意の偏光状態にある入射ビーム（電磁波）２０は、２つの直交する偏光状態、すなわち、ＴＭモードとＴＥモードとに分けることが可能である。そこで、以下、ＴＭモードとＴＥモードとに分けてＦＤＴＤ法のアルゴリズムについて説明する。
【００３７】
[ＴＭモード]
ＴＭモードでは、式（１５）及び式（１６）を中心差分を利用して差分化して解く。
【００３８】
式（１６）中の時間に関する一次導関数の差分形式として式（３２）を利用する。
【数３２】
0005113572.tif 000033
この場合、式（１６）に等価な差分形式の方程式は、以下のように表される。
【数３３】
0005113572.tif 000034
式（３３）中のＤ１は∇に対応する２次の中心差分形式の演算子である。また、式（３３）中のΨｎ＋１−Ψｎは、式（２９）により与えられる。
【００３９】
ここでは、式（１６）に等価な差分形式の方程式を示したが、式（１５）に関しては、従来のＦＤＴＤ法における差分形式である以下の式を用いる。
【数３４】
0005113572.tif 000035
【００４０】
[ＴＥモード]
次にＴＥモードの場合について説明する。
【００４１】
ＴＥモードではマクスウェル方程式である式（１５）及び式（１６）から導かれる次の波動方程式を利用する。
【数３５】
0005113572.tif 000036
本実施形態では∇・Ｅ＝０を仮定するが、高い伝導性を有する金属の場合、それらは自由電荷を長い間保持しないので、この仮定は妥当である。
【００４２】
式（１６）の時間に関する二次導関数の差分形式として式（３６）を利用する。
【数３６】
0005113572.tif 000037
この場合、帰納的畳み込みを適用した式（３５）に等価な差分形式の方程式は、以下のように表される。
【数３７】
0005113572.tif 000038
式（３７）中のＤ１２は∇２に対応する２次の中心差分形式の演算子であり、ｖは位相速度である。また、式（３７）中のΨｎ＋１＋Ψｎ−１−２Ψｎは式（３０）から与えられる。
【００４３】
本実施形態では、入射ビーム２０がＴＭモードで与えられる場合、すなわち、入射ビーム２０がＴＭ波である場合には、電磁界解析用の式（第２の電磁界解析用の式）である式（３３）及び式（３４）と、式（２９）とを利用して電界Ｅ及び磁界Ｈを算出する。入射ビーム２０がＴＥモードで与えられる場合、すなわち入射ビーム２０がＴＥ波である場合には、電磁界解析用の式（第１の電磁界解析用の式）である式（３７）と、式（３０）とを利用して電界Ｅを算出する。
【００４４】
次に、本実施形態におけるＦＤＴＤ領域における空間的な離散化条件の決定方法について説明する。
【００４５】
ＦＤＴＤ法では、図２に示したように、ＦＤＴＤ領域１００に複数のセルからなるグリッドを適用する。本実施形態のシミュレーション法で解くＦＤＴＤ法では、アルゴリズムの安定性に基づいてＦＤＴＤ領域１００の離散化条件を決定する。
【００４６】
式（２９）及び式（３０）より、アルゴリズムの安定性のためにｃ１及びｃ２は次の条件を満たすことが必要である。
【数３８】
0005113572.tif 000039
【数３９】
0005113572.tif 000040
【００４７】
ｃ１，ｃ２は式（２７）及び式（２８）に示すようにωｐ及びνｃを用いて表される。
【００４８】
一方、式（１９）より、誘電率ε（ω）の実部及び虚部をεｒ（ω）及びεｉ（ω）とすると、
【数４０】
0005113572.tif 000041
【数４１】
0005113572.tif 000042
が成り立つ。
【００４９】
よって、ωｐ及びνｃは次のように表される。
【数４２】
0005113572.tif 000043
【数４３】
0005113572.tif 000044
【００５０】
本実施形態では、入射ビーム２０が有する所定波長λに対応する角周波数をωλとしたとき、式（４０）及び式（４１）内のεｒ（ω）及びεｉ（ω）として、角周波数ωλに対する実測値を採用する。
【００５１】
所定波長λに対して誘電率の実測値が得られていない場合には、ドルーデ分散を利用して実測されている誘電率の局所的なフィッティングから上記所定波長λに対する誘電率、すなわち、εｒ（ω）及びεｉ（ω）を得る。上記局所フィッティングは、所定波長λに最も近く誘電率が実測されている波長に対する誘電率の実測値を通るようにドルーデのε−ωプロットを実施することで得ることができる。
【００５２】
例えば、波長５００ｎｍ（角周波数３．７７×１０１５／ｓｅｃ）での誘電率の実測値が得られていない一方、波長４９６ｎｍ（又は角周波数３．８×１０１５／ｓｅｃ）と、波長５５１．１１１ｎｍ（角周波数３．４２×１０１５／ｓｅｃ）に対してそれぞれ誘電率が得られているとする。この場合、波長５００ｎｍでの誘電率は、角周波数３．８×１０１５／ｓｅｃに対する誘電率ε＝−９．５６４＋ｉ０．３０９をとおるドルーデのε−ωプロットから計算する（図３参照）。
【００５３】
より正確な誘電率を得るために、所定波長λに対して高波長側及び短波長側においてそれぞれ所定波長λに最も近く誘電率が実測されている波長に対する誘電率を基に前述したドルーデのε−ωプロットを実施し、所定波長λの誘電率としてそれぞれ得られる２つの推定値の平均値を算出して推定値とすることが好ましい。金属回折格子に入射させる光として例えばＬＥＤから出力された光等を利用される際には、多くの場合、ＬＥＤから出力される光の波長に対して誘電率が実測されていない場合があるため、上記のような算出方法は好ましい。なお、以下の説明では、実測値という場合には、上記ドルーデのε−ωプロットを実施して決まる値を含めるものとする。
【００５４】
また、本実施形態では、式（４２）及び式（４３）における角周波数ωをＦＤＴＤ領域における角周波数ωＦに置き換えて、プラズマ角周波数ωｐ及び衝突周波数νｃをＦＤＴＤ領域１００におけるプラズマ角周波数ωｐ，Ｆ及び衝突周波数νｃ，Ｆに対応させる。この場合、プラズマ角周波数ωｐ，Ｆ及び衝突周波数νｃ，Ｆは以下のように表される。
【数４４】
0005113572.tif 000045
【数４５】
0005113572.tif 000046
式（４４）及び式（４５）中のεｒ及びεｉは角周波数ωλに対する実測値であり定数である。
【００５５】
そして、本実施形態では、実空間での角周波数ωλと、式（４４）及び式（４５）中のＦＤＴＤ領域１００での角周波数（解析用角周波数）ωＦとを次式によって対応づけることによって、ωｐ，Ｆ及びνｃ，Ｆを自由パラメータとしている。
【数４６】
0005113572.tif 000047
式（４６）中のｖは入射ビーム２０の実空間での位相速度である。そして、λＦは、ＦＤＴＤ領域１００において波長λを表すのに必要とする格子点１１０の数であり、ｖＦはＦＤＴＤ領域１００における入射ビーム２０の位相速度（解析用速度）である。
【００５６】
上記のように実空間におけるωｐ及びνｃに対応してＦＤＴＤ領域１００におけるωｐ，Ｆ及びνｃ，Ｆを規定した場合、式（３８），式（３９）に示す安定化条件は、
【数４７】
0005113572.tif 000048
【数４８】
0005113572.tif 000049
と表されることになる。
【００５７】
なお、ｃ１，Ｆ及びｃ２，Ｆは、ｃ１，ｃ２におけるωｐ及びνｃをωｐ，Ｆ及びνｃ，Ｆで置き換えたものであり、以下の式で表される。
【数４９】
0005113572.tif 000050
【数５０】
0005113572.tif 000051
また、アルゴリズムの安定性として、式（４９）及び式（５０）と共に、ＦＤＴＤ法で通常知られている以下の関係式を満たす必要がある。
【数５１】
0005113572.tif 000052
【００５８】
前述したようにεｒ及びεｉとして角周波数ωλに対する所定の値を利用すると共に、式（４６）を採用することによって、ωＰ，Ｆ及びνｃ，Ｆを自由パラメータとした場合、式（４９）及び式（５０）並びに式（５１）を満たすようにλＦ及びｖＦを決定する。λＦは、所定波長λを表すための格子点１１０の数であるため、空間的離散化条件が決定されることになる。
【００５９】
次に、上記アルゴリズムを実施するための解析装置、シミュレーション方法及びシミュレーションプログラムについて説明する。
【００６０】
図４は、本発明に係るシミュレーション方法の一例を実施するための解析装置のブロック図である。
【００６１】
図４に示した解析装置３０は、例えばパーソナルコンピュータ装置やワークステーション等のコンピュータ装置であり、主な構成要素として、入力手段３１、メモリ３２、中央処理装置３３及び出力手段３４を備えている。
【００６２】
入力手段３１は、キーボードなどであり、シミュレーションを実施するためのパラメータを含む入力データを受け付ける。上記パラメータとしては、入射ビーム２０の波長λとそれに対応する角周波数ωλを含む入射ビーム２０に関する諸条件、角周波数ωλに対するεｒ及びεｉの実測値、金属細線１１の条件及び境界条件、及び、シミュレーションでの最大時刻ｔｍａｘ（ｔｍａｘ＝ｎΔｔ）等である。金属細線１１の諸条件としては、金属細線１１の幅ｗ及び厚さｈ並びに周期Λである。また、入射ビーム２０に関する条件としては上記の他に、入射ビーム２０の偏光状態（例えばＴＥモード又はＴＭモード等）、入射ビーム２０の金属回折格子１０への入射角θｉｎ（金属回折格子法線Ｎと入射ビーム２０との為す角）等である。
【００６３】
メモリ３２は、解析装置３０を動作させるための種々のプログラムや、入力手段３１からの入力データや、後述するシミュレーション方法に沿ってＦＤＴＤ解析を実施するためのＦＤＴＤプログラム（シミュレーションプログラム）４０及びＦＤＴＤプログラム４０を実施中又は実施後の計算結果を記憶する。ＦＤＴＤプログラム４０は、式（４４）〜式（５１）に基づいて空間離散化条件、すなわち、λＦを決定せしめる離散化条件決定モジュール４１と、ＴＭモード及びＴＥモードに分けて説明したアルゴリズムに基づいてＦＤＴＤ法を実行するための解析モジュール４２とを含んで構成されている。
【００６４】
中央処理装置３３は、メモリ３２に記憶されているＦＤＴＤプログラム４０を読み込み、シミュレーションのための種々の演算を実行する。また、中央処理装置３３は、メモリ３２に記憶されている各種プログラムを読み込み、解析装置３０の構成要素を制御する機能も有する。
【００６５】
出力手段３４は、例えばディスプレイなどであり、中央処理装置３３で算出されて得られるシミュレーション結果を表示する。
【００６６】
次に、ＦＤＴＤ法に従ったシミュレーション方法について説明する。図５は、シミュレーション方法を示すフローチャートである。シミュレーションは、解析装置３０のメモリ３２に格納されているＦＤＴＤプログラム４０を中央処理装置３３が読み込み実行することで実施される。
【００６７】
図５に示すように、ステップＳ１０で操作者が入力手段３１を介してシミュレーション用のパラメータを入力する（入力工程）。入力手段３１がパラメータを受け付けると、メモリ３２がそのパラメータを記憶する。
【００６８】
次に、ステップＳ２０において、離散化条件を決定する（離散化条件決定工程）。この離散化条件決定工程では、中央処理装置３３が、ステップＳ１０で入力されたεｒ及びεｉの値並びに式（４４）〜式（５１）を利用して、式（４７），式（４８）及び式（５１）を満たすλＦ及びｖＦを決定する。この際、式（４７），式（４８）及び式（５１）を満たすλＦ及びνＦが複数ある場合には、例えば、中央処理装置３３がそれらを比較してλＦが最小となる場合を選択するようにすることが好ましい。これにより、少ない格子点１１０の数で安定したシミュレーションが可能となるからである。
【００６９】
続くステップＳ３０では、ステップＳ２０で決定された空間的な離散化条件で規定されるＦＤＴＤ領域１００に、ステップＳ１０で入力された諸条件を満たす金属回折格子１０を設定し、ＦＤＴＤ領域１００における電磁界を時間ステップ間隔Δｔ毎に算出する（解析工程）。この解析工程（ステップＳ３０）においては、ステップＳ１０で入力された境界条件を必要に応じて使用する。本実施形態では、Δｔ＝１である。
【００７０】
ＴＥモードとＴＭモードとで使用するアルゴリズムが異なるため、ステップＳ３０では、(i)入射ビーム２０がＴＥモードとして与えられている場合、(ii)入射ビーム２０がＴＭモードとして与えられている場合、に分けて解析を実施する。
【００７１】
図６は、図５の解析処理のフローチャートである。解析処理では、先ず、ステップＳ３１において、入射ビーム２０の偏光状態を判定する。以下、ステップＳ１０において入射ビーム２０がＴＥモードで与えられている場合及び入射ビーム２０がＴＭモードで与えられている場合について説明する。
【００７２】
（ｉ）入射ビーム２０がＴＥモードで与えられている場合、ステップＳ３２Ａにおいて、セルの電磁界を初期化する。具体的には、各セルにおける電界Ｅ及び磁界Ｈの各軸方向の成分を０とする。
【００７３】
続いてステップＳ３３Ａに進み、式（３７）及び式（３０）を利用して電界Ｅを算出し、電界Ｅをアップデートすると共に、時刻ｔをアップデートする。なお、式（３０）におけるｃ１，ｃ２の値としては、ステップＳ２０で決定されたλＦ及びｖＦにより規定されるｃ１，Ｆ及びｃ２，Ｆの値を採用すればよい。
【００７４】
ステップＳ３３Ａで全てのセルにおける電界Ｅを算出すると、ステップＳ３４Ａに進み、アップデートされた時刻ｔが最大時刻ｔｍａｘであるか否かを判定する。時刻ｔが最大時刻ｔｍａｘより小さい場合（ステップＳ３４ＡでＹ）には、ステップＳ３３Ａに戻る。
【００７５】
次に、（ｉｉ）入射ビーム２０がＴＭモードで与えられている場合について説明する。この場合、ステップＳ３２Ｂにおいて、セルを初期化する。具体的には、各セルにおける電界Ｅ及び磁界Ｈの各軸方向の成分を０とする。次いで、ステップＳ３３Ｂに進み、式（３３）及び式（３４）並びに式（２９）を利用して各セルの電界Ｅ及び磁界Ｈを算出し、電界Ｅ及び磁界Ｈをアップデートすると共に、時刻ｔをアップデートする。式（２９）におけるｃ１，ｃ２の値に対しては、ＴＥモードの場合と同様である。
【００７６】
ステップＳ３３Ｂで全てのセルにおける電界Ｅ及び磁界Ｈを算出すると、ステップＳ３４Ｂに進み、アップデートされた時刻ｔが最大時刻ｔｍａｘであるか否かを判定する。時刻ｔが最大時刻ｔｍａｘより小さい場合（ステップＳ３４Ｂで「Ｙ」の場合）には、ステップＳ３３Ｂに戻る。なお、図６中では、ＴＥモード及びＴＭモードに分けて説明するために、時刻ｔと最大時刻ｔｍａｘとの大小関係の判定ステップを２つに分けて説明しているが、ステップＳ３４Ａ及びステップＳ３４Ｂは実質的に同じステップである。
【００７７】
ステップＳ３４Ａ又はステップＳ３４Ｂにおいて、時刻ｔが最大時刻tｍａｘ以上である場合（ステップＳ３４Ａ，Ｓ３４Ｂで「Ｎ」の場合）には、図５に記載のステップＳ４０に進み、解析結果を出力手段３４により出力する（出力工程）。
【００７８】
また、本実施形態におけるＦＤＴＤ法の実施では、上述した各式を使用する点以外は、従来のＦＤＴＤ法のシミュレーション方法と同様であることから詳細な説明は省略している。また、金属細線１１以外の部分のセルに関しては、通常の空気中又は真空中を伝搬する電磁波に対するマクスウェル方程式又は波動方程式を扱う場合のＦＤＴＤ法と同様に、ＴＥモードに関しては下記式（５２）を利用し、ＴＭモードに関しては式（５３）及び下記式（５４）を利用すればよい。
【数５２】
0005113572.tif 000053
【数５３】
0005113572.tif 000054
【数５４】
0005113572.tif 000055
式（５２）及び式（５３）中において、ａ＝(σｄΔｔ)／（２εｄ）であり、σｄ及びεｄは隣接する金属細線１１間の媒質の電導度及び誘電率である。また、μは隣接する金属細線１１間の媒質の透磁率である。
【００７９】
以上説明した本実施形態のシミュレーション法及びＦＤＴＤプログラムでは、解析対象としての金属回折格子１０が有する金属細線１１の誘電率に対してドルーデの分散式を適用してアルゴリズムを定式化している。また、金属細線１１の誘電率に対して実測値を採用することでＦＤＴＤ領域１００の空間的な離散化を図っている。そして、空間離散化を図ることでアルゴリズムの安定性も保証している。その結果、金属回折格子１０に入射ビーム２０が入射した場合の電磁界を、より正確に安定してシミュレーションすることが可能である。
【００８０】
また、ＴＥモードに対するシミュレーションでは、マクスウェル方程式から導かれる波動方程式を利用していることから、シミュレーションのアルゴリズムを簡略化することができている。その結果、シミュレーションに要する時間の短縮化も図られている。更に、本実施形態のシミュレーション方法及びＦＤＴＤプログラムでは、ＲＣ法を採用しているため、メモリ３２に蓄えるべきデータ量を低減することが可能となっている。
【００８１】
次に、本実施形態で説明したシミュレーション方法及びＦＤＴＤプログラムを用いた金属回折格子１０のシミュレーション結果について説明する。ここでは、実施例１及び実施例２として２つの金属回折格子１０についてそれぞれシミュレーションを実施した。
【００８２】
実施例１のシミュレーションの条件は以下のとおりである。
・金属細線１１の材料：銀
・金属回折格子１０の周期Λ：１７６．１２５ｎｍ
・金属回折格子１０のフィルファクタ（ｗ／Λ）：０．５
・金属細線１１の厚さｈ：１００ｎｍ
・入射ビーム２０の波長λ：７０４．５ｎｍ（角周波数ωλ：２．６７６×１０１５／ｓｅｃ）
・入射ビーム２０の位相速度ｖ：３×１０８ｍ／ｓｅｃ
・角周波数ωλに対するεｒの値：−２３．４０５
・角周波数ωλに対するεｉの値：０．３８７
【００８３】
実施例２のシミュレーションの条件は以下のとおりである。
・金属細線１１の材料：銀
・金属回折格子１０の周期Λ：４９６ｎｍ
・金属回折格子１０のフィルファクタ（ｗ／Λ）：０．５
・金属細線１１の厚さｈ：４９６ｎｍ
・入射ビーム２０の波長λ：４９６ｎｍ（角周波数ωλ：３．８×１０１５／ｓｅｃ）
入射ビーム２０の位相速度ｖ：３×１０８ｍ／ｓｅｃ
・角周波数ωλに対するεｒの値：−９．５６４
・角周波数ωλに対するεｉの値：０．３０９
【００８４】
実施例１，２では、入射ビーム２０をＴＥモード及びＴＭモードとしてそれぞれ与えた場合に、入射角θｉｎを変えながら図５，６に示したフローチャートに基づいてシミュレーションを実施した。なお、実施例１，２では、金属回折格子１０の大きさはｘ軸方向に無限であることを仮定している。
【００８５】
実施例１の条件の場合、ステップＳ２０において、λが７０４．５ｎｍであり、εｒが−２３．４０５であり、εｉが０．３８７であるとき、λＦは３５．２２７であり、νＦは０．５３８であった。この場合、ｃ１，Ｆは０．２２５であり、ｃ２，Ｆは１．００２である。また、実施例２の条件の場合、ステップＳ２０において、λが４９６ｎｍであり、εｒが−９．５６４であり、εｉが０．３０９であるとき、λＦは３５．２２７であり、νＦは０．５３８であった。この場合、ｃ１，Ｆは０．０９８であり、ｃ２，Ｆは１．００３である。
【００８６】
図７は、実施例１の場合のシミュレーション結果に基づいた透過スペクトル及び反射スペクトルを示している。図７の横軸は入射角θｉｎ（°）を示し、縦軸は透過率（％）及び反射率（％）を示している。透過率及び反射率は、入射ビーム２０の光強度に対する透過ビーム及び反射ビームの光強度の比として定義している。図８は実施例２のシミュレーション結果に基づいた透過スペクトル及び反射スペクトルを示している。図８の横軸及び縦軸は図７の場合と同様である。なお、図７及び図８のシミュレーション結果は、銀から構成される無限の大きさの金属回折格子１０におけるファーフィールドでの透過率及び反射率に基づいている。ファーフィルでの結果は、ＦＤＴＤ法を実行して得られる複素電磁界を金属回折格子１０の周囲を囲む閉曲線に沿って統合することによって計算される。
【００８７】
図７及び図８に示すように、上記シミュレーションにより、ＴＥモード及びＴＭモードに対して透過特性及び反射特性を得ることができる。そして、上記シミュレーションにより、金属回折格子１０の構成の違いによる特性の違いが現れていることがわかる。例えば、図７に示すように、上記シミュレーションにより、実施例１の金属回折格子１０では、ＴＥモードの多くは反射しており、透過が抑制されている。また、実施例１の金属回折格子１０では、ＴＭモードの多くが透過させることができている一方、その一部を反射させることが可能となっている。このように金属回折格子１０の構成では偏光分離ができていることがわかる。また、図８に示すように、実施例２の金属回折格子１０では、ＴＭモード及びＴＥモードが共に透過可能であるが、ＴＥモードの方がより多く透過させることができており、ＴＭモードの透過率を入射角θｉｎが０〜６０°で１０％〜４０％に抑制できていることがわかる。
【００８８】
よって、上記シミュレーション方法及びＦＤＴＤプログラムを利用すれば、金属回折格子１０に関するパラメータを変更することで、所望の透過特性及び反射特性を備えた金属回折格子１０を設計可能であることになる。
【００８９】
図９は、図１に示した金属回折格子が適用された導光板ユニットの構成を模式的に示す側面図である。導光板ユニット５０は、略直方体形状の導光板５１と導光板５１の表面５１ａ上に金属回折格子１０が設けられて構成されている。導光板ユニット５０は、導光板５１の側面５１ｂの側方に光源（不図示）が配置されることで面光源装置として利用することができる。
【００９０】
光源から出力され側面５１ｂから導光板５１内に入射された光（ビーム）は、導光板５１内を側面５１ｃ側に向けて伝搬する。この導光板５１内を伝搬する光が金属回折格子１０に入射すると、金属回折格子１０の回折により、入射された光の一部は導光板ユニット５０から出射され、他の部分は主に導光板５１内に戻される。
【００９１】
このように金属回折格子１０を導光板ユニット５０の一部として利用する場合には、金属回折格子１０が偏光分離素子として設計されていることが好ましい。具体的には、所定の偏光成分（例えば、Ｐ偏光成分としてのＴＭモード）が選択的に透過され、他の光が反射されて導光板５１内に戻されるように、金属回折格子１０が設計されていることが好ましい。この場合、導光板ユニット５０からは所定の偏光成分が支配的なビーム（透過ビーム２１）が出射されることになるので、導光板ユニット５０からの出射ビームを、例えば液晶表示装置におけるバックライトとして好適に使用することができる。そして、このようなバックライトを使用する液晶表示装置では、不要偏光を除去する光学素子が不要となるため、小型化や薄型化を図ることが可能である。
【００９２】
そして、上述したような特性を有する金属回折格子１０を設計するために、前述したシミュレーション方法及びＦＤＴＤプログラムが有効である。金属回折格子１０を、例えば液晶表示装置が有する面光源装置で使用する場合には、所定波長λは、可視領域の波長として、シミュレーションを実施すればよい。
【００９３】
金属回折格子１０を導光板５１上に設ける際には、図１及び図２に示したように、金属回折格子１０のみを解析対象としてシミュレーションを実施してもよいが、例えば、図９に示したように、板状の導光板５１を含む導光板ユニット５０を解析対象としても良い。
【００９４】
この場合、導光板５１に対応するセルでは、ＦＤＴＤ法のアルゴリズムとして以下の式を利用すればよい。すなわち、ＴＥモードの場合には式（５５）を利用し、ＴＭモードの場合には、式（５６）及び式（５７）を利用すればよい。
【数５５】
0005113572.tif 000056
【数５６】
0005113572.tif 000057
【数５７】
0005113572.tif 000058
式（５５）及び式（５６）中においては、ａ＝(σｄΔｔ)／（２εｄ）であり、σｄ及びεｄは導光板５１の電導度及び誘電率である。また、μは導光板５１の透磁率である。
【００９５】
また、上記実施形態では、入射ビーム２０は、ＴＥモード又はＴＭモードとして与えられるとしたが、例えば、ＴＥモード及びＴＭモードを含む非偏光状態の入射ビーム２０を与えても良い。この場合には、解析工程（ステップＳ３０）において、入射ビーム２０中のＴＥモードに対しては、ステップＳ３３Ａで説明した場合と同様にして電界Ｅを算出し、入射ビーム２０中のＴＭモードに対しては、ステップＳ３３Ｂで説明した場合と同様にして電界Ｅ及び磁界Ｈを算出し、それらを組み合わせて各セルの電界Ｅ及び磁界Ｈとすればよい。
【００９６】
ＴＭモード及びＴＥモードに対して算出した電界Ｅ及び磁界Ｈの組み合わせ方としては、各セルに対して算出した電界Ｅ及び磁界Ｈそれぞれの加重平均をとることができる。例えば、一つのセルに対応して、ＴＭモードに対して算出した電界をＥＴＭとし、ＴＥモードに対して算出した電界をＥＴＥとしたとき、そのセルにおける電界を（ＥＴＭ＋ＥＴＥ）／２とすることができる。磁界についても同様である。なお、このように加重平均をとる場合には、例えば、入射ビーム２０をＴＥモード及びＴＭモードからなる光であると仮定して算出した電界Ｅ及ぶ磁界Ｈについてそれぞれ前述のように加重平均を算出することで、非偏光状態の入射ビーム２０が入射した場合の電界Ｅ及び磁界Ｈとすることもできる。
【００９７】
また、金属回折格子１０は、図９を利用して説明したように導光板ユニット５０に利用することができるが、例えば、光学系内の所定の位置に配置される偏光分離素子として利用することもできる。この場合、例えば、ＴＥモード及びＴＭモードを完全分離するように金属回折格子１０を設計することが考えられる。
【００９８】
また、本発明に係るシミュレーションプログラムは、例えば、ＴＥモード及びＴＭモードの何れにも対応可能なプログラムとして説明したが、ＴＥモードに対するプログラム及びＴＭモードに対するプログラムすることもできる。そして、シミュレーション方法における解析工程（ステップＳ３０）では、図６に示したステップＳ３２における判定工程を省略して、図６に示したＴＥモード及びＴＭモードの各工程に沿ってシミュレーションを実施することも可能である。更に、上記実施形態では、２次元のＦＤＴＤ法について説明したが、３次元のＦＤＴＤ法に拡張することもできる。更に、上記実施形態では金属部材を含む構造物として金属回折格子１０を例示して説明したが、上記構造物は金属回折格子に限定されず、本発明に係るシミュレーション方法及びシミュレーションプログラムは、金属部材を含む任意の構造物に対して適用することができる。
【図面の簡単な説明】
【００９９】
【図１】本発明に係るシミュレーション方法の一実施形態適用して設計するための光学素子の模式図である。
【図２】ＦＤＴＤ領域に設定された図１に示した光学素子の模式図である。
【図３】角周波数に対する誘電率の実測値及びフィッティングカーブの一例を示す図面である。
【図４】本発明に係るシミュレーション方法の一実施形態を適用するための解析装置のブロック図である。
【図５】本発明に係るシミュレーション方法の一実施形態のフローチャートである。
【図６】図５の解析処理のフローチャートである。
【図７】実施例１のシミュレーション結果に基づいた透過スペクトル及び反射スペクトルを示す図面である。
【図８】実施例２のシミュレーション結果に基づいた透過スペクトル及び反射スペクトルを示図面である。
【図９】図１に示した光学素子を含む導光板ユニットの構成を模式的に示す側面図である。
【符号の説明】
【０１００】
１０…金属回折格子（構造物）、１１…金属細線（金属部材）、２０…入射ビーム（電磁波）、１００…ＦＤＴＤ領域、１１０…格子点。
【図１】
0005113572.tif 000073
【図２】
0005113572.tif 000074
【図３】
0005113572.tif 000075
【図４】
0005113572.tif 000076
【図５】
0005113572.tif 000077
【図６】
0005113572.tif 000078
【図７】
0005113572.tif 000079
【図８】
0005113572.tif 000080
【図９】
0005113572.tif 000081