COMSOL光子晶体超表面模拟：GH位移计算与避坑指南

做超表面模拟的人看到“GH位移”这个词第一反应多半是这玩意儿不是光学课本里全反射那个几微米的横向偏移吗跟我的光子晶体超表面有什么关系我第一次在透反射相位计算里碰到它也是这个反应。实际做下来才发现只要结构对入射角敏感反射光的等效位置就会偏离几何反射点——这就是广义的Goos-Hänchen位移。超表面的相位梯度可以比普通全反射界面大好几个数量级所以GH位移也能从微米级被放大到几个波长的量级。这篇文章就围绕光子晶体超表面的透反射相位计算把GH位移怎么在COMSOL里算出来、怎么避免踩坑一次性讲清楚。适合谁来参考两类人最有共鸣一类是刚用COMSOL做超表面仿真、对着S参数相位图发懵的同学另一类是已经算出透反射相位但不知道下一步怎么跟光束位移挂上钩的研究者。这里不堆公式也不贴大段理论按我自己跑通的流程来先讲物理图像再讲模型怎么搭然后讲透反射相位的提取细节最后讲怎么从相位梯度得到GH位移以及我踩过的几个隐蔽的坑。1. 先搞懂GH位移到底是什么以及它为什么值得算1.1 从全反射说起相位梯度才是位移的本质经典的Goos-Hänchen位移描述的是全反射时反射光相对几何反射点产生的横向偏移。教科书里常见的解释是入射光束其实不是“一个点”在反射而是一簇不同入射角的平面波分量在竞争每个分量反射时相位变化不同叠加起来等效反射面就后退或前移了。这句话有点绕换成年间的说法相位梯度就是位移的“驱动力”。一束有限宽的光可以拆成一系列朝不同方向传播的平面波每个平面波在界面处经历的反射相位都不一样。反射相位随入射角变化越快这束光叠加之后等效的反射点就偏离得越远。严格公式是[ L_s -\frac{1}{k}\frac{d\phi_r}{d\theta_i} ]这里的 (\phi_r) 是反射相位(\theta_i) 是入射角(k) 是入射空间波矢大小。负号表示相位随角度增加而减小时位移为正。用横向波矢 (k_x k\sin\theta_i) 来表示的话还有更常用的等价形式[ \Delta -\frac{d\phi_r}{dk_x} ]在COMSOL里扫参时我们通常直接扫 (k_x) 或扫入射角再用这个式子换算。1.2 超表面的GH位移为什么值得关注普通全反射界面的GH位移通常只有波长量级几个微米已经算大了。但超表面的结构单元可以提供更强的色散和更大的相位梯度GH位移可以被放大到几个甚至十几个波长。更关键的是超表面还可能存在“负位移”——反射光等效点往入射点前方偏移这在普通反射界面几乎看不到。这个特性直接关系到几个实际应用传感位移对入射角和波长极其敏感微小折射率变化会引起位移显著变化比单纯测反射率锐利得多。光束操控利用GH位移可以把光束在空间上精准搬移做成光开关或光路由器件。模式匹配在波导耦合场景里位移大小决定耦合效率很多片上光子器件设计都需要量化这个参数。所以“透反射相位计算 GH位移”这个组合本质上是在回答一个问题我的超表面把光束“挪”了多远、往哪个方向挪这个值在实验里可以直接用位置灵敏探测器测出来在模拟里则要从相位里抠出来。1.3 模拟中的GH位移不等于全反射公式照搬需要特别提醒一句超表面里的“反射相位”和“透射相位”都不是简单的一个常数而是随入射角变化的函数。因此反射光的GH位移和透射光的GH位移要分开算两套相位曲线导出来的结果可能符号都相反。我习惯在模拟开始前先在纸上把目标明确下来我要算的是反射GH位移、透射GH位移还是两者都要入射侧和出射侧介质是不是同一种我只关心TE横电还是TM横磁还是两个偏振都要扫这三个问题没想清楚就开模型很容易出现“模型跑通了结果却没法解释”的局面。2. COMSOL里光子晶体超表面模型的搭建细节2.1 模型维度与物理场选择2D通常就够了做光子晶体超表面的GH位移计算我绝大多数情况用二维模型就足够了。只要结构在第三个方向上是均匀的比如无限长介质柱组成的周期阵列二维模型和三维模型在物理上是等价的而计算量差了一个数量级。在COMSOL中选择“电磁波频域”接口二维模型下可以分别计算TE和TM模式。TE模式对应电场沿柱体轴向面外方向TM模式对应磁场沿柱体轴向。这里最容易被新手绕晕TE波电场只有 (E_z) 分量边界条件是 (n \times E 0) 在完美电导体上对应切向电场为零TM波磁场只有 (H_z) 分量对应的是 (n \times H 0) 的边界条件。我做仿真时习惯先用TE模式验证流程因为TE模式在介质柱光子晶体里通常带隙更清晰透反射曲线也更平滑方便排查模型设置问题。TM模式等整个流程跑通后再切换。2.2 几何建模与参数化把周期、半径、高度这些变量全部参数化我常用的一套假想超表面参数如下仅作示例你可以按自己课题替换周期 (a 800 \text{ nm})介质柱半径 (r 0.2a 160 \text{ nm})介质柱高度 (h 400 \text{ nm})柱体折射率 (n 3.5)模拟典型半导体材料背景介质折射率 (n_0 1.0)空气衬底折射率 (n_s 1.45)典型玻璃几何建模时最重要的是把“仿真区域”和“结构区域”分开。我的建议是空气层上下各留至少1个波长的厚度作为端口区和模式衰减区介质柱用圆柱体在2D里是圆布尔操作之后与背景空气域合并成计算域左右两侧用周期性边界条件不是PEC/PMC而是“周期性”边界条件里的Floquet周期性用来支持斜入射。所有几何尺寸都设成参数不要用硬编码数值。原因很简单后面对周期、半径做参数扫描时如果尺寸是硬编码的每次改都要重新设一遍几何太容易出错。2.3 端口设置与Floquet周期性的配合端口在COMSOL里可以直接作为激励和S参数提取的载体。在电磁波频域接口下端口类型选“端口”模式类型选“衍射阶次”这样端口会按Floquet模式自动展开各级衍射波。这里要强调一个很容易出错的点端口必须放在周期性边界条件的“同一侧”对齐位置且端口所在的边界必须是平直的不能有结构穿越。也就是说介质柱不能延伸到端口边界上否则端口上的场不是均匀的平面波S参数参考面就乱了。我通常在上下边界各留500 nm以上的均匀空气层再把端口放在最外边界。这样做的意义有两个端口面上的模式场是纯平面波S11/S21的提取值干净计算GH位移需要知道“结构内部”的相位端口位置偏远了也没关系后面可以用传输距离的相位项扣除。关于Floquet周期性需要注意的设置在“周期性端口”的波矢 (k_x) 指定上如果直接扫入射角 (\theta)需要设置 (k_x k_0 \sin\theta)其中 (k_0 2\pi/\lambda)如果直接扫横向波矢 (k_x)要注意它和入射角是一一映射关系但 (k_x) 扫到超过 (k_0) 时对应的是倏逝波入射物理上要谨慎解释。个人经验计算GH位移时扫入射角比扫 (k_x) 更直观因为你最终想看的是位移随角度的曲线。但扫 (k_x) 有两个好处一是数值上更稳定二是对相位梯度公式 (\Delta -d\phi/dk_x) 直接适用。我实际做法是先按入射角扫把物理趋势看清楚要精细求导时再用 (k_x) 的均匀网格重新扫一遍。2.4 网格剖分相位计算对网格的要求比强度更高做光子晶体超表面仿真网格密度是个绕不开的话题。如果只是算反射率/透射率网格稍微粗一点结果也不会差太多但算相位、尤其是相位随角度的导数网格太粗直接带来两个问题相位曲线不光滑数值微分后噪声巨大S参数在不同入射角下网格一致性差导致相位曲线出现“阶梯状”跳动。我的经验是介质柱内部至少保证每波长10个二阶单元空气区域可以适当放宽到每波长6~8个。COMSOL里可以用“波长比例”来控制也可以手动指定最大单元尺寸。这里再补充一个实际心得做参数扫描比如扫描入射角时最好把网格固定下来不要每个角度重新自适应加密。自适应性会让每个角度点的网格略有不同最终相位曲线的“数值噪声”会比固定网格大好几倍。我的做法是先粗网格跑一遍看结果趋势再统一加密固定网格跑正式计算。3. 透反射相位的提取与计算从S参数到连续相位3.1 S参数相位怎么读angle和unwrap的区别COMSOL里端口S参数的结果变量一般是emw.S11、emw.S21这样的复数。很多人直接对复数取角度然后发现相位曲线在 (\pm\pi) 之间来回跳根本没法用。原因很简单atan2得到的角度被限制在 ([-π, π])而真实的相位梯度是连续的跨过边界时被“折叠”了。解决办法是先对角度做“解卷绕”unwrap也就是在相邻数据点之间如果相位跳变超过 (\pi)就自动加上或减去 (2\pi) 的整数倍把相位曲线还原成连续的形式。在COMSOL后处理里直接做unwrap不太方便我的流程是在“派生值”里把angle(emw.S11)和angle(emw.S21)导出为文本或CSV用Python或MATLAB读取数据用numpy.unwrap或MATLAB的unwrap函数解卷绕再做后续处理。这里有个细节unwrap默认是按“相邻两点相位差是否大于π”来判断的但如果扫描步长太大真实相位差真的超过πunwrap就会失效。所以扫描入射角时步长不要超过相位变化 (π) 对应的角间隔。3.2 参考面扣除这一步不做GH位移全是错的端口位置到结构表面之间有一段均匀空气层这段长度会给S参数相位加上一个额外的传播项。假设端口面到结构上表面的距离是 (d)那么[ \phi_{\text{port}} \phi_{\text{structure surface}} k_z d ]其中 (k_z k_0 \cos\theta_i) 是入射方向上的波矢分量。对透射端口来说还要把下表面到出射端口的距离也算进去。如果不扣除这个传播相位相位随入射角的变化里会掺入一个 (k_z d) 的项这个额外项对角度求导之后会给GH位移带来一个系统偏差。扣除参考面具体有两种做法手动扣除从S11和S21的相位里减去 (k_z d_1) 和 (k_z d_2)分别对应入射侧和出射侧的距离再把折射率差异带来的透射相位修正掉在COMSOL里设定参考面位置在端口设置里直接把参考面挪到结构表面位置这样S参数相位就是结构表面的相位。第二种做法看起来简单但挪参考面后COMSOL计算S参数时会在端口面上重新投影场如果结构表面附近网格不够精细反而会产生额外误差。我个人的习惯是用第一种方法手动扣除虽然麻烦一点但每一步都清楚出了问题也容易定位。具体计算时注意斜入射下的传播相位要用 (k_z k_0\sqrt{1-\sin^2\theta_i})不要直接用 (k_0 d)。入射角大了以后(k_z) 逐步减小这个修正量最多能差到百分之几十影响相当可观。3.3 透射相位还得多算一步相位参考基准的统一反射相位本身就有一个自然的参考基准入射面位置。只要规定入射参考面在结构上表面反射相位就是S11修正后的相位。但透射相位要注意一个问题出射介质和入射介质的折射率可能不同透射光的横向波矢 (k_x) 和纵向波矢 (k_z) 与入射侧不一致。此时透射GH位移的公式要改成[ \Delta_t -\frac{d\phi_t}{dk_{x,t}} ]而 (k_{x,t} k_x) 在界面两侧是守恒的连续性的结果所以仍然是 (-\frac{d\phi_t}{dk_x})。这个看似简单的结论在实际的S参数提取里却容易埋坑如果你的透射参考面位置是用“几何距离”扣除的而没有用出射介质的 (k_z) 来算传播相位那相位对 (k_x) 的导数就会被错误平移一个量。我的习惯是不管是反射还是透射统一先用 (k_x) 作为扫描变量保证两侧波矢分量一致再用 (k_z \sqrt{k_0^2 - k_x^2}) 扣除各自的传播项。做完之后把S11相位和S21相位分别对 (k_x) 画出来检查连续性确认没有跳变再继续。3.4 一个实测例子无结构平板的相位基线为了验证流程我建议每位同学先做一个“无结构平板”的基线模型就是均匀介质平板不加任何柱体。理论上这个模型不产生GH位移透射/反射相位对 (k_x) 的导数应该严格对应简单的光学路径差。我自己跑基线模型时用了一个折射率 (n1.5)、厚度 (300\text{ nm}) 的平板波长 (800\text{ nm})从正入射扫到60度。结果是扣除参考面后反射相位随角度基本不变透射相位随角度变化接近 (k_zd) 的贡献这一段恰好对应平板内的传播路径。这个基线一旦跑通说明参考面扣除、端口设置、Floquet边界这几项都正确了后续加结构才有意义。如果你跑基线就发现相位曲线抖动厉害不必急着加结构回头检查网格和端口模式数量准没错。4. 由相位梯度推算GH位移数值处理技巧与误差控制4.1 直接求导 vs 拟合求导噪声控制的胜负手拿到连续的相位曲线之后下一步就是求它对 (k_x)或入射角的导数。这一步最直接的想法是中心差分[ \phi(k_x) \approx \frac{\phi(k_x\Delta k_x) - \phi(k_x-\Delta k_x)}{2\Delta k_x} ]但实际跑下来中心差分对相位噪声很敏感尤其是S21幅值接近0的波长或角度区间——相位噪声会爆炸差分结果简直没法看。我推荐的流程是“先拟合再求导”把相位曲线分成若干段每段用多项式三阶或五阶拟合或者用Savitzky-Golay滤波做局部平滑对拟合多项式解析求导再把导数代进 (\Delta -d\phi/dk_x)得到GH位移曲线。用拟合的好处是可以控制平滑程度和局部偏差在数据本身比较光滑时结果相当稳定。分段拟合时要注意相接处不能有明显跳变可以用重叠窗口取平滑过渡。在实现上Python的SciPy库里scipy.signal.savgol_filter非常方便import numpy as np from scipy.signal import savgol_filter # kx 单位1/mphase 单位rad # window_length 必须是奇数选择 15~25 效果较好 phase_smooth savgol_filter(phase, window_length15, polyorder3) dkx kx[1] - kx[0] dphase_dkx np.gradient(phase_smooth, kx) # GH位移 gh_shift -dphase_dkx # 单位msavgol_filter的窗口长度直接决定平滑力度窗口越大曲线越平滑但可能把真实的小结构细节也磨掉窗口太小噪声压不住。我的经验是先用窗口长度21、阶数3跑一遍再分别用15和25对比如果结果差异小于0.01个波长说明平滑参数选得合适如果差异大说明数据本身噪声严重需要回到网格/端口设置找原因。4.2 扫描范围与步长的选择原则GH位移对入射角的依赖不是单调的尤其在光子晶体带隙边缘附近位移可能出现尖锐的峰或负值。这要求扫描范围必须覆盖感兴趣的物理区间以外一些否则边界效应会影响曲线形状。我常用的做法是从正入射 (\theta0^\circ) 扫到 (\theta70^\circ)步长1度先探趋势在趋势变化剧烈的区域加密到0.2度或0.1度计算GH位移时只用加密后的均匀 (k_x) 网格数据。步长选择有一个参考相邻两点的相位差尽量控制在0.1~0.3 rad以内。如果步长太大unwrap容易出问题太小则计算量激增相位噪声未必能改善。这个准则比单纯“细扫描”更本质。另外扫描角度上限要避开掠入射接近90度此时网格需要急剧加密而且Floquet端口的高阶衍射模式可能开启S参数中多个模式混叠单模S21/H位移计算就失效了。一般做到70度已经足够看见带隙边缘的大部分物理行为。4.3 如何判断GH位移结果的合理性算完GH位移先别急着画图完事有几个基线检查无结构平板位移应该在0附近或等于简单的传播路径差如果偏离明显超过0.05个波长说明参考面扣除有bug正入射时由于对称性GH位移应该严格为0正入射不存在横向方向如果算出来不是0说明网格或边界条件有不对称性位移随波长的积分量或者说在带隙内的形状应该和反射相位奇异点对应如果位移峰位跟反射谷位差太远多半是相位处理出了问题。把这些检查做过一遍再下结论说“这个超表面产生了XX波长的GH位移”才能站得住。4.4 多模波动时的处理高阶衍射端口怎么取舍当斜入射角较大或周期较大时高阶Floquet衍射模式可能变成传播模式不再衰减。此时透射/反射能量会被分到多个衍射阶次单端口S参数不能代表全部物理。对于GH位移计算我建议做两步检查数值上确认扫描范围内只有0级衍射模式是传播的即 (k_0 \sin\theta \pm mG) 的横向波矢满足 (|k_x| k_0) 的只有 (m0)如果高阶模式确实开启了必须为每个传播衍射阶次单独设端口然后分别计算每个端口的相位和位移。总位移在数学上不是简单地按强度加权平均这里涉及多光束干涉的概念比较复杂。对大多数光子晶体超表面设计来说周期通常小于波长斜入射角也不太大所以高阶模式开启的风险不高。但如果你做的是大周期超表面一个周期里有好几个单元就要格外注意这个坑。我自己的习惯是任意一个新结构先把各个衍射阶次的功率输出扫描一遍看有没有明显拐点。一旦功率曲线出现连续跳变就得回头查高阶模。5. 常见问题与避坑实录5.1 相位在±π处反复跳变导致unwrap失效这是新手最容易碰到的问题。相位的原始数据在 ±π 边界来回跳unwrap之后依然不连续往往是因为扫描步长太大真实相位变化超过了π。我遇到过最典型的一次场景在光子晶体带隙边缘反射相位在极小的角度变化内从 (π) 猛跌到 (-π)这个角度间隔内实际相位变化远超π直接unwrap把一条真实的相位悬崖“拧”成了一条平缓曲线后续GH位移计算全部报废。解决办法是先用较密集的扫描步长把相位变化限制在每步小于π以内如果结构对角度过于敏感可以进一步细化。如果细化后数据量太大可以分区域扫描平缓区域用大角度步长陡峭区域用小角度步长各自unwrap后再拼接。5.2 参考面扣除不当引起的系统偏差之前提到端口到结构表面之间的空气层必须扣除传播相位。如果这一步做错GH位移会在全部角度范围内整体偏移一个常数项。我踩过的坑是因为上下两侧空气层厚度不一样我只记得扣一侧忘了扣透射侧。结果透射GH位移整体偏了约0.3个波长但曲线形状和趋势看起来没问题。这种“看起来合理”的错误最坑人因为不仔细对照基线模型很难发现。所以基线模型的价值再强调一遍任何新流程先跑无结构平板确认反射和透射的GH位移都接近0再上真实结构。5.3 网格不对称带来正入射位移不为零理论上正入射不产生横向位移但我第一次跑就发现正入射时GH位移有0.05个波长的残差。排查了很久最后发现是网格剖分不对称导致的。COMSOL默认的网格生成有时会根据几何位置产生非对称分布尤其当介质柱不在计算域正中心时。解决办法是手动在结构周围添加对称的网格控制比如在介质柱左右两侧各加一个同样尺寸的“网格尺寸区域”确保每个角度下结构的左右两侧网格密度一致。做完这个操作正入射的GH位移残差降到了10的负四次方波长量级才算可以接受。5.4 S21接近零时相位噪声放大透射谷带隙位置的S21幅值很小相位对噪声极其敏感。这时候从S21算出来的GH位移会发飘噪声值远超物理信号。一个有效的方法是对复数S21直接做平滑而不是对角度做平滑。因为复数平滑考虑了幅值和相位的联合变化可以有效抑制“幅值很小时的相位噪声”。具体做法是S21_smooth savgol_filter(S21, window_length15, polyorder3) # 对复数数组平滑 phase_smooth np.angle(S21_smooth)注意这里需要先对复数数组做平滑再取相位角。如果S21数组跨越了实部或虚部变号的那一段复数平滑比单对相位角平滑稳健得多。这也是很多做过雷达或通信信号处理的人会自然想到的做法但在光学超表面仿真里提的人不多。5.5 表格速查GH位移计算的常见问题定位现象可能原因检查与解决相位曲线在±π处乱跳扫描步长太大加密扫描或用复数平滑正入射GH位移不为0网格不对称或边界不对称检查网格密度分布手动对称控制全角度范围位移整体偏移参考面未扣除干净跑无结构平板基线验证S21相位噪声大幅值接近0对复数做savgol滤波再取相位位移曲线出现尖刺数值微分噪声放大先拟合/滤波再求导调整平滑窗口入射角较大时结果异常高阶衍射端口未处理检查各衍射阶次功率增设对应端口扫描结果不收敛单个角度网格重剖固定网格禁用自适应剖分这张表我打印出来贴在了显示器边上每次仿真出问题下意识先扫一眼多数情况能立刻定位方向。遇到实在难以判定的就从“无结构平板”基线开始一步步加复杂度这是最笨也最稳的排查方式。写在最后的一点体会把透反射相位计算和GH位移彻底跑通之后我对超表面模拟的理解有了一个很大的转变以前我只盯着反射率、透射率这些标量结果觉得只要幅值对了就万事大吉。现在更多关注相位以及相位对角度的响应因为很多有趣的现象——比如光束位移、负折射、波前调控——全都藏在这条相位曲线的斜率里。过程中最深的体会是COMSOL本身只是工具真正费时间的是对物理图像的理解和数值处理的细节。尤其对于GH位移这种由相位导数决定的量一点点参考面设置偏差、一次unwrap处理失误、网格剖分不对称都会让结果完全失真。但只要把无结构基线模型先跑通、把每一步的相位检查做扎实后面再复杂的超表面结构也能稳当地推进。最后分享一个日常小习惯做任何参数扫描前先导出一组“正入射S参数”存成独立文件以后想验证模型设置有没有被改动过直接跑正入射对比一下就行不用每次都重头核对边界条件和端口设置。这个习惯帮我省了非常多重复性排查的时间也希望能帮到你。