Meep(FDTDシミュレータ)

前置き

原文ページはこちら

Meepは計算電磁気学の有限差分時間領域法(FDTD法)を実装します。
この手法は離散空間・離散時間での場の解析において、より細かく実際の連続な方程式に近いシミュレーションを行う際に広く用いられているものであり、
多くの実用的な問題を本質的かつ正確にシミュレートすることができます。

本項では、MeepのFDTD法で用いられる方程式・電磁気単位およびMeepにおけるFDTD法へのアプローチを紹介します。
電磁気学においては、様々な特殊な用途に適した多くの手法が存在しますが、FDTD法はそのうちの一つに過ぎません。
本項ではまた、いくつかの他の手法を説明し、FDTD法が適している場合・他の手法を考慮すべき場合を判断するヒントを提供します。

本項では、Meepのインターフェースを紹介しません。
本項の焦点はシミュレーションの対象の概要です。
インターフェースについてはこちらを参照してください。

マクスウェル方程式
- Meepにおける単位の取り扱い
境界条件と対称性
FDTD法
- Meepにおける離散性の隠蔽
その他の計算機的手法
FDTD法の応用

マクスウェル方程式

Meepは電磁界の相互作用を記述した方程式であるマクスウェル方程式をシミュレートします。
特に、場の展開の方程式は以下のように記述されます。

$\frac{d\vec{B}}{dt}=-\nabla\times\vec{E}-\vec{J}_B-\sigma_B\vec{B}$
$\frac{d\vec{D}}{dt}= \nabla\times\vec{H}-\vec{J}-\sigma_D\vec{D}$
$\vec{B}=\mu\vec{H}$
$\vec{D}=\varepsilon\vec{E}$

Dは変位場、εは誘電率、Jは電荷電流密度、J_Bは磁荷電流密度、Bは磁束密度、μは透磁率、Hは磁界を表します。
σ_Bおよびσ_Dはそれぞれ磁気伝導率及び電気伝導率に相当します。
発散方程式は暗黙的に以下の式で与えられます。

$\nabla\cdot\vec{B}=-\int^t\nabla\cdot\{\vec{J}_B(t')+\sigma_B\vec{B}\}dt'$
$\nabla\cdot\vec{D}=-\int^t\nabla\cdot\{\vec{J}(t')+\sigma_D\vec{D}\}dt'\equiv\rho$

一般に、εは位置だけでなく周波数(物質の分散)や自身内部の電界B(非線形)にも依存し、損失や利得を含むこともあります。
Meepはこれらの効果をサポートしており、Meepにおける材質の取り扱いにて解説されています。

Meepは円筒座標系におけるシミュレーションをサポートしています。
詳しくはMeepにおける円筒座標系を参照してください。

Meepにおける単位の取り扱い

Meepにおいて、ε₀, μ₀, cといった煩雑な定数が欠如していることにお気づきでしょうか。
これは、Meepが無次元単位系を利用しておりすべての定数が統一されているためです。
実際問題として、計算すべきほとんどすべてのものは比によって表現されるものであり、最終的に単位は解消されます。

特に、マクスウェル方程式がスケール不変であるために、スケール不変の単位を選択することは電磁気の問題において都合がいいです。
ゆえに、我々は本システムにおいていくつかの特徴的な長さスケールaを選択し、それらを距離の単位として利用しました。

さらに、Meepの単位系においては光速c=1であるため、a(正確にはa/c)は時間の単位でもあります。
特に、Meepにおける周波数fは常にc/aの単位で明示されており、これは単位a/cで表される光周期の逆数Tでありf=1/Tと等価です。
同様に真空中の波長λを用いてf=a/λとも表せます。

例えば、ミクロンサイズの赤外域ナノフォトニック構造を考えてみてください。
a=1μmであるとします。
光源波長を1.55μmとしたい場合、周波数fをf=1/1.55=0.6452と設定します。
時間軸100点に渡りシミュレーションを行いたい場合は時間幅155(=100/f)に渡り実行すればいいことが分かります。

透過スペクトルは入射光に対する透過光の強度比になり、ゆえに電界Eの単位は非線形でない限り不適切なものとなるでしょう。

境界条件と対称性

コンピュータ上でのシミュレーションは有限の空間領域でしか実行できず、ゆえにいくつかの境界条件を設定してシミュレーションを打ち切らなければなりません。
Meepでは3タイプの標準的な終端形式をサポートしています。

Bloch周期境界
金属壁
PML吸収層

また、計算要件を大幅に軽減するために問題の対称性を利用することができます。

サイズLのセルにおける通常の周期的境界について、電界成分はf(x+L)=f(x)を満たします。
Bloch周期は、Blochの波数ベクトルkに対して

$f(x+L)=e^{ik_xL}f(x)$

で一般化したものです。
これを利用することで、MPBのように光学結晶や導波路などのモードについて解析することができます。

平坦で単純な境界条件が金属壁であり、セルが完全に金属で覆われているかのように境界面での電界が強制的に0になります。
より一般的に言えば、例えば任意形状の金属製空洞をシミュレートするために、評価セル中の任意の位置に完全な金属を設置することができます。

境界条件が開放である場合をシミュレートするために、入射波の全てを吸収する無反射な境界を設定したいことがあります。
これは完全適合層(Perfect Matched Layers/PML)によって実装されます。
厳密に言えば、PMLは境界条件ではなく、むしろ境界に配置された特殊吸収性材料であると言えます。
実際にはPMLは仮想的な材料であり、その境界面において無反射であると設定されています。
理論的連続系においてPMLは無反射ですが、実際の離散系においては不完全であり微小な反射を持ちます。
ゆえに、吸収が徐々に有効となるようPMLに有限の厚みを設定する必要があります。
より多くの情報が必要な場合は完全適合層を参照してください。

評価セルのサイズを低減する別の方法として、対称性を利用する手法があります。
例えば、系が(構造及び光源の両者に対して)鏡対称面を有していることが既知である場合、構造の半分のみをシミュレートしもう半分を鏡面反射によって取得することで計算量を半減させることができます。
Meepではいくつかの鏡面対称性・回転対称性を利用することができます。
これらサポートされた対称性は完全な最適化を施され、計算に対称性を指定する他は全く同じ方法で設定されています。
詳しくはMeepにおける対称性の利用を参照してください。

FDTD法

FDTD法は空間および時間を有限の長方形格子に分割します。
後述のように、Meepは可能な限り離散性をユーザーから隠蔽しますが、離散化の影響がいくらか存在しており、それらに精通しておくことが好ましいと言えます。

恐らく、理解されるべき最大の要点は、格子が空間分解能Δxを持つ場合、離散時間ステップΔtはΔt=SΔxで与えられるということです。
ただし、Sは

$S<n_{min}/\sqrt{\#\mathrm{dimensions}}$

を満たす定数であり、n_minは最小の屈折率(通常1)です。
これは、FDTD法が安定である(発散しない)ために必要です。
(Meepにおいては[1～3次元に十分なように]S=0.5がデフォルトで与えられていますが、ユーザーが自由に設定可能です)
ゆえに格子分解能を2倍にした場合は(同じシミュレーション時間幅に対して)時間ステップも同様に2倍になります。
そして3次元空間においては、分解能を2倍にするとメモリ使用量が8倍になり、計算時間に至っては16倍にもなってしまいます。

理解を要する次の点は、方程式を二次の正確性で離散化するために、FDTD法は各格子それぞれに別々の電界成分を割り当てるという点です。
この離散化はYee格子として知られるものです。
結果として、結合・比較・電界成分の出力を行う際(例えばエネルギー密度と流を計算するとき)にMeepは共通の格子点に電界成分を補完します。
この処理は自動であるため大抵は補完処理に気を使う必要はありません。
しかし、E,Dは誘電境界で不連続かもしれないにも拘らず単純な線形補完が施されるため、補完されたE,Dの精度は誘電境界付近で期待されるものに劣る可能性があります。

電磁気学のためのFDTD法に関する多くの参考文献が利用できます。
例えば、以下の文献を参照してください。

Allen Taflove and Susan C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method (Artech: Norwood, MA, 2000).

Meepにおける離散性の隠蔽

FDTDは本質的に離散化空間・時間を利用していますが、Meepは可能な限りあたかも連続的なシステムであるかのように振る舞います。
シミュレーションを始めるにあたって、ユーザーは空間分解能を指定しますが、一般的にこれ以降は自身で選択した単位に従い連続座標系で作業することになります。

例えば、誘電率を連続変数xで定義される関数ε(x)であったり、あるいは球や円筒などのオブジェクトとして設定した場合、Meepはそれらの誘電率分布が離散的なグリッド上でどのように表現されるか見つけ出す必要があります。
もしくは、単に点xの位置に点光源を設定したい場合、Meepはxの至近点を探し、それらにxからの距離に応じて電流を重み付けして配分します。
xを連続的に変化させるならば、Meep内で処理される電流も(重み付けを変化させることによって)連続的に変化します。
特定の矩形を通る電磁束を求める場合、Meepはグリッド上の場から矩形中の電磁界を線形補完します。

一般に、Meepにおける界面の原理は広汎的な補完であり、入力を連続的に変化させたとしても、Meepでのシミュレーション結果も同じく連続的に変化し、また分解能を向上するに従いシミュレーション結果は可能な限り高速かつ滑らかに連続解へと収束します。

例えば、Meep内部で利用される誘電率の関数は、ユーザーが設定した関数を単純に離散サンプリングしたものではありません。
むしろ、各格子点は周囲のピクセルが持つ誘電率の平均に近い値を持ちます。
Meepにおけるサブピクセル平均は"ジャギー"やその他シャープな界面に起因する誤差を最小化するよう特別に設計されており、これは過去存在したFDTD法の実装に対し有意な改善であると言えるでしょう。
詳しくはMeepの引用に掲載されているFarjadpourらの論文を参照してください。

その他の計算機的手法

当然ながら、FDTD法は計算電磁気学における唯一の方法という訳でも、常に最善の選択肢という訳でもありません。
一般に、引き出しは多いほうがいいし、目的に沿った最善の策を選ぶことも重要です。
(詳しくは私たちのオンラインテキスト付録Dを参照してください)

例えば、FDTD法で電磁固有モード(後述)を計算することができますが、無損失構造においては、一般的にMPBパッケージのような固有モードに特化したソルバーの方がより高速かつ容易であり、さらには信頼性も高いです。
MPBマニュアルの周波数領域と時間領域や後述の共振モードも合わせて参照してください。

単一周波数における構造の電磁界分布・応答を計算するに当たり、時間の反復よりも対応する線型方程式から直接求める方がむしろ効率的かもしれません。
実際に、私たちはMeepにこの手法(有限差分周波数領域法/FDFD法)を試験的に直接実装してみました。
詳しくはMeepの周波数領域法を参照してください。
とはいえ、特に(例えば金属薄膜を伴っているような)大きなスケール差があるような場合、有限要素法や境界要素法など、空間領域における可変解像度を利用できる別の手法を用いる方がいいかもしれません。
境界要素法は、巨大容積内における微小物体による散乱を計算するときなど、大きな体積対面積比を持つときに特に強力です。

時間領域法の強みは、パルス応答に対するフーリエ変換やHarminvなどのより洗練された単一プロセスにより、ただ一度のシミュレーションで周波数スペクトル全体(あるいは固有周波数)を取得できる点にあります。
有限要素法もまた時間項を持つ電磁界に対し有効ですが、有限差分法と比べると、一般的に安定性のため暗黙の時間ステップを利用しなければならず、また全てのステップにおいて(線形系の解析のために)行列を反転しなければならないという重大な欠点を持っています。

最後の手法として、断面が一定な導波路の連続、円筒の集合、或いは多層膜など、解析が容易な少数の部品からなる系については、転送/散乱行列法は特に魅力的であると言えるでしょう。
これらの手法は個々の単純な要素をいくつかの分析的または半分析的な様式として扱い、これにより全体の構造に対し非常に高速かつ正確なシミュレーションが可能となります。
ここに簡潔に記述するにはこのような手法は多すぎますが、広汎な構造を扱える有用なフリーツールを一つ紹介するとすればCAMFRが挙げられるでしょう。

(その他の標準的手法として光伝搬法[BPM]がありますが、これは構造が一次元的にゆっくりと変化する場合にのみ適しています)

FDTD法の応用

本節では、FDTD法によって電磁気問題を分析できるいくつかの基本的手法を概説します。
これらの手法をMeepで使うためのより詳細な例についてはチュートリアルを参照してください。

電磁界の分布とグリーンの関数

時間領域シミュレーションが可能な最も明白なことは、与えられた光源に対する電磁界分布をただ図示することであり、もしくは界の時間進行を動画化することでしょう。

特定周波数ωの局所光源による電磁界分布は系の持つグリーンの関数により与えられます。
より限定的には、電流源分布が点x'に対し密度Jで与えられる場合に、点xにおける電界Eの一次の成分を与える典型的な『動的』グリーン関数のひとつは次の式で表されます。

$G_{ij}(\omega;\mathbf{x},\mathbf{x}')$

ただし、Jは次の式で与えられているものとします。

$\mathbf{J}(\mathbf{x})=\hat{\mathbf{e}}_j\cdot\exp(-i\omega t)\cdot\sigma(\mathbf{x}-\mathbf{x}')$

これをFDTD法で得るためにすべきことは、点x'に必要な点光源を置き、他の全ての周波数成分が消えるまで十分な時間待つだけです(単にt=0に電流源をオンにすることは周波数スペクトルを導くことに注意)。
或いは、Meep周波数領域ソルバーを用いて(関連する線型方程式を解くことにより)直接的に周波数応答を得ることも可能です。

グリーンの関数について考えると、放射流束から(Σ_iG_iiに比例する)局所的な状態密度、小さな散乱によるボーン近似に至るまで有用なことを幅広く計算することができます。
さらに強力なことに、後述のパルス応答をフーリエ変換することで多くの多重周波数を計算できます。

透過・反射スペクトル

FDTD法を最も多く利用する目的は、おそらく、励起に応答する共振空洞などの有限構造の透過・散乱スペクトルを計算することでしょう。
もちろん、前述のとおり各周波数ωに対して別個に電磁界(および透過流束)を計算することもできます。
しかし、ショートパルスに対する単一のフーリエ変換演算を通して伝搬スペクトル応答を計算することはより効果的であると言えます。

例えば、ある構造の透過パワーを求める場合を考えます。
与えられた周波数ωに対する電磁界について、透過パワーは(法線n方向の)ポインティングベクトルの、構造の向こう側にある平面上での積分で与えられます。

$P(\omega)=\mathrm{Re}\left[\hat{\mathbf{n}}\cdot\int\mathbf{E}_\omega(\mathbf{x})\times\mathbf{H}_\omega(\mathbf{x})d^2\mathbf{x}\right]$

この時、ショートパルスを入力することを考えるとき、各時刻tに対してポインティングベクトルの積分P(t)を計算し、P(ω)を求めるためにフーリエ変換するという操作は魅力的です。
しかし、求めたいものがフーリエ変換された電界Eおよび磁界Hであり、(電磁束は電磁界の線形な関数ではないため)電磁束の変換とは一致しないことを考えると、これが誤りであると分かります。

代わりにすべきことは、(離散)時間ステップnに渡り総和を取ることで、流束平面の全点に対しフーリエ変換E_ω(X)およびH_ω(x)を累算することです。

$\tilde{f}(\omega)=\frac{1}{\sqrt{2\pi}}\sum_ne^{i\omega n\Delta t}f(n\Delta t)\Delta t \approx \frac{1}{\sqrt{2\pi}}\int e^{i\omega t}f(t)dt$

そして、時間ステップの終了後、フーリエ変換された電磁界から流束を求め、P(ω)を計算すればば透過パワーを得ることができます。
もちろんMeepはユーザーのためにこれらの問題を自動的に全て取り扱うので、ユーザーはただ流束を積分すべき範囲と計算すべき周波数を設定するだけで問題を解決できます。

(もちろん、他の時系列分析も使えるでしょう。
時に有効である一つの方法は、複数の点で界の時系列のPadé近似を構築することです。
そこからは多くの場合比較的短い時系列から鋭いピークと他の共振特徴を含む非常に正確な離散フーリエ変換を推定することができます。
MeepはPadé近似法を提供しませんが、任意の位置について電磁界の時間変化を出力することができ、標準的な手法を用いることでPadé近似を計算することができます。
ただし、理想的には、ある点を経由した波による透過スペクトルを求めるため、シングルモード導波路に対する解析が望ましいです)

パワーP(ω)は余り便利であるとはいえません。
透過スペクトルを得るためには、各周波数における入射パワーで除算して正規化する必要があります。
典型的に、シミュレーションの実行によりこれは二度行われます。
正規化を利用した計算を行う必要が出た際、一度目は無散乱構造と入射波のみで、二度目は散乱構造も用いて行われます。

これは、透過スペクトル同様に反射スペクトルを求めようとした際より扱いにくいものとなります。
パワーP(ω)は入射波と反射波の合計であるため、単純に逆方向流束を計算することはできません。
更には、(入射波と反射波の間に生じる)除外されない干渉効果を含んでいるため、透過強度を求めるにあたり、逆方向りゅそくから入射波を単純減算することもできません。
反射・散乱強度を得るには入射電磁界E_ω⁽⁰⁾(x),H_ω⁽⁰⁾(x)を減算する必要があります。

$P_r(\omega) = \mathrm{Re }\left[ \hat{\mathbf{n}} \cdot \int\{\mathrm{E}_\omega(\mathbf{x})-\mathbf{E}_\omega^{(0)}(\mathbf{x})\}^*\times\{\mathbf{H}_\omega(\mathbf{x})-\mathbf{H}_\omega^{(0)}(\mathbf{x})\}d^2\mathbf{x} \right]$

繰り返しになりますが、実際には散乱がある場合とない場合の二度シミュレーションを行い、流束計算前に反射面からフーリエ変換を減算するよう設定することで簡単にこれを行うことができます。
また、反射強度を求めたのち、反射スペクトルを得るために入射強度で正規化することになるでしょう。

何が起きているかの見当がついている限り、Meepはこれらの計算を容易に行えるよう設計されています。

共振モード

FDTD法のもう一つの用途は、与えられた構造の共振モードまたは固有モードを計算することです。
例えば、光学結晶(断続的誘電構造)や導波路があり、波数kに対する共振モード(ωで定義される)を求めたいとします。
もしくは、光を長時間捉える共振空洞があり、共振周波数ωと崩壊寿命(線質係数)Qを求めたいとします。
こういった場合には往々にしてそれらモードにおける電磁界分布も求めたいこともあるでしょう。

FDTD法で周波数や崩壊寿命を引き出すにあたり、標準的な方策は単純です。
閉鎖系か開放系かに応じてBloch周期境界または吸収層を持った構造を設定します。
そして、空洞・導波路・その他の構造の内部に直接配置した電流源から短パルス(伝搬帯域)でモードを励起します。
最後に、系内部で振動する電磁界が得られるので、周波数と減衰率を抽出するためにこれを分析します。

共振分析のもっとも単純な形式は、複数の点で電磁界のフーリエ変換を計算することでしょう。
共振モードはスペクトル中に鋭いピークを引き起こしているはずです。
しかしながら、この方法は、高周波解析に膨大な時間を要するという深刻な欠陥を抱えており、更には減衰率を抽出する問題は不完全条件非線形フィッティング問題につながります。
そのかわりMeepは、NMR分光分析法から借用したより洗練されている信号処理アルゴリズムを提供します。
これはフィルタ対角化と呼ばれるアルゴリズムで、私たちが提供しているHarminvパッケージで実装されています。
Harminvは、周波数の全てとそれらの減衰率(と強度)を短時間で正確に引き出します。
例えば、わずか数百回の演算で10⁹の点から寿命Qを探し出します。

モードの周波数が分かったのち、電磁界分布が欲しい場合は、問題とする周波数のみを励起するための狭帯域(長時間)パルスで再度シミュレーションを行う必要があります。
(他のモードが十分減衰するだけ長大なシミュレーションを行ったうえで、最長寿命モードが欲しい場合はこの限りではありませんが)
電磁界分布が与えられた場合にはほかの分析法を利用することもできます。
(例えば、電磁束計算を介して異なる方向の減衰率にQを分解する、モーダルボリュームを見つける、など)

モード計算の際、なぜMPBの代わりにMeepを使うべきなのでしょうか。
MPBと違い、Meepは金属や吸収性材料をサポートしており、損失のある(共振)モードを計算でき、単一の短パルスから膨大な数の周波数について高速かつただ一度で計算でき、スペクトル(例えばバンドギャップ)の中から効率的にモードを抽出することができます。
では、なぜMPBも使うべきなのでしょうか。
MPBはMeepよりも最低周波数モードを高速に計算でき、周波数と電磁界を同時に出力でき、密集した周波数領域を解決するにあたり問題ありません。
さらには、時間領域でモードを計算することは幾分捉え難いものであり、十分な注意が必要です。
例えば、光源がモードにほぼ直行するように配置されているか、他のモードにあまりにも近い場合、モードを見失ってしまうことがあります。
逆に、信号処理はときに誤ったピーク周波数を特定することもあります。
また、非矩形格子状のセルを持った周期系を勉強することは、MPB上でするよりもMeep上でのほうがより煩雑です。
MPBは、いくつかの点で、より限定されたとはいえ、はるかに簡単で信頼性が高くなります。

タグ：

+ タグ編集

「前置き」をウィキ内検索

最終更新：2013年10月25日 14:13