前几天有个师弟跑来问我说用 Comsol 算光栅复合波导的反射相位想看看准BIC能不能把古斯汉森位移拉起来。这个问题我刚好做过一轮当时把结构调到准BIC附近时反射光束的横向位移从几个波长一下跳到了几十个波长后来继续精调参数甚至到了百微米量级。这么明显的增强效应对传感和光束偏转来说很有价值但仿真里要算准、算稳并不简单尤其是相位提取和端口参考面这两步稍不留神结果就偏得离谱。这篇就以我实际跑过的复合波导光栅模型为例把物理机制、Comsol 建模流程、准BIC 的寻找方法、古斯汉森位移的计算步骤和常见坑位一次性讲清楚。适合正在做微纳光学仿真、想复现高 Q 共振增强现象的研究生和工程师也适合刚接触 Comsol 波光学模块、想找个完整案例练手的朋友。1. 先从物理上把增强机制捋清楚1.1 古斯汉森位移从一条相位梯度说起古斯汉森位移可以理解成一束光在界面发生全反射时反射光并不在几何反射点立刻弹走而是稍微“滑”进第二介质再弹出来所以反射光束的实际横向位置比几何位置偏移了一段距离。教科书里最常见的公式是[ \Delta -\frac{\lambda}{2\pi}\frac{d\phi_r}{d\theta} ]其中 (\phi_r) 是反射系数的相位(\theta) 是入射角。也就是说反射相位随角度变化越剧烈横向位移越大。很多人在算位移时只关注反射率峰值忘了相位才是真正的关键。共振结构之所以能增强 GH 位移本质上就是它在角度谱上产生了一个极陡的相位跳变。普通的介质界面反射相位变化很慢位移只有波长量级而高 Q 共振可以在非常窄的角度范围内把相位从 (-\pi) 快速卷到 (\pi)导数数值可以非常大位移自然就被放大了。所以要做增强型 GH 位移核心思路很清晰找一个角度附近的相位突变点同时保证该角度附近有较高的反射率否则光束会大量透射或吸收位移再大也没有实际信号可测。准BIC正好提供了这样一个低辐射损耗、高相位灵敏度的通道。1.2 BIC与准BIC为什么相位能突变连续谱中的束缚态英文叫 Bound States in the Continuum简称 BIC。简单说它就是一个被“困”在连续辐射通道里的模式理论上既满足波动方程又完全不会往远场辐射能量。这个听起来很反直觉就像一个在水流湍急的河道里既不随波逐流也不被冲走的漩涡看起来不可能但在特定对称性或参数点上是真实存在的。严格意义上的 BIC 因为辐射损耗为零反射谱或透射谱上看不到共振峰它不会和自由空间平面波直接耦合所以无法被外部光激发。但只要把结构做一点扰动破坏掉产生 BIC 的对称条件原本“隐身”的暗态就变成了可以耦合的亮态。这时候它不再严格无辐射而是变成一个辐射损耗很小的共振模式也就是准BIC。准BIC 最吸引人的地方是它的 Q 值非常高。共振的相位响应大致可以写成[ \phi(\omega) \approx \phi_{\text{背景}} \arg\left(\frac{1}{\omega-\omega_0i\gamma}\right) ]当辐射损耗 (\gamma) 很小的时候相位在共振频率附近的变化极其陡峭。对于角度谱也有类似的规律共振位置随入射角改变所以在准BIC 对应的工作点附近反射相位对角度的导数可以做到非常大。这就是 GH 位移增强的根本原因。1.3 复合波导光栅把准BIC“引”到自由空间光栅加波导的结构很常见本质上是用周期折射率调制把自由空间光耦合进波导模式。当满足相位匹配条件时入射波可以激发波导模式再通过同一个光栅耦合出来反射谱或透射谱上会出现一个窄的共振峰这叫导模共振。导模共振本身就能产生 Fano 线型Q 值可以做到几千甚至几万。复合波导在这里的价值是可以提供更多的模式调控自由度和更灵活的带隙位置。比如在高折射率波导层下方再叠一层低折射率介质层两个波导层之间的耦合会形成对称和反对称模式通过调整层厚可以构造不同横向模式分布从而使某些模式在特定动量点与自由空间通道解耦形成 BIC。再加一维光栅就可以方便地用平面波从上方或下方激励这个模式。我用的是“衬底/低折射率间隔层/高折射率波导层/光栅层/空气”的经典五层结构光栅周期控制在亚波长范围这样在目标入射角附近只有 0 级衍射传播其他衍射阶次都是倏逝波计算时只需要保留 0 阶反射和透射通道。这个设计能有效避免高阶衍射带来的能量分配问题也方便后续提取反射系数。2. 用Comsol搭模型的完整流程2.1 几何与材料先把参数定好我习惯在开始建模前先画一张层结构示意图把每一层材料和厚度都标出来然后确定要扫描哪个参数。以我跑过的近红外波段案例为例工作波长设在 1550 nm 附近这样材料折射率数据比较容易找。一套可复现的初始参数如下结构层材料厚度折射率覆盖层空气半无限1.0光栅层TiO₂ 或 SiN100 ~ 150 nm2.0 ~ 2.4高折射率波导层SiN150 ~ 250 nm2.0间隔层SiO₂200 ~ 400 nm1.45衬底SiO₂ or 玻璃半无限1.45光栅在 x 方向周期为 (P)脊宽为 (w)占空比 (fw/P)。我先取 (P800,\text{nm})、(f0.5) 作为初始设计点。如果做对称破缺引入准BIC可以从改变光栅两侧的空气隙宽度入手把占空比从 0.5 变成 0.45 或 0.55也可以上下两层光栅错开一个偏移量 (\Delta)。材料折射率如果有损耗最好设置成复折射率。不过前期为了方便找共振和准BIC我会先设成无吸收材料把损耗放在后边单独讨论。这样计算得到的 Q 值反映的是辐射损耗物理图像更干净。仿真层面任何折射率变化都可以用参数化扫描实现不必每次手动改几何。2.2 物理场、端口与周期性边界关键一步Comsol 里做这种周期结构推荐直接选用“电磁波频域”接口把模型空间设为二维x 方向为周期方向z 方向为法线方向。二维模型对应实际结构在 y 方向无限长的情况这也是光栅导模共振最标准的近似。左右两侧边界用周期性边界条件端口类型选择“周期性端口”这样可以自动处理 Floquet 周期的相位关系。如果入射光是斜入射需要设定端口上的波矢分量 (k_x k_0\sin\theta)。端口阶次只需要保留 0 阶因为亚波长光栅的高阶衍射都是倏逝模不会携带净能量。端口的上方和下方还要加 PML 吸波层避免反射波从截断边界传回来干扰结果。这里有个容易被忽略的细节端口参考面的位置。端口模式在离结构有一定距离的边界上定义所以算出来的反射相位包含了端口到结构之间这段传播路径的相位积累。后面在算 GH 位移时必须扣除否则会混入一个与角度相关的几何项。我在建模时习惯把端口边界面到结构表面的距离记成一个参数 (L_\text{ref})后处理时统一补偿。材料参数直接按常数折射率设定即可。如果要用更真实的色散模型可以从材料库导入但前期建议先用常数因为准BIC 调节和共振位置主要靠几何参数材料色散对相对变化的影响在窄光谱范围内不大。2.3 研究步骤与网格策略算得快还得算得准这个模型至少需要两类研究一是特征频率研究用于找模式的实频和虚部从而算 Q 值二是频域扫描用于扫反射谱 S11 或者固定角度扫波长。具体做法是先在特征频率研究中设置搜索区间比如在 190 THz 附近找模式然后对每种对称破缺参数计算模式频率和 Q 值判断是否处于准BIC 区。找到合适的结构后再切换到频域研究扫反射谱和相位。网格是这类仿真最容易翻车的地方。光栅层的拐角会带来很强的近场局域网格粗了会出现明显频偏和虚部膨胀。我通常先用最大单元尺寸 (\lambda/(6n_\text{max})) 做过渡然后对光栅层和波导层做局部细化最大单元尺寸压到 (\lambda/(10n_\text{max}))光栅拐角处再加角细化。PML 区域的网格用默认映射网格就行不用过度加密。求解器方面频域扫描我推荐用直接求解器 PARDISO虽然内存占用比迭代法高一些但稳定性好尤其是扫相位曲线的时候不会因为预条件器收敛不好而出现莫名其妙的尖刺。固定角度扫描频率时扫描步长要足够细。准BIC 模式线宽可能只有几 GHz对应到频率扫描至少要跨越线宽多倍并且在线宽内至少设置 10 ~ 20 个采样点不然后面做相位差分根本不够用。3. 准BIC的寻找与参数调控3.1 打破对称把暗态变成亮态严格的 BIC 形成通常需要某个对称操作下的模式不辐射比如一个沿 z 方向偶对称的模式在结构上下对称的条件下无法与平面波耦合。要让这个暗态变成准BIC最简单的做法是打破结构的对称性。一维光栅中常用的办法是让光栅脊不再居中也就是把占空比改成不对称。还有一个办法是双层光栅让上下两层脊之间加一个横向偏移 (\Delta)。这种偏移会把原本偶对称的模式变成奇偶混合模式从而获得净辐射。在 Comsol 里实现起来并不复杂就是把上层光栅的相对位置绑定到辅助参数 delta然后做参数扫描。从理论上讲远离严格 BIC 点之后辐射损耗大致随破缺参量的平方变化所以 Q 值的倒数 (1/Q) 通常和 (\Delta^2) 呈线性关系。这个规律可以反过来用来判断哪些模式是真正的 BIC 衍生模式如果你看到 Q 值随 (\Delta) 增大而快速下降并且 (1/Q) 对 (\Delta^2) 的曲线接近一条直线说明这个模式大概率是准BIC。3.2 从特征频率提取Q值准BIC 的 Q 值不需要直接看反射峰半高宽因为反射峰在背景相位叠加下不一定是对称洛伦兹型。更可靠的方法是用特征频率研究直接算出共振模式的复频率 (\omega \omega_r i\omega_i)。Comsol 的特征频率结果默认是频率实部虚部在结果表格里能看到符号约定各个版本略有差异。我一般取[ Q \frac{|\omega_r|}{2|\omega_i|} ]也就是用实频绝对值除以两倍的虚频绝对值先避开符号问题。如果模式虚部很小接近零说明这个模式离严格 BIC 很近。不过这时候特征频率求解器可能很难找到模式因为它和连续谱的辐射模混在一起。技巧是指在求解器设置中把搜索频率区间缩小或者先在大破缺参数下算出模式实频再慢慢缩小破缺参数用上一步结果当预热初值。还有一个需要注意的点特征频率算准BIC 时PML 的覆盖范围要足够大。如果 PML 太薄漏泄模式会被人工反射影响虚部会偏大。我通常会留至少两个中心波长的 PML 厚度缩放系数保持默认基本可以接受。3.3 参数扫描找GH位移最大的工作点仅仅 Q 值高还不够GH 位移要看反射相位对角度的变化率。所以找到准BIC 之后要扫的是“反射相位 vs 入射角”。具体流程是先固定频率在准BIC 共振频率附近然后做一个以入射角为参数的频域扫描提取反射系数 S11 的相位。这里的角度扫描范围并不需要很大通常只需要在共振角附近扫 1~2 度。但步长一定要细比如设置 0.01°甚至更小。因为准BIC 反射相位的变化集中在极窄的角度窗口内步长太大直接会把相位跳变平滑掉。我在工程上有个习惯先用较粗的步长比如 0.05°快速看反射率曲线确定共振角大概位置然后再用 0.002° 或 0.005° 的步长在那个角附近精细扫描。这样既能避免计算量浪费又能保证后处理的相位差分足够可信。扫描结果里如果反射率出现接近 100% 的峰同时相位在该角度附近有一个接近 (2\pi) 的连续翻转那么 GH 位移的峰值一定不小。4. 古斯汉森位移的计算与验证4.1 反射相位提取与解包裹Comsol 里提取反射系数相位有很多种方式。最简单的方法是在“全局计算”里选择端口 S 参数再取它的相位。但直接拿到的相位可能被包裹在 ([-180^\circ,180^\circ])在共振附近会出现从 (180^\circ) 跳到 (-180^\circ) 的假跳变必须做相位解包裹。Comsol 的“相位”函数本身不会自动解包裹需要手动处理。我的做法是把 S11 的实部和虚部分别导出来再用 Python 脚本处理。Comsol 可以很方便地在结果中生成二维数据表把入射角、S11 实部、S11 虚部导出为文本或 CSV然后在 Python 里用numpy.angle加numpy.unwrap得到连续相位谱。这比在 Comsol 里面搞复杂后处理要省事得多也更容易排查问题。这里还要特别注意端口参考面补偿。反射的电磁波在端口到结构表面之间传播时积累了一个额外的传播相位。对于厚度为 (d) 的参考空气层这个附加相位是 (2k_0 d\cos\theta)。如果角度变化这个附加相位也会变化它的导数贡献恰好是 (2d\sin\theta)。假如不扣除最终算出来的 GH 位移里就会多一项和端口距离线性相关的伪位移。我第一次跑这个模型时就吃了这个亏结构参数怎么调都不对后来把参考面补偿项加上之后一切都对上了。4.2 数值差分算角导数得到解包裹后的反射相位 (\phi_r(\theta)) 之后GH 位移就是相位对角度的导数再乘上负波长因子。数值差分我习惯用中心差分[ \frac{d\phi_r}{d\theta}\bigg|_{\theta_0} \approx \frac{\phi_r(\theta_0\Delta\theta)-\phi_r(\theta_0-\Delta\theta)}{2\Delta\theta} ]角度步长 (\Delta\theta) 的选择要谨慎。步长太大共振峰附近的导数会被明显低估步长太小S11 数值噪声会主导导数出现毛刺。一个比较实用的经验是让角度步长小于共振线宽的三分之一并且保证在共振角附近连续扫描时相位变化单调。如果差分后曲线仍然有抖振可以用窗宽适中的 Savitzky-Golay 滤波先平滑相位再求导。差分结果乘上 (-\lambda/(2\pi))就能得到每一角度下的 GH 位移。在准BIC 共振角附近我会看到一条极为尖锐的峰值峰值半高宽通常远小于 0.1°。这个尖锐峰意味着实验上要做高精度角度控制才能测到否则一束光的角度抖动就能让位移掉下来。4.3 高斯光束模拟做交叉验证相位导数法算的是无限平面波极限下的 GH 位移。实际实验中入射光总是有限束宽的高斯光所以我不建议只看相位导数法最好在模型快要定稿时再做一次高斯光束模拟来交叉验证。严格来说周期边界条件默认结构无限大没办法直接模拟有限大小的光束。一个可行的折中办法是取一个足够大的超胞结构四周用 PML 包围然后在超胞顶部用散射边界条件加上一个角度偏转的高斯背景场。这种模型的计算量比原周期模型大不少但能够直观看到反射光束在空间上的横向偏移也可以检验准BIC 的角谱响应会不会把高斯光束的空间轮廓撕裂。我跑过的结果显示当高斯光束束腰半径大于几百个波长时数值模拟得到的 GH 位移和相位导数法的结果差异在 5% 以内。如果束腰太窄共振引起的角谱滤波效应会让反射光束失真甚至出现多峰这时候直接套用平面波相位导数公式就会高估位移。所以如果需要设计实验束腰半径尽量留大一些。5. 仿真中踩过的坑和排查方法5.1 共振峰频偏、Q值不收敛这类问题十有八九出在网格上。准BIC 模式的近场往往集中在光栅拐角或波导层边界那里的电场变化非常剧烈。如果网格不够密相当于人为引入了一个变形的结构共振频率自然偏移虚部也会变大Q 值被压低。排查时可以先做一个网格收敛测试把最大单元尺寸从 (\lambda/(6n)) 细化到 (\lambda/(12n))看共振频率和 Q 值的变化是否在 1% 以内。如果变化很大继续加密直到稳定。还有一个常见原因是 PML 厚度不够。特征频率计算中 PML 的作用是把向外辐射的能量吸收掉如果它太薄模式会感受到一个虚假的反射边界这就等于额外引入损耗通道虚部偏大。应对方法很简单加厚 PML通常两个波长以上基本够了。5.2 相位曲线跳变和端口参考面问题如果后处理得到的相位谱总是不连续即使做了 unwrap 还是看起来很奇怪优先检查角度扫描步长。准BIC 的相位变化窗口窄步长稍微大一点就会漏掉中间连续区域unwarp 后就会在共振处留下一个不太自然的台阶。把角度间隔缩小到共振线宽的三分之一以内这个问题基本会消失。在提取 S11 时我还遇到过端口定义了多个高阶模式的情况。如果端口阶次设置不当S11 对应的可能是多个衍射模式的叠加相位自然乱成一团。亚波长光栅下应该只有 0 阶衍射传播可以在端口属性里检查所有衍射阶次的传播常数把虚数或负数的阶次排除掉只保留 0 阶。另外一个容易被忽视的是入射角符号和 Floquet 周期边界条件的一致性。在斜入射设定中端口波矢和周期性边界的相位差必须设置成同一个值。如果用“自动”默认有时会因为方向定义不一致导致共振角偏移看起来像相位曲线整体平移了。排查方法很简单先在无光栅的平板结构上做一次相同设置看镜面反射相位是否符合解析结果。5.3 内存、收敛与版本兼容性参数化扫描量大时Comsol 的内存占用上升很快尤其是频域研究配合细网格直接求解器可能吃掉几十 GB 内存。如果不想换机型可以分段跑每次只扫一个参数的局部区间把结果存成独立文件再用脚本汇总。在 Linux 服务器上跑也是常见做法用命令行参数即可和 Windows 版相比文件格式完全兼容。新版 Comsol 6.4 在周期性端口和参数化扫描的底层存储上做了一些优化对于这种大型扫描任务体验会比旧版本流畅不少。还有一个小坑是特征频率求解器可能搜到不需要的 PML 模式。PML 区域本身也有伪模式解决方法是把求解器的搜索区域限定在目标频率附近很小的范围内并在结果中通过电场分布判断模式是否主要集中在波导和光栅区域。如果电场最大点出现在 PML 内部那是伪模式直接丢弃即可。6. 从准BIC到工程应用我的一些实用心得这个仿真链路跑完后最大的收获是准BIC 增强 GH 位移并不是简单地提高 Q 值就行。Q 值高意味着共振线宽窄相位变化确实快但可用角度窗口也窄。如果把它用到位移传感上需要在峰值灵敏度和工作动态范围之间做取舍。比如 Q 值太高环境温度和入射角稍微漂移输出信号就会剧烈波动反而对检测系统稳定性要求苛刻。还有一个心得是复合波导层的参数调控能力比单层波导强很多。单层波导光栅也能做出导模共振但模式频率和辐射通道之间的解耦往往只依赖一个几何参数复合波导多了一个间隔层厚度相当于多了第九根琴弦你可以通过调节两层模式之间的杂化来把准BIC 的残余辐射损耗压得更低。我做参数扫描时发现间隔层厚度从 300 nm 调到 350 nmQ 值能变化一个数量级以上这是单层结构很难实现的。后处理流程建议尽早脚本化。我现在的标准操作是把 Comsol 作为求解核心几何参数和扫描范围通过模型方法或命令行传进去计算完成后直接导出一组 CSV再用 Pyhton 做 unwrap、补相位、差分和画图。这个流程一旦跑通后期替换材料和波长就非常快。刚开始把所有后处理都放在 Comsol 界面的结果节点里每改一个参数就要重新设置一次效率很低。如果你也想快速复现这一类结果我个人的建议是先跑一个完全对称的结构把光栅周期调宽一些尽量在可见光或近红外波段找一个明显的导模共振峰并用自己的代码算出 GH 位移。这个基础流程理顺之后再加一层间隔层、引入对称破缺参数慢慢逼近准BIC。等 S11 相位曲线在共振角附近出现连续大斜率变化时你会非常直观地理解为什么说准BIC 是相位调控的利器。最后再补一句不算题外话的题外话仿真中算出来的百微米量级位移实验上不一定能直接看到因为普通反射镜的光斑移动通常只有微米甚至纳米级。可如果你把这种结构集成到芯片波导端口或者用近场探针在像平面上扫描结果就会非常明显。我自己的下一步计划是把电磁场功率流和时间平均坡印廷矢量导出来看一看位移增强时能量到底走了哪条路径这也能反过来验证相位导数法的可靠性。
基于Comsol的准BIC增强古斯汉森位移仿真：复合波导光栅建模与计算
前几天有个师弟跑来问我说用 Comsol 算光栅复合波导的反射相位想看看准BIC能不能把古斯汉森位移拉起来。这个问题我刚好做过一轮当时把结构调到准BIC附近时反射光束的横向位移从几个波长一下跳到了几十个波长后来继续精调参数甚至到了百微米量级。这么明显的增强效应对传感和光束偏转来说很有价值但仿真里要算准、算稳并不简单尤其是相位提取和端口参考面这两步稍不留神结果就偏得离谱。这篇就以我实际跑过的复合波导光栅模型为例把物理机制、Comsol 建模流程、准BIC 的寻找方法、古斯汉森位移的计算步骤和常见坑位一次性讲清楚。适合正在做微纳光学仿真、想复现高 Q 共振增强现象的研究生和工程师也适合刚接触 Comsol 波光学模块、想找个完整案例练手的朋友。1. 先从物理上把增强机制捋清楚1.1 古斯汉森位移从一条相位梯度说起古斯汉森位移可以理解成一束光在界面发生全反射时反射光并不在几何反射点立刻弹走而是稍微“滑”进第二介质再弹出来所以反射光束的实际横向位置比几何位置偏移了一段距离。教科书里最常见的公式是[ \Delta -\frac{\lambda}{2\pi}\frac{d\phi_r}{d\theta} ]其中 (\phi_r) 是反射系数的相位(\theta) 是入射角。也就是说反射相位随角度变化越剧烈横向位移越大。很多人在算位移时只关注反射率峰值忘了相位才是真正的关键。共振结构之所以能增强 GH 位移本质上就是它在角度谱上产生了一个极陡的相位跳变。普通的介质界面反射相位变化很慢位移只有波长量级而高 Q 共振可以在非常窄的角度范围内把相位从 (-\pi) 快速卷到 (\pi)导数数值可以非常大位移自然就被放大了。所以要做增强型 GH 位移核心思路很清晰找一个角度附近的相位突变点同时保证该角度附近有较高的反射率否则光束会大量透射或吸收位移再大也没有实际信号可测。准BIC正好提供了这样一个低辐射损耗、高相位灵敏度的通道。1.2 BIC与准BIC为什么相位能突变连续谱中的束缚态英文叫 Bound States in the Continuum简称 BIC。简单说它就是一个被“困”在连续辐射通道里的模式理论上既满足波动方程又完全不会往远场辐射能量。这个听起来很反直觉就像一个在水流湍急的河道里既不随波逐流也不被冲走的漩涡看起来不可能但在特定对称性或参数点上是真实存在的。严格意义上的 BIC 因为辐射损耗为零反射谱或透射谱上看不到共振峰它不会和自由空间平面波直接耦合所以无法被外部光激发。但只要把结构做一点扰动破坏掉产生 BIC 的对称条件原本“隐身”的暗态就变成了可以耦合的亮态。这时候它不再严格无辐射而是变成一个辐射损耗很小的共振模式也就是准BIC。准BIC 最吸引人的地方是它的 Q 值非常高。共振的相位响应大致可以写成[ \phi(\omega) \approx \phi_{\text{背景}} \arg\left(\frac{1}{\omega-\omega_0i\gamma}\right) ]当辐射损耗 (\gamma) 很小的时候相位在共振频率附近的变化极其陡峭。对于角度谱也有类似的规律共振位置随入射角改变所以在准BIC 对应的工作点附近反射相位对角度的导数可以做到非常大。这就是 GH 位移增强的根本原因。1.3 复合波导光栅把准BIC“引”到自由空间光栅加波导的结构很常见本质上是用周期折射率调制把自由空间光耦合进波导模式。当满足相位匹配条件时入射波可以激发波导模式再通过同一个光栅耦合出来反射谱或透射谱上会出现一个窄的共振峰这叫导模共振。导模共振本身就能产生 Fano 线型Q 值可以做到几千甚至几万。复合波导在这里的价值是可以提供更多的模式调控自由度和更灵活的带隙位置。比如在高折射率波导层下方再叠一层低折射率介质层两个波导层之间的耦合会形成对称和反对称模式通过调整层厚可以构造不同横向模式分布从而使某些模式在特定动量点与自由空间通道解耦形成 BIC。再加一维光栅就可以方便地用平面波从上方或下方激励这个模式。我用的是“衬底/低折射率间隔层/高折射率波导层/光栅层/空气”的经典五层结构光栅周期控制在亚波长范围这样在目标入射角附近只有 0 级衍射传播其他衍射阶次都是倏逝波计算时只需要保留 0 阶反射和透射通道。这个设计能有效避免高阶衍射带来的能量分配问题也方便后续提取反射系数。2. 用Comsol搭模型的完整流程2.1 几何与材料先把参数定好我习惯在开始建模前先画一张层结构示意图把每一层材料和厚度都标出来然后确定要扫描哪个参数。以我跑过的近红外波段案例为例工作波长设在 1550 nm 附近这样材料折射率数据比较容易找。一套可复现的初始参数如下结构层材料厚度折射率覆盖层空气半无限1.0光栅层TiO₂ 或 SiN100 ~ 150 nm2.0 ~ 2.4高折射率波导层SiN150 ~ 250 nm2.0间隔层SiO₂200 ~ 400 nm1.45衬底SiO₂ or 玻璃半无限1.45光栅在 x 方向周期为 (P)脊宽为 (w)占空比 (fw/P)。我先取 (P800,\text{nm})、(f0.5) 作为初始设计点。如果做对称破缺引入准BIC可以从改变光栅两侧的空气隙宽度入手把占空比从 0.5 变成 0.45 或 0.55也可以上下两层光栅错开一个偏移量 (\Delta)。材料折射率如果有损耗最好设置成复折射率。不过前期为了方便找共振和准BIC我会先设成无吸收材料把损耗放在后边单独讨论。这样计算得到的 Q 值反映的是辐射损耗物理图像更干净。仿真层面任何折射率变化都可以用参数化扫描实现不必每次手动改几何。2.2 物理场、端口与周期性边界关键一步Comsol 里做这种周期结构推荐直接选用“电磁波频域”接口把模型空间设为二维x 方向为周期方向z 方向为法线方向。二维模型对应实际结构在 y 方向无限长的情况这也是光栅导模共振最标准的近似。左右两侧边界用周期性边界条件端口类型选择“周期性端口”这样可以自动处理 Floquet 周期的相位关系。如果入射光是斜入射需要设定端口上的波矢分量 (k_x k_0\sin\theta)。端口阶次只需要保留 0 阶因为亚波长光栅的高阶衍射都是倏逝模不会携带净能量。端口的上方和下方还要加 PML 吸波层避免反射波从截断边界传回来干扰结果。这里有个容易被忽略的细节端口参考面的位置。端口模式在离结构有一定距离的边界上定义所以算出来的反射相位包含了端口到结构之间这段传播路径的相位积累。后面在算 GH 位移时必须扣除否则会混入一个与角度相关的几何项。我在建模时习惯把端口边界面到结构表面的距离记成一个参数 (L_\text{ref})后处理时统一补偿。材料参数直接按常数折射率设定即可。如果要用更真实的色散模型可以从材料库导入但前期建议先用常数因为准BIC 调节和共振位置主要靠几何参数材料色散对相对变化的影响在窄光谱范围内不大。2.3 研究步骤与网格策略算得快还得算得准这个模型至少需要两类研究一是特征频率研究用于找模式的实频和虚部从而算 Q 值二是频域扫描用于扫反射谱 S11 或者固定角度扫波长。具体做法是先在特征频率研究中设置搜索区间比如在 190 THz 附近找模式然后对每种对称破缺参数计算模式频率和 Q 值判断是否处于准BIC 区。找到合适的结构后再切换到频域研究扫反射谱和相位。网格是这类仿真最容易翻车的地方。光栅层的拐角会带来很强的近场局域网格粗了会出现明显频偏和虚部膨胀。我通常先用最大单元尺寸 (\lambda/(6n_\text{max})) 做过渡然后对光栅层和波导层做局部细化最大单元尺寸压到 (\lambda/(10n_\text{max}))光栅拐角处再加角细化。PML 区域的网格用默认映射网格就行不用过度加密。求解器方面频域扫描我推荐用直接求解器 PARDISO虽然内存占用比迭代法高一些但稳定性好尤其是扫相位曲线的时候不会因为预条件器收敛不好而出现莫名其妙的尖刺。固定角度扫描频率时扫描步长要足够细。准BIC 模式线宽可能只有几 GHz对应到频率扫描至少要跨越线宽多倍并且在线宽内至少设置 10 ~ 20 个采样点不然后面做相位差分根本不够用。3. 准BIC的寻找与参数调控3.1 打破对称把暗态变成亮态严格的 BIC 形成通常需要某个对称操作下的模式不辐射比如一个沿 z 方向偶对称的模式在结构上下对称的条件下无法与平面波耦合。要让这个暗态变成准BIC最简单的做法是打破结构的对称性。一维光栅中常用的办法是让光栅脊不再居中也就是把占空比改成不对称。还有一个办法是双层光栅让上下两层脊之间加一个横向偏移 (\Delta)。这种偏移会把原本偶对称的模式变成奇偶混合模式从而获得净辐射。在 Comsol 里实现起来并不复杂就是把上层光栅的相对位置绑定到辅助参数 delta然后做参数扫描。从理论上讲远离严格 BIC 点之后辐射损耗大致随破缺参量的平方变化所以 Q 值的倒数 (1/Q) 通常和 (\Delta^2) 呈线性关系。这个规律可以反过来用来判断哪些模式是真正的 BIC 衍生模式如果你看到 Q 值随 (\Delta) 增大而快速下降并且 (1/Q) 对 (\Delta^2) 的曲线接近一条直线说明这个模式大概率是准BIC。3.2 从特征频率提取Q值准BIC 的 Q 值不需要直接看反射峰半高宽因为反射峰在背景相位叠加下不一定是对称洛伦兹型。更可靠的方法是用特征频率研究直接算出共振模式的复频率 (\omega \omega_r i\omega_i)。Comsol 的特征频率结果默认是频率实部虚部在结果表格里能看到符号约定各个版本略有差异。我一般取[ Q \frac{|\omega_r|}{2|\omega_i|} ]也就是用实频绝对值除以两倍的虚频绝对值先避开符号问题。如果模式虚部很小接近零说明这个模式离严格 BIC 很近。不过这时候特征频率求解器可能很难找到模式因为它和连续谱的辐射模混在一起。技巧是指在求解器设置中把搜索频率区间缩小或者先在大破缺参数下算出模式实频再慢慢缩小破缺参数用上一步结果当预热初值。还有一个需要注意的点特征频率算准BIC 时PML 的覆盖范围要足够大。如果 PML 太薄漏泄模式会被人工反射影响虚部会偏大。我通常会留至少两个中心波长的 PML 厚度缩放系数保持默认基本可以接受。3.3 参数扫描找GH位移最大的工作点仅仅 Q 值高还不够GH 位移要看反射相位对角度的变化率。所以找到准BIC 之后要扫的是“反射相位 vs 入射角”。具体流程是先固定频率在准BIC 共振频率附近然后做一个以入射角为参数的频域扫描提取反射系数 S11 的相位。这里的角度扫描范围并不需要很大通常只需要在共振角附近扫 1~2 度。但步长一定要细比如设置 0.01°甚至更小。因为准BIC 反射相位的变化集中在极窄的角度窗口内步长太大直接会把相位跳变平滑掉。我在工程上有个习惯先用较粗的步长比如 0.05°快速看反射率曲线确定共振角大概位置然后再用 0.002° 或 0.005° 的步长在那个角附近精细扫描。这样既能避免计算量浪费又能保证后处理的相位差分足够可信。扫描结果里如果反射率出现接近 100% 的峰同时相位在该角度附近有一个接近 (2\pi) 的连续翻转那么 GH 位移的峰值一定不小。4. 古斯汉森位移的计算与验证4.1 反射相位提取与解包裹Comsol 里提取反射系数相位有很多种方式。最简单的方法是在“全局计算”里选择端口 S 参数再取它的相位。但直接拿到的相位可能被包裹在 ([-180^\circ,180^\circ])在共振附近会出现从 (180^\circ) 跳到 (-180^\circ) 的假跳变必须做相位解包裹。Comsol 的“相位”函数本身不会自动解包裹需要手动处理。我的做法是把 S11 的实部和虚部分别导出来再用 Python 脚本处理。Comsol 可以很方便地在结果中生成二维数据表把入射角、S11 实部、S11 虚部导出为文本或 CSV然后在 Python 里用numpy.angle加numpy.unwrap得到连续相位谱。这比在 Comsol 里面搞复杂后处理要省事得多也更容易排查问题。这里还要特别注意端口参考面补偿。反射的电磁波在端口到结构表面之间传播时积累了一个额外的传播相位。对于厚度为 (d) 的参考空气层这个附加相位是 (2k_0 d\cos\theta)。如果角度变化这个附加相位也会变化它的导数贡献恰好是 (2d\sin\theta)。假如不扣除最终算出来的 GH 位移里就会多一项和端口距离线性相关的伪位移。我第一次跑这个模型时就吃了这个亏结构参数怎么调都不对后来把参考面补偿项加上之后一切都对上了。4.2 数值差分算角导数得到解包裹后的反射相位 (\phi_r(\theta)) 之后GH 位移就是相位对角度的导数再乘上负波长因子。数值差分我习惯用中心差分[ \frac{d\phi_r}{d\theta}\bigg|_{\theta_0} \approx \frac{\phi_r(\theta_0\Delta\theta)-\phi_r(\theta_0-\Delta\theta)}{2\Delta\theta} ]角度步长 (\Delta\theta) 的选择要谨慎。步长太大共振峰附近的导数会被明显低估步长太小S11 数值噪声会主导导数出现毛刺。一个比较实用的经验是让角度步长小于共振线宽的三分之一并且保证在共振角附近连续扫描时相位变化单调。如果差分后曲线仍然有抖振可以用窗宽适中的 Savitzky-Golay 滤波先平滑相位再求导。差分结果乘上 (-\lambda/(2\pi))就能得到每一角度下的 GH 位移。在准BIC 共振角附近我会看到一条极为尖锐的峰值峰值半高宽通常远小于 0.1°。这个尖锐峰意味着实验上要做高精度角度控制才能测到否则一束光的角度抖动就能让位移掉下来。4.3 高斯光束模拟做交叉验证相位导数法算的是无限平面波极限下的 GH 位移。实际实验中入射光总是有限束宽的高斯光所以我不建议只看相位导数法最好在模型快要定稿时再做一次高斯光束模拟来交叉验证。严格来说周期边界条件默认结构无限大没办法直接模拟有限大小的光束。一个可行的折中办法是取一个足够大的超胞结构四周用 PML 包围然后在超胞顶部用散射边界条件加上一个角度偏转的高斯背景场。这种模型的计算量比原周期模型大不少但能够直观看到反射光束在空间上的横向偏移也可以检验准BIC 的角谱响应会不会把高斯光束的空间轮廓撕裂。我跑过的结果显示当高斯光束束腰半径大于几百个波长时数值模拟得到的 GH 位移和相位导数法的结果差异在 5% 以内。如果束腰太窄共振引起的角谱滤波效应会让反射光束失真甚至出现多峰这时候直接套用平面波相位导数公式就会高估位移。所以如果需要设计实验束腰半径尽量留大一些。5. 仿真中踩过的坑和排查方法5.1 共振峰频偏、Q值不收敛这类问题十有八九出在网格上。准BIC 模式的近场往往集中在光栅拐角或波导层边界那里的电场变化非常剧烈。如果网格不够密相当于人为引入了一个变形的结构共振频率自然偏移虚部也会变大Q 值被压低。排查时可以先做一个网格收敛测试把最大单元尺寸从 (\lambda/(6n)) 细化到 (\lambda/(12n))看共振频率和 Q 值的变化是否在 1% 以内。如果变化很大继续加密直到稳定。还有一个常见原因是 PML 厚度不够。特征频率计算中 PML 的作用是把向外辐射的能量吸收掉如果它太薄模式会感受到一个虚假的反射边界这就等于额外引入损耗通道虚部偏大。应对方法很简单加厚 PML通常两个波长以上基本够了。5.2 相位曲线跳变和端口参考面问题如果后处理得到的相位谱总是不连续即使做了 unwrap 还是看起来很奇怪优先检查角度扫描步长。准BIC 的相位变化窗口窄步长稍微大一点就会漏掉中间连续区域unwarp 后就会在共振处留下一个不太自然的台阶。把角度间隔缩小到共振线宽的三分之一以内这个问题基本会消失。在提取 S11 时我还遇到过端口定义了多个高阶模式的情况。如果端口阶次设置不当S11 对应的可能是多个衍射模式的叠加相位自然乱成一团。亚波长光栅下应该只有 0 阶衍射传播可以在端口属性里检查所有衍射阶次的传播常数把虚数或负数的阶次排除掉只保留 0 阶。另外一个容易被忽视的是入射角符号和 Floquet 周期边界条件的一致性。在斜入射设定中端口波矢和周期性边界的相位差必须设置成同一个值。如果用“自动”默认有时会因为方向定义不一致导致共振角偏移看起来像相位曲线整体平移了。排查方法很简单先在无光栅的平板结构上做一次相同设置看镜面反射相位是否符合解析结果。5.3 内存、收敛与版本兼容性参数化扫描量大时Comsol 的内存占用上升很快尤其是频域研究配合细网格直接求解器可能吃掉几十 GB 内存。如果不想换机型可以分段跑每次只扫一个参数的局部区间把结果存成独立文件再用脚本汇总。在 Linux 服务器上跑也是常见做法用命令行参数即可和 Windows 版相比文件格式完全兼容。新版 Comsol 6.4 在周期性端口和参数化扫描的底层存储上做了一些优化对于这种大型扫描任务体验会比旧版本流畅不少。还有一个小坑是特征频率求解器可能搜到不需要的 PML 模式。PML 区域本身也有伪模式解决方法是把求解器的搜索区域限定在目标频率附近很小的范围内并在结果中通过电场分布判断模式是否主要集中在波导和光栅区域。如果电场最大点出现在 PML 内部那是伪模式直接丢弃即可。6. 从准BIC到工程应用我的一些实用心得这个仿真链路跑完后最大的收获是准BIC 增强 GH 位移并不是简单地提高 Q 值就行。Q 值高意味着共振线宽窄相位变化确实快但可用角度窗口也窄。如果把它用到位移传感上需要在峰值灵敏度和工作动态范围之间做取舍。比如 Q 值太高环境温度和入射角稍微漂移输出信号就会剧烈波动反而对检测系统稳定性要求苛刻。还有一个心得是复合波导层的参数调控能力比单层波导强很多。单层波导光栅也能做出导模共振但模式频率和辐射通道之间的解耦往往只依赖一个几何参数复合波导多了一个间隔层厚度相当于多了第九根琴弦你可以通过调节两层模式之间的杂化来把准BIC 的残余辐射损耗压得更低。我做参数扫描时发现间隔层厚度从 300 nm 调到 350 nmQ 值能变化一个数量级以上这是单层结构很难实现的。后处理流程建议尽早脚本化。我现在的标准操作是把 Comsol 作为求解核心几何参数和扫描范围通过模型方法或命令行传进去计算完成后直接导出一组 CSV再用 Pyhton 做 unwrap、补相位、差分和画图。这个流程一旦跑通后期替换材料和波长就非常快。刚开始把所有后处理都放在 Comsol 界面的结果节点里每改一个参数就要重新设置一次效率很低。如果你也想快速复现这一类结果我个人的建议是先跑一个完全对称的结构把光栅周期调宽一些尽量在可见光或近红外波段找一个明显的导模共振峰并用自己的代码算出 GH 位移。这个基础流程理顺之后再加一层间隔层、引入对称破缺参数慢慢逼近准BIC。等 S11 相位曲线在共振角附近出现连续大斜率变化时你会非常直观地理解为什么说准BIC 是相位调控的利器。最后再补一句不算题外话的题外话仿真中算出来的百微米量级位移实验上不一定能直接看到因为普通反射镜的光斑移动通常只有微米甚至纳米级。可如果你把这种结构集成到芯片波导端口或者用近场探针在像平面上扫描结果就会非常明显。我自己的下一步计划是把电磁场功率流和时间平均坡印廷矢量导出来看一看位移增强时能量到底走了哪条路径这也能反过来验证相位导数法的可靠性。
