在几何光学中，我们追踪光线经过透镜后向哪里传播；在波动光学中，我们追踪振幅和相位如何衍射与干涉。可当一片透镜被替换成密密麻麻的亚波长纳米柱时，光还会在柱内形成模式、在界面上反射、在相邻结构之间耦合。此时，仅给每个位置指定一个理想相位，往往不足以回答“真实器件究竟能聚焦多少光”。

时域有限差分法（finite-difference time-domain，FDTD）直接离散 Maxwell 方程，在空间网格上逐步推进电场和磁场。 它既能描述传播，也能描述反射、干涉、衍射和材料响应。漂亮的场图只是其中一种输出；更有价值的是在明确误差范围内获得透射率、相位、焦斑、效率和光谱。

FDTD 与介质超构透镜的概念封面

封面｜由生图工具生成的概念插画。网格、纳米柱和会聚波前仅用于引入主题，不代表精确结构、实际入射配置或仿真场分布。下文图 1–5 由代码按给定坐标、公式或求解器结果绘制。

本文面向具备大学电磁学基础的读者。第 1–5 节讲算法，第 6–9 节讲材料、边界和输出，第 10–13 节结合超构透镜建立验证流程，第 14 节给出可运行代码。文中的器件参数是教学假设，没有把示意图当作已经完成的三维超构透镜仿真。

1. FDTD 到底在求什么？

1.1 从旋度方程出发

采用 SI 单位，记电场、磁场、电位移和磁感应强度为 。没有磁流源时，

$$ \frac{\partial\mathbf B}{\partial t}=-\nabla\times\mathbf E,\qquad \frac{\partial\mathbf D}{\partial t}=\nabla\times\mathbf H-\mathbf J. $$

先考虑各向同性、线性、非色散介质：

$$ \mathbf D=\varepsilon\mathbf E,\qquad \mathbf B=\mu\mathbf H,\qquad \mathbf J=\sigma\mathbf E+\mathbf J_{\mathrm{src}}. $$

这里 、，电导率 表示欧姆损耗， 表示外加激励。于是，

$$ \mu\frac{\partial\mathbf H}{\partial t}=-\nabla\times\mathbf E,\qquad \varepsilon\frac{\partial\mathbf E}{\partial t} =\nabla\times\mathbf H-\sigma\mathbf E-\mathbf J_{\mathrm{src}}. $$

这组方程天然适合“知道现在，计算下一刻”：已知电场的空间变化，就能更新磁场；已知磁场的空间变化，就能更新电场。FDTD 所做的，是把导数改写成附近几个采样值之间的差分。

全文以标准显式、二阶 Yee 格式为主。高阶格式、隐式格式和非结构网格方法不在本文推导范围内。

1.2 散度方程去哪了？

Maxwell 方程还包含

$$ \nabla\cdot\mathbf D=\rho,\qquad \nabla\cdot\mathbf B=0. $$

它们是必须满足的约束。连续数学中，“旋度的散度为零”意味着：若初始场满足约束，而且电流与电荷满足连续性方程，后续旋度推进就会保持约束。

Yee 网格具有相应的离散结构。不过，这并不意味着任意初始化、任意电流注入和任意插值都能自动正确。带电粒子耦合、复杂材料界面或自定义源仍需要检查离散电荷守恒。

1.3 “全波”意味着什么？

全波表示求解电磁场的波动方程组，而不是仅追踪光线或套用近轴传播公式。它并不意味着没有近似：材料模型、几何离散、有限计算区域、有限时间和源的理想化仍然限制结果。

对一个已定义的模型，FDTD 回答的是“这个离散模型会如何响应”，然后通过收敛与基准验证，逐步说明它与目标物理问题有多接近。

2. Yee 网格：为什么要把电场和磁场错开放？

2.1 中心差分的几何意义

将一个光滑函数在中心两侧展开相减，可以得到

$$ \left.\frac{\partial F}{\partial x}\right|_{x} =\frac{F(x+\Delta x/2)-F(x-\Delta x/2)}{\Delta x} +O(\Delta x^2). $$

导数位于两个采样点的中间。如果电场存储在一组位置，更新磁场时需要的电场差分自然位于两点之间。因此，把磁场放到半格偏移的位置，能够让空间差分和场分量的位置对应起来。

在三维中，一种常用坐标约定如下。表中坐标分别以 为单位：

场分量 存储坐标 时间层

电场可以看作位于单元边中点，磁场位于相应面的中心。整体平移半格，或交换整数与半整数时间层，也可以得到等价约定；关键是差分必须在正确的位置相遇。

2.2 空间交错，时间也交错

推进顺序为

$$ \left(\mathbf E^n,\mathbf H^{n-1/2}\right) \longrightarrow\mathbf H^{n+1/2} \longrightarrow\mathbf E^{n+1}. $$

这称为蛙跳（leapfrog）更新。电场和磁场不是在同一时间点同时更新，而是交替跨过对方。

二维 Yee 网格与时间蛙跳推进

图 1｜左图采用二维 约定：圆点为出平面的电场，橙色和蓝色短线分别表示沿 x、y 方向的磁场。虚线框提示围绕中心 的磁场差分位置。右图给出更新时间顺序，箭头表示依赖关系。

空间与时间的中心差分在光滑区域具有二阶精度，即误差的主导项通常随 和 缩小。界面、几何台阶和局部奇异场可能改变实际观测到的收敛阶，不能机械地把整个复杂器件都视为二阶收敛。

Yee 的原始工作发表于 1966 年；交错离散的基础推导也可参阅 Schneider《Understanding the FDTD Method》第 3 章。

3. 先从一维推导出完整更新方程

令波沿 z 方向传播，仅保留 和 ，并暂时忽略损耗和源：

$$ \frac{\partial H_y}{\partial t} =-\frac1\mu\frac{\partial E_x}{\partial z},\qquad \frac{\partial E_x}{\partial t} =-\frac1\varepsilon\frac{\partial H_y}{\partial z}. $$

这里两个负号来自所选场方向和传播坐标。改用 时符号会不同，不能只凭记忆抄写。

设

$$ E_i^n=E_x(i\Delta z,n\Delta t),\qquad H_{i+1/2}^{n+1/2} =H_y((i+\tfrac12)\Delta z,(n+\tfrac12)\Delta t). $$

在磁场位置做中心差分：

$$ \frac{H_{i+1/2}^{n+1/2}-H_{i+1/2}^{n-1/2}}{\Delta t} =-\frac1{\mu_{i+1/2}}\frac{E_{i+1}^n-E_i^n}{\Delta z}. $$

移项后得到

$$ \boxed{ H_{i+1/2}^{n+1/2} =H_{i+1/2}^{n-1/2} -\frac{\Delta t}{\mu_{i+1/2}\Delta z}(E_{i+1}^n-E_i^n) }. $$

随后在电场位置更新：

$$ \boxed{ E_i^{n+1} =E_i^n-\frac{\Delta t}{\varepsilon_i\Delta z} \left(H_{i+1/2}^{n+1/2}-H_{i-1/2}^{n+1/2}\right) }. $$

这就是一个最小 FDTD 求解器的核心。每个场值只依赖相邻位置，既不需要在每一步求解全局线性方程组，也不需要提前假设波的形状。不过，边界和初始条件必须另行指定。

若取归一化单位 ，定义 ，核心循环只有两行：

1
2
h -= S * (np.roll(e, -1) - e)
e -= S * (h - np.roll(h, 1))

这里 np.roll 把右端接回左端，实现的是周期边界。它并不吸收波；把这段代码用于自由空间脉冲时，波会从另一端重新进入。第 14 节会给它配上真正适合周期边界的解析基准。

4. 从一维扩展到二维与三维

若结构和场沿 z 不变，可取 三个非零分量，本文将其称为 TM。不同软件对 TE/TM 的参照轴可能不同，交流模型时应同时给出非零分量。

无损情况下，

$$ \frac{\partial H_x}{\partial t}=-\frac1\mu\frac{\partial E_z}{\partial y},\qquad \frac{\partial H_y}{\partial t}=+\frac1\mu\frac{\partial E_z}{\partial x}, $$
$$ \frac{\partial E_z}{\partial t} =\frac1\varepsilon\left( \frac{\partial H_y}{\partial x}-\frac{\partial H_x}{\partial y} \right). $$

图 1 的空间约定给出

$$ H_{x,i,j+1/2}^{n+1/2} =H_{x,i,j+1/2}^{n-1/2} -\frac{\Delta t}{\mu\Delta y}(E_{z,i,j+1}^n-E_{z,i,j}^n), $$
$$ H_{y,i+1/2,j}^{n+1/2} =H_{y,i+1/2,j}^{n-1/2} +\frac{\Delta t}{\mu\Delta x}(E_{z,i+1,j}^n-E_{z,i,j}^n), $$
$$ \begin{aligned} E_{z,i,j}^{n+1}=E_{z,i,j}^n+\frac{\Delta t}{\varepsilon} \bigg[& \frac{H_{y,i+1/2,j}^{n+1/2}-H_{y,i-1/2,j}^{n+1/2}}{\Delta x}\\ &-\frac{H_{x,i,j+1/2}^{n+1/2}-H_{x,i,j-1/2}^{n+1/2}}{\Delta y} \bigg]. \end{aligned} $$

为简化书写，这里没有在 上标出空间下标；实际程序必须读取各场分量所在位置的材料系数。

三维算法沿同一逻辑更新六个分量。例如，

$$ \frac{\partial E_x}{\partial t} =\frac1\varepsilon\left( \frac{\partial H_z}{\partial y}-\frac{\partial H_y}{\partial z} \right). $$

其余分量按坐标循环置换得到。计算量增加主要来自网格体积和分量数量，而不是某个新的数学原理。

二维模型的适用条件是物理结构在省略方向上不变，或采用另一个经过说明的降维模型。 把一排纳米柱的截面做成二维模型，通常等价于沿省略方向无限延伸的条带，不能直接代表由有限圆柱构成的三维超构透镜。

5. CFL 条件保证稳定，网格收敛决定精度

5.1 时间步不能随意增大

对于均匀正交网格上的标准显式 Yee 格式、正的非色散介电常数和磁导率，常用 CFL 上界为

$$ \Delta t\le \frac{1}{ v_{\max}\sqrt{\Delta x^{-2}+\Delta y^{-2}+\Delta z^{-2}} }, $$

其中 是模型中最大的传播速度。二维和一维分别去掉不存在的空间项。三维等距网格、最快速度为真空光速 c 时，

$$ \frac{c\Delta t}{\Delta}\le\frac1{\sqrt3}. $$

实际常在理论上界内留一定余量。该公式不是所有材料模型的通用稳定性保证：色散辅助方程、增益、特殊网格或边界实现可能再施加约束。Meep 的算法说明也强调时间步与空间分辨率共同限制稳定性。

例如，三维均匀网格 ，取真空 CFL 上限的 0.95 倍，则

$$ \Delta t=0.95\frac{20\,\mathrm{nm}}{c\sqrt3} \approx36.59\,\mathrm{as}. $$

模拟 200 fs 至少约需 5466 步。这里 。这是给定网格的算术估计，不表示 200 fs 对所有超构透镜都足够。

5.2 为什么“没有发散”仍然可能算错？

在均匀无损介质中，将离散平面波代入更新式，可得到数值色散关系：

$$ \left[\frac{\sin(\omega\Delta t/2)}{v\Delta t}\right]^2 = \left[\frac{\sin(k_x\Delta x/2)}{\Delta x}\right]^2+ \left[\frac{\sin(k_y\Delta y/2)}{\Delta y}\right]^2+ \left[\frac{\sin(k_z\Delta z/2)}{\Delta z}\right]^2. $$

网格无限细时，，它才还原成连续关系 。有限网格下，数值相速度通常随波长和传播方向变化：即使物理材料没有色散，算法也可能产生数值色散。

网格分辨率与传播方向导致的数值色散

图 2｜由离散色散关系直接求得。左图是一维真空情形，横轴按连续介质中的物理波长计数；右图是二维方格网格。二维曲线显示了网格引入的方向依赖。左图的 只用于一维，不能直接搬到二维或三维。

对聚焦器件，相位差特别敏感。即使局部相速度偏差很小，传播距离较长时也会积累相位误差，改变焦点、焦斑和干涉结构。相关推导可参阅 Schneider 第 7 章。

5.3 网格需要同时看波长和几何

一个粗略起点是按高折射率区域内的最短波长估计：

$$ \lambda_{\mathrm{mat}}\approx \frac{\lambda_0}{\operatorname{Re}n},\qquad \Delta\lesssim\frac{\lambda_{\mathrm{mat,min}}}{N_\lambda}. $$

对于弱损耗介质，可以先尝试每个介质波长 10–20 个网格，再做收敛研究。这只是初值经验。若最小缝隙 20 nm，而网格也是 20 nm，即使“每波长网格数达标”，几何仍然几乎没有被解析。

金属皮肤深度、纳米间隙、尖角附近的强场和高 Q 共振也可能提出更细要求。圆柱被台阶化后，等效直径与真实直径存在偏差；共形网格或亚像素材料平均可以改善误差，但其适用性应结合材料与求解器实现检查。

6. 材料如何进入更新时间循环？

6.1 电导损耗：对耗散项做时间中心平均

将电场方程写为

$$ \varepsilon\frac{E^{n+1}-E^n}{\Delta t} =C^{n+1/2} -\sigma\frac{E^{n+1}+E^n}{2}, $$

其中 表示该分量的磁场旋度减去外加电流。整理可得

$$ E^{n+1}=C_aE^n+C_bC^{n+1/2}, $$
$$ C_a=\frac{1-\sigma\Delta t/(2\varepsilon)} {1+\sigma\Delta t/(2\varepsilon)},\qquad C_b=\frac{\Delta t/\varepsilon} {1+\sigma\Delta t/(2\varepsilon)}. $$

损耗通过系数进入更新，仍可以预计算系数并显式推进。不过，常数电导率对应特定频率依赖，不能任意代替整个宽波段的实测吸收。

6.2 色散材料：需要记住材料的历史

真实材料的极化不会总是瞬时跟随电场。可写为

$$ \mathbf D=\varepsilon_0\varepsilon_\infty\mathbf E+ \sum_m\mathbf P_m. $$

一种 Lorentz 极化模型是

$$ \frac{\mathrm d^2\mathbf P_m}{\mathrm dt^2} +\gamma_m\frac{\mathrm d\mathbf P_m}{\mathrm dt} +\omega_m^2\mathbf P_m =\varepsilon_0F_m\mathbf E. $$

这里 具有角频率平方的量纲。全文频域采用 约定，对应

$$ \varepsilon_r(\omega)=\varepsilon_\infty+ \sum_m\frac{F_m}{\omega_m^2-\omega^2-i\gamma_m\omega}. $$

设置 且保留适当的 ，可得到 Drude 型项。一个中心差分的极化递推为

$$ \mathbf P_m^{n+1}= \frac{ (2-\omega_m^2\Delta t^2)\mathbf P_m^n -(1-\gamma_m\Delta t/2)\mathbf P_m^{n-1} +\varepsilon_0F_m\Delta t^2\mathbf E^n }{ 1+\gamma_m\Delta t/2 }. $$

这称为辅助微分方程（ADE）思路。推进极化和 后，再由本构关系求出 。程序需要额外保存极化历史，所以色散材料会增加内存和计算量。具体耦合格式的稳定范围必须另外验证。

6.3 从实测 n、k 到求解器材料

对非磁材料，令复折射率为 ，则

$$ \varepsilon_r=(n+i\kappa)^2 =(n^2-\kappa^2)+i\,2n\kappa. $$

在本文时间约定下，被动吸收材料通常具有正的 。将实验表格用于时域仿真时，需要拟合为求解器支持的因果色散模型，而不是在一个宽带运行中把不同波长的 n 任意填进一个常数材料。

实际应检查三件事：拟合频带是否覆盖源与监视器；实部和虚部是否都拟合合理；结构对象最终绑定的材料、各向异性轴和求解器读回曲线是否正确。材料名字正确并不足以证明物理参数正确。Meep 材料文档给出了用 Lorentz 项拟合 n、k 数据及其稳定性限制。

7. 边界条件决定你究竟模拟了哪个物理问题

7.1 PML：让离开的波进入吸收层

有限网格的最外层若直接设为零，通常会产生反射。开放辐射问题常在计算区外侧加入完全匹配层（PML）：它通过特殊的场方程或坐标伸缩，使出射波进入后衰减。

连续理论中的理想匹配不等于离散实现绝对无反射。有限厚度、网格误差和大角度入射都可能留下残余反射。应分别检查 PML 厚度与结构到 PML 的距离，并把物理监视器放在 PML 之外；靠得过近的吸收层还可能干扰倏逝场或共振。

基底等背景材料通常需要连续延伸到吸收边界，以免人为制造一个材料终止面；这与“让 PML 避开器件局域近场”并不矛盾。具体介质和传播方向还可能存在 PML 失效情形，参阅 Meep PML 说明。

7.2 周期与 Bloch 边界：模拟无限重复

普通周期条件为

$$ \mathbf F(\mathbf r+\mathbf a)=\mathbf F(\mathbf r). $$

斜入射或一般 Bloch 问题则满足

$$ \mathbf F(\mathbf r+\mathbf a) =e^{i\mathbf k_\parallel\cdot\mathbf a}\mathbf F(\mathbf r). $$

一个周期单元的透射响应描述“该单元在无限重复阵列中的表现”。有限口径超构透镜的边缘、非均匀邻居和整体散射不会自动包含在里面。给整片透镜的侧面加周期条件，得到的也会是透镜阵列，而不是孤立透镜。

宽带斜入射还需特别小心：固定切向波矢与固定物理入射角并非等价。软件如何实现角度随频率的关系，必须在源和边界配置中确认。

7.3 对称边界：结构对称还不够

只有几何、材料、激励及场分量的奇偶性同时满足要求，才能用对称面缩小计算区域。电场是极向量，磁场是轴向量，镜像变换下各分量的奇偶性并不相同。

圆形口径不自动意味着任意偏振、斜入射和纳米柱排布都可采用四分之一计算域。用小规模完整模型比较一次，通常比仅看几何图案更可靠。

8. 光源与监视器：从时间波形得到频谱

8.1 宽带脉冲和单频激励

高斯包络脉冲可以写成

$$ s(t)=\exp\left[-\frac{(t-t_0)^2}{2\tau^2}\right] \cos[\omega_0(t-t_0)]. $$

一次线性时不变系统的时域运行，可以在有效激励频带内提取多个频率响应。不过，源频谱尾部很弱时，除以源谱会放大噪声；宽带源不意味着所有频率都同样可靠。

单频连续波激励适合少量目标频点，但要等待瞬态充分衰减。高 Q 共振的衰减时间可能远大于普通传播时间。停止条件应结合目标频率、监视位置和结果变化，而不是只看源已经关闭。

源实现也有区别：硬源直接覆盖场值，可能干扰返回波；软源或电流源将激励加入更新方程；总场/散射场（TFSF）方法通过边界修正分离入射场与散射场。选择哪种方式，应与反射或散射的提取方法配套。

8.2 在线傅里叶积累

采用 约定，频域场可由

$$ \widetilde{\mathbf E}(\omega)\approx \Delta t\sum_n w_n\mathbf E^n e^{+i\omega n\Delta t} $$

提取。 为时间窗；不必保存所有时刻，只需维护所需频点的复数累加器。若只有少数频点，这比保存整段四维场数据节省得多。上式给出傅里叶累积量，换成单频稳态相量或频谱功率时，需要采用相应的时间窗与源谱归一化；下文的功率公式使用已按一致约定归一化的场。

有限观测时长 T 带来的直接傅里叶频率分辨尺度约为

$$ \Delta f\sim\frac1T. $$

窗函数、谱线形状和参数拟合会影响实际分辨能力；简单增加 FFT 补零只能让曲线插值更密，不能创造新的观测信息。

磁场的真实采样时刻是 ，做傅里叶积累时应使用对应时间相位。计算功率前还需要把交错的场插值到共同空间位置。直接拿不同位置、不同时间层的 E 和 H 相乘，可能引入系统偏差。

8.3 透射相位与功率不是同一个量

在相同参考面和同一输出模式下，可以定义归一化复振幅

$$ t(\omega)=\frac{a_{\mathrm{out}}(\omega)} {a_{\mathrm{ref}}(\omega)},\qquad \phi(\omega)=\arg t(\omega). $$

对于周期单元，若关心零级透射，应提取相应衍射级次或模式振幅。结构附近某一点的电场同时含有近场和多个级次，不宜直接当作单元的远场透射系数。

一般情况下，功率应由时间平均 Poynting 矢量积分：

$$ P(\omega)=\frac12\operatorname{Re} \int_A \left(\widetilde{\mathbf E}\times\widetilde{\mathbf H}^{*}\right) \cdot\hat{\mathbf n}\,\mathrm dA. $$

只有在相应归一化条件满足时， 才直接等于功率透射率。例如，法向入射、两侧无损非磁介质、t 按电场振幅定义时，

$$ T=\frac{n_{\mathrm{out}}}{n_{\mathrm{in}}}|t|^2. $$

不同软件的模式振幅可能已经按功率归一化，使用时应确认其定义。

9. 先算算：一次三维 FDTD 要付出多少代价？

设均匀网格数为 ，固定物理时间内的步数为 ，则基本代价近似为

$$ \text{内存}\sim O(N_xN_yN_z),\qquad \text{运算量}\sim O(N_xN_yN_zN_t). $$

在固定物理区域和总仿真时间下，将三维网格间距全部减半，空间单元数约乘 8，CFL 时间步约减半，时间步数约乘 2，因此运算量约乘 16。实际耗时还受内存带宽、并行通信、GPU 数据搬运和监视器影响。

例如， 的均匀网格区域，若三个方向都采用 20 nm 间距，则有

$$ 600\times600\times400=1.44\times10^8 $$

个单元。若仅为六个场分量各存一个 32 位浮点数，就需要约 3.456 GB，约合 3.22 GiB。这只是场数组的简化下限估计：材料系数、色散变量、PML、复数频域监视器和运行开销还会继续增加内存。

因此，设计超构透镜时常将“单元电磁响应”“有限区域验证”“器件外部传播”分层计算。能用受验证的传播算法完成的均匀空间，不必全部填满极细 FDTD 网格。

10. 超构透镜设计：把目标波前变成纳米结构

10.1 理想聚焦相位从哪里来？

设平面波从下方照亮位于 的薄相位面，目标焦点为 。输出侧为均匀介质，波数

$$ k_{\mathrm{out}}(\omega)= \frac{n_{\mathrm{out}}(\omega)\omega}{c}. $$

半径 r 处到焦点的距离为 。为了让各位置到达焦点时同相，透镜提供的相位必须补偿额外传播相位：

$$ \boxed{ \phi(r,\omega)= -k_{\mathrm{out}}(\omega) \left(\sqrt{f^2+r^2}-f\right)+\phi_0(\omega) }. $$

这里沿 +z 传播的相量带有 ，因此补偿项取负号。 为不影响单频聚焦的整体相位。若入射本身带有非均匀相位，还要减去入射波前在该平面上的相位分布。

超构透镜的目标相位与等光程示意

图 3｜教学参数：真空波长 532 nm、空气中直径 10 μm、焦距 10 μm。左图由理想相位公式计算；右图仅画出通向目标焦点的几何路径，没有计算实际场强或聚焦效率。

此例边缘比中心多传播约 ，所以需要约 rad 的相对补偿，约为 个周期。实际选型通常使用包裹到一个 区间内的相位。

空气中的数值孔径为

$$ \mathrm{NA}=\frac{D/2}{\sqrt{f^2+(D/2)^2}} \approx0.4472. $$

该双曲相位公式本身没有使用近轴展开；但将其实现为局域薄相位面，以及后续采用标量传播，仍各自包含建模假设。

10.2 单元库：求的是复数响应

以介质纳米柱为例，固定周期 p、高度 h、材料和基底，扫描直径 d 或其他几何参数。每个周期单元计算

$$ t(d,\lambda)=|t(d,\lambda)|e^{i\phi(d,\lambda)}. $$

单元库应同时保存相位、透射功率、偏振转换和参数依赖。相位覆盖 到 只是基本条件；若某一段相位只能由低透射的强共振单元提供，最终器件可能损失很多能量。

两段接近 0 和 的相位实际上很接近，因此匹配时可使用圆周相位误差

$$ \delta\phi= \arg\!\left[e^{i(\phi_{\mathrm{cell}}-\phi_{\mathrm{target}})}\right]. $$

一种示例选型目标为

$$ d_*=\arg\min_d \left[ w_\phi\,\delta\phi^2+ w_T(1-T)^2+ w_{\mathrm{fab}}C_{\mathrm{fab}}(d) \right]. $$

权重需要结合设计任务归一化和选择；制造项可表示最小间隙、尺寸灵敏度等约束。这是设计策略示例，并非一套保证最优的通用公式。

对双折射单元或偏振转换器件，应使用 Jones 矩阵或相应模式耦合矩阵。一个标量 t 无法描述完整偏振响应。

10.3 单元库的局域周期近似

周期仿真假设每根纳米柱的周围都是相同单元。真实超构透镜中的相邻直径则不断变化，因此把周期库直接拼成器件，依赖局域周期近似（LPA）。

这种近似可能在相位梯度大、相邻单元突变、强共振或高 NA 情况下变差。它忽略的部分包括非均匀邻域耦合及口径边缘影响。Ansys 的超构透镜工作流介绍区分了完整器件 FDTD、单元场拼接和光线传播等不同层级。

一个有用的采样提醒来自空气中的边缘相位梯度：

$$ \left|\frac{\mathrm d\phi}{\mathrm dr}\right|_{\mathrm{edge}} =k_0\mathrm{NA}. $$

若希望相邻采样点的相位变化小于 ，可得到粗略条件 。它只约束目标相位的采样，不能替代对基底衍射级次、实际散射和单元耦合的检查。

11. 从单元到整片透镜：边界、近场与传播

周期单元与有限口径透镜的计算区域对照

图 4｜横截面示意，不按实际尺度绘制。左侧侧边为周期边界，右侧侧边为 PML；两者求解的是不同物理问题。图中的监视面位于器件上方、PML 内缘以下。

11.1 一个可复现的工作流

  1. 明确设计输入。 记录真空波长或频带、焦距、口径、入射方向与偏振、基底、材料来源及制造约束。
  2. 生成并验证单元库。 固定参考面和归一化方法；对候选尺寸检查网格、频带和材料拟合。
  3. 生成一次器件排布。 将每个位置、尺寸、旋转角和材料写成可读取的几何清单，保存版本或哈希。
  4. 做局部与完整器件验证。 可先模拟非均匀相邻单元的小簇或代表性子区域，再对可承受尺寸的有限口径器件做三维全波计算。
  5. 提取复数近场并传播。 在合适的均匀输出区域记录场，通过经过验证的近远场变换或传播算法得到焦域。
  6. 计算预先定义的指标。 保存原始功率、归一化分母、焦平面位置和积分半径，并完成数值收敛。

Ansys 的小尺度超构透镜示例使用了直接器件仿真与重建场传播之间的比较。这类交叉检查能帮助界定简化模型可用的范围。

11.2 为什么可以不把焦点放进 FDTD 区域？

器件附近需要细网格来解析纳米结构；离开器件后的均匀空间主要发生传播。若已获得足够范围、足够采样的复数场，就可以把传播计算移交给其他方法。

以均匀、各向同性、无损输出介质中的角谱法为例，

$$ \widehat{\mathbf E}(k_x,k_y;z) =\widehat{\mathbf E}(k_x,k_y;z_0) e^{ik_z(z-z_0)},\qquad k_z=\sqrt{k_{\mathrm{out}}^2-k_x^2-k_y^2}. $$

选取向 +z 出射且倏逝分量不增长的平方根分支。矢量传播还应满足横向性条件，并正确处理纵向分量及 E、H 的关系。

这一步需要注意监视面截断、空间采样、传播时的 FFT 周期延拓以及倏逝分量。监视面过小会丢失大角度能量；简单补零可以减轻周期卷绕和改善插值，却不能补回已被截掉的场。

若监视面和目标面之间还有分层介质或其他结构，不能直接套用单一均匀介质传播核。软件的“远场”命令也可能针对特定背景假设，应核对其适用条件。

11.3 对称计算能证明到什么程度？

在严格满足对称条件时，半域或四分之一域计算能够恢复相应完整解。但它并不是独立的完整口径交叉验证。若要研究偏心、制造不对称、斜入射或偏振变化，必须重新判断对称性，而不是沿用旧配置。

12. 怎样定义“聚焦得好”？

12.1 焦距、焦斑与旁瓣

可通过轴上强度最大值定义某种轴上焦点，也可按平面包围功率或其他任务指标寻找最佳焦平面；这些定义在存在多峰、像差时可能不一致，应明确说明。

焦斑常报告 x、y 两方向的强度 FWHM。对于均匀照明圆孔径的标量衍射参考，中央 Airy 斑的 FWHM 全宽近似为

$$ \mathrm{FWHM}\approx0.514\frac{\lambda_0}{\mathrm{NA}}. $$

本例给出约 。这是理想参考尺度，不是本文器件的 FDTD 预测值。高 NA 矢量效应、非均匀振幅、偏振和离散单元都会改变实际焦斑。

只报很小的 FWHM 也不充分：能量可能转移到了强旁瓣。应同时报告二维分布、旁瓣和包围能量。

12.2 聚焦效率必须写清分母与积分区域

例如，可定义

$$ \eta_{\mathrm{focus}}(a,z_f) = \frac{ \displaystyle\int_{(x-x_f)^2+(y-y_f)^2\le a^2} \langle S_z(x,y,z_f)\rangle\,\mathrm dx\,\mathrm dy }{ P_{\mathrm{inc,aperture}} }. $$

这里 a 是明确指定的积分半径， 是所选焦点，分母为参考入射场穿过设计口径的功率。若采用“半径为若干倍 FWHM”，需说明 FWHM 是哪一个方向的全宽，以及它如何换算成圆盘半径。

另一个指标是

$$ \eta_{\mathrm{enc}} =\frac{P_{\mathrm{focus}}}{P_{\mathrm{transmitted}}}. $$

它衡量透射光中有多少集中到焦区，不能与以前述入射功率为分母的聚焦效率混用。若透射率和积分面也一致，两者满足 。

对高斯照明，还应记录束腰、口径截断和器件位置。仿真源覆盖了多大区域，并不自动定义效率分母。

12.3 能量检查要覆盖所有出路

对被动稳态问题，在一致参考与完整功率通道下，应满足反射、透射和吸收之和接近入射功率：

$$ R+T+A\approx1. $$

但有限口径结构可能有大量侧向散射。如果 T、R 只在有限平面上积分，上式可能遗漏侧向通量；不能把所有未收集到的能量都叫作材料吸收。更可靠的检查是使用闭合通量面或明确列出所有功率通道。瞬态计算还要考虑计算区内部的储能变化。

13. 多波长、消色差与数值收敛

13.1 扫描一片固定器件，才是在研究它的色差

对一片设计好的透镜进行波长扫描时，应保持几何排布不变，仅改变光源、监视频点以及由同一物理材料模型给出的频率响应。

如果每到一个波长就重新匹配一套纳米柱，那么得到的是多片分别优化的器件，无法由此说明原透镜的宽带性能。因此，应保存几何版本，同时读取各波长的焦点、FWHM、效率和旁瓣。

13.2 单频相位覆盖不等于消色差

保持固定焦距 f 时，需要同时满足

$$ \phi(r,\omega) =-k_{\mathrm{out}}(\omega)\Delta L(r)+\phi_0(\omega), \qquad \Delta L(r)=\sqrt{f^2+r^2}-f. $$

对频率求导，

$$ \frac{\partial\phi}{\partial\omega} =-\frac{\partial k_{\mathrm{out}}}{\partial\omega} \Delta L(r)+\frac{\mathrm d\phi_0}{\mathrm d\omega}. $$

这表明宽带同焦还要求不同位置具有合适的相对群时延，更宽的频带可能继续约束高阶色散。在本文相位约定下， 对应延迟；上式描述的是相对分布，公共相位项可改变整体延迟。

计算相位斜率时，需要在透射振幅可靠的频段展开相位，并扣除一致参考面的传播贡献。透射零点附近的相位导数可能非常敏感。仅在若干离散波长发现一个亮斑，不足以宣称连续频带消色差。

13.3 收敛研究应该检查哪些量？

不要只问“场图有没有变”，应预先指定可比较的指标。一个适合超构透镜的记录表如下：

检查项 保持不变 改变什么 比较什么
空间网格 同一物理几何和材料 网格间距、界面处理 焦距、相位、效率、FWHM
仿真时间 网格与频带 运行时长、终止阈值 频谱及焦域指标
PML 器件与源 厚度、参数 反射、谱线、效率
边界距离 PML 设置 空气缓冲区、侧向余量 焦域与功率通量
监视与传播 同一器件响应 监视面尺寸、采样、传播设置 积分功率、焦斑与旁瓣
对称性 相同物理激励 对称域与小型完整域 场的奇偶性和器件指标
材料 相同实验数据来源 拟合阶数与频带 n、k 误差和最终光学响应

对非零指标 Q，可记录

$$ \delta_Q= \frac{|Q_{\mathrm{fine}}-Q_{\mathrm{coarse}}|} {|Q_{\mathrm{fine}}|}. $$

Q 接近零时用相对误差会失去意义，应改用绝对容差；相位使用圆周差或一致展开后的差。比如将“效率相差小于 1%”作为目标时，要区分相对百分比与效率的百分点。

至少三档网格更有助于判断趋势。只有确认进入平滑的渐近收敛区，才适合做 Richardson 外推。界面重划分导致非单调变化时，不应为了得到一个外推数值而强行拟合。

14. 可运行示例：用解析驻波检验一维 FDTD

先验证基本更新时间循环，再叠加 PML、材料和复杂几何，通常更容易定位错误。这里采用长度为 1 的周期区域、归一化单位 ，选取解析解

$$ E_x(z,t)=\sin(kz)\cos(kt),\qquad H_y(z,t)=-\cos(kz)\sin(kt),\qquad k=4\pi. $$

电场从 开始，磁场必须初始化在 ，并位于半格偏移位置。这同时检查时间和空间交错是否被正确实现。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
import numpy as np

def periodic_mode(n, final_time=0.37, max_courant=0.5):
    dz = 1.0 / n
    steps = int(np.ceil(final_time / (max_courant * dz)))
    dt = final_time / steps       # 各网格严格比较同一个终止时刻
    S = dt / dz
    z = np.arange(n) * dz
    k = 4 * np.pi
    e = np.sin(k * z)             # E 位于 z_i，t = 0
    h = np.cos(k * (z + dz/2)) * np.sin(k * dt/2)
                                  # H 位于 z_(i+1/2)，t = -dt/2
    for _ in range(steps):
        h -= S * (np.roll(e, -1) - e)
        e -= S * (h - np.roll(h, 1))
    exact = np.sin(k * z) * np.cos(k * final_time)
    return np.sqrt(np.mean((e - exact)**2))

errors = [periodic_mode(n) for n in (40, 80, 160, 320)]
orders = np.log2(np.array(errors[:-1]) / np.array(errors[1:]))
print(errors)
print(orders)

本文配套脚本的实际运行结果为：

周期区域网格数 电场 RMS 误差 与上一档比较的收敛阶
40 —
80 2.0129
160 2.0140
320 2.0071

实际运行的一维 Yee 求解器与二阶收敛验证

图 5｜左图比较数值场和连续方程解析解；右图显示网格加密后误差近似按二阶下降。它验证的是无损、均匀、周期一维基准，不包含 PML、色散材料或三维超构透镜。

这组结果的意义是：代码实现了预期的交错更新，并在这个光滑基准上呈现二阶收敛。更完整的求解器仍应加入其他独立基准，例如介质平面界面的 Fresnel 反射、均匀层的传输矩阵解、散射体的解析解，以及 PML 残余反射测试。

下载 完整绘图与验证脚本 和 本次运行数值记录。安装 NumPy、Matplotlib 后，在脚本所在目录运行：

1
python figures-and-checks.py --output-dir figures

脚本会生成本文五张教学图和 validation.json。AI 概念封面独立于数值脚本生成；其生成提示词与来源记录也一并提供。

15. 从“会运行”走到“结果可信”

FDTD 的核心并不复杂：用中心差分替代导数，将电场与磁场错开放置，再交替更新时间层。难点在于让材料、几何、激励、边界和测量定义共同对应同一个物理问题。

在超构透镜中，目标相位给出设计方向，周期单元库提供局部响应，有限器件计算揭示耦合与边缘效应，传播与功率积分连接到焦斑和效率。每一步都有可以检验的假设，也有各自的误差来源。

读一份仿真报告时，可以沿着五个问题检查它：实际求解的几何和材料是什么；电磁边界代表什么环境；输出量怎样归一化；结果是否随网格、时间和边界收敛；简化模型是否经过独立基准或更完整模型的比较。能够回答这些问题，场图中的亮点才有清楚的物理含义。

参考资料与进一步阅读

  1. K. S. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Transactions on Antennas and Propagation, 14(3), 302–307 (1966). DOI: 10.1109/TAP.1966.1138693。交错网格 FDTD 的原始论文。
  2. John B. Schneider, Understanding the FDTD Method。开放教材；可重点阅读一维更新、数值色散、二维与三维算法等章节。
  3. Meep 文档：Introduction、Materials、Perfectly Matched Layers。用于理解实际求解器中的网格、色散材料与吸收层。
  4. Ansys Optics：Introduction to metalens workflows。比较不同尺度下的超构透镜建模方法。
  5. Ansys Optics：Small-Scale Metalens – Field Propagation。单元响应、场重建、有限器件仿真及传播的应用示例。