做光子晶体超表面绕不开透反射相位计算。最近我在COMSOL里模拟一个硅基反射式超表面本来只是想算透射和反射相位覆盖范围顺手把入射角扫了一圈结果发现相位曲线在共振波长附近变得非常陡进一步算Goos-Hänchen位移时反射光束的位置能偏出好几十微米。很多初做超表面的人只关心“相位能不能覆盖2π”其实相位对角度或者空间频率的斜率同样重要GH位移就是它最直观的表现。这篇文章把我在COMSOL里做这套仿真时踩过的坑、整理出来的流程和几个常用脚本一起分享出来适合正在做超表面单元设计、导模共振器件或者打算处理COMSOL透反射相位数据的朋友参考。1. 项目拆解相位、超表面和GH位移之间是什么关系1.1 透反射相位计算解决什么问题超表面这类人工周期结构相比传统透镜和波片的核心优势在于能把空间每个位置的相位单独设计。我们用介质柱、孔洞或者金属天线组成周期性排布的单元每个单元相当于一个微型谐振器。波穿过或者被它反射时振幅基本不变但相位可以依据几何尺寸和材料被大幅改变。只要阵列里每个单元的透反射相位都能独立控制就能拼出任意波前这也是超构透镜、异常反射、涡旋光束等器件的基本原理。所以透反射相位计算是所有超表面设计的起点。实际工程中通常先建立单个单元的“单元库”固定材料和周期扫柱子半径或者长度记录每个变体对应的反射/透射系数其中就包括相位。最理想状态是一条从0到2π的平滑覆盖曲线。很多器件还要求反射率或透射率在目标频段内尽量高这就要用参数扫描找出振幅和相位的平衡点。如果只靠手算或者查理想模型很难准确反映真实结构中的法布里-珀罗干涉、材料损耗和邻近耦合效应这也是必须用COMSOL这类全波仿真工具的原因。1.2 GH位移到底意味着什么Goos-Hänchen位移用大白话说就是光打在界面上被反射时反射光并不是严格从“几何接触点”返回而是在界面附近滑行了一段后才离开肉眼可能看不见但仪器能测出来。传统全反射棱镜中这个位移通常在波长量级。对于很多宏观光学系统来说可以忽略但在精密测量、传感器设计或者光束对准任务里就会变成系统误差。超表面和光子晶体中GH位移被放大的机制很有意思。因为单元结构内部往往存在导模共振、连续域束缚态或表面等离激元共振这些模式让反射相位和透射相位对入射角非常敏感。相位变化越陡GH位移越大。用严格的电磁场语言解释反射光束本质上是一组不同平面波分量的叠加每个分量因为入射角不同反射相位也不同。当相位在某个角度范围内发生急剧变化时不同平面波分量的相干叠加会使光束包络中心发生横向移动。这些内容都直接和透反射相位计算关联所以我更喜欢把GH位移看作“相位对角度的导数”而不是额外引入的玄学现象。1.3 我给自己定的仿真目标这次项目我主要想回答三个问题一是确定在1550 nm通信波段下硅基超表面反射/透射系数和相位随几何参数的变化规律二是把一个足够薄的周期性结构放到COMSOL里准确提取S参数相位并扣除参考面影响三是利用反射相位对入射角的斜率估算GH位移看看在共振角度附近位移能达到多大。整个过程从建模到后处理核心还是落在透反射相位计算上。无论你是做反射式超表面还是透射式超表面只要核心参数是相位下面这套思路都能复用到你自己的结构上。2. COMSOL建模准备几何、端口和材料怎么设置2.1 模型简化先算单胞再组装阵列直接在COMSOL里画几百几千个纳米柱阵列虽然能跑但计算量会很大而且未必更准确。绝大部分透反射问题通过单胞模型加周期性边界就能描述无限阵列的行为。单胞建模的好处是不光节省网格还能天然引入Bloch周期性条件后续波矢或入射角变化都是很方便的连续参数。我这次用三维模型一个硅柱放在二氧化硅基底顶部四周是空气或者匹配介质整个模型的横向外边界为一个晶格周期。支柱高度设为220 nm周期900 nm半径从50 nm到300 nm做参数扫描。在层结构上顶部空气层留出约一个波长的距离底部衬底也留出几百纳米避免端口边界离结构太近产生数值倏逝波影响。因为要研究的只是单胞的反射透射并不需要把大尺寸阵列全部画出来。很多初学者到了这里会犹豫要不要在2D里简化。如果你的结构沿某个方向完全不变比如无限长光栅沟道用二维模型没问题。但如果结构是圆柱、方块这类三维单元用二维强行等效往往会丢掉偏振耦合的细节算出的相位不可靠。建议至少在阶段验证用三维单胞模型至于超大阵列的多次散射效应之后可以再做全结构验证。2.2 周期性边界、端口和Floquet周期模式COMSOL的“电磁波频域”接口是处理这类问题的主要入口。建模时要在横向相对的边界上加周期性条件。很多人会在这里混淆x方向和y方向是周期边界z方向上下两端不要加周期性条件而是加端口。具体设置步骤大致如下在“电磁波频域”接口下选择两对相对边界。应用“周期性条件”类型选“周期性”。如果是有倾斜入射还要勾选“Floquet周期性”并输入面内Bloch波矢的k_x、k_y值。在模型顶部空气边界添加端口“Port 1”作为激励端口。在模型底部衬底边界添加端口“Port 2”作为透射波接收端口。两侧的边界保持为周期性不要在端口外再加PML因为周期边界已经描述了无限周期单元的横向特性。如果直接用两个端口上下夹住结构端口本身会提供入射平面波激发并在后处理给出完整的S参数矩阵。为了让COMSOL正确计算端口模式通常会指定端口模式类型为“Floquet周期性”填入入射角。很多人第一次跑出来的相位完全不对多数是因为把θ角当成与z轴的夹角或与x轴的夹角填错了。这里我建议先固定一个已知例子比如无结构的平面介质板对比理论反射相位和仿真值确定角度定义正确后再计算自己的超表面。2.3 材料参数与扫描变量的取舍材料参数出错后续相位计算没有任何意义。硅在1550 nm附近折射率约3.45到3.48几乎无吸收但如果你研究的是高Q共振即使是极小的消光系数也会明显压低峰值反射率所以要不要加虚部取决于你想研究什么。我在仿真里给了一个主参数表方便后面看相位和GH位移。参数数值说明中心波长 λ01550 nm光通信C波段周期 a900 nm需保证工作波长不发生高阶衍射柱高 h220 nm常见SOI工艺厚度柱半径 r50~300 nm扫描得到相位库背景与填充材料空气/硅n_air1.0, n_Si3.45基底SiO2n_sub1.45入射角 θ0°~60°用于反射相位和GH位移材料色散可暂时在1550 nm附近取常数如果要扫1500到1600 nm的宽谱就需要导入硅和二氧化硅的色散曲线。COMSOL材料库自带不少数据但使用前最好看一眼损耗虚部的趋势。很多现成材料库里默认数据偏保守实际工艺制备的样品可能和理想模型差别不小。最保险的做法是把材料折射率做成可调参数先算损耗为零的理想情况再逐步加入实损耗看相位覆盖和共振深度如何变化。这样可以判断结果对你的体系敏感还是不敏感。3. 相位提取与参考面校准S参数不等于最终相位3.1 从S参数到相位arg(S11)只是一个起点COMSOL完成频域计算后可以直接在“派生值”里选“端口-S参数”得到S11和S21。S参数是复数所以相位理论上就是arg(S11)和arg(S21)。大多数超表面设计教程到这里就会画相位曲线初学者也喜欢直接看结果。但实际操作时你会发现直接从端口S参数得到的相位里混入了端口到结构之间均匀介质段的传播相位。这段均匀区的长度虽然不大但在纳米级结构里足以让相位曲线整体偏离物理结果。透反射相位计算追求的往往是“超表面本身带来的额外相位”也就是入射面到结构表面、结构本身、再到出射面的总相位变化中剔除掉没有结构时纯背景波导的传播相位。如果不做这一步后续扫描几何参数时得到的“相位覆盖范围”虽然也能看相对变化但和文献里常说的单元相位库对不上也很难用来反推GH位移。另一个常见坑是S参数在某个频点或角度会变得很小比如透射率低到接近零此时相位完全由数值噪声主导画出来的相位图会出现一些毫无规律的波动。处理方法是检查对应的振幅确认是否真的进入暗模或禁带区域。如果是该区域的相位本身没有明确物理意义不必强求光滑。3.2 去掉均匀段的相位累加参考面校准怎么做端口参考面默认位于端口边界。如果端口距离超表面表面有额外距离L那么入射波和反射波都会在均匀介质中额外传播一段距离。反射相位需要扣除双程长度透射相位需要扣除从入射端口到出射端口的均匀区传播相位。比如我在模型顶部留了500 nm空气层顶部端口在这个空气层外面底部端口在衬底下面。计算反射相位时空气层和衬底的均匀段都需要扣除。扣除公式很简单[ \phi_r^{corr} \arg(S_{11}) - 2 k_0 n_{air} L_{air} - 2 k_0 n_{sub} L_{sub} ][ \phi_t^{corr} \arg(S_{21}) - k_0 n_{air} L_{air} - k_0 n_{sub} L_{sub} ]其中k0是真空波数。如果端口斜入射还要把沿法线方向的投影距离换算进去通常是除以cosθ。后处理时直接在COMSOL表达式里写arg(S11) - 2*(2*pi/1550[nm])*1.0*500[nm] - 2*(2*pi/1550[nm])*1.45*300[nm]写成参数化形式会更方便。我习惯把基底层厚度设为参数校准表达式用参数表示以后扫参数不用每次改公式。如果不想做扣除也可以直接把端口边界紧贴结构表面。这个方法看似简洁但实际不好因为结构表面的倏逝场可能干扰端口模式计算尤其在强共振结构附近仿真误差会变大。留白再扣除反而是更稳的做法。3.3 参数扫描之后用unwrap处理2π跳变相位作为反正切结果天然被包裹在(-π, π]区间里。超表面结构在共振点附近相位常常会快速从接近π掉到接近-π画图时看起来就像一条地震波断裂带。要算GH位移或做波前设计必须判断这是真实相位快速变化还是单纯的数值包裹。判断方法很简单加密波长或角度扫描步长后观察那块“断裂带”如果断点前后相位差接近2π且曲线连续就是包裹如果相位斜率依然很陡但不再存在固定2π跳变那就是真实物理。处理包裹的常见办法是把相位数据导出来放到Python或MATLAB里用unwrap函数。COMSOL本身没有封装好的unwrap功能所以我都是做成数据导出后处理。用Python时要注意COMSOL导出的相位单位是弧度还是度。如果因为勾选习惯导出的是度数要先统一换算到弧度再unwrap否则unwrap会失效。4. GH位移计算的两种路线与后处理实现4.1 路线一Artmann公式平面波扫角度快速估算严格计算GH位移最稳妥的方法是仿真有限尺寸高斯光束然后对比入射与出射光束的质心位置。但这类模型计算量大也不便于扫参。更常用的工程近似是Artmann思路把平面波反射系数相位对横向波矢k_x求导得到反射光束中心的横向位移。简化到角度扫描数据时可以写成[ \Delta \approx \frac{\partial \Phi}{\partial k_x} \frac{1}{k_0\cos\theta}\frac{\partial \Phi}{\partial \theta} ]需要提醒的是不同文献对符号和方向的定义不完全一致你算出来的正负号取决于相位定义和坐标方向。实际使用前最好先用一个已知的介质平板反射算例验证一下位移方向是否符合预期。由于这个公式本质上是相位斜率所以它和透反射相位计算的精度强相关。只要你的相位曲线有抖动求导后会出现很大噪声。解决办法是角度扫描步长不要太大我一般扫0.1°或0.2°然后用局部线性拟合代替直接相邻差分。若共振非常尖锐还需要在共振角度附近把步长进一步加密到0.02°否则会直接跨过共振点测得的最大GH位移被低估。4.2 路线二高斯光束直接仿真Artmann公式是线性近似它假设反射过程中光束横截面形状不变只是整体平移。如果相位斜率本身在角度或频率上变化剧烈光束会同时发生畸变这时一阶导数近似就不再够用。更严谨的做法是让高斯光束从模型一端入射整个结构作为散射体然后在反射端或透射端提取场分布用质心法计算位移。这种全光束仿真可以放在2D截面里比如把三维超表面等效成二维光栅或者如果单胞本身有轴对称性就选取沿入射平面的剖面引入高斯光束作为背景场。好处是能看到光斑变形和分裂坏处是计算量明显增加而且对边界条件要求更高因为高斯光束在横向不是周期函数不能再使用无限周期单胞边界。我的做法是两步走先用Artmann公式快速扫描角度找到GH位移峰值大概出现在哪个入射角范围然后把角度固定在那几个候选值附近用高斯光束做一次全波验证。如果两者量级一致那结论就比较可靠。如果差距很大通常是因为结构共振让反射波产生了高阶衍射或者光束畸变这时候再回头检查单胞设计中有没有遗漏的衍射级。4.3 后处理导出与Python脚本计算GH位移透反射相位计算之后的GH位移数据处理我通常这样导数据在COMSOL“派生值”里使用“全局评估”表达式填S11的相位或校准后的相位勾选“参数扫描中所有参数”最终导出成一张逗号分隔的表格。第一列是入射角第二列是反射相位。下面是处理脚本示例读者可以根据自己导出的文件结构微调import numpy as np import matplotlib.pyplot as plt k0 2 * np.pi / 1.55e-6 # 假设导出的CSV有两列角度(deg)和相位(rad) data np.loadtxt(s11_phase_scan.csv, delimiter,, skiprows1) theta_deg data[:, 0] phase_wrapped data[:, 1] # 1. 相位解包裹 phase np.unwrap(phase_wrapped) # 2. 转为弧度角度 theta np.deg2rad(theta_deg) # 3. 中心差分求相位对角度的导数 dphi_dtheta np.zeros_like(phase) dphi_dtheta[1:-1] (phase[2:] - phase[:-2]) / (theta[2:] - theta[:-2]) dphi_dtheta[0] (phase[1] - phase[0]) / (theta[1] - theta[0]) dphi_dtheta[-1] (phase[-1] - phase[-2]) / (theta[-1] - theta[-2]) # 4. GH位移近似 GH dphi_dtheta / (k0 * np.cos(theta)) # 符号取决于坐标定义 # 5. 画图 plt.figure(figsize(6, 4)) plt.plot(theta_deg, GH * 1e6, linewidth1.5) plt.xlabel(Incident angle (deg)) plt.ylabel(GH shift (um)) plt.title(Reflection GH shift from phase derivative) plt.grid(True) plt.show()脚本跑出来如果发现GH位移曲线全是毛刺多半不是物理问题而是角度步长太粗或相位没有正确解包裹。还有一种可能是共振点附近S11振幅太小相位噪声被放大。此时改用复数S参数本身来算梯度而不是先取相位[ \frac{d\phi}{d\theta} \operatorname{Im}\left(\frac{1}{r}\frac{dr}{d\theta}\right) ]这样不会出现unwrap带来的跳变问题但表达式更复杂一些适合用COMSOL直接后处理或者在脚本里对复数S11做差分。5. 从仿真到设计闭环相位覆盖扫描与共振增强5.1 单元库扫描实现0到2π的相位覆盖超表面设计最关心的一个指标是相位覆盖范围。我通常把周期、柱高固定扫描半径或边长的离散值导出每个半径下的校准后反射相位。画出来就是一条类似“S”型曲线理想情况下能从0铺到接近2π并且反射率保持在0.8以上。实际项目中如果中心波长落在共振区附近反射率可能会掉得很厉害这时相位覆盖倒是够了但振幅不够均匀仍然无法做高效器件。要解决这个问题可以引入多变量扫描。比如同时扫半径和柱子高度在结果里筛出“相位每变化10°且振幅最高”的候选几何组合。COMSOL参数扫描虽然能扫多个维度但三维模型多参数扫描速度会比较慢。建议先固定一组材料与周期只看半径和高度两条扫描轴找到规律后再加密细化。如果计算资源紧还可以利用结构的对称性把全空间模型缩小到1/4区域在切割面加对称/反对称条件也能大幅提速。5.2 共振增强GH位移在传感器里的用途回到GH位移本身最大的工程价值之一就是做角度传感或折射率传感。结构在某个入射角附近发生导模共振时反射相位对角度几乎成垂直下降相位导数极大因此GH位移可以比普通介质表面大几个数量级。这种机制非常灵敏样品表面覆盖少量待测物共振角度偏移同一入射角下的GH位移就会显著改变。不过这类结构对损耗极其敏感尤其是金属超表面欧姆损耗会压低共振Q值相位斜率下降GH位移随之减小。所以在COMSOL里仿真传感器时不能只算理想无损金属必须给金的介电常数加实际复数值。金的Drude模型参数在不同来源差别不小建议至少用两种已知参数对比一下看共振波长变化是否在可接受范围内。5.3 超表面设计和异常反射不要忽略空间相位梯度还有一个容易混淆的地方传统超表面设计里说的“空间相位梯度”和这里的GH位移本质上同源都是相位对波矢或者说对横向位置的响应。当一个超表面被设计成偏折角很大的异常折射器件时整个界面上会引入很大的相位梯度。这种梯度不仅决定了出射方向也会使光束在界面附近等效产生位移。如果你在后续系统里要精确对焦或配准这部分位移会表现为系统误差。因此透反射相位计算不只是单元阶段的表面功夫它甚至可以反过来解释整个器件的宏观光学行为。6. 常见问题与排查技巧实录6.1 扫参时相位曲线在共振点附近“乱跳”这个问题我一开始也遇到过。有三种主要原因一是相位包裹没有处理曲线在图上是-π到π的硬切二是共振频率附近S参数振幅极小数值噪声让相位不稳定三是扫描步长太大跨过了完整的相位快速变化区域。处理方法分别是unwrap、检查振幅、加密步长。判断是否属于正常共振可以在共振点附近加一次极小步长扫描如果相位曲线逐步过渡且无2π跳变那物理上是正常的。不要看到尖峰就觉得仿真出错了。6.2 端口模式方向或Floquet波矢方向填错斜入射时端口里的θ、φ设置必须和几何坐标严格对应。我常用的办法是先建一个简单的裸衬底模型不加任何超表面结构然后用菲涅尔公式算出理论反射相位角或反射率再与COMSOL结果对照。如果吻合说明端口、Floquet周期性、材料折射率这三个基础环节都没有问题。这个方法虽然多花几分钟但能避免你后排查真正的超表面结构时被一堆耦合因素干扰。6.3 S参数对网格敏感需要提前收敛性检查高频电磁仿真里如果网格过粗相位曲线会有波动尤其是共振结构需要很细的网格。建议至少取两种不同最大网格尺寸比较目标频段内S11的模和相位是否变化小于设定容差。介质柱的曲面边界尽量用边界层网格或“较细化”网格不要偷懒用均匀粗网格。另一个实用技巧是观察共振点反射率最低点是否随网格加密继续加深。如果每加密一次反射率就变低一次说明还没有收敛因为真实介质共振在无损耗下反射率应接近确定值不会一直随网格漂移。6.4 几何报错与CAD内核拓扑问题有人会碰到COMSOL提示CAD内核转换时不支持某些拓扑操作通常是因为把所有几何体用布尔并集困在了一起产生退化边或零厚度面。解决方法是把基底、柱子和空气域建为独立几何域只在物理边界条件中通过“连续性”完成连接或者使用“分区”功能让重合边界自动对齐。简单说能用装配体就尽量装配体不要盲目做并集。6.5 时谐约定和材料符号问题COMSOL电磁模块默认时谐约定是e^{-iωt}许多光学教材习惯用e^{iωt}。当你从文献导入折射率虚部或者色散拟合参数前要看清楚对方用的是一哪种约定否则材料增益和损耗会颠倒相位也会莫名其妙差个负号。这里没有更轻松的捷径只能多验证多对比。真要深究可以用一块已知折射率的均匀介质板做透射相位对比测试把约定问题一次性试出来。我在实际做这个项目时还有一个体会透反射相位计算和GH位移不是两个割裂的任务而是一条完整数据链。相位算得准位移趋势才靠谱相位曲线有一点系统偏移可能相位覆盖图仍然好看但位移数值就会偏得离谱。建议每一步都保留好校准后的相位结果再做任何求导或空间处理别所有数据堆到最后一口气算完。后续如果你想把同一套流程扩展到大尺寸超表面透镜或者可调谐器件只要把这里学到的“相位参考面校准unwrap”基本功打扎实基本就能顺下来了。