电脑里的光：FDTD 如何把麦克斯韦方程算成光学仿真
引言：仿真软件里的”完美响应”是怎么来的？
最近几篇超表面文章反复出现同一个词：仿真。单元库扫描要预计算几百个单元的响应，伴随方法要在每一步迭代里跑正向和反向两次求解，加工误差分析要统计不同尺寸下的性能波动——它们的前提都是同一个东西：一个能算的电磁仿真器。
仿真软件对很多读者来说是个黑箱：填好结构、点运行、出来一条透射率曲线。但”能算”本身是 1966 年之后才逐渐解决的问题。那一年，Kane Yee 在 IEEE 期刊上发表了后来被称为 FDTD（Finite-Difference Time-Domain，有限时域差分）的方法，把连续空间中的麦克斯韦方程搬到了矩形网格上，用最朴素的时间步进方式直接”播放”电磁场的演化（S1、S8）。此后半个多世纪，FDTD 相关出版物数量一路增长（S1）。
图 1：FDTD 相关出版物数量逐年增长，反映该方法在计算电磁学中的普及。来源：Wikimedia Commons（公有领域）。
这篇文章不教你操作具体软件，而是回答一个更基本的问题：FDTD 到底是怎么把麦克斯韦方程”算”出来的？ 理解了它的网格、精度和稳定性，你才能看懂超表面仿真里那些最常见的建议——“网格要加密””时间步长要收敛””边界要用 PML”——到底是为什么。
核心思想：把连续的场变成网格上的积木
真空中电磁场的演化由两个旋度方程支配：
$$-\mu\frac{\partial \mathbf{H}}{\partial t} = \nabla\times\mathbf{E}, \qquad \varepsilon\frac{\partial \mathbf{E}}{\partial t} = \nabla\times\mathbf{H}$$
左边是时间导数，右边是场的空间变化（旋度）。如果能把空间和时间都切成小步长，用有限差分近似导数，方程就变成了”已知此刻的场，推下一刻的场”的递推式——这就是 FDTD 的基本思路（S5）。
为了让差分公式保持对称和稳定，Yee 的巧妙之处在于把电场和磁场分量交错排列：电场的每个分量周围环绕着磁场的分量，反之亦然，就像两块互相咬合的积木（S8、S1）。
图 2：Yee 网格：电场与磁场分量在空间上错开半个网格单元交错排列，使旋度项可以用中心差分近似。来源：FDominec，Wikimedia Commons（CC BY-SA 4.0）。
以一维的横电磁波为例（场只有 $E_x$ 和 $H_y$，沿 $z$ 传播），Yee 网格上的更新公式可以写成：
$$H_y^{n+1/2}(k+1/2) = H_y^{n-1/2}(k+1/2) - \frac{\Delta t}{\mu\Delta z}\left[E_x^{n}(k+1) - E_x^{n}(k)\right]$$
$$E_x^{n+1}(k) = E_x^{n}(k) - \frac{\Delta t}{\varepsilon\Delta z}\left[H_y^{n+1/2}(k+1/2) - H_y^{n-1/2}(k-1/2)\right]$$
其中上标 $n$ 是时间步序号，括号里的 $k$ 是空间网格序号。电场和磁场在时间上也错开半步，形成”蛙跳”式（leapfrog）交替更新：先用旧电场推新磁场，再用新磁场推新电场（S1、S4）。
图 3：Yee 单元局部示意：黑点代表电场分量，白点代表磁场分量，各分量错开半格排布。来源：Zohar0729，Wikimedia Commons（CC BY-SA 4.0）。
这套递推有两个漂亮的性质。第一，它不需要解线性方程组：每一步只用到邻近网格点的旧值，显式推进即可（S5）。第二，它对几何几乎没有限制：介质、金属、非线性材料都可以通过逐点的介电常数 $\varepsilon$ 和磁导率 $\mu$ 放进网格里，这就是为什么 FDTD 能处理超表面里任意形状的柱体（S2、S3）。
精度的代价：网格里的光，速度不再严格等于 c
离散化不是免费的。在连续空间中，平面波的频率 $\omega$ 与波矢 $k$ 满足 $\omega = ck$；到了网格上，有限差分会把这条色散关系悄悄改掉。一维 FDTD 的数值色散关系是：
$$\left[\frac{1}{c\Delta t}\sin\left(\frac{\omega\Delta t}{2}\right)\right]^2 = \left[\frac{1}{\Delta z}\sin\left(\frac{\tilde{k}\Delta z}{2}\right)\right]^2$$
当 $\Delta t$ 和 $\Delta z$ 都趋于零时，$\sin(x)\approx x$，上式退化为 $\omega = c\tilde{k}$，与物理一致；但网格有限时，数值波矢 $\tilde{k}$ 与频率之间不再是线性关系，数值相速度会依赖频率和传播方向——这就是数值色散（S4）。后果很具体：一个宽谱脉冲在网格里传播一段距离后会被”拉宽”、峰值位置偏移，仿真的谐振峰也可能偏离真实位置。
控制误差的办法是加密网格：每波长采样点数越多，色散误差越小。常见经验范围是每波长 10–20 个网格点起步，但正确的做法是把网格加密一倍做收敛性测试——如果结果变化明显，说明网格还不够细（S3、S4）。
代价也很直接：三维 FDTD 的网格点数和时间步数都随分辨率上升，分辨率翻倍意味着计算量大幅增长（S2）。这就是为什么”仿真精度”永远是一场与算力的交易。
稳定的红线：Courant 条件
显式时间步进有一个致命弱点：时间步长不能随意取。如果 $\Delta t$ 太大，数值解会指数发散。保证稳定的条件叫 Courant 条件（CFL 条件），三维均匀网格下写作：
$$\Delta t \le \frac{1}{c\sqrt{\frac{1}{\Delta x^2}+\frac{1}{\Delta y^2}+\frac{1}{\Delta z^2}}}$$
一维时退化为 $\Delta t \le \Delta z/c$（S1、S4）。物理直觉是：一个时间步内，信息最多只能跨过一个网格单元，所以时间步长必须小于”光走过一个网格所需的时间”。这也是为什么 FDTD 里网格加密后，时间步长必须跟着缩小——两者是绑在一起的（S4）。
边界与激励：PML 吸收边界和”制造”入射波
仿真域必须是有限的，但光会一直往外跑。粗暴截断边界会产生反射、污染结果。1994 年 Berenger 提出的完美匹配层（PML）解决了这个问题：在计算域外围铺一层特殊介质，让入射波”走进去就出不来”，几乎不产生反射（S1、S4）。今天的 FDTD 软件（包括 Meep）普遍用 PML 包裹计算域（S2、S3）。
入射波怎么放进去？一种常用方案是总场/散射场（TFSF）技术：把计算域分成总场区和散射场区，在分界面上注入入射波，这样散射体之外记录的是散射场，方便提取吸收截面、散射截面这类物理量（S1）。
图 4：FDTD 计算域示意图：中部为总场/散射场（TFSF）源注入区域，外围深色区域为 PML 吸收边界。来源：Dominator9000 等，Wikimedia Commons（CC BY-SA 4.0）。
时域的红利：一个源，算出一整段光谱
FDTD 是时域方法，这个”缺点”（必须逐步推进）也带来一个独特红利：一次仿真可以覆盖宽频带。注入一个高斯脉冲而不是单频连续波，让它在结构里传播足够长时间，再用傅里叶变换把时域响应变成频谱，就能一次得到一段频率范围内的透射率、反射率曲线（S1、S2、S3）。相比每算一个频率都要重新求解的频域方法，宽带问题在 FDTD 里天然高效。
真实应用也是这样做的。Sur 等人 2026 年的工作就用 FDTD 全波仿真研究了超均匀无序纳米孔图案对钙钛矿太阳能电池的宽带吸收增强——这类结构既没有周期性可以利用，又要求宽谱响应，正是 FDTD 的舒适区（S9）。
什么时候不该用 FDTD：与 RCWA 的取舍
FDTD 并非万能。超表面设计里最常见的做法——预计算单元库——用的是另一种方法：严格耦合波分析（RCWA）。RCWA 把周期性结构中的电磁场展开为平面波/傅里叶级数，在频域求解，天然假设结构无限周期（S7）。对二维周期单元做光学仿真时，RCWA 通常比 FDTD 这类全波方法快约一个数量级，且精度足够用于光学设计（S6）。两种方法的取舍可以概括成一张表：
| 对比项 | FDTD | RCWA |
|---|---|---|
| 求解思路 | 时域显式步进 | 频域傅里叶展开 |
| 周期性假设 | 不需要，任意几何 | 需要（无限周期） |
| 宽谱问题 | 一次算完（傅里叶变换） | 逐频率求解 |
| 典型开销 | 随网格数与时间步数增长 | 随展开阶数增长，周期单元快 |
| 适合场景 | 任意几何、宽带、非线性、无序结构 | 周期单元库、光栅、严格周期超表面 |
表 1：FDTD 与 RCWA 的典型取舍。对比依据：S5（时域/频域方法综述）、S6（RCWA 与 FDTD 速度对比）、S7（RCWA 原理）。
所以工程上的常见分工是：用 RCWA 快速扫描周期单元、建立设计库，再用 FDTD 对最终的非周期超表面器件做全波验证。上一篇文章讲的伴随方法也依赖正向求解器——它要求正向和反向各跑一次，FDTD 的宽带特性在这里同样有用。
结论与展望
FDTD 用四个简单的设计选择换来了极大的通用性：交错网格（Yee 网格）、中心差分、显式时间步进、PML 截断边界。理解这四个选择，就能看懂仿真软件里几乎所有关键参数：
- 网格分辨率决定数值色散误差，用收敛性测试而不是经验值确认；
- 时间步长受 Courant 条件约束，网格加密后必须同步缩小；
- PML 厚度与距离决定边界反射水平，太薄或太近都会污染结果；
- 源与监视器决定你提取的是什么物理量，宽带问题优先用时域脉冲。
它也有明确的使用边界：三维仿真对内存和时间的需求随分辨率急剧增长，不适合超大面积或超精细问题；这类问题要么用频域方法（RCWA/FEM）缩小问题规模，要么用近似方法，这也是 FDTD 与频域方法互补共存的场景。下次再看到”仿真结果很漂亮，但实验对不上”，除了工艺误差（参见《超表面为什么怕加工误差》），也可以回头检查一下：网格够细吗？时间步稳定吗？PML 够远吗？（S1、S3、S4）
参考资料
- S1: Finite-difference time-domain method（Wikipedia）— https://en.wikipedia.org/wiki/Finite-difference_time-domain_method
- S2: Oskooi et al., Meep: A flexible free-software package for electromagnetic simulations by the FDTD method（Computer Physics Communications 181, 2010；MIT DSpace 开放副本）— https://dspace.mit.edu/bitstream/handle/1721.1/60946/Johnson_MEEP%20A.pdf
- S3: Meep 官方文档 Introduction — https://meep.readthedocs.io/en/latest/Introduction/
- S4: J. B. Schneider, Understanding the FDTD Method（华盛顿州立大学公开教材）— https://eecs.wsu.edu/~schneidj/ufdtd/
- S5: H. De Raedt et al., Numerical methods for solving the time-dependent Maxwell equations（arXiv:physics/0210035）— https://arxiv.org/abs/physics/0210035
- S6: B. Gao, H. Gersen, S. Hanna, On the suitability of rigorous coupled-wave analysis for fast optical force simulations（arXiv:2306.17016）— https://arxiv.org/abs/2306.17016
- S7: Rigorous coupled-wave analysis（Wikipedia）— https://en.wikipedia.org/wiki/Rigorous_coupled-wave_analysis
- S8: K. S. Yee, Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media（IEEE Trans. Antennas Propag. 14, 1966）— https://ieeexplore.ieee.org/document/1138693/
- S9: A. Sur, K. Nath, A. Zubair, Nature-Inspired Hyperuniform Nanohole Patterning for Robust Broadband Absorption Enhancement in Perovskite Solar Cells（arXiv:2604.11264）— https://arxiv.org/abs/2604.11264