基于FETI-DP和麦克斯韦应力张量的纳米光学捕获精准高效模拟

摘要

数值模拟在基于等离子体纳米光学镊子的光学捕获预测中起重要作用。然而，复杂结构和等离子体效应的剧烈局域场增强给传统数值方法带来巨大挑战。本文提出一种基于对偶-原始有限元撕裂互联（FETI-DP）和麦克斯韦应力张量的精准高效数值模拟方法，用于计算捕获纳米颗粒的光学力和势。引入低秩稀疏化方法进一步提升FETI-DP模拟性能。该方法利用非重叠域分解和灵活网格离散化，将大规模复杂问题分解为小规模简单问题，展现高效率和并行性。数值结果证明了该方法在纳米尺度光学捕获预测和分析中的有效性。

引言

基于表面等离子体（SPs）的等离子体光学镊子备受关注，广泛用于捕获纳米颗粒。SP是特定波长入射光与金属-介质界面自由电子耦合引起的共振现象。SPs使光学镊子突破衍射极限。此外，SPs的剧烈局域场增强可降低入射光强度需求。然而，SPs与物体材料、尺寸及入射光波长密切相关，实际中需大量实验确定SP光学镊子最佳参数。因此，模拟方法作为SP光学镊子设计优化的辅助手段日益重要。这些模拟中需要计算光学力以预测稳定捕获。对于球体等规则物体，光学力可从广义洛伦兹-米理论解析导出。然而，对于复杂构型物体，求解麦克斯韦方程的严格数值方法对于建模电磁场及后续光学力和势是必要的。

这些数值方法主要分为微分方程方法（DEMs）和积分方程方法（IEMs）。相比IEMs，DEMs在处理复杂几何和组件方面表现优越。DEMs还具有直接计算近场分布的优势，这在SPs分析中起重要作用。作为代表性DEM，有限差分时域（FDTD）法在时域实现，易获取宽带信息和瞬态响应。然而，FDTD需精确色散模型描述SPs中频率相关材料性质，其解精度高度依赖该色散模型近似精度。此外，FDTD依赖结构化网格，常导致曲面阶梯误差。作为另一代表性DEM，有限元方法（FEM）因易于在频域处理色散材料并用非结构化网格消除阶梯误差而被广泛采用。相比FDTD，FEM可直接采用实测材料参数无需近似色散模型。然而，SPs中的剧烈局域场增强要求FEM离散化使用精细网格。大尺寸和多物体物体将显著增加未知数，导致病态矩阵系统和巨大计算消耗，给传统FEM分析SP增强光学捕获带来巨大挑战。

本文引入高效对偶-原始有限元撕裂互联（FETI-DP）方法模拟纳米尺度光学捕获。FETI-DP采用非重叠域分解方案，将原始大规模复杂问题分解为一系列小规模简单问题求解。它在子域界面强制传输条件确保场连续性，并引入对偶变量通过拉格朗日乘子将原始三维（3D）问题降为二维（2D）问题。在子域角提取原始变量以加速对偶问题迭代解收敛。开发低秩稀疏化方法提升FETI-DP性能。它利用数据稀疏算法提高求解子域问题和对偶问题的效率。该方法提供完全解耦子域，使纳米颗粒捕获光学力并行模拟成为可能。采用揭示电磁场与机械动量关系的麦克斯韦应力张量（MST）评估光学力。基于获得的光学力，可进一步计算光学势以分析稳定捕获。相比IEMs，该方法在处理复合材料和求解SP基光学捕获近场方面更强大。相比FDTD，该方法可准确处理SP基光学捕获系统中的色散金属材料并消除曲面边界物体的阶梯误差。相比FEM，该方法适合光学捕获大规模计算。通过若干例子分析，数值结果证明了该方法在纳米尺度光学捕获预测和分析中的精度和效率。

方法

FETI-DP公式化

对于FETI-DP实现，原始计算域Ω首先撕裂为一系列非重叠子域Ωⁱ（i=1,2,3⋯, Nₛ），如图1所示。在每个子域Ωⁱ中，可从矢量波方程导出子域有限元系统：

∇ × μᵣ⁻¹ ∇ × Eⁱ - k₀² εᵣ Eⁱ = j k₀ η₀ J_impⁱ in Ωⁱ (1)

n̂ × ∇ × Eⁱ + j k₀ n̂ × (n̂ × Eⁱ) = 0 on Γ_ABCⁱ (2)

其中Eⁱ为Ωⁱ中待求未知电场，k₀和η₀分别为自由空间波数和本征阻抗，J_impⁱ为激励电流。Γ_ABCⁱ为吸收边界条件（ABC）以截断无限开放区域。注意若周围介质非自由空间（光学捕获常见），k₀应替换为介质波阻抗。在子域界面Γⁱ，需假设边界条件以在Ωⁱ生成完整边界值问题。此处施加带未知辅助变量Λⁱ的Robin型传输条件：

n̂ⁱ × (μᵣ⁻¹ ∇ × Eⁱ) + αⁱ n̂ⁱ × (n̂ⁱ × Eⁱ) = Λⁱ on Γⁱ (3)

其中n̂ⁱ为子域界面Γⁱ处单位外法向量，αⁱ为复参数常选为j k₀。所有子域随后用四面体单元离散化。在每个单元中，用矢量基函数N和未知电场系数E展开E：

E = Σ_{p=1}^s E_p N_p (4)

s为每个四面体单元中矢量基函数数量。基于边的传统低阶基函数s=6，高阶矢量基函数s更大，因引入基于面或体的额外基函数。

应用Galerkin法，可得Ωⁱ中关于未知电场系数Eⁱ的FEM矩阵方程：

[Kᵣᵣⁱ Kᵣcⁱ; K_cᵣⁱ K_ccⁱ] [Eᵣⁱ; E_cⁱ] = [fᵣⁱ - Bᵣⁱᵀ λⁱ; f_cⁱ] (5)

下标c和r区分角自由度（DOFs）和其余DOFs，提取角DOFs作为原始变量构建对偶-原始（DP）方案。K为FEM系统矩阵，f为激励向量。B为提取子域界面DOFs的布尔矩阵。λ为从展开Λⁱ生成的对偶变量，也称拉格朗日乘子。

然后通过在界面强制切向电场和磁场连续性互联相邻子域。组装所有子域界面并消除所有子域内部未知数Eⁱ。可得关于对偶变量λ的约化全局界面方程：

[K̃ᵣᵣ + K̃ᵣc K̃_cc⁻¹ K̃_cᵣ] λ = f̃ᵣ - K̃ᵣc K̃_cc⁻¹ f̃_c (6)

方程(6)可用迭代法求解，如广义最小残量（GMRES）法。K̃_cc为全局角系统，可在原始空间加速迭代收敛。可采用合适预条件子提高迭代收敛率，如近似逆和不完全LU分解。一旦求解对偶变量λ，各子域内电场可由(5)独立求值。构建全局矩阵K̃ᵣᵣ需构建子域数值格林函数Zᵣᵣⁱ，形式为：

Zᵣᵣⁱ = Bᵣⁱ (Kᵣᵣⁱ)⁻¹ Bᵣⁱᵀ (7)

其中包含子域FEM矩阵(Kᵣᵣⁱ)⁻¹的逆。此外，矩阵K̃ᵣc、K̃_cᵣ、K̃_cc和向量f̃ᵣ、f̃_c也需要计算(Kᵣᵣⁱ)⁻¹。在预处理阶段构建(Kᵣᵣⁱ)⁻¹及其迭代求解阶段的矩阵-向量乘积（MVPs）计算量大。虽然Kᵣᵣⁱ稀疏，但(Kᵣᵣⁱ)⁻¹稠密，意味着高计算成本。接下来引入低秩稀疏化方法加速(Kᵣᵣⁱ)⁻¹计算。因全局界面系统中某些子矩阵可用低秩矩阵形式表示，其计算可用低秩算法执行，提升FETI-DP性能。FETI-DP提供独立子域操作，可利用并行计算提高效率。高效并行方案的域分解原则是使各子域DOFs数量尽可能平衡。因此，子域尺寸应与网格离散化密度相关。通常精细网格区域采用小子域，粗糙网格区域采用大子域。

低秩稀疏化

提出低秩稀疏化方法以数据稀疏方式提升FETI-DP效率。此处“数据稀疏”指这些矩阵虽非稀疏，但某些子块可用低秩分解矩阵形式表示：

M = X Yᵀ (M ∈ ℂ^{m×n}, X ∈ ℂ^{m×k}, Y ∈ ℂ^{n×k}) (8)

其中X和Y为全矩阵形式，秩k远小于m和n。矩阵(Kᵣᵣⁱ)⁻¹具积分算子矩阵性质，可用数据稀疏矩阵形式表示。因此，若(Kᵣᵣⁱ)⁻¹在给定子域具低秩性质，可用低秩稀疏化方法高效以数据稀疏形式计算和存储，加速迭代求解中的MVPs。

低秩稀疏化过程分步骤：(1)通过细分各子域基函数集构建簇树，(2)通过两个簇树交互构建块簇树，(3)通过可接受性条件生成Kᵣᵣⁱ的数据稀疏形式，(4)执行低秩格式化算法获取(Kᵣᵣⁱ)⁻¹的数据稀疏表示，(5)用数据稀疏算法进入FETI-DP系统求解。可采用合适预条件子加速求解。注意采用数据稀疏LU分解(Kᵣᵣⁱ)_DS = (Lᵣᵣⁱ)_DS (Uᵣᵣⁱ)_DS替代矩阵求逆(Kᵣᵣⁱ)_DS⁻¹。采用嵌套剖分技术进一步提升低秩稀疏化效率。嵌套剖分用分隔符产生大非对角零子块，在LU分解中保持零从而显著减少填入。

为生成(Kᵣᵣⁱ)_DS，首先用包围盒递归细分子域基于边的基函数集I={1,2,……N}构建簇树T_I。采用嵌套剖分，包围盒内簇t分为三个后继{s₁, s_sep, s₂}，s₁和s₂为两个不相连包围盒索引集，s_sep为分隔符索引集。图2a展示此过程简单例子。然后通过两个簇树T_I交互构建块簇树T_{I×I}，如图2b所示，可选为原始基于边基函数集的簇树和Galerkin法中测试基函数集的簇树。接着基于嵌套剖分引入可接受性条件区分T_{I×I}中的全块、低秩分解块和非对角零块。从而通过用Kᵣᵣⁱ非零条目填充相应块生成(Kᵣᵣⁱ)_DS。最后，可递归计算数据稀疏LU分解(Kᵣᵣⁱ)_DS = (Lᵣᵣⁱ)_DS (Uᵣᵣⁱ)_DS：

Kᵣᵣⁱ = [K₁₁ 0 K₁₃; 0 K₂₂ K₂₃; K₃₁ K₃₂ K₃₃] = [L₁₁ 0 0; 0 L₂₂ 0; L₃₁ L₃₂ L₃₃] [U₁₁ 0 U₁₃; 0 U₂₂ U₂₃; 0 0 U₃₃] (9)

其中常规全矩阵运算替换为其数据稀疏对应项。自适应截断误差ε_t控制低秩近似精度。获得的LU因子(Lᵣᵣⁱ)_DS和(Uᵣᵣⁱ)_DS存储并用于通过以下构建Zᵣᵣⁱ：

Zᵣᵣⁱ = Bᵣⁱ (Uᵣᵣⁱ)_DS⁻¹ (Lᵣᵣⁱ)_DS⁻¹ Bᵣⁱ (10)

其中Bᵣⁱ (Uᵣᵣⁱ)_DS⁻¹和(Lᵣᵣⁱ)_DS⁻¹ Bᵣⁱ可通过数据稀疏上/下三角求解器计算。(Lᵣᵣⁱ)_DS、(Uᵣᵣⁱ)_DS和Zᵣᵣⁱ随FETI-DP计算进入数据稀疏前向后向替换（FBSs）和数据稀疏MVPs。

光学力与势

根据电动力学理论，光学力可通过揭示电磁场与机械动量关系的麦克斯韦应力张量（MST）评估。一旦获得物体周围电磁场分布，可通过在包围物体的闭合面上积分MST计算光学力。基于获得的电场分布，任意坐标处MST可由下式构建：

T↔ = ½ Re[ε EE* + μ HH* - ½(ε|E|² + μ|H|²) I↔] (11)

上标星号表示电场或磁场共轭，ε和μ为介电常数和磁导率，I↔为3×3单位矩阵。通过矢量外积，T↔的张量形式可写为：

T↔ = [Tₓₓ Tₓᵧ Tₓ_z; T_yₓ T_yy T_yz; T_zₓ T_zy T_zz] = [εEₓEₓ*+μHₓHₓ*-(ε|E|²+μ|H|²)/2 εEₓEᵧ*+μHₓHᵧ* εEₓE_z*+μHₓH_z*; εEᵧEₓ*+μHᵧHₓ* εEᵧEᵧ*+μHᵧHᵧ*-(ε|E|²+μ|H|²)/2 εEᵧE_z*+μHᵧH_z*; εE_zEₓ*+μH_zHₓ* εE_zEᵧ*+μH_zHᵧ* εE_zE_z*+μH_zH_z*-(ε|E|²+μ|H|²)/2] (12)

下标x、y、z表示三个方向分量。根据(4)中E的展开，MST条目T_mn（m,n=x,y,z）在FETI-DP计算中可转换为展开形式：

T_mn = Σ_{p,q=1}^s E_p E_q {ε(N_p)_m (N_q*)_n - 1/(ω²μ) (∇×N_p)_m (∇×N_q*)_n - ½[ε(N_p)(N_q*) - 1/(ω²μ)(∇×N_p)(∇×N_q*)]} if m=n (13)

T_mn = Σ_{p,q=1}^s E_p E_q {ε(N_p)_m (N_q*)_n - 1/(ω²μ) (∇×N_p)_m (∇×N_q*)_n} if m≠n (14)

其中ω为角频率；N和s已在方程(4)描述。

最终，物体受光学力可通过在包围它的任意闭合面上积分MST计算：

F = ∮_S (T↔ · n̂) dS (15)

注意光学力计算也可并行实现，因MST积分分配给对应子域。稳定光学捕获主要条件之一是梯度力大于散射力。换言之，总力方向应与梯度力方向一致，后者总指向电场强度最强位置。

光学势是揭示光学捕获稳定性的另一重要参数。基于获得的光学力，位置r₀处光学势深度U可计算为：

U(r₀) = -∫_∞^{r₀} F(r) · r (16)

下标∞指无限远定义为零势参考点。U值可用k_B T表示，k_B为玻尔兹曼常数1.3806488×10⁻²³ J/K，T为环境温度。通常当U > 1 k_B T满足时，颗粒可克服溶液中布朗运动稳定捕获。否则无法稳定捕获。因总光学力包含保守梯度力分量和非保守散射力分量，(15)得总光学力F为非保守力。然而，若纳米颗粒运动限制在一维，(16)仍给出明确光学势定义，即使总光学力非保守。

结果与讨论

呈现三个例子演证所提方法有效性。因贵金属常用于激发表面等离子体，我们选择代表性金和银材料分析。第一个例子计算银纳米颗粒光学力验证所提方法精度。第二、三个例子模拟讨论金纳米颗粒光学捕获。所有例子中，无限域用ABC截断，ABC与物体距离设为一个波长，足以达到可接受精度。所有计算在配备3.6 GHz Intel Xeon处理器的Dell工作站上进行。

银纳米胶囊

首先考虑银纳米胶囊物体测试所提FETI-DP方法预测光学力的精度和效率。图3a、b给出其构型和尺寸。银本构参数均取自[32]实测值。实施FETI-DP方案，整个分析域首先分为24个子域。金属表面附近需更密网格以模拟等离子体局域场增强效应。采用四面体单元离散化，共6.9×10⁵未知数，含4.1×10⁴对偶未知数和313角未知数。入射光沿+z方向照射，电场偏振方向为-x。

首先，将入射光波长λ从200 nm变至400 nm模拟纳米胶囊受光学力。因FETI-DP在频域工作，在15个采样频率点计算光学力。图4显示银纳米胶囊受光学力随波长λ变化曲线。为指示FETI-DP精度，将FETI-DP光学力结果与商业软件Lumerical FDTD Solutions[33]结果比较，观察到良好一致性。

然后，测试不同子域数量下FETI-DP性能。保持离散化密度，将子域数从4增加到24。每个处理器处理一个子域。表1报告构建全局界面方程(6)和总求解时间。可见FETI-DP可充分利用并行计算资源显著提高求解效率。此外，子域数量增加时FETI-DP精度也经检验并报告于表1。此处精度定义为光学力2-范数相对误差δ_OF = ‖OF_i - OF_ref‖/‖OF_ref‖，OF_i为i个子域下光学力，OF_ref为两子域下参考光学力。可见精度随子域数量增加基本保持不变。

金纳米球二聚体

第二个例子利用金纳米球二聚体分析金纳米球光学捕获。二聚体间隙处等离子体效应可有效增强捕获纳米颗粒光学力。图5a、b给出该系统构型和尺寸。金本构参数均取自[32]实测值。周围介质为水，相对折射率n=1.33。入射光为平面波，功率10 mW/μm²，电场偏振方向+x，入射方向-z。FETI-DP方法计算物体纳米球受光学力。FETI-DP实现中，整个计算域分为32个子域，四面体网格离散化，产生3.5×10⁶未知数，含1.6×10⁵对偶未知数和1738角未知数。

首先，用不同处理器数测试所提FETI-DP并行性能。表2报告方程(6)求解时间及总求解时间，并提供并行计算加速比。加速比定义为Speedup = T₁ / T_{N_p}，T_{N_p}为使用N_p处理器总墙钟时间。可见FETI-DP显著提高求解效率并展现良好并行加速比。对此大量未知数，所有处理器总内存使用仅57.2 GB。

然后，考察低秩稀疏化方法有效性。采用低秩稀疏化后，子域矩阵可用数据稀疏算法分解并存储为数据稀疏矩阵。构建时间和内存使用仅18 s和0.5 GB，而常规矩阵算法为67 s和1.7 GB。可见获得72%时间节省和70%内存压缩。相关内存使用，后续MVPs也获70%时间节省。

接下来，FETI-DP测试波长λ从277 nm到818 nm变化的光学力计算。实践中，不同波长入射光下光学力分析常用于寻找等离子体共振波长，剧烈场增强发生处可获最强光学力。考虑纳米球位于(0,0,20 nm)和(0,0,-20 nm)两种情况。图6a、b绘制不同λ下纳米球受光学力计算结果。可见最大光学力出现在λ=472 nm，即等离子体共振波长。该共振波长下光学力较非共振波长增强近40倍。此外，光学力始终指向二聚体间隙（图6所示），即电场强度最强处，也是捕获物体的梯度力方向。图7a、b分别显示非共振波长λ=300 nm和共振波长λ=472 nm下系统在xoz平面电场增强分布。可见电场强度因等离子体共振效应增加近250倍。

此外，计算纳米球从(0,0,-30 nm)沿z轴移至(0,0,-17 nm)时的光学力和光学势。因沿z方向作用的捕获力和势最典型有趣，我们考虑沿z轴积分的轴向捕获势。因纳米颗粒运动限制在一维，(16)给出明确光学势定义，即使(15)总光学力非保守。如图8a、b所示，随纳米球移向二聚体间隙，光学力和光学势深度明显增加。在(0,0,-17 nm)位置，产生4.6 k_B T光学势深度，足以克服水中布朗运动实现稳定光学捕获。

最后，测试介电基底对此例影响。分别计算有无基底时光学力。两种情况纳米球均位于(0,0,-20 nm)，入射波长选为共振波长。无基底情况计算光学力|F₀|=0.769 pN。有基底情况，金纳米球二聚体置于厚度60 nm、相对介电常数ε_r=2.25介电基底上。计算光学力|F₁|=0.761 pN。两光学力结果相对误差约1.0×10⁻²，定义为|F₁-F₀|/|F₀|。

金截断锥二聚体

第三个例子处理利用金截断锥二聚体捕获金纳米球。图9给出该系统构型和尺寸。金本构参数取自[32]。介电基底相对介电常数ε_r=2.25。周围介质为水，相对折射率n=1.33。入射光为平面波，功率10 mW/μm²，电场偏振方向+x，入射方向-z。整个计算域分为32个子域，四面体网格离散化，产生3.1×10⁶未知数，含1.3×10⁵对偶未知数和1227角未知数。

首先，分析λ从277 nm到818 nm变化的光学力。图10绘制不同λ下纳米球受光学力计算结果。纳米球位于(0,0,35 nm)。可见最大光学力出现在λ=464 nm（等离子体共振波长），此处光学力较非共振波长增强近30倍。此外，总光学力始终指向-z（图10所示），即梯度力方向。这证实梯度力大于散射力，是纳米球可稳定捕获条件之一。图11a、b分别展示非共振波长λ=300 nm和共振波长λ=464 nm下系统在xoz平面电场分布。可见电场强度因局域表面等离子体共振增加近500倍。

然后，将纳米球位置改为0、5、35 nm观察光学力。图12给出计算光学力，可见明显y分量光学力，但存在更大z分量光学力。总光学力仍指向电场最强位置捕获纳米球。

进一步，分析纳米球从(0,0,50 nm)沿z轴移至(0,0,20 nm)时光学势。此处考虑沿z方向轴向捕获势，将纳米颗粒运动限制在一维从而给出明确光学势定义。计算光学力和势均。如图13a、b观察，随纳米球移向二聚体间隙，光学力和光学势深度明显增加。在(0,0,20 nm)获得3.8 k_B T光学势深度，足以克服水中布朗运动实现稳定光学捕获。

最后，基于不同网格尺寸将未知数从100万变至320万测试FETI-DP计算成本。实践中，不同网格密度测试常需满足不同精度要求。如此大规模复杂问题给常规数值方法带来巨大挑战。然而FETI-DP可轻松处理。采用32个处理器进行FETI-DP模拟，每个处理器处理一个子域。表3报告FETI-DP计算成本。可见FETI-DP展现高模拟效率和低内存需求。

结论

提出结合低秩稀疏化的FETI-DP方法用于金属纳米颗粒光学捕获预测和分析。该方法提供完全解耦子域问题，将大规模复杂问题转化为一系列小规模简单问题。适合并行计算，可显著提高数值模拟效率。例子证明该方法展现优异大规模计算性能，适合纳米尺度光学捕获快速精准模拟。

数据与材料可用性

本研究期间生成或分析的所有数据均包含在本文中。

缩略语

ABC

吸收边界条件

DOF

自由度

FDTD

有限差分时域

FEM

有限元方法

FETI-DP

对偶-原始有限元撕裂互联

MST

麦克斯韦应力张量

MVP

矩阵-向量乘积