去年年中的时候我给自己定了一个有点“硬核”的任务复现一篇发表在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中心波长附近需满足亚波长条件避免高级衍射纳米柱高度600 nm固定高度通过直径调相位纳米柱直径80 nm ~ 320 nm扫描区间步长设为 5 nm基底材料熔融石英折射率约 1.46纳米柱材料TiO2可见光波段折射率约 2.4为什么把周期定在 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π相位覆盖在多个波长下同时成立”这件事上。如果你也正在复现类似的光学超表面论文我的建议是先把单元库地基打好因为后面所有看似华丽的结果归根到底都是那些不起眼的扫描数据支撑起来的。