(57)【要約】【課題】本発明は、ＦＤＴＤ法によるシミュレーションをより安定して実施可能なシミュレーション方法及びシミュレーションプログラムを提供する。【解決手段】  一実施形態のシミュレーション方法は、ＦＤＴＤ領域１８のセル２２の大きさ及び時間ステップ単位を設定し、電磁界をＦＤＴＤ法に従って算出する。解析モデル２０のうち物体モデル２２Ａ及び媒質モデル２２Ｂに対して、第１及び第２の電磁界解析用の式を用いて電磁界を算出し、第１の電磁界解析用の式は、マクスウェル方程式に帰納的畳込み法を組み込むことで導かれる。少なくとも物体モデルに対するセルの第１〜第３の辺の長さをΔｓ１、Δｓ２及びΔｓ３とし、時間ステップ単位をΔｔとしたとき、Δｓ１、Δｓ２、Δｓ３及びΔｔは、安定化条件：Δｔ／Δｓｊ≦０．３５４／ｖに基づいて設定される。ただし、ｊ＝１，２，３であり、ｖは真空における光速である。【選択図】  図５
2014026524.tif 000002
【特許請求の範囲】
【請求項１】
誘電率の実部が正である媒質内に誘電率の実部が０以下であり虚部が正の値である物体が配置された解析対象物に電磁波が入射された場合の電磁界を、ＦＤＴＤ法に従って３次元的に解析するためのシミュレーション方法であって、
前記解析対象物の解析モデルを配置するＦＤＴＤ領域であって、複数のセルから構成される前記ＦＤＴＤ領域の前記セルの大きさ及び前記ＦＤＴＤ法において前記電磁界を更新するための時間ステップ単位を、コンピュータによって設定する解析条件設定工程と、
前記解析モデルに前記電磁波が入射した場合の電磁界を、前記ＦＤＴＤ法に従って前記コンピュータにより算出する解析工程と、
を備え、
前記解析工程において、前記解析モデルのうち前記物体に対応する物体モデルに対しては、マクスウェル方程式から導かれる第１の電磁界解析用の式を用いて電磁界を前記コンピュータにより算出し、前記解析モデルのうち前記媒質に対応する媒質モデルに対しては、マクスウェル方程式から導かれる第２の電磁界解析用の式を用いて電磁界を前記コンピュータによって算出し、
前記第１の電磁界解析用の式は、前記マクスウェル方程式に帰納的畳込み法を組み込むことによって導かれており、
真空における光速をｖとし、少なくとも前記物体モデル内の前記セルが有する互いに直交する第１〜第３の辺の長さをそれぞれΔｓ１、Δｓ２及びΔｓ３とし、時間ステップ単位をΔｔとしたとき、前記解析条件設定工程において、前記Δｓ１、Δｓ２及びΔｓ３と及び前記Δｔは、以下の安定化条件を示す式（１）に基づいて設定される、
シミュレーション方法。
【数１】
2014026524.tif 000053
ただし、ｊ＝１，２，３
【請求項２】
前記電磁波は単色の電磁波である、請求項１記載のシミュレーション方法。
【請求項３】
前記物体及び前記媒質における前記誘電率の実部及び虚部の値には、実測値又は実測値を補間して得られる補間値が使用される、請求項１又は２記載のシミュレーション方法。
【請求項４】
前記物体の材料は、金、銀、アルミニウム及び銅のうちの少なくとも一つである、請求項１〜３の何れか一項記載のシミュレーション方法。
【請求項５】
コンピュータに、請求項１〜４のいずれか一項に記載のシミュレーション方法を実行させる、シミュレーションプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、ＦＤＴＤ法を利用したシミュレーション方法及びシミュレーションプログラムに関する。
【背景技術】
【０００２】
所定物体が媒質に囲まれた系の中を進行する電磁波の挙動を、コンピュータを用いて算出する電磁波のシミュレーション方法は、以下のような種々の設計に有用である。
・拡散板やカラーフィルターなどのディスプレイ用部材の設計
・高周波回路、レーダー及びアンテナなどのエレクトロニクス部材の設計
・インキ、塗料、プラスチック、及び、染色などの着色工業における材料の設計
・リモートセンシングや気象科学や医療の分野での各種測定装置の設計
・ＬＥＤデバイス、特に、ナノ微粒子又はナノ構造を含むＬＥＤデバイスの設計に必要とされる、ＬＥＤの光取出し効率の計算
・太陽電池における光―電気変換効率の改善のための設計
【０００３】
電磁波のシミュレーションにおいては、２次元平面内のシミュレーションよりも、実際の電磁波の進行に近い３次元空間内のシミュレーションが求められている。そこで、物体が媒質に囲まれた系の中で、電磁波が進行する仮想空間内に格子状のセルを設定し、ＦＤＴＤ（Finite  Difference  Time  Domain）法を用いて計算を行う電磁波のシミュレーション方法が提案されている（特許文献１参照。）。ＦＤＴＤ法は、上記仮想空間であるＦＤＴＤ領域を複数のセルに分割して、マクスウェル方程式内の空間および時間微分を有限差分によって近似することでマクスウェル方程式を直接解くものである。
【０００４】
しかしながら、系の中の物体が金属を含む場合、物体の誘電率の実部が、負になるので、ＦＤＴＤ法の計算において場の更新が進むにつれて、すなわち、時間ステップの更新が進むにつれて、計算が発散する傾向にある。
【０００５】
このような発散を抑制するために、例えば、デバイ（Debye）、ドルーデ(Drude)、またはローレンツ（Lorentz）など、誘電率についての様々な解析的モデルが、周波数領域畳込みを使用することによって、ＦＤＴＤに組み込むことができる。このような畳込みを利用した方法は、帰納的畳込み（Recursive Convolution）法（以下、単にＲＣ法ともいう）として知られている。
【先行技術文献】
【特許文献】
【０００６】
【特許文献１】特開２０１０−２０４８５９号公報
【発明の概要】
【発明が解決しようとする課題】
【０００７】
このような帰納的畳込み法を利用しても、空間的離散化の値等の条件設定によっては、シミュレーションが不安定になる傾向にあった。
【０００８】
そこで、本発明は、ＦＤＴＤ法によるシミュレーションをより安定して実施可能なシミュレーション方法及びシミュレーションプログラムを提供すること目的とする。
【課題を解決するための手段】
【０００９】
本発明の一側面に係るシミュレーション方法は、誘電率の実部が正である媒質内に誘電率の実部が０以下であり虚部が正の値である物体が配置された解析対象物に電磁波が照射された場合の電磁界を、ＦＤＴＤ法に従って３次元的に解析するためのシミュレーション方法である。このシミュレーション方法では、解析対象物の解析モデルを配置するＦＤＴＤ領域であって、複数のセルから構成されるＦＤＴＤ領域のセルの大きさ及びＦＤＴＤ法において電磁界を更新するための時間ステップ単位を、コンピュータによって設定する解析条件設定工程と、解析モデルに電磁波が入射した場合の電磁界を、ＦＤＴＤ法に従ってコンピュータにより算出する解析工程と、を備える。上記解析工程において、解析モデルのうち物体に対応する物体モデルに対しては、マクスウェル方程式から導かれる第１の電磁界解析用の式を用いて電磁界をコンピュータにより算出し、解析モデルのうち媒質に対応する媒質モデルに対しては、マクスウェル方程式から導かれる第２の電磁界解析用の式を用いて電磁界をコンピュータによる算出する。第１の電磁界解析用の式は、マクスウェル方程式に帰納的畳込み法を組み込むことによって導かれており、真空における光速をｖとし、少なくとも物体モデル内のセルが有する互いに直交する第１〜第３の辺の長さをそれぞれΔｓ１、Δｓ２及びΔｓ３とし、時間ステップ単位をΔｔとしたとき、ＦＤＴＤ領域設定工程において、Δｓ１、Δｓ２、Δｓ３及びΔｔは、以下の安定化条件を示す式（１）に基づいて設定される。
【数１】
2014026524.tif 000003
ただし、ｊ＝１，２，３
【００１０】
この方法では、式（１）で表される安定化条件に基づいて、少なくとも物体モデルのセルのサイズであるΔｓ１、Δｓ２、Δｓ３、及び時間ステップ単位Δｔが設定される。その結果、シミュレーションを安定して実施可能である。
【００１１】
一実施形態において、上記電磁波は単色の電磁波であり得る。
【００１２】
一実施形態において、物体及び媒質における誘電率の実部及び虚部の値には、実測値又は実測値を補間して得られる補間値が使用されてもよい。
【００１３】
このように実測値又は補間値を利用することで、より正確にシミュレーションを実施することが可能である。
【００１４】
一実施形態において、上記物体の材料は、金、銀、アルミニウム及び銅のうちの少なくとも一つであり得る。
【００１５】
本発明の他の側面は、コンピュータに、本発明の一側面に係る上記シミュレーション方法を実行させるためのシミュレーションプログラムに係る。
【発明の効果】
【００１６】
本発明によれば、ＦＤＴＤ法に従ってシミュレーションを安定して実施することができる。
【図面の簡単な説明】
【００１７】
【図１】本発明に係るシミュレーション方法の一実施形態を適用する解析対象物の模式図である。
【図２】図１に示した解析対象物に対応する解析モデルを配置するＦＤＴＤ領域の模式図である。
【図３】図２に示したセルの拡大図である。
【図４】本発明に係るシミュレーション方法の一実施形態を適用するための解析装置のブロック図である。
【図５】本発明に係るシミュレーション方法の一実施形態のフローチャートである。
【図６】補間値と、最適化されたパラメータを使用したドルーデの分散式によって算出された値との比較を示す図面である。
【図７】シミュレーションの一例の結果を示す図面である。
【発明を実施するための形態】
【００１８】
以下、図面を参照して本発明に係るシミュレーション方法及びシミュレーションプログラムの実施形態について説明する。図面の説明においては同一要素には同一符号を付し、重複する説明を省略する。図面の寸法比率等は、説明のものと必ずしも一致していない。
【００１９】
本発明に係るシミュレーション方法及びシミュレーションプログラムは、誘電率の実部が正である媒質中に、誘電率の実部が０以下である物体を解析対象物に電磁波が入射した場合における電磁界を、ＦＤＴＤ法（Finite Difference Time Domain Method）に従って３次元的に解析するためのものである。上記物体の例は、金属微粒子や、複数の金属微粒子が連結された部材が例示され得る。上記媒質の例は、樹脂などの誘電体、空気及び真空を含み得る。以下の説明では、特に断らない限り、媒質中の物体は、金属微粒子であり、媒質は誘電体である。本発明は、例えば、金属微粒子などが拡散剤として利用された光学シートの設計などに適用され得る。
【００２０】
図１は、一実施形態に係る解析対象物の一例を示す模式図である。解析対象物１０は、金属微粒子１２と、その金属微粒子１２周囲において光を伝搬させるための伝搬媒質１４とから構成され得る。金属微粒子１２の一例は、金属から構成されるナノ微粒子である。金属微粒子１２は、例えば、銀、金、アルミニウム及び銅のうち少なくとも一つから構成され得る。誘電体としての伝搬媒質１４の例は、樹脂である。
【００２１】
本実施形態に係るシミュレーション方法及びシミュレーションプログラムは、図１に示したような解析対象物１０に単色の光（以下、入射光とも称す）１６が入射した場合における電磁界をＦＤＴＤ法（Finite Difference Time Domain Method）に従って解析するためのものである。入射光１６の波長λは、紫外線領域、可視光領域、近赤外光領域のうちの範囲から選択され得る。波長λは、好ましくは、３００ｎｍ〜９００ｎｍの範囲から選択され、より好ましくは、３００ｎｍ〜７００ｎｍの範囲から選択される。
【００２２】
ＦＤＴＤ法について説明する。ＦＤＴＤ法では、コンピュータで設定する仮想空間としてのＦＤＴＤ領域１８内に、解析対象物１０に対応する解析モデル２０を配置する。図２は、解析モデルを示す模式図である。解析モデル２０において、金属微粒子１２に対応する金属微粒子モデルを金属微粒子モデル２０Ａと称し、媒質１４に対応する媒質モデルを媒質モデル２０Ｂと称する。図２では、金属微粒子モデル２０Ａを明確に示すために破線で示している。図２では、金属微粒子モデル２０Ａの外形を曲線で表しているが、後述するように、ＦＤＴＤ領域１８は空間的に離散化される。従って、金属微粒子モデル２０Ａの外形形状は、微細な直線（或いはセルの辺）の組み合わせで表される。このような物体の表現方法は、階段近似(staircace approximation)として知られている。空間の離散化が微細になれば、物体の境界面はより実際の曲線に近づく。そのため、後述するセル（又はグリッド）を用いて物体を表現する際の正確性を確保するためには、微細な空間離散化が好ましい。図２において、解析モデル２０の配置領域の大きさと、ＦＤＴＤ領域１８の大きさとは一致している。
【００２３】
ＦＤＴＤ法は、式（２ａ）及び式（２ｂ）に示すマクスウェル方程式内の空間および時間微分を有限差分によって近似することでマクスウェル方程式を直接解く方法である。
【数２】
2014026524.tif 000004
ｒは、ＦＤＴＤ領域１８において設定した原点からの位置ベクトルである。μは透磁率である。本実施形態の解析対象物１０において、透磁率μは、通常、定数である。Ｅは電界ベクトル、Ｈは磁界ベクトルである。「∇」は、ｘｕ、ｙｕ及びｚｕをｘ軸方向、ｙ軸方向及びｚ軸方向の単位ベクトルとしたとき、ｘｕ(∂／∂ｘ)＋ｙｕ(∂／∂ｙ)＋ｚｕ（∂／∂ｚ)で定義されるベクトル微分演算子である。「×」は、回転演算子であり、いわゆる外積を表す。
【００２４】
式（２ｂ）中のＤは電束密度ベクトルであり、周波数領域において、Ｄは次式で与えられる。
【数３】
2014026524.tif 000005
式（３）中のε（ｒ，ω）は誘電率である。
【００２５】
マクスウェル方程式を空間的及び時間的な有限差分で近似して解くために、ＦＤＴＤ法では、ＦＤＴＤ領域１８を、図２に示したように、複数のセル２２に分割する。セル２２の個数は、解析モデル２０の大きさ、計算時間及び計算に要するメモリの容量などに応じて設定すればよい。ＦＤＴＤ領域１８を構成する各セル２２の形状の例は立方体である。従って、セル２２において、互いに直交するｘ軸方向の辺（第１の辺）、ｙ軸方向の辺（第２の辺）及びｚ軸方向の辺（第３の辺）の長さをそれぞれΔｘ（＝Δｓ１）、Δｙ（＝Δｓ２）及びΔｚ（＝Δｓ３）としたとき、Δｘ＝Δｙ＝Δｚである。以下の説明では、立方体のセル２２の一辺の長さをΔｓとする。この場合、Δｘ＝Δｙ＝Δｚ＝Δｓである。
【００２６】
図３は、図２に示したセル２２の拡大図である。図３において、白抜き矢印は電界を示し、黒色矢印は磁界を示す。図３に示すように、セル２２を構成する各辺の中心で電界強度が計算され、セル２２を構成する面の中心で磁界強度が計算される。
【００２７】
次に、電界及び磁界を計算するための式について、金属微粒子モデル２０Ａに対する計算式と、媒質モデル２０Ｂに対する計算式とに分けて説明する。以下の説明において、誘電率をεとする。誘電率εの複素数表示は、実部をεｒ及び虚部をεｉとした場合、ε＝εｒ−ｉεｉである。まず、金属微粒子モデル２０Ａに対する計算式について説明する。
【００２８】
電磁界計算のための時間ステップ単位をΔｔと表し、ｎ（ｎは０以上の整数）を時間ステップ数とする。この場合、ある時刻ｔはｔ＝ｎΔｔである。時刻ｔにおける電界ベクトル及び磁界ベクトルをＥｎ及びＨｎと表す。
【００２９】
電磁波の伝搬はＦＤＴＤ法を利用して実時間で計算されるので、式（３）で表される電束密度ベクトルを時間領域において取り扱う必要がある。時間領域における電束密度ベクトルＤは、式（３）の両辺でフーリエ変換をとり、畳込み定理を適用することによって次式で表される。
【数４】
2014026524.tif 000006
式（４）において、ε（ｔ）＝Ｆ｛ε（ω）｝であり、Ｆはフーリエ変換を表す。
【００３０】
誘電率ε（ω）が角周波数ωのある特殊な関数で表されるとき、式（４）の畳込み積分は帰納的に計算可能である。これは帰納的畳込み法（ＲＣ法）として知られている。
【００３１】
本実施形態では、金属微粒子１２の誘電率ε（ω）を表す関数として、ドルーデの一次近似式（ドルーデの分散式）を採用する。この場合、誘電率ε（ω）は式（５）のように表される。
【数５】
2014026524.tif 000007
式（５）において、ωｐはプラズマ周波数、νｃは衝突周波数である。ε∞は、周波数が無限大のときの誘電率であり、１に等しい。χ（ω）は電気比感受率である。
【００３２】
χ（ω）をフーリエ変換したもの、すなわち、Ｆ｛χ（ω）｝をχ（ｔ）と表した場合、χ（ｔ）は次式で表される。
【数６】
2014026524.tif 000008
Ｕ（ｔ）は、時間の単位ステップ関数であって、ｔ＝０の場合、Ｕ（ｔ）＝０であり、ｔ＞０の場合、Ｕ（ｔ）＝１である。
【００３３】
離散化された時間領域中の時刻ｔ＝ｎΔｔ（ｎは０以上の整数）において、区間［ｍΔｔ，（ｍ＋１）Δｔ］（ｍは０以上ｎ以下の整数）にわたって電界ベクトルＥを定数ベクトルとして、式（４）の畳込み積分を評価した場合、次式が得られる。
【数７】
2014026524.tif 000009
式（７）中において、ε０は、真空の誘電率である。
【００３４】
式（７）は、以下のように書き表すことができる。
【数８】
2014026524.tif 000010
Ψは、累積変数ベクトルであり下記式（９）が成り立つ。
【数９】
2014026524.tif 000011
【００３５】
式（９）は、以下の式（１０）のように表される。
【数１０】
2014026524.tif 000012
この式（１０）において、
【数１１】
2014026524.tif 000013
である。式（１１）より、
【数１２】
2014026524.tif 000014
が成立する。式（１２）において、
【数１３】
2014026524.tif 000015
である。式（１２）において、Ψ０＝０とすると共に、Ψ１＝Ｅ１χ０とする。
【００３６】
式（１０）〜式（１３ｂ）より、累積変数Ψに関する以下の帰納的関係式が得られる。
【数１４】
2014026524.tif 000016
【００３７】
式（１４）は、χ０＝０とすると共に、ｎ≦０においてΨを０とすることによって解かれ得る。
【００３８】
ＦＤＴＤ法において、式（２ａ）及び式（２ｂ）は、中心差分法を利用して解かれる。式（２ａ）中の一次時間導関数を近似するために、二次精度有限差分表現を利用したとき、式（２ａ）は、次式のように表される。
【数１５】
2014026524.tif 000017
【００３９】
磁界ベクトルＨの更新方程式は、式（２ａ）に等価な差分形式の方程式から導かれて、以下のように表される。
【数１６】
2014026524.tif 000018
式（１６）中において、ｄはベクトル差分演算子である。ｄは下記式で定義される。
【数１７】
2014026524.tif 000019
ｄｘ、ｄｙ及びｄｚは、一次微分である∂／∂ｘ、∂／∂ｙ及び∂／∂ｚの二次精度有限差分近似である。任意の関数をｆ（ｘ，ｙ，ｚ）で表したとき、ｄｘｆ（ｘ，ｙ，ｚ）、ｄｙｆ（ｘ，ｙ，ｚ）及びｄｚｆ（ｘ，ｙ，ｚ）は、式（１８ａ）〜式（１８ｃ）で与えられる。
【数１８】
2014026524.tif 000020
【００４０】
ＦＤＴＤ法で電磁界を３次元的に計算する場合、磁界ベクトルＨは、ｘ軸方向、ｙ軸方向及びｚ軸方向にそれぞれ成分を有する。磁界ベクトルＨの各成分の更新方程式は、式（１６）、式（１７）及び式（１８ａ）〜式（１８ｃ）を考慮すれば、次のように表される。
【数１９】
2014026524.tif 000021
【００４１】
次に、電界ベクトルＥについて説明する。式（２ｂ）中の一次時間導関数を近似するために、二次精度有限差分表現を利用したとき、式（２ｂ）は、次式のように表される。
【数２０】
2014026524.tif 000022
【００４２】
電界ベクトルＥの更新方程式は、式（２ｂ）に等価な差分形式の方程式から導かれて、以下のように表される。
【数２１】
2014026524.tif 000023
【００４３】
式（２１）中のΨｎ＋１−Ψｎは、式（１４）により与えられる。式（２１）中のｄは、式（１７）で定義される。ｄｘｆ（ｘ，ｙ，ｚ）、ｄｙｆ（ｘ，ｙ，ｚ）及びｄｚｆ（ｘ，ｙ，ｚ）は、式（２２ａ）〜式（２２ｃ）で与えられる。
【数２２】
2014026524.tif 000024
【００４４】
磁界ベクトルＨの場合と同様に、電界ベクトルＥは、ｘ軸方向、ｙ軸方向及びｚ軸方向にそれぞれ成分を有する。電界ベクトルの各成分の更新方程式は、式（２１）、式（１７）及び式（２２ａ）〜式（２２ｃ）を考慮すれば、次のように表される。
【数２３】
2014026524.tif 000025
【００４５】
次に、媒質１４に対する計算式について説明する。透磁率μが、正の値で定数である場合、磁界ベクトルＨの更新方程式は、式（１９ａ）、式（１９ｂ）及び式（１９ｃ）と同じである。
【００４６】
一方、誘電体内では、誘電率εの実部εｒ及び虚部εｉがいずれも正の値である。媒質モデルに対する電界ベクトルＥの各成分の更新方程式は、累積変数Ψを用いずに、以下のように表される。
【数２４】
2014026524.tif 000026
【００４７】
式（２４ａ）〜式（２４ｃ）内のａは、σを伝導率とした場合、（σΔｔ）／（２εｒ）＝（εｉωΔｔ）／（２εｒ）で定義され得る。ただし、εｉ≧０であり、εｒ＞εｉである。
【００４８】
［安定化条件］
次に、ＦＤＴＤ法の計算を発散させないための安定化条件について説明する。帰納的畳込み法を組み込んでいないアルゴリズムに対しては、すなわち、媒質モデル２０Ｂに対しては、いわゆるクーラン条件に基づいて安定条件を算出すればよい。ただし、帰納的畳込み法を組み込んだアルゴリズムに対する安定化条件の方がより厳しい条件であるため、帰納的畳込み法を利用したアルゴリズムに対する安定化条件について説明する。
【００４９】
帰納的畳込み法を組み込んだＦＤＴＤ法のアルゴリズムの安定化させるためには、式（１４）におけるｃ１及びｃ２は、以下の条件を満たす必要がある。
【数２５】
2014026524.tif 000027
ｃ１及びｃ２は、計算に使用する波長に対する金属微粒子の誘電率の値を用いて計算され得る。ここで、入射光１６の波長λ（又は角周波数ω）に対応する金属微粒子１２の誘電率ε（＝εｒ−ｉεｉ）に対して、プラズマ周波数ωｐ及び衝突周波数νｃは、式（５）より、以下の式（２６）及び式（２７）で表される。
【数２６】
2014026524.tif 000028
【数２７】
2014026524.tif 000029
従って、式（１３ａ）、式（１３ｂ）、式（２６）及び式（２７）より以下の式が成り立つ。
【数２８】
2014026524.tif 000030
式（２８ａ）及び式（２８ｂ）を考慮すれば、εｉが０より大きく、εｒが０以下である限り、ｃ２は、正であり１より小さい。すなわち、εｉが０より大きく、εｒが０以下である金属微粒子に対しては、（式２５ｂ）は成立する。更に、式（２５ａ）を満たすように、ωΔｔが選択され得る。
【００５０】
次に、金属微粒子モデル２０Ａに対して適用する式（２３ａ）〜式（２３ｃ）を利用する場合の安定化条件について説明する
【００５１】
有限差分時間微分演算子をｄｔとしたとき、式（２１）及び式（１４）は次のように表され得る。
【数２９】
2014026524.tif 000031
ただし、ｄｔは、ｄｔΨｎ＋１／２＝Ψｎ＋１−Ψｎで定義される演算子である。
【００５２】
式（２９ａ）、（２９ｂ）及び式（１９ａ）〜（１９ｃ）より、電磁界の更新方程式は、行列を利用すると以下の式に表される。
【数３０】
2014026524.tif 000032
行列Ｍ及びＮは、係数行列である。式（３０）より、次式が得られる。
【数３１】
2014026524.tif 000033
式（３１）において、行列Ｉは、恒等行列である。行列Ｍ及びＮの成分は、式（２９ａ）、（２９ｂ）及び式（１９ａ）〜（１９ｃ）を考慮すると以下のように得られる。
【数３２】
2014026524.tif 000034
ここで、ｐＥは、（１／μ）（Δｔ／Δｓ）（ｄ×）で表され、電界ベクトルＥに作用する有限差分演算子である。ｐＨは、（１／ε０）（Δｔ／Δｓ）（ｄ×）で表され、磁界ベクトルＨに作用する有限差分演算子である。
【００５３】
ここで、行列Ｃを以下のように定義する。
【数３３】
2014026524.tif 000035
式（３３）を利用すると式（３１）は、次式で表される。
【数３４】
2014026524.tif 000036
Ｌｊ及びｑｊ（ｊ＝１，２，３）を行列Ｃの固有ベクトル及び固有値として、線形システム理論を利用することによって、式（３５）及び式（３６）が得られる。
【数３５】
2014026524.tif 000037
【数３６】
2014026524.tif 000038
ただし、ａｊは定数である。
【００５４】
アルゴリズムの安定化のためには、|ｑｊ|≦１であることが必要である。ｑｊは、行列Ｃの固有値である。固有値は、式（３７）に示す固有方程式の解である。
【数３７】
2014026524.tif 000039
式（３７）においてベクトルＱは、式（３８）で定義される対角行列である。
【数３８】
2014026524.tif 000040
式（３７）より式（３９）が成り立つ。
【数３９】
2014026524.tif 000041
式（３９）の解の一つは０である。式（３９）の他の２つの解をｑ±として表すと、ｑ±は、式（４０）で与えられる。
【数４０】
2014026524.tif 000042
アルゴリズムの安定化のためには、|ｑ±|≦１であることが必要である。ｃ１及びｃ２は共に１以下であるので、ｐＥｐＨ≦２であれば、|ｑ±|≦１となる。
【００５５】
ベクトル差分演算子ｐＥｐＨは、式（４１）で与えられる。
【数４１】
2014026524.tif 000043
従って、アルゴリズムの安定性のためには、式（４２）が成立すればよい。
【数４２】
2014026524.tif 000044
【００５６】
ここで、真空の光速（電磁波の速度）をｖ＝１／（με０）１／２とすると、式（４２）は式（４３）で表される。
【数４３】
2014026524.tif 000045
【００５７】
演算子”ｄ×ｄ×”を調和ベクトル場Ｅに適用することによって演算子”ｄ×ｄ×”の最大値を評価する。調和ベクトル場Ｅは、式（４４）で表される。
【数４４】
2014026524.tif 000046
ただし、Ｅ０は、Ｅ０＝ｘｕＥ０ｘ＋ｙｕＥ０ｙ＋ｚｕＥ０ｚで表されるベクトルである。ｘｕ、ｙｕ、及びｚｕは、前述したように、それぞれｘ軸方向、ｙ軸方向及びｚ軸方向の単位ベクトルである。式（４４）中において、ｋは、波数ベクトルであり、ｒは、位置ベクトルである。
【００５８】
３次元において、ｄ×Ｅ及びｄ×ｄ×Ｅは、それぞれ式（４５ａ）及び式（４５ｂ）で与えられる。
【数４５】
2014026524.tif 000047
【００５９】
一次の差分演算子の定義から次式が成り立つ。
【数４６】
2014026524.tif 000048
式（４６）中において、α及びβは独立した添え字である。α及びβは、それぞれｘ，ｙ又はｚである。
【００６０】
ｄ×ｄ×Ｅのｘ成分を（ｄ×ｄ×Ｅ）ｘと表したとき、（ｄ×ｄ×Ｅ）ｘは、次式で与えられる。
【数４７】
2014026524.tif 000049
【００６１】
着目しているベクトル場におけるベクトルＥ０の方向に対する極角及び方位角をθ及びφとすると、Ｅ０ｘ、Ｅ０ｙ及びＥ０ｚはそれぞれ式（４８ａ）、式（４８ｂ）及び式（４８ｃ）で表される。
【数４８】
2014026524.tif 000050
【００６２】
|ｓｉｎθｃｏｓφ|≦１、|ｓｉｎθｓｉｎφ|≦１及び|ｃｏｓθ|≦１であると共に、|ｓｉｎ（ｋｍΔｓ／２）|≦１であることから、式（４７）の右辺の最大値は１６、すなわち、ｄ×ｄ×Ｅのｘ成分を（ｄ×ｄ×Ｅ）ｘの最大値は１６と評価され得る。同様に、ｄ×ｄ×Ｅのｙ成分及びｚ成分の最大値も１６と評価され得る。
【００６３】
従って、式（４３）より、以下の式が成り立つ。
【数４９】
2014026524.tif 000051
式（４９）より、アルゴリズムの安定化のために、時間と空間の離散化の比の最大値は、０．３５４／ｖ（ただし、ｖは真空の光速）であることがわかる。
【００６４】
よって、金属微粒子モデル２０Ａに対して適用する式（２３ａ）〜式（２３ｃ）を利用する場合は、式（４９）を満たすようにΔｔ及びΔｓを設定すればよい。前述したように、金属微粒子モデル２０Ａに対する安定化条件は、媒質モデル２０Ｂ安定化条件の方がより厳しい。従って、式（２３ａ）〜式（２３ｃ）に対する安定化条件を満たすように、Δｔ及びΔｓを設定すれば媒質モデル２０Ｂに対しても安定して計算可能である。
【００６５】
次に、ＦＤＴＤ法において使用する誘電率の値について説明する。例示した金属微粒子モデル２０Ａに対する電界ベクトルＥ及び磁界ベクトルＨの更新方程式（第１の電磁界解析用の式）では、式（１３ａ）及び式（１３ｂ）中のプラズマ周波数ωｐ及び衝突周波数νｃに金属微粒子１２の誘電率εが含まれている。
【００６６】
より具体的には、入射光１６の波長λ（又は振動数）に対応する金属微粒子１２の誘電率ε（＝εｒ−ｉεｉ）に対して、プラズマ周波数ωｐ及び衝突周波数νｃは、式（２６）及び式（２７）で表される。
【００６７】
媒質モデルに対する電界ベクトルＥ及び磁界ベクトルＨの更新方程式（第２の電磁界解析用の式）では、式（２４ａ）〜（２４ｃ）中の“ａ”に媒質１４の誘電率εの実部εｒ及び虚部εｉが含まれる。
【００６８】
本実施形態では、このような誘電率εの実部εｒ及び虚部εｉの値として実測値又はそれら実測値を補間して得られる値（以下、補間値と称す）を使用する。本実施形態において、実測値とは、光学特性など一覧として掲載されている書籍記載のデータ値など公知の実測値、又は、シミュレーションのために実際に使用する金属微粒子１２に対して個別に実測した値であり得る。光学特性など一覧として掲載されている書籍の例は、「C.L.Foiles, ‘Optical properties of pure metals and binary alloys’,Chapter 4 of Landolt-Bornstein Numerical Data and Functional Relationships inScience and Technology New Series, Vol.15, Subvolume b, K. H.Hellwege and J.L.Olsen, Ed.Springer-Verlag, Berlin 1985,pp. 228」である。補間値を得るための補間手法の例は、三次スプライン補間である。
【００６９】
次に、上記アルゴリズムを実施するための解析装置２４、シミュレーション方法及びシミュレーションプログラムについて説明する。
【００７０】
図４は、本発明に係るシミュレーション方法の一例を実施するための解析装置２４のブロック図である。
【００７１】
図４に示した解析装置２４は、例えばパーソナルコンピュータ装置やワークステーション等のコンピュータ装置であり、主な構成要素として、入力手段２６、メモリ２８、中央処理装置３０及び出力手段３２を備えている。
【００７２】
入力手段２６は、キーボードなどであり、シミュレーションを実施するためのパラメータを含む入力データを受け付ける。上記パラメータとしては、入射光１６の波長λを含む入射光１６に関する諸条件、波長λに対する金属微粒子１２及び媒質１４それぞれのεｒ及びεｉの値、解析対象物１０の諸条件並びに境界条件、及び、シミュレーションでの最大時刻ｔｍａｘ（ｔｍａｘ＝ｎΔｔ）等である。
【００７３】
解析対象物１０の諸条件は、解析対象物１０に対応する解析モデル２０を解析装置２４内で仮想的に構築できる条件であればよい。解析対象物１０の諸条件の例は、媒質１４及び金属微粒子１２の大きさ（サイズ）及び媒質１４内での金属微粒子１２の位置である。誘電率εの実部εｒ及び虚部εｉの値は、前述したように、実測値（既知の値）又は補間値である。境界条件としては、種々の吸収境界条件又は種々の周期的境界条件が使用され得る。吸収境界条件の例は、二次精度Ｍｕｒ境界条件、ＰＭＬ（Perfectly matched layer）及びＵＰＭＬ（Uniaxial PML）境界条件を含み得る。
【００７４】
メモリ２８は、解析装置２４を動作させるための種々のプログラムや、入力手段２６からの入力データや、後述するシミュレーション方法に沿ってＦＤＴＤ解析を実施するためのＦＤＴＤプログラム（シミュレーションプログラム）３４及びＦＤＴＤプログラム３４を実施中又は実施後の計算結果を記憶する。ＦＤＴＤプログラム３４は、金属微粒子モデル及び媒質モデルにそれぞれ適用するための更新方程式を利用して安定的にＦＤＴＤ方法を実施するための解析条件を設定せしめる解析条件設定モジュール３４Ａと、ＦＤＴＤ法を実行するための解析モジュール３４Ｂとを含む。
【００７５】
中央処理装置３０は、メモリ２８に記憶されているＦＤＴＤプログラム３４を読み込み、シミュレーションのための種々の演算を実行する。中央処理装置３０は、メモリ２８に記憶されている各種プログラムを読み込み、解析装置２４の構成要素を制御する機能も有する。
【００７６】
出力手段３２は、例えばディスプレイなどであり、中央処理装置３０で算出されて得られるシミュレーション結果を表示する。
【００７７】
次に、ＦＤＴＤ法に従ったシミュレーション方法について説明する。図５は、シミュレーション方法を示すフローチャートである。シミュレーションは、解析装置２４のメモリ２８に格納されているＦＤＴＤプログラム３４を中央処理装置３０が読み込み実行することで実施される。
【００７８】
図５に示すように、ステップＳ１０で操作者が入力手段２６を介してシミュレーション用のパラメータを入力する（入力工程）。入力手段２６がパラメータを受け付けると、メモリ２８がそのパラメータを記憶する。
【００７９】
次に、ステップＳ１２において、解析条件を設定する（解析条件設定工程）。この解析条件決定工程では、中央処理装置３０が、安定化条件に基づいて解析条件を設定する。解析条件は、例えば、式（４９）に基づいて設定され得る。具体的には、時間ステップ単位Δｔが一定の場合、式（４９）からΔｓ、すなわち、解析対象物１０に対する解析モデル２０を複数のセル２２に分ける際のセル２２の一辺の大きさが設定され得る。なお、式（４９）を満たすΔｓを設定する場合は、計算時間が短縮出来るようなΔｓをより優先的に設定する。
【００８０】
他の実施形態では、Δｔ及びΔｓを、式（４９）を満たすように任意に設定してもよい。
【００８１】
更に他の実施形態では、次のステップによってΔｓ及びΔｔを設定してもよい。
【００８２】
・ステップI：式（４９）のｖを２．９９８×１０８ｍ／ｓにセットする。
・ステップII：式（２８ａ）及び式（２８ｂ）で表されるｃ１及びｃ２がそれぞれ１以下となるようにΔｔを選択する。
・ステップIII:式（４９）を満たすように、Δｓを選択する。
【００８３】
続くステップＳ１４では、ステップＳ１２で設定された空間的な離散化条件で規定されるＦＤＴＤ領域１８を設定し、ＦＤＴＤ領域１８に解析モデル２０を配置する。
【００８４】
その後、ステップＳ１６において、ＦＤＴＤ領域１８を構成する全てのセル２２の電磁界を初期化する。具体的には、各セル２２における電界及び磁界の各軸方向の成分を０とする。
【００８５】
ステップＳ１８において、時刻ｔにおける磁界及び電界を、全てのセル２２において算出し、時刻ｔをアップデートする（解析工程）。具体的には、媒質モデル２０Ｂに対しては、式（１９ａ）〜式（１９ｃ）及び式（２４ａ）〜式（２４ｃ）に基づいて磁界及び電界を算出し、金属微粒子モデル２０Ａに対しては、式（１９ａ）〜式（１９ｃ）及び式（２３ａ）〜式（２３ｃ）に基づいて磁界及び電界を算出する。式（１９ａ）〜式（１９ｃ）、式（２３ａ）〜式（２３ｃ）及び式（２４ａ）〜式（２４ｃ）に示すように、ＦＤＴＤ法では、電界を算出する場合には、電界を算出する位置において、１時間ステップ単位前に算出された電界と、その電界を算出する位置の近傍（その位置を囲む位置）で半時間ステップ単位前に算出された磁界とに基づいて電界が算出される。一方、磁界を算出する場合には、磁界を算出する位置において、１時間ステップ単位前に算出された磁界と、その磁界を算出する位置の近傍（その位置を囲む位置）で半時間ステップ単位前に算出された電界とに基づいて磁界が算出される。
【００８６】
ステップＳ２０において、アップデートされた時刻ｔが最大時刻ｔｍａｘであるか否かを判定する。時刻ｔが最大時刻ｔｍａｘより小さい場合（ステップＳ２０で「ＹＥＳ」）には、ステップＳ１８に戻る。
【００８７】
ステップＳ２０において、時刻ｔが最大時刻ｔｍａｘ以上である場合（ステップＳ２０で「ＮＯ」の場合）には、ステップＳ２２に進み、解析結果を出力手段３２により出力する（出力工程）。
【００８８】
以上説明した本実施形態のシミュレーション法及びＦＤＴＤプログラム３４では、安定化条件に基づいて、Δｓ及びΔｔを設定している。そのため、誘電率εの実部εｒの値が負である金属微粒子１２を解析対象に含んでいても安定してＦＤＴＤ法を実行することが可能である。
【００８９】
式（２６）及び式（２７）に含まれる金属微粒子１２の誘電率εの実部εｒ及び虚部εｉの値として、解析に使用する波長λに対する実測値（既知の値）を利用することによって、より正確にシミュレーションすることが可能である。
【００９０】
ただし、誘電率εの実測値は、通常、光のエネルギー又は周波数において離散的に取得されている。波長は、周波数に反比例するので、仮に、エネルギー領域において一定間隔で誘電率εの実測値が存在しても、波長領域において一定間隔で実測値は存在しない場合があり得る。従って、解析に使用する波長λに対する誘電率εの実測値が存在しない場合があり得る。この場合、前述したように、実測値から補間法を使用して得られた補間値を用いてシミュレーションを実施する。これにより、本実施形態では、より正確にシミュレーションすることが可能である。補間法の例は、三次スプライン(Cubic spline)補間である。
【００９１】
誘電率εの実部εｒ及び虚部εｉの値として、ある波長λに対する実測値が存在しない場合、例えば、パラメータを最適化したドルーデの分散式を利用して誘電率εを算出することも考えられる。
【００９２】
しかしながら、図６に示すように、補間値を利用する方が、より正確に誘電率εを算出可能である。
【００９３】
図６は、銀の誘電率の実部の三次スプライン補間値を示す図面である。図６には、比較のため、最適化されたパラメータを使用したドルーデの分散式によって算出された値も示している。図６中の横軸は波長（ｎｍ）を示し、縦軸は、誘電率の実部（Ｆ／ｍ）を示す。図６中のマーク「×」は、実測値を示す。具体的には、「C.L.Foiles, “Optical properties of pure metals and binary alloys”,Chapter 4 of Landolt-Bornstein Numerical Data and Functional Relationships inScience and Technology New Series, Vol.15, Subvolume b, K. H.Hellwege and J.L.Olsen, Ed.Springer-Verlag, Berlin 1985,pp. 228」記載の値である。図６中のマーク「○」は、三次スプライン補間を利用し、実測値から直接算出した補間値である。図６中のマーク「●」は、最適化したドルーデパラメータ（プラズマ周波数ωｐ及び衝突周波数νｃ）を利用してドルーデの一次分散式から得た値である。
【００９４】
ドルーデモデルにおける最適化パラメータは、遺伝的アルゴリズム（例えば、「S. Banerjee and L. N. Hazra, “Experiments with a genetic algorithmfor structural design of cemented doublets with prespecified aberration targets”,App. Opt. Vol. 40, No. 34, 2001, pp. 6265.」参照）を利用して算出した。具体的には、ドルーデモデルにおけるパラメータであるプラズマ周波数ωｐ及び衝突周波数νｃで表される誘電率の実部及び虚部を含むメリット関数の最適化を、遺伝的アルゴリズムを用いて行う。図６のマーク「●」で示した値は、プラズマ周波数ωｐを１．２８８×１０１６ｓ−１とし、衝突周波数νｃを６．４７１３×１０１３ｓ−１として算出した値である。
【００９５】
図６から理解されるように、三次スプライン補間を利用した方が、最適化パラメータを含むドルーデの分散式で算出される値より、実測値から得られる傾向（例えば、実測値のフィティング曲線）により近い値であることが理解され得る。特に、短波長側で、三次スプライン補間を利用した値の方が、実測値から得られる傾向により近い値である傾向にある。
【００９６】
従って、波長λに対応する実測値がある場合は、その実測値を採用し、波長λに対応する実測値がない場合には、三次スプライン補間を利用した補間値を採用することで、シミュレーションをより正確に行うことができる。
【００９７】
図７は、シミュレーション結果の一例を示す図面である。図７では、金属微粒子モデル２０Ａ近傍における金属微粒子モデル２０Ａによる光の散乱をシミュレーションした結果であり、ｘｚ面のシミュレーション結果を示している。図７において、濃淡は光の強度を示しており明るいほど光の強度が強いことを示す。図７に示したシミュレーションの条件は次の通りである。
・金属微粒子１２の材料：銀
・金属微粒子１２の粒子径：２０ｎｍ
・金属微粒子１２の透磁率μ：４π×１０−７［Ｈ／ｍ］
・媒質１４：真空（すなわち、金属微粒子１２は真空中に配置されているとした）。
・媒質１４の誘電率ε：ε＝ε０＝１／μｖ２  なお、本シミュレーションでは、ｖは、真空における光速（２．９９８×１０８ｍ／ｓ）である。
・媒質１４の透磁率μ：４π×１０−７［Ｈ／ｍ］
・セル２２の形状：立方体
・セル２２の一辺の長さΔｓ：４ｎｍ
・入射光１６の波長λ：４４０ｎｍ
・ｘ軸方向のセル数：１０１
・ｙ軸方向のセル数：１０１
・ｚ軸方向のセル数：１３１
・時間ステップ単位Δｔ：４．７２×１０−１６ｓ
・ステップ数ｎ：１８６６０
【００９８】
銀の波長４４０ｎｍに対する誘電率は、少なくとも、上記誘電率に対する説明で挙げた文献に記載されていないので、三次スプライン補間により誘電率を算出した。算出した誘電率は、−６．５１１＋０．１９３ｉであった。上記シミュレーション条件では、安定化条件である式（３４）を満たす。
【００９９】
図７に示すように、金属微粒子１２の近傍における電磁界を計算できている。すなわち、安定してＦＤＴＤ法を実行できていることが理解され得る。
【０１００】
以上、本発明の種々の実施形態について説明したが、本発明は、上記実施形態に限定されず、本発明の趣旨を逸脱しない範囲で種々の変形が可能である。
【０１０１】
例えば、図２に示した形態では、セル２２の大きさは、金属微粒子モデル及び媒質モデルにおいて同じであり、式（４９）を満たしていた。しかしながら、金属微粒子モデルに対するセル２２及び媒質モデルに対するセル２２それぞれの大きさが異なっても良い。具体的には、金属微粒子モデルに対するセル２２に対しては、式（４９）を満たすようにΔｓを設定し、媒質モデルに対するセル２２については、媒質モデルに対する式（１９ａ）〜式（１９ｃ）及び式（２４ａ）〜式（２４ｃ）に対する安定化条件（例えば、クーラン条件）を満たすように、Δｓを設定してもよい。このように、金属微粒子モデル及び媒質モデルそれぞれに対してセル２２の大きさを設定する場合、通常、媒質モデルに対するセル２２のΔｓが金属微粒子モデルに対するセル２２のΔｓより大きい。これは、金属微粒子モデルに対するセル２２の安定化条件がより厳しいからである。従って、金属微粒子モデル及び媒質モデルそれぞれに対してセル２２の大きさを設定する場合、金属微粒子モデルに対するセル２２のΔｓで解析モデル２０のＦＤＴＤ領域１８を分割する場合より、セル２２の個数が低減される。すなわち、ＦＤＴＤ法によって電磁界を計算するための処理時間の短縮を図ることが可能である。
【０１０２】
セル２２の形状は立方体としたが、セル２２の形状は直方体形状であればよい。セル２２の形状が直方体形状である場合、セル２２を構成する互いに直交する第１〜第３の辺（ｘ軸、ｙ軸及びｚ軸方向に延在する辺）の長さであるΔｘ＝Δｓ１、Δｙ＝Δｓ２及びΔｚ＝Δｓ３は、下記の式（５０ａ）〜式（５０ｃ）を満たすように設定されればよい。
【数５０】
2014026524.tif 000052
【０１０３】
更に、Δｔは、予め設定した値としたり、式（４９）又は式（５０ａ）〜式（５０ｃ）を満たすように、Δｔ及びΔｓを変更しながら決定さればよい。また、入射光（電磁波）は単色としたが、単色に限定されない。更に、媒質内に配置される物体として、金属微粒子を例示したが、誘電率の実部が０以下であり、虚部が正の値である物体であって、３次元的な形状を有するものであればよい。
【符号の説明】
【０１０４】
１０…解析対象物、１２…金属微粒子（物体）、１４…伝搬媒質（媒質）、１６…入射光（光）、１８…ＦＤＴＤ領域、２０…解析モデル、２２…セル、３４…シミュレーションプログラム。
【図１】
2014026524.tif 000054
【図２】
2014026524.tif 000055
【図３】
2014026524.tif 000056
【図４】
2014026524.tif 000057
【図５】
2014026524.tif 000058
【図６】
2014026524.tif 000059
【図７】
2014026524.tif 000060