OghmaNano における FDTD

1. 導入

有限差分時間領域法（FDTD） は、離散空間グリッド上で時間領域において Maxwell 方程式を直接解くための全波数値手法です。レイトレーシング、モードソルバー、 transfer-matrix 法とは対照的に、FDTD は幾何学的仮定やモード仮定を置きません。 その代わりに、電場と磁場は curl 方程式、すなわち Ampère の法則と Faraday の法則を 離散化することで時間的に陽的に前進されます。その結果、電磁波はシミュレーション領域内を 自然に伝搬、散乱、干渉、回折します。

\[ \nabla \times \boldsymbol{H} = \sigma \boldsymbol{E} + \epsilon \frac{\partial \boldsymbol{E}}{\partial t}, \qquad \nabla \times \boldsymbol{E} = - \mu \frac{\partial \boldsymbol{H}}{\partial t}. \]

有限差分時間領域法（FDTD）では、Maxwell 方程式は、空間微分を千鳥配置グリッド （通常は Yee 格子）上の有限差分で置き換え、電磁場を時間発展させることで、時間領域で 直接解かれます。この手法は伝搬方向、モード構造、あるいは電磁場の挙動について本質的に 仮定を置かないため、単一の FDTD シミュレーションで広帯域応答、回折と干渉、 近接場効果、サブ波長幾何構造、過渡電磁現象を一つの枠組みで扱うことができます。 この単純化仮定の欠如こそが FDTD を非常に汎用的にしている理由ですが、同時に計算コストを 高くする理由でもあります。空間グリッドは最小の物理長さスケールを分解できなければならず、 時間刻みは数値安定性によって厳しく制約され、三次元シミュレーションは見かけ上は控えめな サイズの形状であっても急速にメモリと計算時間を消費します。

2. FDTD を使うべき場合と使うべきでない場合

実用的な光学問題の大半では、FDTD は最初に選ぶべき手法ではありません。 大規模光学系、長距離伝搬、あるいは幾何光学が支配的な状況では、 FDTD は極めて非効率な選択です。イメージング系、マクロなレンズ、 開口、およびミリメートルからセンチメートルスケールの光学アセンブリは レイトレーシングで扱う方がはるかに適切であり、薄膜太陽電池、 OLED スタック、コーティングのような平面多層構造は transfer-matrix 法により より正確かつ効率的にモデル化できます。これらの領域では、FDTD は 計算上の大槌のようなものであり、実行時間とメモリ使用量を大幅に増加させる一方で、 追加される物理的洞察はほとんど、あるいは全くありません。

実際には、FDTD は最後の手段となる手法とみなすべきであり、 より単純なモデルが本当に破綻する場合に限定して用いるべきです。その強みは、 サブ波長格子による回折、フォトニック結晶、メタサーフェス、 強く無秩序または粗い界面、そしてレイ光学や層状媒質モデルでは捉えられない 近接場結合効果のように、完全な電磁場を明示的に解く必要がある状況にあります。 そのような場合であっても、FDTD はその数値コスト、安定性制約、 メモリ要求を明確に理解した上で、限定的かつ意図的に適用すべきです。 したがって、以下の節では OghmaNano で使用される離散更新方程式を Ampère の法則と Faraday の法則から直接導出しますが、それは無差別な利用を 勧めるためではなく、FDTD に何ができるのか、そしてなぜ慎重に使うべきなのかを 明示するためです。

3. 理論的背景

マニュアルのこの節では、理解を助け、誤りを見つけやすくするために、 詳細な導出を含めて FDTD コードを完全に説明することを目的としています。

Ampere の法則は \[\sigma \boldsymbol{E} + \epsilon \frac{\partial \boldsymbol{E}}{\partial t} = \nabla \times \boldsymbol{H} = \begin{vmatrix} \hat{\boldsymbol{x}} & \hat{\boldsymbol{y}} & \hat{\boldsymbol{z}} \\ \frac{\partial}{\partial x} & \frac{\partial}{\partial y} & \frac{\partial}{\partial z} \\ H_{x} & H_{y} & H_{z} \end{vmatrix}\] で与えられます。

これを展開すると

\[\sigma E_{x} + \epsilon \frac{\partial E_{x}}{\partial t} = \frac{\partial H_{z}}{\partial y}-\frac{\partial H_{y}}{\partial z}\]

\[\sigma E_{y} + \epsilon \frac{\partial E_{y}}{\partial t} = -\frac{\partial H_{z}}{\partial x}+\frac{\partial H_{x}}{\partial z}\]

\[\sigma E_{z} + \epsilon \frac{\partial E_{z}}{\partial t} = \frac{\partial H_{y}}{\partial x}-\frac{\partial H_{x}}{\partial y}\]

\(\frac{\partial}{\partial y}=0\) の場合

\[\begin{split} &\sigma E_{x} + \epsilon \frac{\partial E_{x}}{\partial t} =-\frac{\partial H_{y}}{\partial z}\\ &\sigma E_{y} + \epsilon \frac{\partial E_{y}}{\partial t} = -\frac{\partial H_{z}}{\partial x}+\frac{\partial H_{x}}{\partial z}\\ &\sigma E_{z} + \epsilon \frac{\partial E_{z}}{\partial t} = \frac{\partial H_{y}}{\partial x} \end{split}\]

\(E_{x}\) について \[\begin{split} &\sigma E_{x} + \epsilon \frac{\partial E_{x}}{\partial t} =-\frac{\partial H_{y}}{\partial z}\\ &\sigma \frac{E_{x}^{t+1}[]+E_{x}^{t}[]}{2} + \epsilon \frac{E_{x}^{t+1}[]-E_{x}^{t}[]}{\Delta t} = -\frac{H_{y}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{y}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta z}\\ &\sigma \frac{E_{x}^{t+1}[]}{2} + \epsilon \frac{E_{x}^{t+1}[]}{\Delta t} = -\frac{H_{y}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{y}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta z}-\sigma \frac{E_{x}^{t}[]}{2}+\epsilon \frac{E_{x}^{t}[]}{\Delta t}\\ &\sigma \frac{E_{x}^{t+1}[]}{2} + \epsilon \frac{E_{x}^{t+1}[]}{\Delta t} = -\frac{H_{y}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{y}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta z}-\sigma \frac{E_{x}^{t}[]}{2}+\epsilon \frac{E_{x}^{t}[]}{\Delta t}\\ & \frac{\sigma \Delta t + 2 \epsilon }{ 2 \Delta t}E_{x}^{t+1}[] = -\frac{H_{y}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{y}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta z}-\sigma \frac{E_{x}^{t}[]}{2}+\epsilon \frac{E_{x}^{t}[]}{\Delta t}\\ & E_{x}^{t+1}[] = \left ( -\frac{H_{y}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{y}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta z}-\sigma \frac{E_{x}^{t}[]}{2}+\epsilon \frac{E_{x}^{t}[]}{\Delta t} \right ) \frac{2 \Delta t}{\sigma \Delta t + 2 \epsilon} \end{split}\]

\(E_{y}\) について \[\begin{split} &\sigma E_{y} + \epsilon \frac{\partial E_{y}}{\partial t} = -\frac{\partial H_{z}}{\partial x}+\frac{\partial H_{x}}{\partial z}\\ &\sigma \frac{E_{y}^{t+1}[]+E_{y}^{t}[]}{2} + \epsilon \frac{E_{y}^{t+1}[]-E_{y}^{t}[]}{\Delta t} = -\frac{H_{z}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{z}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta x}+\frac{H_{x}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{x}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta z}\\ &\sigma \frac{E_{y}^{t+1}[]}{2} + \epsilon \frac{E_{y}^{t+1}[]}{\Delta t} = -\frac{H_{z}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{z}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta x}+\frac{H_{x}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{x}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta z}-\sigma \frac{E_{y}^{t}[]}{2} + \epsilon \frac{E_{y}^{t}[]}{\Delta t}\\ &E_{y}^{t+1}[] = \left ( -\frac{H_{z}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{z}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta x}+\frac{H_{x}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{x}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta z}-\sigma \frac{E_{y}^{t}[]}{2} + \epsilon \frac{E_{y}^{t}[]}{\Delta t} \right ) \frac{2 \Delta t}{\sigma \Delta t + 2 \epsilon} \end{split}\]

\(E_{z}\) について \[\begin{split} &\sigma E_{z} + \epsilon \frac{\partial E_{z}}{\partial t} = \frac{\partial H_{y}}{\partial x}\\ &\sigma \frac{E_{z}^{t+1}[]+E_{z}^{t}[]}{2} + \epsilon \frac{E_{z}^{t+1}[]-E_{z}^{t}[]}{\Delta t} = \frac{H_{y}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{y}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta x}\\ &\sigma \frac{E_{z}^{t+1}[]}{2} + \epsilon \frac{E_{z}^{t+1}[]}{\Delta t} = \frac{H_{y}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{y}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta x}-\sigma \frac{E_{z}^{t}[]}{2} + \epsilon \frac{E_{z}^{t}[]}{\Delta t}\\ &E_{z}^{t+1}[]= \left ( \frac{H_{y}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{y}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta x}-\sigma \frac{E_{z}^{t}[]}{2} + \epsilon \frac{E_{z}^{t}[]}{\Delta t} \right ) \frac{2 \Delta t}{\sigma \Delta t + 2 \epsilon}\\ \end{split}\]

Faraday の法則は \[-\sigma_{m} \boldsymbol{H} - \mu \frac{\partial \boldsymbol{H}}{\partial t} = \nabla \times \boldsymbol{E} = \begin{vmatrix} \hat{\boldsymbol{x}} & \hat{\boldsymbol{y}} & \hat{\boldsymbol{z}} \\ \frac{\partial}{\partial x} & \frac{\partial}{\partial y} & \frac{\partial}{\partial z} \\ E_{x} & E_{y} & E_{z} \end{vmatrix}\] で与えられます。

これを展開すると次を得ます：

\[-\sigma_{m} H_{x} - \mu \frac{\partial H_{x}}{\partial t} = \frac{\partial E_{z}}{\partial y}-\frac{\partial E_{y}}{\partial z}\]

\[-\sigma_{m} H_{y} - \mu \frac{\partial H_{y}}{\partial t} = -\frac{\partial E_{z}}{\partial x}+\frac{\partial E_{x}}{\partial z}\]

\[-\sigma_{m} H_{z} - \mu \frac{\partial H_{z}}{\partial t} = \frac{\partial E_{y}}{\partial x}-\frac{\partial E_{x}}{\partial y}\]

\(\sigma_m=0\) かつ \(\frac{\partial}{\partial y}=0\) のとき

\[\begin{split} &\frac{\partial H_{x}}{\partial t} = \frac{1}{\mu} \left ( \frac{\partial E_{y}}{\partial z} \right )\\ &\frac{\partial H_{y}}{\partial t} = \frac{1}{\mu} \left ( \frac{\partial E_{z}}{\partial x}-\frac{\partial E_{x}}{\partial z} \right )\\ &\frac{\partial H_{z}}{\partial t} = - \frac{1}{\mu} \left ( \frac{\partial E_{y}}{\partial x} \right ) \end{split}\]

これを離散化すると

\[\begin{split} & H_{x}^{t+1} = \frac{1}{\mu} \left ( \frac{E_{y}^{t+\frac{1}{2}}[\frac{1}{2}]-E_{y}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta z} \right ) \Delta t + H_{x}^{t}[]\\ & H_{y}^{t+1} = \frac{1}{\mu} \left ( \frac{E_{z}^{t+\frac{1}{2}}[\frac{1}{2}]-E_{z}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta x}-\frac{E_{x}^{t+\frac{1}{2}}[\frac{1}{2}]-E_{x}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta z} \right ) \Delta t+ H_{y}^{t}[]\\ & H_{z}^{t+1} = \frac{1}{\mu} \left ( - \frac{E_{y}^{t+\frac{1}{2}}[\frac{1}{2}]-E_{y}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta x} \right ) \Delta t + H_{x}^{z}[] \end{split}\]

4. コンピュータ上での数値解法

上で導出した方程式は、電場と磁場の連続時間・連続空間における発展を記述しています。 これをコンピュータで解くには、空間と時間の両方を離散化する必要があります。 シミュレーション領域は規則的なセル格子に分割され、電場成分は離散的な空間位置と 離散的な時間ステップに格納されます。標準的な FDTD 形式では、電場と磁場は 空間と時間の両方で千鳥配置され、各場成分が curl 方程式において自然に必要となる 位置で評価されるようになっています。

実際には、これは電場成分が整数時間ステップ （\(t, t+\Delta t, t+2\Delta t\)）に格納され、 磁場成分が半整数時間ステップ （\(t+\tfrac{1}{2}\Delta t, t+\tfrac{3}{2}\Delta t\)）に 格納されることを意味します。 \(\partial / \partial x\) や \(\partial / \partial z\) のような空間微分は、 隣接する格子点間の中心差分で置き換えられます。この千鳥配置により、 電磁場が時間方向に互いを飛び越える単純な陽的更新スキームが得られます。

各時間ステップにおいて、まず磁場が Faraday の法則の離散化形式を用いて更新され、 電場の curl は現在の電場値から評価されます。磁場が半時間ステップ進んだ後、 電場は Ampère の法則の離散化形式を用いて更新されます。このとき、誘電率、透磁率、 導電率といった材料パラメータも含まれます。全体行列を解く必要はなく、 各格子セルにおける各場成分は、前の時間ステップにおける局所的な近傍値だけを使って 更新されます。

この陽的時間積分手順は、所望のシミュレーション時間に達するまで繰り返されます。 この方法は時間領域法であるため、単一のシミュレーションに本質的に広い周波数範囲の情報が 含まれます。スペクトルや定常場分布のような周波数領域量は、通常、時間の関数として場の値を 記録し、それらを Fourier 変換で後処理することで得られます。ソース、境界条件、 そして吸収層（完全整合層など）は、適切な格子位置で更新方程式に直接組み込まれます。

更新方程式が単純であることにより、FDTD は概念的に分かりやすく、高い並列化適性を持ちますが、 同時に厳しい数値制約も課されます。時間刻みは空間格子間隔と材料特性に依存する安定条件を 満たさなければならず、メモリ使用量は格子セル数と保存される場成分数に直接比例します。 これらの制約は、実際に実行可能なシミュレーションのサイズと時間長を最終的に制限し、 格子分解能、領域サイズ、モデル複雑性の慎重な選択を必要とします。