没用过 COMSOL 去碰周期性波导的人,可能很难理解我第一次在计算结果里看到本征频率虚部趋于极小数时的感觉:一条本应该带损耗的模式,在布里渊区中心点上,虚部直接掉到了浮点误差级别,Q 值趋近于无穷。这就是波导里的 BIC——连续谱束缚态。这类模式因为对称性保护而“锁”在波导里,既不向外辐射,也不和入射光耦合,在整个光子学设计里几乎是教科书级的理想状态。我这次就用 COMSOL Multiphysics 完整跑了一遍受陷波导结构的 BIC 求解,从几何搭建、模式扫参到 Q 因子追踪都有,把这套流程和踩过的坑一起整理出来,给正在做微纳光学仿真、特别是对周期性波导和准 BIC 感兴趣的朋友一个可直接参考的方案。
先说清楚这次要做什么:在一个硅基光栅波导结构里,用 COMSOL 的电磁波频域接口计算 Floquet 周期边界下的特征频率,扫描面内波矢 kx,找出色散关系中虚拟性消失的那个点。BIC 不是靠肉眼在模场图上“看”出来的,它的本质是模式在辐射连续谱内部却没有能量泄漏,所以在数值上最直接的判定标准就是特征频率的虚部接近零、Q 因子发散。这个过程听起来不难,真正做起来涉及几何对称性的设置、网格策略、求解器配置和结果后处理,每一个环节都会直接影响你最终能不能看到这个“无穷大 Q”。
1. BIC 是什么:为什么波导能锁住本该漏走的光
1.1 从波导模式说起——连续谱里为什么会有束缚态
普通波导里的束缚模式很好理解:光被限制在高折射率芯层里,通过全反射沿着波导传播,理论上不存在辐射损耗。但周期性波导不一样,光栅把原本只存在于导模内部的模式“衍射”出去,色散曲线被折叠,模式跟自由空间或者衬底里的辐射模式形成耦合,能量就会沿着泄漏通道跑掉。这是大部分周期性波导设计的默认情况:导模变成了真正意义上的“泄漏模”,Q 因子有限,器件表现为宽带辐射器。
但 BIC 恰恰是那个例外。它对应的模式在频率上落在辐射连续谱的范围内,也就是说它满足波矢匹配和能量守恒,按道理完全有资格向外辐射,但因为它跟辐射通道之间的耦合矩阵元恰好为零,能量被原先的对称性“保护”住了。你可以把它想成两个同频率的弹簧摆在同一个桌面上,其中一个的运动方式让桌子沿特定方向完全不动,另一个再使劲也推不动它——不是没力,是方向不匹配。BIC 就是那个方向不匹配的模式。
这个特性放在器件设计里太值钱了。如果能设计一个结构让某个波导模式工作在 BIC 附近,那就是在原本泄漏的结构里获得极高 Q,而结构本身还是开放辐射型,不需要像微环那样靠一圈圈绕来锁能量。更妙的是,只要轻微打破对称性,BIC 就退化成准 BIC,Q 虽然不是无穷大,但可以轻松做到 10⁴、10⁵ 甚至更高,同时还能人为控制耦合强度,做出可调谐的高 Q 谐振器。所以这个方向几乎所有工作都是同一个思路:先找到对称性保护的 BIC,再引入可控的扰动,把 Q 调到想要的量级。
1.2 对称性保护与“消失的耦合矩阵元”
BIC 的成因有很多种,比如调谐类、对称保护类、奇异点类等等。波导光栅里最经典、也最容易在仿真里稳定复现的就是对称性保护的 BIC。这类 BIC 出现在布里渊区的高对称点,比如 Γ 点或者 X 点,原因是结构本身具有某个反射或者旋转对称性,而 BIC 模式的辐射分量跟外部连续谱模式在对称性操作下呈现相反的变换行为,二者之间没有公共的辐射通道。
举一个最典型的例子:一个在 y 方向上完全对称的平板波导光栅,TE 模式具有偶对称性,电场分量在上下两个半平面里呈现相同相位;而可辐射的连续谱模式在这个对称面上是奇对称的,或者说具有相反的对称性。偶模和奇模之间没有矩阵元联系,偶模就不会衰减。COMSOL 里做模式分析的时候,如果你把模型的对称性保留得足够干净,算出来的本征频率虚部就应当是一个极小量——注意不是恰好为零,因为数值离散会让理想 BIC 变成一个很小的数值尾巴,但尾巴应该随着网格加密不断消失。
几何扰动破坏这个对称性之后,BIC 模式开始跟辐射通道有微弱交叠,虚部从零变成有限值,Q 因子从发散变成有限值。理论上 Q 跟扰动参数 δ 之间呈倒数平方关系,Q ~ δ⁻²,这个二次方标度也经常被用来验证你说的模式是不是真的源于对称保护 BIC。我在扫描光栅槽的横向偏移量时,把 Q 取对数画出来,斜率确实是 −2,跟预期完全吻合。这一步验证很重要,很多初学仿真的人算出了一个高 Q 模式就以为是 BIC,实际上那可能只是个低泄漏的普通模式,斜率检查能帮你排除大量假象。
1.3 在 COMSOL 里算 BIC 和常规光子晶体仿真有什么不同
做过光子晶体能带计算的人都知道,COMSOL 里最常见套路是画一个单胞,加 Floquet 周期边界,扫描 kx 得到色散关系。算 BIC 的流程看起来一样,但有几个关键差异必须重视。第一,BIC 对几何对称性极其敏感,你画图的时候哪怕只是在深宽比上差了 1 nm,模式也会从理想的无穷 Q 变成有限高 Q,所以几何必须参数化、精确取整。第二,普通能带计算关心的是实部频率随 kx 变化的色散关系,虚部通常是顺带看两眼,但 BIC 判定完全依赖虚部,虚部计算对网格和求解器容差要求都苛刻得多。第三,BIC 经常出现在连续谱内部,附近往往还挤着其他泄漏模式,特征值搜索范围设置不好很容易漏掉想要的那个值。
还有一点容易被忽略的是研究类型的设置。COMSOL 里算模场和 Q 因子最直接的方法是用“特征频率”研究结合 Floquet 周期条件,不是用频域扫描。频域扫描得到的是给定频率下的场分布,你得自己花 Q 因子或者拟合 Fano 线形才能间接推断 BIC;特征频率研究则直接把频率的实部和虚部都算出来,Q = Re(f)/(2|Im(f)|),一秒搞定。当然,后面如果要研究准 BIC 在入射光下的透射光谱,还得回到频域扫描,但那是第二步的事,第一步的 BIC 定位用特征频率就够了。
2. 建模前的决策:几何结构、材料参数和求解策略
2.1 选结构:为什么用“顶着浅光栅的埋氧硅波导”
BIC 要能在波导里出现,结构至少要满足两个条件:存在一个低损耗的导波层,同时导波层上方又有能提供辐射连续谱的周期结构。我这次用的是绝缘体上硅平台上做的浅光栅波导,最底层是硅衬底,中间是埋氧层 SiO₂,最上面一层是做光栅的硅薄膜。光栅槽不刻穿整个芯层,只在下层留一个连续平板部分,这样既保持了波导的导波能力,又引入了周期微扰。
选择这个结构有几个现实理由。硅在通信波段折射率 3.48,跟空气的折射率差很大,导波能力很强,模式限制得好,BIC 的 Q 值发散特性也会更稳定;埋氧层和衬底组成一个半开放系统,能够提供真实的辐射泄漏通道;从工艺角度讲,浅刻蚀光栅比全刻蚀要容易做,只是对后续实验验证友好的。当然你可以换成氮化硅、聚合物或者全介电双层光栅,物理本质一样,仿真流程不需要改。
几何参数我用一套近红外的典型配置:光栅周期 Λ = 1 μm,光栅占空比 f = 0.5,光栅槽深 d = 200 nm,硅芯厚度 h = 400 nm,埋氧层厚度 h_box = 1 μm,衬底厚度取 1 μm 就够了。注意这里衬底不能无限厚,你需要截断并加 PML 把它处理成开放式边界,别为了省事直接用默认的完美电导体边界,那会让泄漏模式被金属壁反射回来,虚部算出来完全是错的。PML 的问题后面单独说,这里先记住结论。
2.2 材料参数:折射率设不对,BIC 位置会漂
COMSOL 里材料模型的复杂度是个陷阱。算 BIC 这种本征模问题,一般情况下不需要用色散模型,直接在材料属性里填折射率常数就够了。硅在 1550 nm 用 3.48,SiO₂ 用 1.45,空气 1.0。如果你非要严谨地给硅加一个波长相关的折射率,其实也没问题,只是特征频率扫描迭代的时候每次都要重新计算材料插值,速度会慢不少,而且对 BIC 位置的判定基本没有实质影响。
真正需要注意的是材料损耗。默认情况下,COMSOL 材料库里的硅可能带了一个很小的虚部折射率,这会导致所有模式的虚部本征值整体偏大,Q 因子被“压”下来。一个 Q 值 10⁶ 的理想准 BIC 模式,如果材料损耗虚部设置不当,算出来可能只剩 10⁴,你甚至会误以为收敛有问题。对于纯数值探索 BIC 的阶段,强烈建议把硅和 SiO₂ 的损耗都设为 0,把问题干净的理想化,等确认几何和网格方案没问题以后再去加材料吸收。这样你可以清楚地区分什么损耗是几何泄漏引起的、什么损耗是材料吸收引起的,两种物理机制不应该混在一起。
2.3 求解策略:特征频率扫参数是主力,频域扫描做验证
COMSOL 里算 BIC 的完整策略可以拆成三步走。第一步,固定某个 kx,跑特征频率研究,找出虚部最小的几个候选模式;第二步,对 kx 做参数化扫描,追踪这个候选模式的实频和虚频随 kx 的变化,找到虚部掉到极小的那个点;第三步,在准 BIC 点附近用频域研究加端口激发或者偶极子源,看透射峰或者近场增强峰,验证特征频率法得到的 Q 值。
特征频率研究里,你要选择“电磁波,频域”接口,物理场设置为“模式分析”或者“特征频率”都可以,本质上是求解电场本征方程。扫描参数的时候把 kx 设成辅助扫描参数,从布里渊区边界一路扫到 Γ 点。因为 BIC 常出现在高对称点,你会在扫到边界附近时看到虚部开始快速下降,那个趋势非常漂亮,像是有一条模式准备“潜水”进入实轴。
频域验证的时候可以在波导上方设置一个入射端口,监测透射率或者反射率谱线,准 BIC 会在光谱里留下一个不对称的 Fano 线形。这里提醒一下,特征频率模式下算出来的 Q 是模式本身的本征 Q,频域谱拟合出来的 Q 还包含了入射耦合展宽,两者在弱耦合极限下应该接近但不完全相等,别拿两边的数字直接对不上就怀疑模型。
3. COMSOL 实操:从几何搭建到特征值提取的完整流程
3.1 参数化几何搭建:一个变量改遍全结构
进入 COMSOL 后先建一个二维模型,把求解域设成 y > 0 的上半空间,x 方向是一个周期单元,两侧用 Floquet 周期条件连起来。为了保持几何对称性,所有尺寸都通过全局参数定义,不要直接在几何里填死数字。全局参数表这样建:
| 参数 | 数值 | 物理意义 |
|---|---|---|
| Lambda | 1e-6 m | 光栅周期 |
| duty | 0.5 | 光栅占空比 |
| d_depth | 2e-7 m | 光栅刻蚀深度 |
| h_core | 4e-7 m | 硅芯层厚度 |
| h_box | 1e-6 m | 埋氧层厚度 |
| h_sub | 1e-6 m | 截断衬底厚度 |
| t_pml | 5e-7 m | PML 厚度 |
| n_si | 3.48 | 硅折射率 |
| n_ox | 1.45 | 二氧化硅折射率 |
| kx_param | 0 rad/m | 面内波矢分量 |
几何画法就是几个矩形拼接:从下往上依次是 PML 域、衬底、埋氧层、硅芯层,然后硅芯层顶部再加两个矩形做光栅的齿和槽。齿的宽度是 duty*Lambda,槽的宽度是 (1-duty)*Lambda,直接用参数表达式驱动。这样当你后面想扫占空比或者槽深研究扰动时,只需要把几何重建一遍,所有物理场自动跟着变,不需要重设边界条件。
一个容易被忽略的细节是 PML 域和衬底之间必须留出足够的间隔,一般至少 1 μm,否则 PML 会像一面镜子一样把倏逝场反射回波导区域,特征频率虚部被污染。同样,PML 的背边绝对不能设置什么特殊边界条件,保持默认的连续性,让吸收层工作在纯吸收模式。
3.2 物理场与边界条件的正确配置
物理场接口选择“电磁波,频域”,研究类型用“特征频率”。在“电磁波”节点下,把求解域里的所有材料都分配好,然后把上下边界条件配置清楚。核心操作有三个:第一个是给模型的左右两侧添加“周期性边界条件”,周期类型选 Floquet,波矢 k 的 x 分量设为 kx_param,y 分量设为 0。这个很重要,Floquet 边界是模拟无限周期结构的关键,如果没有它你就要建整个阵列,计算量直接膨胀几百倍。
第二个是给最高处的边界加“散射边界条件”或者默认开放边界,然后在它背后加上 PML 域。我一直强调这里别用完美电导体或者完美磁导体,因为那些条件会把辐射模式反射回计算域,只有当模型本身是一个闭合腔时它们才正确,开放波导结构用它们就是灾难。
第三个是根据你要算的模式类型设置恰当的“面内波数”搜索方向。TE 模式在二维 COMSOL 里表现为面外磁场分量 H_z 为主,TM 模式表现为面外电场 E_z 为主,这两个模式的对称性行为差别很大,BIC 的对称保护机制也不同。通常我们关注 TE 类模式,因为硅波导的 TE 基模限制更强,你也可以两个都算,反正求解器给的是一堆特征值,自己按场分布筛选即可。
还有一个小技巧:在使用特征频率研究前,在“特征频率”节点里把“所需模式数”预设为 8 到 12,然后把“搜索范围”设置在目标频率附近一个相对窄的区间,比如 180 THz 到 220 THz。这个做法能避免特征值求解器一股脑把所有高频模式都算出来,把计算时间浪费在表面上。不过搜索范围也不能太窄,因为随着 kx 扫描,模式的色散会移动,固定窄范围会漏掉目标模式。
3.3 网格划分:BIC 模拟成败的隐形决定因素
网格策略这方面我得说重点:BIC 数值模拟里,网格尺寸直接影响你能否观察到虚部消失。一个普通能带计算对网格不太敏感,网格粗一点,带的位置偏了 1%,也能接受;但 BIC 的虚部本质上是本应精确抵消的两个大数相减之后剩下的残差,网格一旦不够密,数值不对称性就会造出虚假的辐射损耗,模式看起来就不是 BIC 而是普通高 Q 模式。
我的建议是先在光栅层和波导芯层使用最大单元尺寸 Lambda/20 左右,测试粗网格,然后收敛到 Lambda/40、Lambda/80,观察虚部的变化趋势。真正的 BIC 应该在网格加密时虚部持续下降,向零逼近;如果某个高 Q 模式的虚部在网格加密后保持一个稳定值不再变化,那它就不是 BIC,只是结构里真实存在的一个低损耗模式。这个判据比什么后处理技巧都好用,可以说是唯一的金标准。
光栅转角和槽底处的网格需要额外细化,因为场在这些位置变化剧烈,粗糙的三角形会引入不真实的能量散射。我在这些区域设置尺寸为 Lambda/100 的局部细化,计算量还在可接受范围。Floquet 边界两侧的网格必须是一一对应的,否则周期性条件在离散层面出现错配,也会产生虚假泄漏。检查方法很简单:把网格显示打开,看左边界和右边界上的网格节点分布是否完全一致,如果不一致就调整边界细化方式。
3.4 求解器设置:让虚部算得足够干净
设置好上述一切后,点“研究”里的“特征频率”开始求解。求解器默认用 MUMPS 直接法,对这个规模的二维问题没什么压力,但特征值求解器在寻找模式时会对初始估计很敏感。我在计算中发现,即使设置了“所需模式数”,有时候虚部最小的那个 BIC 候选不会在第一次扫描中出现,因为它的辐射损耗太小、虚部跟其他模式差了好几个量级,数值上像被淹没了一样。
解决办法是在“特征频率”设置里调整搜索方法,从“区域”改为“点”附近,给定一个接近目标模式的起始频率。另一个办法是分两步:先在 kx = 0 处跑一个较宽的频率范围,找到大致频率;然后固定这个频率,缩小搜索区间重新求解。掌握这个技巧后,追踪模式变得很顺滑,特别是做 kx 参数扫描时,每一轮都用上一轮的解作为初始估计,能大幅减少漏解和模式跳跃。
关于特征值精度还有一个细节:COMSOL 默认的相对容差可能不够。在求解器配置里把“特征值容差”从默认的 1e-6 调到 1e-8 或者更小,虚部的数值噪声会明显降低,BIC 附近的判定会更干净。代价是每个特征点要多个几次迭代,但为了看到真实虚部,这点开销完全值得。
4. 结果解读:从本征虚部到 Q 因子标度律
4.1 模场分布与“对称性对不上”现象
算完特征频率后,先看模场再谈数字。选择虚部最小的那个模式,画出它的电场模分布。BIC 的场分布有一个标志性特征:能量高度集中在波导芯层里,上方空气区域几乎看不到辐射场,哪怕你把颜色范围放到很小,场也像是被隔离在芯层内。作为对比,旁边的普通泄漏模会在空气区域留下明显的传播波条纹,这种条纹就是辐射通道存在的直观证据。
然后点开“全局计算”,查看该模式的电场分量在反射对称面上的表现。以 y = 0 为镜面,BIC 模式的某个电场分量在镜像操作下要么是对称要么是反对称,但绝不会同时具备与辐射连续谱模式相同的对称性。数值上你怎么检查这一点?最实用的方式是对称和反对称分量做内积,如果内积极小,说明耦合为零,对应的就是 BIC。COMSOL 里可以直接用“集成”算子在不同域上计算,当然更直接的办法是肉眼判断场图加上前面说的网格收敛测试。
还有一个易混淆的点:不要把所有看起来像局域模式的都叫 BIC。在周期性结构中,有些模式在 Γ 点处因为处于带隙里、压根没有可耦合的辐射通道,它也成为束缚态,但这不是连续谱束缚态,它只是普通带隙模式。两者的区别要看它在布里渊区里的位置:带隙模式处于频率带隙内,旁边没有连续谱;BIC 位于连续谱内部,旁边有大量的辐射模式。判别方式很简单,看这个频率处是否存在泄漏模色散曲线经过,如果存在且模式仍是束缚的,那才是 BIC。
4.2 扫描 kx 追踪虚部:寻找那个“跳水点”
特征频率扫描完成后,我做了一张 kx 从 0.2×(2π/Λ) 到 0 的扫描曲线图,横坐标是归一化的 kx,纵坐标是 Im(ω)(或者直接画 Q)。从扫描图可以非常清楚看到,随着 kx 接近 Γ 点,模式的虚部先是缓慢下降,然后在靠近 kx = 0 时像坐滑梯一样急剧掉下去几个数量级,Q 因子一路上扬。这个跳水点就是 BIC 的数值指纹。
画图的时候建议纵坐标用对数尺度。因为虚部的动态范围可能从 1e-4 一直到 1e-10,用线性轴根本看不出下降趋势,对数轴才能把这个跨越七个量级的过程完整展示出来。同时把实频也画在同一张图上,确认频率在这个扫描范围内没有跨越其他模式或者出现带交叉,避免你追踪的模式在某个 kx 处突然“换轨”,跑到另一条模式曲线上。
关于扫描步长,在远离 Γ 点的区域用 Δkx = 0.02×(2π/Λ) 就够了,靠近 Γ 点时步长要缩小到 0.002×(2π/Λ),因为虚部的变化在 BIC 附近是非线性的,指数式下降,步长太大会漏掉最深的那个点。每次算完后把特征值按虚部排序,挑最小的记录,这就是当前 kx 下最接近 BIC 的候选模式。
4.3 扰动与准 BIC:Q 是否按 δ⁻² 标度
找到精确对称结构里的 BIC 后,下一步通常是研究它的准 BIC 行为,因为真实器件不可避免地有加工误差和形状偏移,模拟纯粹理想的 BIC 用处有限。我在光栅几何上引入一个扰动参数 δ,让原本左右完全对称的光栅齿在 x 方向偏移一段距离,从 0 逐渐增加到 50 nm,每步重跑一次特征频率计算,记录最低虚部模式的 Q 因子。
对一组计算结果做对数拟合,我发现 Q 和 δ 之间呈很好的直线关系,斜率约为 −2,跟理论预期一致。这个标度律检验非常值得做,它不仅验证了你的模式确实是对称保护的,还能反过来确认仿真过程中没有引入其他不受保护的限制条件。如果斜率明显偏离 −2,比如变成 −1 或者更随机的数值,那大概率是几何设置里还有其他破坏对称性的因素,例如材料折射率的小数位数不对导致镜像对称被破坏,或者网格不对称造成数值耦合。
在实际器件设计中,这个 δ 与 Q 的关系是很有用的设计曲线。你想要 Q = 10⁵,反推需要的结构偏移是 δ ≈ 十几纳米,这个值落在电子束光刻的可控范围内;想要 Q = 10⁶,偏移就得控制到几纳米,接近工艺极限。有了这条曲线,做实验之前就能先决定工艺窗口和容差范围,这是仿真对项目最直接的贡献。
5. 踩坑记录:BIC 模拟里最常见的五个问题
5.1 假性 BIC:连续边界条件导致的隐形偏差
第一次做周期性波导模拟的人,很容易在 Floquet 边界设置上踩坑。默认的“周期性边界条件”会要求边界两侧场满足 E_dst = E_src × exp(-ik·r),如果你把这个波矢的符号设置反了,表面看色散关系没受太大影响,但模式对称性被错误地判断,某些模式的虚部会出现异常小的值,看起来像 BIC 实际却不是。检查方法是在求解前先做一次 kx = 0 的解析检查:在该点,Floquet 边界退化为普通周期性边界,如果你得到的模式分布不是严格的周期对称,就说明边界条件设置有问题。
还有一种更隐蔽的情况是使用了“连续周期边界条件”而不是 Floquet。很多人为了方便直接右键点“周期”选连续性,这在低频或者零波矢下没问题,但一旦要扫描非零 kx,连续边界没有相位偏移,算出来的能带是完全错误的,BIC 自然无从谈起。强烈建议所有周期性波导问题统一使用 Floquet,手动输入复相位参数。
5.2 收敛假象与特征值漏解
有一次我在一个看似已经收敛的网格上算出了 Q = 10⁷ 的高 Q 模式,心里还美滋滋,结果把网格加密一倍后,这个模式的虚部直接跳了两个数量级,完全不是 BIC 该有的行为。后来排查才知道,那其实是网格太粗导致的人为局域态——粗网格让某些高频模式原本存在的泄漏通道被数值离散“堵住”了,模型变相成了一个腔体,模式被虚假地锁住。
这个问题最有效的排查手段就是做系统的网格收敛性测试,不要只在一个网格密度下做判断。正确流程是固定所有物理参数,分别用 Lambda/20、Lambda/40、Lambda/80 三套网格计算目标模式的 Q 因子,观察它是否单调变化。如果 Q 随网格加密而持续上升,那就是潜在的 BIC;如果 Q 上升到一个平台后就不再变化,那是真实的低损耗模式;如果 Q 反而下降,那你之前看到的根本是网格导致的假象。
还有一个常见问题:特征值求解器漏解。特别是扫描 kx 时,有时候你明明知道某个频率附近应该存在一个 BIC 模式,但求解结果里就是找不到。这时候检查一下“所需模式数”和“搜索频率区间”,把区间往目标频率周围放大一圈;再不行就手动给一个“特征值估计”作为求解起点。做参数扫描时建议用“上一步解”作为初始值,能有效避免模式跳变和漏解。
5.3 PML 参数选择:厚一点还是薄一点
关于 PML,我踩过的坑是在算准 BIC 的时候,把 PML 的吸收系数调得过大。COMSOL 默认的 PML 参数在大多数问题里没问题,但对于 BIC 这种虚部信号极其微弱的计算,PML 的吸收特性会引起色散误差,导致本应完全消失的辐射通道被“部分吸收”掉,虚部被污染成一个有限值。解决办法是用两层结构:PML 区域保持默认吸收强度,但把它和求解域之间的边界距离拉大,减少倏逝波尾的泄漏。
理想情况下,你算出的 BIC 虚部应该不受 PML 厚度和吸收强度的影响。做一次双参数扫描测试:分别改变 PML 厚度和吸收系数,看目标模式的虚部是否保持稳定。如果虚部随 PML 参数变化了,说明你的计算域截断没有真正起到开放边界的作用,或者 PML 和物理域之间的耦合接口设置有问题。这种稳定化测试虽然浪费几个求解时长,但比调试半天找不到虚部异常的原因要高效得多。
5.4 性能优化:别再让求解器等一晚上
BIC 模拟最痛苦的不是求解器报错,而是你明明只扫了 30 个 kx 点,结果等了一晚上没跑完。二维模型本来求解很快,但如果你把整个模型域包括厚衬底和厚 PML 全部用细网格铺满,计算量就爆炸了。我建议用映射网格处理衬底和埋氧层,用三角形网格处理光栅区域,再在两个区域交界处使用“连接”功能保证不连续网格正常传递。这样大部分计算资源都集中在最需要分辨率的波导芯层和光栅附近,背景介质区域用粗网格即可。
“自适应网格细化”功能对特征频率研究是有用的,但别无脑开启。细化过程会在每次迭代中重新划分网格,虽然提升了模式精度,但耗时成倍增加。我的经验是先关闭自适应,跑完一遍找到候选模式,再局部加密特定区域做最终的精度验证。还有一点,扫描 kx 时可以先用较少模式数计算粗扫,找到 BIC 附近的频率范围后,再对缩小后的范围精细扫描。把“所需模式数”从 12 减少到 6,速度能提升一半以上。
5.5 后处理里的虚部陷阱:别把数值噪声当物理
最后提醒一个后处理层面的坑。COMSOL 里查看本征频率时,有时虚部会显示为负值,有时又显示为正值,这取决于时谐约定。在 e^{jωt} 约定下,模式的虚部为负表示衰减;在 e^{-iωt} 约定下,虚部为正表示衰减。如果你从不同版本或者不同接口得到的结果符号不一致,千万别慌,先确认全局定义里用的时谐约定,再把虚部取绝对值算 Q 因子。我见过不止一个同事因为符号问题把衰减模式当成增益模式,绕了大半天。
另外在画 Q 因子曲线时,如果某个 kx 点远处虚部小到跟浮点相对误差一个量级,比如 1e-12 以下,已经没法区分是真实 BIC 还是数值残留。这时候别过度解读那条 Q 曲线,它已经趋近发散,数值上达到了设备能表达的下限。你在文章里描述的应该是“虚部低于可分辨阈值”,而不是具体写出一个貌似精确的 Q 值,否则实验上根本无法复现。
6. 一点个人体会
跑完这一整套流程,我最大的感受是 BIC 模拟的门槛不在物理方程,而在数值细节的把控。方程本身在 COMSOL 里就是现成的电磁波接口,真正决定成败的是你对对称性、网格和求解容差这三个变量的处理。对称性决定了 BIC 是否存在,网格决定了虚部能不能算准,求解容差决定了你是否能稳定地追踪模式。这三个环节任何一个稍微马虎,结果就会从“理想的无穷 Q”退化成一堆到处漏光的普通泄漏模。
也是在做这个课题之后,我才真正理解为什么那么多人说 BIC 是“藏起来的完美模式”——它明明处在连续谱里,周围能量都在向外跑,只因为一个对称性就能纹丝不动地待在那里。算出来的那一刻,你会有一种“物理定律真的如此精巧”的实感。对于刚上手 COMSOL 波导仿真的朋友,我的建议是从最简单的对称结构开始,先复现别人文章里的 BIC 位置和 Q 因子标度律,再逐步加入扰动和复杂度。这个过程虽然要花些时间,但你会对整个仿真流程的每一步都建立起直觉,后续做任何周期性光子器件设计都能少走很多弯路。