复现任意阶宽带贝塞尔光束超表面：FDTD建模全流程解析

去年年中的时候，我给自己定了一个有点“硬核”的任务：复现一篇发表在Light: Science & Applications上的超表面论文，题目方向是宽带、任意阶贝塞尔光束。当时我手头的工具是Lumerical FDTD，目标很明确，就是从零搭一个模型，把论文里的现象在仿真里跑出来。说实话，这种复现工作最折磨人的不是FDTD软件本身，而是那些论文不会写进去的细节：单元结构怎么扫、相位数据库怎么建、宽带条件怎么加、布局生成器怎么写、仿真结果怎么判定“我复现对了”。这篇博客就把我当时从读论文到仿真收尾的完整过程捋一遍，重点放在贝塞尔光束超表面模型从理论到FDTD仿真的落地细节，尤其是任意阶数和宽带这两个看起来轻描淡写、做起来全是坑的地方。适合正在复现超表面论文、或者准备用FDTD做微纳光学仿真的朋友参考，哪怕你还不熟悉贝塞尔光束，跟着我的思路走一遍也能建立起完整的模型框架。

1. 这篇Light论文究竟做了什么：贝塞尔光束的超表面实现逻辑

1.1 贝塞尔光束的“无衍射”到底是怎么来的

贝塞尔光束是亥姆霍兹方程在柱坐标下的一组特殊解，它的理想形式是横向电场分布随时间演化保持不变。换句话说，一个理想的贝塞尔光束沿传播方向走多远，它的中心光斑尺寸都不会扩散。这和我们常见的普通圆光斑完全不同——高斯光束发散，平面波也会因为衍射扩展，唯独贝塞尔光束理论上可以一直保持“细而直”的状态。

不过理想贝塞尔光束携带无穷大能量，现实中永远无法严格生成，实际使用中用有限孔径去截取它的一段近似无衍射区域。这段区域的长度大致由孔径半径 R 和圆锥角 β 决定，近似关系是 Z_max ≈ R / tanβ。我刚开始复现时总想着用尽量大的孔径去提升无衍射距离，后来发现孔径一旦大了，FDTD的计算量就会呈指数级上涨，这在后面章节会专门说。

贝塞尔光束还有一个“阶数”的概念。零阶贝塞尔光束中心是一个亮核，这是最常见的光镊、成像应用场景；高阶贝塞尔光束在中心存在相位奇点，中心是暗的，横向强度分布呈环形，而且携带轨道角动量（OAM），阶数 l 就是拓扑荷，一束 l 阶贝塞尔光束意味着波前在绕传播轴一圈时累计积累了 2πl 的相位变化。论文里的“任意阶”，指的就是可以按需生成 l = 0、1、2、3 任意拓扑荷的贝塞尔光束。

1.2 超表面怎么把轴棱锥和涡旋板合成一个平面结构的

传统贝塞尔光束生成器件，从物理上看需要两个东西：一个是轴向的锥形相位，负责把入射平面波“掰弯”到同一圆锥面上，形成无衍射干涉；另一个是方位角方向的螺旋相位，负责给光束加上轨道角动量。轴向锥形相位用玻璃轴棱锥就能做到，螺旋相位用螺旋相位板也能做到，把两者叠加起来就成了任意阶贝塞尔光束发生器。

问题在于这类传统器件体积大、加工精度粗糙，而且锥角和螺旋相位都是固定死的，换个波长就不匹配。超表面方案的核心思路，是用一层亚波长纳米柱阵列，在同一平面上同时编码这两部分相位。每个纳米柱相当于一个微小的相位延迟元件，通过调整它的直径、高度或者旋转角，可以让经过它的局部波前获得目标相位延迟。把一个区域的相位分布按照锥形加螺旋的解析公式排布出来，就能在一层平面上实现传统两个光学元件串联才能完成的功能。论文里用的应该就是介质纳米柱方案，这种方案的优势是透过率高、相位调控灵活，而且可以做得非常轻薄。

1.3 宽带为什么难，论文的分量就在这

只做单一波长的任意阶贝塞尔超表面，其实难度低不少——只要在一个波长下找到一组纳米柱几何参数覆盖 0 到 2π 的相位就行。但“宽带”两个字把难度拉到另一个量级。

原因是超表面单元的相位响应是色散的。同一个纳米柱，在不同波长下的相位延迟并不相同，而且这个变化不是线性的。如果你的目标相位分布是 φ(ρ, φ) = -(2π/λ)·ρ·sinβ + l·φ，那么问题就来了：不同波长下，锥形相位项的斜率本身就和 λ 相关，如果单元色散不能把这个差异补偿掉，那么你会看到不同颜色的光从同一个超表面出来时，具有不同的圆锥角——这就不是宽带贝塞尔光束了，而是把白光打散成了一个个不同锥角的单色环。

论文给出的解决思路，本质是通过合理设计纳米柱的几何色散，让整个超表面单元在目标波长范围内等效相位响应满足“色散匹配”条件。复现时这一部分最难落地，因为你不能只在一个波长下扫描单元库，而要在多个波长下同时评估相位误差，再做组合优化。

2. 从零搭单元库：材料选型、结构扫描与相位数据库构建

2.1 我复现时选用的材料与几何参数

复现的第一步是确定单元结构。论文工作波长在可见光范围，我选了 TiO2 纳米圆柱，基底为熔融石英。之所以选 TiO2，是因为它在可见光波段折射率较高、透明性好，用不太高的柱子就能积累足够的相位延迟。如果你工作在近红外，可以考虑非晶硅方案，折射率更高，能进一步缩小器件尺寸。

我当时用的具体参数范围如下表，注意这些数值是基于常规复现实践的合理选择，不同论文的具体尺寸可能有出入，但参数扫描的流程完全一样：

为什么把周期定在 500 nm 而不是更小？因为圆柱直径最大有 320 nm，如果周期太小，相邻柱间距会被压缩，加工上不好实现；如果周期太大，超表面会出现非亚波长衍射级，影响无衍射光束的质量。这个平衡在做参数扫描之前就要心里有数。

2.2 用参数扫描脚本建立0到2π相位库

单元库扫描是这个复现工程里最耗时、最机械、但也最基础的一步。我当时用 Lumerical FDTD 的脚本模式批量跑，核心流程是：每个直径对应一个 500 nm 周期的单单元仿真，用周期边界条件模拟无限阵列，平面波正入射，提取透射场的复振幅，从而得到该直径下的透射相位。

写扫描脚本时有一个容易踩的细节：相位提取要扣掉参考面带来的背景相位累积。我第一期没注意，在基底内部放了一个监视器，结果直径和相位曲线整体带上了一个随波长变化的偏置，导致后续布局映射全部出错。后来改成在纳米柱上方空气区域取一个固定参考面，并且每个直径相对同一个“无结构参考”做相位差，才算真正拿到干净的相位响应。

扫描完后的典型结果是：随着直径从 80 nm 增到 320 nm，相位从接近 0 一路爬到 2π 附近，呈现一条近似单调的 S 形曲线。但单调并不意味着可以直接用——关键是看这条曲线的斜率分布，尤其在宽带约束下，不同直径的相位色散特性差异很大。数据库建完后，不能只存中心波长的相位，要把多个采样波长下的透射相位全部存下来。后面做宽带布局时，这个多波长数据库就是选柱的基本依据。

2.3 多波长约束下的单元选型思路

单波长设计，你只要找一个直径使得单元相位尽可能接近目标相位即可；多波长宽带设计，要求是同一个直径在所有设计波长下的相位误差总和最小。比如我设计波段取 450 nm、550 nm、650 nm、750 nm 四个采样波长，那么对于布局中某个位置，目标相位在不同波长下分别是 φ1、φ2、φ3、φ4，我需要找到数据库中直径 D 使得加权误差 ε(D) = Σ wλ·|φλ(D) - φtargetλ|² 最小。

注意，如果目标相位公式里含有 1/λ 项，不同波长的目标相位方向完全一致，相对值随波长差异很大，这时仅用直径一个几何自由度往往不够。更稳妥的单元库要做“直径 + 高度”二维扫描，同时记录多个波长的相位响应。我在前期只扫描了直径一维，后来发现宽带约束下无法同时满足四个波长的相位匹配，于是补充了一组高度参数（500 nm、600 nm、700 nm），构建了二维查找表。这一步虽然增加了扫描时间，但对宽带复现来说是必须的。

3. 布局生成器：把理想相位掩膜变成纳米柱坐标

3.1 目标相位掩膜的数学形式与物理含义

超表面布局生成的核心，是把目标相位掩膜离散化到每一个单元位置。任意阶宽带贝塞尔光束的理想相位掩膜可以写成：

φ(ρ, φ) = -(2π/λ₀)·ρ·sinβ + l·φ

其中 ρ 是径向坐标，φ 是方位角，β 是设计圆锥角，l 是拓扑荷，λ₀ 是设计中心波长。物理上，第一项是锥形相位，第二项是螺旋相位，二者线性叠加。因此超表面只需一个相位层就能同时完成“轴棱锥 + 涡旋相位板”两个功能。

有一点要特别注意，就是相位掩膜的符号约定。不同文献对符号的取法不同，可能影响涡旋的旋转方向。复现时以你自己设定的坐标方向为准，保持整个仿真流程中符号一致即可，不要一边用正的螺旋项、一边用负的锥形项，否则结果会自相矛盾。

3.2 Python生成超表面布局的完整流程

布局生成器我直接用 Python 写的，核心就是遍历极坐标网格，计算每个位置的目标相位，然后在单元数据库中查找最优几何参数。拿 l = 2、β = 1.3°、中心波长 600 nm 举例，生成过程如下：

import numpy as np # 设计参数 lam0 = 0.6 # 中心波长，单位微米 beta = np.radians(1.3) # 圆锥角 l = 2 # 拓扑荷 pitch = 0.5 # 单元周期，单位微米 R = 15.0 # 超表面半径，单位微米 # 遍历网格坐标 coords = [] for ix in range(-int(R / pitch), int(R / pitch) + 1): for iy in range(-int(R / pitch), int(R / pitch) + 1): x, y = ix * pitch, iy * pitch rho = np.hypot(x, y) if rho > R: continue # 目标相位掩膜 phi = -2 * np.pi * rho * np.sin(beta) / lam0 + l * np.arctan2(y, x) phi = np.mod(phi, 2 * np.pi) coords.append((x, y, phi)) # 根据phi查数据库，得到对应的纳米柱直径 # 这里省略了数据库查找函数，实际流程中这一步会做多波长加权寻优

这段代码只是很基础的流程原型，但复现时它的作用很关键。你会发现目标相位中螺旋项在中心处存在奇点，即 ρ = 0 时 arctan2(y, x) 没有定义。处理方式我后面专门讲，这个地方如果不处理，布局生成出来的中心区域会非常难调整。

实际工程版的布局脚本还应该包括：把坐标和直径输出为 CSV 文件、生成 Lumerical 可直接导入的脚本、检查相邻柱的最小间距是否满足加工能力。我当时没有单独画 GDS，直接在 Lumerical 里用脚本循环添加圆柱体结构，对纯仿真验证来说足够；如果后续准备送加工，则需要再写一层 GDS 导出逻辑。

3.3 采样密度、相位跳变与工艺约束的取舍

布局生成过程中有三组矛盾需要平衡。

第一，采样密度与仿真内存的矛盾。周期越小，单元分辨率越高，无衍射光束质量越好，但同样的器件面积单元数量成倍增加，FDTD 计算量暴涨。亚波长原则下 500 nm 周期已经能保证不错的效果，实际复现时不必一味压缩。

第二，相位跳变与数据库覆盖范围的矛盾。相位掩膜在 2π 边界处会产生跳变，数据库在边界附近的相位覆盖如果存在空洞，生成的布局会产生局部误差。一个常见做法是把目标相位按 2π 取模，同时确保数据库中 0 到 2π 区间都有密集覆盖，必要时对数据库做插值，而不是只保留扫描过的离散点。

第三，单元尺寸的工艺约束与光学性能的矛盾。扫描数据库时理论上直径越细越好，但实际加工线宽有限制。我当时在数据库里保留了 80 nm 到 320 nm 的区间，布局查找时如果某个位置的最佳直径落在边界附近，我就会检查相邻位置是否同样落在边界区域，避免出现大面积“顶格”填充，因为工艺偏差下这些位置的相位误差最敏感。

4. FDTD仿真配置与结果提取

4.1 光源、监视器、边界条件的设置

布局生成之后，进入正式的 FDTD 仿真环节。光源我用的平面波，从基底一侧垂直入射，覆盖整个超表面区域。偏振方面，介质圆柱单元本身没有很强的偏振依赖性，用线偏振即可，不影响贝塞尔光束的基本形态。

边界条件一律使用 PML，但注意 PML 距离超表面的侧向边缘要留至少一个波长的缓冲区域，否则光在器件边缘会有非物理反射，污染轴向场分布。我一开始把 PML 贴得很近，仿真出来无衍射条纹在尾部出现明显的驻波调制，排查了很久才发现是边界反射导致的。

监视器核心放三个：

• 一个横向 XY 剖面监视器，放在超表面后方的某个距离处，记录横截面光强和相位分布；
• 一个纵向 YZ 剖面监视器，完整记录光从超表面出来后沿传播方向的演化；
• 一组轴向方向的多点监视器或线监视器，用于直接读取中心光强随传播距离的变化。

4.2 从场剖面提取贝塞尔光束的阶数与锥角

跑完仿真后，第一个看的是横向剖面。零阶贝塞尔光束应该看到一个中心亮核加一圈圈旁瓣，亮度分布沿半径方向满足零阶贝塞尔函数平方的形式；高阶 l 阶贝塞尔光束的中心应呈暗核，横向强度呈环形结构。观察这些特征时，图得用强度显示，而判断拓扑荷是否准确需要提取相位分布：在横向剖面上看相位等值线绕中心一圈累计变化多少倍的 2π。我复现时为了确认这一点，直接把相位剖面导出来，在极坐标下数相位曲线绕了几圈，绕一圈就是拓扑荷 1，绕三圈就是拓扑荷 3——这个方法比肉眼观察强度分布可靠得多。

纵向剖面则是验证无衍射特性的关键。理想状态下，你会看到一束“不发散”的光沿着传播方向走很长的距离。把不同传播距离处的横向剖面拉出来对比，如果中心亮核的半径保持稳定，就说明这一段距离内无衍射特性成立。测量得到的无衍射长度可以换算回 Z_max = R / tanβ，来验证器件口径范围下的理论极限是否对得上。

锥角的验证有好几种方式，最直观的是做空间傅里叶变换，在频域出现一个圆环，圆环的半径对应贝塞尔光束的圆锥角。Lumerical 的远场投影功能可以直接把横向近场转换到角谱空间，环的半径与设计值 sinβ 对比即可。我第一次做的时候没注意这个环节，直到后来用角谱验证才发现单元库的相位色散导致锥角偏移了大约 0.1°，问题出在宽带布局的多波长权重分配对中心波长不够偏重。

4.3 宽带性能怎么量化评估

宽带复现不是跑完一个波长就完事。我在设计波段内取了 450、550、650、750 nm 四个典型波长分别跑仿真，再评估三个指标。

第一，每个波长下的横向强度剖面是否还保持贝塞尔形态。这个判断很直观，但也很容易被忽略——因为宽带单元色散会导致某些波长下旁瓣结构明显畸变，强度剖面已经不是理想的贝塞尔函数形式。

第二，每个波长下提取得到的锥角是否保持一致。理想的宽带设计是各波长拥有相同圆锥角，如果某个波长的锥角偏差过大，说明单元库中该波长的相位匹配没有做好，需要回头调整加权权重，优先压低偏差最大的波长的误差贡献。

第三，透过率或能量效率的趋势。宽带超表面设计往往存在效率随波长起伏的问题，典型原因是材料折射率色散导致纳米柱相位覆盖范围在不同波长下不一致。通过记录透过率随波长的变化曲线，可以判断这个设计到底是在全波段“平均兼顾”还是在某个波段“局部最优”。

5. 复现中踩过的坑与针对性调优

5.1 高阶贝塞尔的核心奇点被网格抹掉的伪影问题

高阶贝塞尔光束的核心在理论上是绝对暗的——那里是相位奇点，强度严格为零。但我在 FDTD 仿真结果里看 l = 2 的中心区域，竟然出现了一个小小的亮斑，一度让我以为布局生成错了。

后来排查发现，这是数值网格在相位奇点附近的固有缺陷。当仿真网格划分没有让采样点恰好避开中心奇点时，离散化过程会把几个网格点的相位平均掉，形成强度非零的伪影。解决办法有两个方向：一是在布局生成时把中心位置留一个专门的空白单元或让中心落在网格间隙，这比较敏感；二是把中心区域的网格局部细化，并且结果提取时多取几个不同的中心位置来统计强度最低值。我实测下来，把网格从 10 nm 细化到 5 nm 能明显抑制伪影强度，但计算时间几乎翻倍。如果不是盯着中心区域做量化研究，通常不需要对整个模型做这种细网格，只对中心局部区域开细化网格就够了。

5.2 大尺寸模型的内存开销与降规模验证策略

这是复现过程中最影响效率的“大坑”。贝塞尔光束器件通常设计成几十微米口径，抄着论文里的全尺寸参数直接跑，FDTD 网格会达到数千万甚至上亿个，内存直接爆掉。哪怕有 GPU 版本，全尺寸三维模型的单次运行也可能在几十个小时量级，根本没有办法做参数调优。

我的策略是分阶段降规模验证。先跑一个口径只有 10 微米左右的原型模型，确认整个流程——布局生成、单元库匹配、仿真设置、结果提取——没有基本功问题；再在原型模型基础上做参数扫描调优，比如调整多波长权重、修正锥角偏差；最后才根据计算资源跑一个中等尺寸的模型验证最终的场分布特征。论文里的无衍射长度和效率数据，是在他们那个口径条件下成立的，复现时并不是必须全尺寸重现才算“做出来”，验证设计逻辑正确、场分布符合预期，就已经达到复现目标。

另外需要提醒一句，贝塞尔超表面布局是螺旋结构，不具备镜像对称性，所以 FDTD 的对称/反对称边界条件在这里基本用不上，能用的只有减小来源的“单次小模型策略”和利用 GPU 加速。想靠对称边界省内存的同学可以省省心了。

5.3 相位拟合误差与透过率的矛盾权衡

宽带布局的核心就是找“最优直径”，但这个“最优”本身不是唯一的，不同评价标准会得到不同的布局。如果只看中心波长的相位误差，其他波长偏差大，结果就是宽带性能差；如果对四个波长等权加权，中心波长的锥角反而不一定最准。

我最终发现，宽带复现中相位拟合误差和透过率有一个隐蔽的权衡关系。有些直径组合虽然相位匹配得很好，但对应的透射振幅很低，光功率大量被反射或吸收，实际生成的贝塞尔光束看似漂亮，强度却弱得离谱。更合理的目标函数应该同时评估复振幅误差——相位差和透射振幅差二者合并计算，而不只是相位差。我在后期把单元库的选型标准从最小相位误差改成了最小复振幅误差后，仿真出来的透过率提升非常明显，而贝塞尔形态几乎没有劣化，这对后续想送加工把器件做出来的朋友尤其重要。

5.4 周期性边界反噬与宽谱光源的效率陷阱

单元库扫描阶段会用周期边界理论模拟无限阵列，这个阶段提取的相位数据和实际布局中纳米柱之间的耦合效应存在偏差。尤其当相邻柱直径差异很大时（比如直径 80 nm 旁边放了一个 300 nm 的柱子），近场耦合会带来额外的相位偏移。我在第一阶段布局仿真时直接用数据库相位匹配，结果实际场分布和理论剖面有明显错位。

解决这个误差的思路是迭代修正：先用当前布局跑一次有限阵列仿真，提取实际波前相位剖面，和理想剖面做差，然后把误差补偿回布局公式，重新生成一版布局。这个过程迭代两三轮后，相位误差可以降到很低的水平。如果计算资源紧张，也可以先只修正布局中误差最集中的几个环带区域，不必全盘重来。

宽谱光源方面，FDTD 确实支持宽频脉冲一次跑多个波长，看起来比逐波长扫描快很多。但要注意，宽谱光源的仿真结果受材料色散精度影响极大，TiO2 在可见光波段的折射率曲线如果不准确，多个波长下的单元相位数据库会整体失真。稳妥的做法是先逐波长扫描建立单元库，确认各波长下单元相位响应和实际结构匹配，再在最后的整体仿真阶段用宽谱光源提速验证趋势，两套结果互相印证。我在复现时为了省时间直接用宽谱光源建库，数据库和最终模型的相位响应对不上，返工了一轮，这个教训还是挺贵的。

做完整套复现之后再回头看，最核心的经验就是：复现论文不是照着摘要做一遍仿真，而是要把论文每一条结果背后的设计约束拆出来，变成自己可执行的单元库、布局脚本和仿真流程。贝塞尔超表面这个方向，理论公式就一行，真正的工程量全在“如何在一个实际制造的单元结构上，让0到2π相位覆盖在多个波长下同时成立”这件事上。如果你也正在复现类似的光学超表面论文，我的建议是先把单元库地基打好，因为后面所有看似华丽的结果，归根到底都是那些不起眼的扫描数据支撑起来的。