FETI-DP法とマクスウェル応力テンソルを用いたナノスケール光学トラッピングの精密効率的シミュレーション

概要

数値シミュレーションは、プラズモニックナノ光学ピンセットベースの光学トラッピング予測に重要な役割を果たす。しかし、複雑な構造とプラズモニック効果による劇的な局所場増強が従来の数値手法に大きな課題をもたらす。本稿では、デュアル・プライマル有限要素ティアリング相互接続法（FETI-DP）とマクスウェル応力テンソルに基づく精密かつ効率的な数値シミュレーション手法を提案し、ナノ粒子トラッピングの光学力とポテンシャルを計算する。FETI-DPシミュレーション性能をさらに向上させるため、低ランク疎化アプローチを導入。提案手法は、非重複ドメイン分割と柔軟なメッシュ離散化により、大規模複雑問題を小規模単純問題に分解し、高効率と並列化を実現。数値結果は、ナノスケール光学トラッピングの予測と解析における提案手法の有効性を示す。

導入

表面プラズモン（SP）ベースのプラズモニック光学ピンセットは多くの注目を集め、ナノ粒子捕獲に広く応用されている[1-6]。SPは、特定波長の入射光と金属-誘電体界面の自由電子との結合により生じる共鳴現象[7]。SPにより光学ピンセットは回折限界を突破可能。さらに、SPの劇的な局所場増強により入射光強度要求を低減できる[7, 8]。しかし、SPは物質と物体寸法、入射光波長に密接に関連し、実践ではSP光学ピンセットの最適パラメータ決定に多数の実験を要する。このため、シミュレーション手法がSP光学ピンセットの設計最適化の補助手段としてますます重要になっている[9]。これらのシミュレーションでは、安定トラッピング予測のため光学力計算が必要。球体など規則的物体では、一般化ローレンツ・ミー理論から光学力を解析的に導出可能[10, 11]。しかし、複雑形状物体では、マクスウェル方程式を厳密に解く数値手法が電磁場モデリングとその後の光学力・ポテンシャル算出に必要。

これら数値手法は主に微分方程式法（DEM）と積分方程式法（IEM）に大別[12-15]。IEMと比較して、DEMは複雑幾何形状と成分の処理に優れる。DEMは近接場分布の直接計算にも優れ、SP解析に重要。代表的DEMとして、有限差分時間領域（FDTD）法は時間領域で実装され、広帯域情報と過渡応答を容易に取得可能[16, 17]。しかし、FDTDはSPにおける周波数依存材料特性を記述する精密分散モデルを必要とし、FDTD解精度はこの分散モデル近似精度に強く依存[18]。また、FDTDは構造化メッシュに依存し、曲面で階段誤差を生じやすい。もう一つの代表的DEMである有限要素法（FEM）は、周波数領域で分散材料を容易に扱え、非構造化メッシュで階段誤差を排除できるため広く採用[19-22]。FDTDと比較して、FEMは測定材料パラメータを近似分散モデルなしで直接採用可能。しかし、SPにおける劇的局所場増強はFEM離散化で細かいメッシュを必要とし、大型寸法・多数物体では未知数が劇的に増大。これら因子は悪条件行列系と巨大計算コストを引き起こし、SP増強光学トラッピング解析における従来FEMに大きな課題をもたらす。

本稿では、ナノスケール光学トラッピングシミュレーションのため、効率的なデュアル・プライマル有限要素ティアリング相互接続（FETI-DP）法を導入。FETI-DPは非重複ドメイン分解スキームを採用し、元の大規模複雑問題を一連の小規模単純問題に分割して克服。サブドメイン界面で伝達条件を課し場の連続性を保証し、ラグランジュ乗数により双対変数を導入して元の三次元（3D）問題を二次元（2D）問題に還元。サブドメインコーナーのプライマル変数を抽出し、双対問題の反復解収束率を加速[23-26]。低ランク疎化アプローチを開発し、FETI-DP性能を改善。データ疎アルゴリズムでサブドメイン問題と双対問題の解法効率を向上[27, 28]。提案手法は完全分離サブドメインを提供し、ナノ粒子トラッピング光学力の並列シミュレーションを可能にする。電磁場と力学的運動量の関係を明らかにするマクスウェル応力テンソル（MST）を採用し光学力を評価[29]。得られた光学力から光学ポテンシャルをさらに計算し、安定トラッピング解析を行う。IEMと比較して、提案手法は複合材料処理とSPベース光学トラッピングの近接場解析に強力。FDTDと比較して、SPベース光学トラッピングシステムでの分散金属材料を精密に扱え、曲面境界物体の階段誤差を排除。FEMと比較して、提案手法は光学トラッピング大規模計算に適する。複数例を分析し、数値結果がナノスケール光学トラッピング予測解析における提案手法の精度と効率を実証。

手法

FETI-DP定式化

FETI-DP実装では、元計算領域Ωをまず一連の非重複サブドメインΩⁱ（i=1,2,3⋯, N_s）に分割（図1参照）。各サブドメインΩⁱで、ベクトル波動方程式からサブドメイン有限要素系を導出：

$$ \nabla \times {\mu}_r^{-1}\nabla \times {\mathbf{E}}^i-{k}_0^2{\varepsilon}_r{\mathbf{E}}^i={jk}_0{\eta}_0{\mathbf{J}}_{\mathrm{imp}}^i \quad \mathrm{in} \quad {\Omega}^i \quad (1) $$

$$ \hat{n}\times \nabla \times {\mathbf{E}}^i+{jk}_0\hat{n}\times \left(\hat{n}\times {\mathbf{E}}^i\right)=0 \quad \mathrm{on} \quad {\Gamma}_{\mathrm{ABC}}^i \quad (2) $$

ここで、EⁱはΩⁱで解く未知電場、k₀とη₀は自由空間波数と固有インピーダンス、J_i^impは印加電流。Γ_ABCⁱは無限開領域を切断する吸収境界条件（ABC）。k₀は周囲媒質が自由空間でない場合（光学トラッピングで一般的）媒質波インピーダンスで置換。サブドメイン界面Γⁱで、Ωⁱで完全境界値問題を生成するため仮定境界条件が必要。ここでは未知補助変数Λ_iを用いたRobin型伝達条件を課す：

$$ {\hat{n}}^i\times \left({\mu}_r^{-1}\nabla \times {\mathbf{E}}^i\right)+{\alpha}^i{\hat{n}}^i\times \left({\hat{n}}^i\times {\mathbf{E}}^i\right)={\boldsymbol{\Lambda}}^i \quad \mathrm{on} \quad {\Gamma}^i \quad (3) $$

ここで、n̂ⁱはサブドメイン界面Γⁱでの外向き単位法線ベクトル、αⁱはしばしばjk₀と選べる複素パラメータ。全サブドメインを四面体要素で離散化。各要素で、ベクトル基底関数Nと未知電場係数EでEを展開：

$$ \mathbf{E}=\sum \limits_{p=1}^s{E}_p{\mathbf{N}}_p \quad (4) $$

ここで、sは各四面体要素のベクトル基底関数数。伝統的低次エッジベース基底関数ではs=6、高次ベクトル基底関数では面や体積ベース基底関数が追加されるためより大きい。

ガルキン法を適用し、Ωⁱでの未知電場係数EⁱについてのFEM行列方程式を取得：

$$ \left(\begin{array}{cc}{\mathbf{K}}_{rr}^i & {\mathbf{K}}_{rc}^i\\ {\mathbf{K}}_{cr}^i & {\mathbf{K}}_{cc}^i\end{array}\right)\left(\begin{array}{l}{E}_r^i\\ {E}_c^i\end{array}\right)=\left(\begin{array}{l}{f}_r^i-{\mathbf{B}}_r^{i^T}{\lambda}^i\\ {f}_c^i\end{array}\right) \quad (5) $$

添え字cとrでコーナー自由度（DOF）と残りDOFを区別し、コーナーDOFをプライマル変数として抽出しデュアル・プライマル（DP）スキームを構築。KはFEMシステム行列、fは励起ベクトル。Bはサブドメイン界面DOFを抽出するブール行列。λはΛ_i展開から生成される双対変数（ラグランジュ乗数）。

隣接サブドメインを、界面での接線電場と磁場連続性強制で相互接続。全サブドメイン界面をアセンブルし、全サブドメイン内部未知数Eⁱを消去。双対変数λについての簡約大域界面方程式を取得：

$$ \left[{\tilde{\mathbf{K}}}_{rr}+{\tilde{\mathbf{K}}}_{rc}{\tilde{\mathbf{K}}}_{cc}^{-1}{\tilde{\mathbf{K}}}_{cr}\right]\lambda ={\tilde{f}}_r-{\tilde{\mathbf{K}}}_{rc}{\tilde{\mathbf{K}}}_{cc}^{-1}{\tilde{f}}_c \quad (6) $$

式(6)は一般化最小残差（GMRES）法など反復法で解ける。K̃_ccは大域コーナー系で、プライマル空間での反復収束を加速。適切な前処理子（近似逆行列、不完全LU分解など）で反復収束率改善可能。双対変数λ解得後、各サブドメイン内電場を式(5)で独立評価可能。大域行列K̃_rr構築には、サブドメイン数値グリーン関数Z_rrⁱ構築が必要で、その形式は：

$$ {\mathbf{Z}}_{rr}^i={\mathbf{B}}_r^i{\left({\mathbf{K}}_{rr}^i\right)}^{-1}{{\mathbf{B}}_r^i}^T \quad (7) $$

ここでサブドメインFEM行列の逆行列<{(K_rrⁱ)}^{-1}>が含まれる。また、行列K̃_rc, K̃_cr, K̃_ccとベクトルf̃_r, f̃_cでも<{(K_rrⁱ)}^{-1}>計算が必要。前処理段階での<{(K_rrⁱ)}^{-1}>構築と反復解段階での行列ベクトル積（MVP）は計算コスト大。K_rrⁱは疎だが、<{(K_rrⁱ)}^{-1}>は密行列で高コスト。次に、<{(K_rrⁱ)}^{-1}>計算加速のため低ランク疎化手法を導入。大域界面系の一部部分配列が低ランク行列形式で表現可能なため、低ランクアルゴリズムで計算可能でFETI-DP性能改善。FETI-DPは独立サブドメイン操作を提供し、並列計算活用で効率改善可能。効率的並列スキームでは、各サブドメインDOF数を可能な限りバランスよくする原則。したがって、サブドメインサイズはメッシュ離散化密度に関連。通常、細かいメッシュ領域では小サブドメイン、粗いメッシュ領域では大サブドメインを採用。

低ランク疎化

FETI-DP効率改善のため、データ疎方式の低ランク疎化アプローチを提案。ここでは、これら行列が実際には疎ではないが、特定部分配列が低ランク分解行列形式で表現可能な意味で疎であることを指す：

$$ \mathbf{M}={\mathbf{XY}}^{\mathrm{T}} \quad \left(\mathbf{M}\in {\mathrm{\mathbb{C}}}^{m\times n},\mathbf{X}\in {\mathrm{\mathbb{C}}}^{m\times k},\mathbf{Y}\in {\mathrm{\mathbb{C}}}^{n\times k}\right) \quad (8) $$

ここでXとYはフル行列形式、ランクkはm,nよりはるかに小さい。行列<{(K_rrⁱ)}^{-1}>は積分作用素の行列特性を持つためデータ疎行列形式で表現可能。したがって、与えられたサブドメインで<{(K_rrⁱ)}^{-1}>が低ランク特性を持てば、低ランク疎化アプローチでデータ疎形式で効率的に計算・保存可能で、反復解MVPを加速。

低ランク疎化アプローチプロセスを以下ステップに分割：(1) 各サブドメインで基底関数セットを細分化しクラスタツリー構築、(2) 2クラスタツリー相互作用でブロッククラスタツリー構築、(3) 許容条件でK_rrⁱのデータ疎形式生成、(4) 低ランクフォーマットアルゴリズムで<{(K_rrⁱ)}_{\mathrm{DS}}^{-1}>のデータ疎表現取得、(5) データ疎アルゴリズムでFETI-DP系解法に移行。適切な前処理子で解法高速化可能。データ疎LU分解<{(K_rrⁱ)}_{\mathrm{DS}}={L_rrⁱ}_{\mathrm{DS}}{U_rrⁱ}_{\mathrm{DS}}>を採用し、行列逆計算<{(K_rrⁱ)}_{\mathrm{DS}}^{-1}}に置換。ネステッド分解技術で低ランク疎化効率をさらに改善。ネステッド分解はセパレータを用いて大きな非対角ゼロ部分配列を生成し、LU分解中にゼロを保持しフィルインを大幅削減。

<{(K_rrⁱ)}_{\mathrm{DS}}>生成のため、まず境界ボックスを用いてサブドメインエッジベース基底関数セットI={1,2,……N}を再帰細分化しクラスタツリーT_I構築。ネステッド分解で、対応境界ボックス内クラスタtを3後継{s₁, s_sep, s₂}に分割。s₁, s₂は切断された2境界ボックスのインデックスセット、s_sepはセパレータインデックスセット。図2aにこのプロセス簡単例示。次に、2クラスタツリーT_I相互作用でブロッククラスタツリーT_I×I構築（図2b）、これはガルキン法での元エッジベース基底関数セットクラスタツリーとテスト基底関数セットクラスタツリーとして選択可能。次に、ネステッド分解ベース許容条件を導入し、T_I×Iでフルブロック、低ランク分解ブロック、非対角ゼロブロックを区別[23]。これにより、K_rrⁱの非ゼロエントリで対応ブロック埋め、<{(K_rrⁱ)}_{\mathrm{DS}}>生成可能。最終的に、データ疎LU分解<{(K_rrⁱ)}_{\mathrm{DS}}={L_rrⁱ}_{\mathrm{DS}}{U_rrⁱ}_{\mathrm{DS}}>を以下から再帰計算：

$$ {\mathbf{K}}_{rr}^i=\left[\begin{array}{ccc}{\mathbf{K}}_{11}& & {\mathbf{K}}_{13}\\ {}& {\mathbf{K}}_{22}& {\mathbf{K}}_{23}\\ {\mathbf{K}}_{31}& {\mathbf{K}}_{32}& {\mathbf{K}}_{33}\end{array}\right]=\left[\begin{array}{ccc}{\mathbf{L}}_{11}& & \\ {}& {\mathbf{L}}_{22}& \\ {\mathbf{L}}_{31}& {\mathbf{L}}_{32}& {\mathbf{L}}_{33}\end{array}\right]\left[\begin{array}{ccc}{\mathbf{U}}_{11}& & {\mathbf{U}}_{13}\\ {}& {\mathbf{U}}_{22}& {\mathbf{U}}_{23}\\ {}& & {\mathbf{U}}_{33}\end{array}\right] \quad (9) $$

ここで従来フル行列演算をデータ疎対応物で置換[28]。適応切断誤差ε_tで低ランク近似精度制御。得られたLU因子<{(L_rrⁱ)}_{\mathrm{DS}}>と<{(U_rrⁱ)}_{\mathrm{DS}}>を保存し、以下でZ_rrⁱ構築：

$$ {\mathbf{Z}}_{rr}^i={\mathbf{B}}_r^i{\left({\mathbf{U}}_{rr}^i\right)}_{\mathrm{DS}}^{-1}{\left({\mathbf{L}}_{rr}^i\right)}_{\mathrm{DS}}^{-1}{\mathbf{B}}_r^i \quad (10) $$

ここで、B_rⁱ<{(U_rrⁱ)}_{\mathrm{DS}}^{-1}>と<{(L_rrⁱ)}_{\mathrm{DS}}^{-1}>B_rⁱはデータ疎上下三角ソルバーで計算可能。<{(L_rrⁱ)}_{\mathrm{DS}}>, <{(U_rrⁱ)}_{\mathrm{DS}}>, Z_rrⁱがデータ疎前進後退代入（FBS）とデータ疎MVPでFETI-DP計算に投入。

光学力とポテンシャル

電磁気学理論によれば、光学力は電磁場と力学的運動量の関係を明らかにするマクスウェル応力テンソル（MST）で評価可能[29]。物体周囲電磁場分布取得後、物体を囲む任意閉曲面上でMST積分により光学力計算可能。得られた電場分布ベースで、任意座標でのMSTを以下で構築：

$$ \overleftrightarrow{\mathbf{T}}=\frac{1}{2}\operatorname{Re}\left[\varepsilon {\mathbf{EE}}^{\ast }+\mu {\mathbf{HH}}^{\ast }-\frac{1}{2}\left(\varepsilon {\left|\mathbf{E}\right|}^2+\mu {\left|\mathbf{H}\right|}^2\right)\overleftrightarrow{\mathbf{I}}\right] \quad (11) $$

上付きアスタリスクは電場または磁場共役、εとμは誘電率と透磁率、Iは3×3単位行列。ベクトル外積により、Tのテンソル形式を以下で記述：

$$ \overleftrightarrow{\mathrm{T}}=\left[\begin{array}{lll}{T}_{xx}& {T}_{xy}& {T}_{xz}\\ {T}_{yx}& {T}_{yy}& {T}_{yz}\\ {T}_{zx}& {T}_{zy}& {T}_{zz}\end{array}\right]=\left[\begin{array}{ccc}\varepsilon {E}_x{E}_x^{\ast }+\mu {H}_x{H}_x^{\ast }-\frac{\varepsilon {\left|\mathbf{E}\right|}^2+\mu {\left|\mathbf{H}\right|}^2}{2}& \varepsilon {E}_x{E}_y^{\ast }+\mu {H}_x{H}_y^{\ast }}& \varepsilon {E}_x{E}_z^{\ast }+\mu {H}_x{H}_z^{\ast }\\ \varepsilon {E}_y{E}_x^{\ast }+\mu {H}_y{H}_x^{\ast }}& \varepsilon {E}_y{E}_y^{\ast }+\mu {H}_y{H}_y^{\ast }-\frac{\varepsilon {\left|\mathbf{E}\right|}^2+\mu {\left|\mathbf{H}\right|}^2}{2}& \varepsilon {E}_y{E}_z^{\ast }+\mu {H}_y{H}_z^{\ast }\\ \varepsilon {E}_z{E}_x^{\ast }+\mu {H}_z{H}_x^{\ast }}& \varepsilon {E}_z{E}_y^{\ast }+\mu {H}_z{H}_y^{\ast }}& \varepsilon {E}_z{E}_z^{\ast }+\mu {H}_z{H}_z^{\ast }-\frac{\varepsilon {\left|\mathbf{E}\right|}^2+\mu {\left|\mathbf{H}\right|}^2}{2}\end{array}\right] \quad (12) $$

添え字x,y,zは三方向成分。式(4)のE展開により、MSTエントリT_mn（m,n=x,y,z）をFETI-DP計算で展開形式に変換：

$$ {T}_{mn}=\sum \limits_{p,q=1}^s{E}_p{E}_q\left\{\varepsilon {\left({\mathbf{N}}_p\right)}_m{\left({\mathbf{N}}_q^{\ast}\right)}_n-\frac{1}{\omega^2\mu }{\left(\nabla \times {\mathbf{N}}_p\right)}_m{\left(\nabla \times {\mathbf{N}}_q^{\ast}\right)}_n\right. $$

$$ \left.-\frac{1}{2}\left[\varepsilon \left({\mathbf{N}}_p\right)\left({\mathbf{N}}_q^{\ast}\right)-\frac{1}{\omega^2\mu}\left(\nabla \times {\mathbf{N}}_p\right)\left(\nabla \times {\mathbf{N}}_q^{\ast}\right)\right]\right\} \quad \mathrm{if}\ m=n. \quad (13) $$

$$ {T}_{mn}=\sum \limits_{p,q=1}^s{E}_p{E}_q\left\{\varepsilon {\left({\mathbf{N}}_p\right)}_m{\left({\mathbf{N}}_q^{\ast}\right)}_n-\frac{1}{\omega^2\mu }{\left(\nabla \times {\mathbf{N}}_p\right)}_m{\left(\nabla \times {\mathbf{N}}_q^{\ast}\right)}_n\right\} \quad \mathrm{if}\ m\ne n. \quad (14) $$

ここでωは角周波数；Nとsは式(4)で説明済み。

最終的に、物体を囲む任意閉曲面上でMST積分により物体に働く光学力を計算：

$$ \mathbf{F}={\oint}_S\left(\overleftrightarrow{\mathbf{T}}\cdot \hat{n}\right)\ dS. \quad (15) $$

光学力計算も並列実装可能（MST積分が対応サブドメインに割り当てられるため）。安定光学トラッピング主要条件の一つは、勾配力が散乱力より大きいこと。言い換えれば、総力方向が勾配力方向と同一で、常に電場強度最強位置を指すこと。

光学ポテンシャルは光学トラッピング安定性を示すもう一つの魅力的パラメータ。得られた光学力ベースで、位置r₀での光学ポテンシャル深さUを以下で計算：

$$ \mathbf{U}\left({r}_0\right)=-{\int}_{\infty}^{r_0}\mathbf{F}\left(\mathbf{r}\right)\cdot \mathbf{r}, \quad (16) $$

添え字∞はゼロポテンシャル基準点として定義される無限大。U値はk_BTで表現可能（k_Bはボルツマン定数1.3806488×10⁻²³J/K、Tは環境温度）。一般に、U > 1 k_BT満たせば溶液中ブラウン運動克服し安定トラッピング可能。さもなければ安定トラッピング不可。総光学力は保存的勾配力成分と非保存的散乱力成分を含むため、式(15)の総光学力Fは非保存的[30, 31]。しかし、ナノ粒子運動が一次元制限されれば、式(16)から全光学力が非保存的でも光学ポテンシャルの明確定義が得られる。

結果と考察

提案手法有効性実証のため3例提示。貴金属が表面プラズモン励起に常用されるため、代表的金と銀材料を選択解析。第1例で銀ナノ粒子光学力計算し提案手法精度検証。第2,3例で金ナノ粒子光学トラッピングシミュレーションと考察。全例で無限領域をABCで切断、ABCと物体距離を1波長に設定し許容精度達成。全計算を3.6GHz Intel Xeonプロセッサ搭載Dellワークステーションで実施。

銀ナノカプセル

銀ナノカプセル物体をまず検討し、提案FETI-DP法の光学力予測精度と効率をテスト。図3a,bに構成と寸法示す。銀構成パラメータはすべて[32]実測値。FETI-DPスキーム実装で、全解析領域をまず24サブドメインに分割。金属表面近くでプラズモニック局所場増強効果モデル化のため密メッシュ必要。四面体要素採用で離散化、総計6.9×10⁵未知数（双対未知数4.1×10⁴、コーナー未知数313含む）。入射光は+z方向照射、電場偏光方向は-x。まず、入射光波長λを200nmから400nmまで変化させ、ナノカプセルに働く光学力シミュレーション。FETI-DPは周波数領域動作のため、15サンプリング周波数点で光学力計算。図4に銀ナノカプセルに働く光学力計算曲線示す。FETI-DP精度示すため、市販ソフトLumerical FDTD Solutions[33]結果と比較、良好な一致観察。次に、サブドメイン数違いでFETI-DP性能テスト。離散化密度維持しサブドメイン数4から24へ増加。各プロセッサに1サブドメイン割当。表1に大域界面式(6)構築時間と総解時間報告。FETI-DPが並列計算リソース完全活用し解効率大幅改善。また、サブドメイン数増加時のFETI-DP精度も検証表1報告。精度は光学力2ノルム相対誤差δ_OF = ‖OF_i − OF_ref‖/‖OF_ref‖で定義（OF_iはiサブドメイン使用光学力、OF_refは2サブドメイン使用基準光学力）。サブドメイン数増加で精度ほぼ一定維持。

金ナノスフィアダイマー

第2例で金ナノスフィアダイマー使用し金ナノスフィア光学トラッピング解析。ダイマーギャッププラズモニック効果がナノ粒子トラッピング光学力効果的増強。図5a,bに系構成と寸法示す。金構成パラメータすべて[32]実測値。周囲媒質は水（相対屈折率n=1.33）。入射光は平面波（パワー10mW/μm²）、電場偏光方向+x、入射方向-z。物体ナノスフィアに働く光学力をFETI-DP法で計算。FETI-DP実装で全計算領域32サブドメインに分割、四面体メッシュ離散化、3.5×10⁶未知数（双対未知数1.6×10⁵、コーナー未知数1738含む）。まず、様々プロセッサ数で提案FETI-DP並列性能テスト。表2に式(6)解時間と総解時間報告。並列計算スピードアップも表2提供。スピードアップ定義：

$$ \mathrm{Speed\ up}=\frac{T_1}{T_{N_p}} \quad (17) $$

T_{N_p}はN_pプロセッサ使用総ウォールクロック時間。FETI-DPが解効率大幅改善、良好並列スピードアップ示す。この大量未知数で全プロセッサ総メモリ使用量わずか57.2GB。次に、低ランク疎化アプローチ有効性検証。低ランク疎化でサブドメイン行列をデータ疎アルゴリズムで因子化しデータ疎行列として保存。構築時間とメモリ使用量はわずか18秒と0.5GB、従来行列アルゴリズムでは67秒と1.7GB。72%時間短縮と70%メモリ圧縮達成。メモリ使用量関連で後続MVPも70%時間短縮。次に、波長λを277nmから818nmまで変化させFETI-DP光学力計算テスト。実践では異波長入射光下光学力解析がプラズモニック共鳴波長探索（劇的場増強発生、最強光学力取得）にしばしば必要。ナノスフィア位置(0,0,20nm)と(0,0,-20nm)2ケース検討。図6a,bに異λでのナノスフィア受光学力計算結果プロット。最大光学力はλ=472nmで発生、プラズモニック共鳴波長。この共鳴波長光学力は非共鳴波長比ほぼ40倍増強。さらに、光学力は常にダイマーギャップ方向を指す（図6参照）、電場強度最強位置。また勾配力方向で物体トラップ。図7a,bに非共鳴波長λ=300nmと共鳴波長λ=472nmでの電場増強分布（xoz平面）示す。プラズモニック共鳴効果で電場強度ほぼ250倍増加。さらに、ナノスフィアを(0,0,-30nm)から(0,0,-17nm)へz軸沿い移動させ光学力と光学ポテンシャル計算。トラッピング力ポテンシャルで最も典型的興味深い挙動はz方向作用もの、ここではz軸沿い積分で軸方向トラッピングポテンシャル考慮。ナノ粒子運動が一次元制限されるため、式(16)から全光学力非保存的でも光学ポテンシャル定義明確。図8a,b示す通り、ナノスフィアダイマーギャップ移動に伴い光学力と光学ポテンシャル深度明確増加。(0,0,-17nm)位置で4.6 k_BT光学ポテンシャル深度生成、水中ブラウン運動克服し安定光学トラッピング達成十分。最後に、誘電基板効果テスト。基板あり/なしで光学力計算。両ケースでナノスフィア位置(0,0,-20nm)、入射波長共鳴波長選択。基板なしケース光学力|F₀|=0.769pN。基板ありケースで金ナノスフィアダイマーを厚さ60nm、相対誘電率ε_r=2.25誘電基板上配置。光学力|F₁|=0.761pN。2光学力結果相対誤差約1.0×10⁻²（|F₁ − F₀|/|F₀|で定義）。

金切り欠き円錐ダイマー

第3例で金切り欠き円錐ダイマー使用し金ナノスフィア光学トラッピング扱う。図9に系構成と寸法示す。金構成パラメータ[32]から採用。誘電基板相対誘電率ε_r=2.25。周囲媒質水（相対屈折率n=1.33）。入射光平面波（パワー10mW/μm²）、電場偏光方向+x、入射方向-z。全計算領域32サブドメイン分割、四面体メッシュ離散化、3.1×10⁶未知数（双対未知数1.3×10⁵、コーナー未知数1227含む）。まず、λを277nmから818nmまで変化させ光学力解析。図10に異λでのナノスフィア受光学力計算結果プロット。ナノスフィア位置(0,0,35nm)。最大光学力はλ=464nmで発生、プラズモニック共鳴波長、非共鳴波長比ほぼ30倍増強。さらに総光学力常に-z方向指す（図10参照）、勾配力方向。勾配力>散乱力確認、ナノスフィア安定トラッピング条件一つ満たす。図11a,bに非共鳴波長λ=300nmと共鳴波長λ=464nmでの電場分布（xoz平面）示す。局在表面プラズモン共鳴で電場強度ほぼ500倍増加。次に、ナノスフィア位置を0,5,35nmに変更し光学力観測。図12に明確なy成分光学力観測、より大きなz成分光学力存在。総光学力依然最強電場位置指しナノスフィアトラップ。さらに、ナノスフィアを(0,0,50nm)から(0,0,20nm)へz軸沿い移動させ光学ポテンシャル解析。ここでもz方向軸方向トラッピングポテンシャル考慮、ナノ粒子運動一次元制限で光学ポテンシャル明確定義。光学力とポテンシャル両方計算。図13a,b観察通り、ナノスフィアダイマーギャップ移動に伴い光学力と光学ポテンシャル深度明確増加。(0,0,20nm)で3.8 k_BT光学ポテンシャル深度取得、水中ブラウン運動克服し安定光学トラッピング達成十分。最後に、メッシュサイズ違いで未知数100万から320万まで変化させFETI-DP計算コストテスト。実践では異メッシュ密度テストが異精度要求満足に通常必要。大規模複雑問題は従来数値手法に大き課題。しかしFETI-DP容易対処。32プロセッサでFETI-DPシミュレーション、各プロセッサ1サブドメイン担当。表3にFETI-DP計算コスト報告。FETI-DPが高シミュレーション効率と低メモリ要求示す。

結論

金属ナノ粒子光学トラッピング予測解析のため、低ランク疎化組み合わせFETI-DP法を提案。提案手法は完全分離サブドメイン問題提供し、大規模複雑問題を一連の小規模単純問題に変換。並列計算に適し、数値シミュレーション効率大幅改善。例が示す通り、提案手法は大規模計算優秀性能示し、ナノスケール光学トラッピング高速精密シミュレーションに適する。