COMSOL实现非偏振等离激元BIC超表面建模与Q因子分析

1. 等离激元BIC是个什么东西先通俗过一遍1.1 BIC在“连续谱”里赖着不走的模式先讲物理。光学超表面或等离激元阵列中电磁波在自由空间、衬底、周期结构之间来回耦合大多数模式都会朝外辐射能量能量流出去叫“泄漏”。但有一类特殊模式频率明明落在辐射连续谱的范围里却因为自身场分布与所有外部辐射通道的对称性互相排斥能量怎么也漏不出去只在结构局部维持振荡——这就是“连续区束缚态”简称BIC。可以这样理解这事你站在嘈杂的大厅里周围全是说话声这就是连续谱你想跟旁边的人说一句悄悄话理论上声音会被喧闹淹没。但如果你恰好站在声学死角后背贴着墙嘴对着墙角你那个频率的声音虽然落在所有人说话的频率范围之内却没有向外传播只有贴着墙的这个人能听见。光学BIC也是靠“波函数节点”和“远场通道对称性失配”实现困光的本质上不是一个多神秘的机制而是一个线性系统里特征模式解不出的漏波解。在周期性微纳结构中BIC最常见的实现是“对称性保护BIC”模式在单元胞内的场分布相对于某个对称操作是奇函数而平面波辐射通道对应的场分布是偶函数两者积分结果为零耦合矩阵这一项直接消掉。于是这一模式既不会被平面波激发也不会自发辐射。你从实验上观察它它始终是“暗”的只有打破结构对称性它才被“释放”出来变成可观测的高Q模式。1.2 等离激元BIC和纯介质BIC的差别介质超表面的BIC通常用高折射率纳米柱硅、二氧化钛来实现Q值动辄上万。等离激元BIC换成金属/介质复合结构后局域表面等离激元共振参与进来给模式带来了高度集中的电场增强能力。金属引入欧姆损耗纯材料的Q值被压低但换个角度想局域场增强、表面灵敏度反而提升了像做传感、增强非线性、荧光增强这类应用恰恰需要这种“又把光锁住、又把场压到极小空间”的特点。从COMSOL建模的角度看这两类BIC的差别很实际。纯介质BIC往往在长方体柱阵列中实现建模简单材料折射率直接给常数就行。等离激元BIC则必须处理金属的复数介电常数色散还要保证网格在金属表面足够细仿真代价和收敛难度都高一截。哪个更难收敛我在实际仿真里觉得最麻烦的不是金属本身而是“如何让一个本不该辐射的模式在数值上确实保持不辐射”——这是后面要重点聊的。1.3 为什么要纠结“非偏振”多数的BIC结构设计要从一个对称性完整的单元出发再人为破坏对称性BIC才能转成准BIC被外部平面波激发。问题是很多人破坏对称性时只破坏一个方向圆盘改成单方向椭圆孔只往一个方向偏于是结构变成了对x偏振敏感、对y偏振不敏感的单偏振器件。现实中的光谱检测、荧光增强、偏振无关光学滤波等场景都要求器件对入射偏振不敏感。也就是说结构既要“破缺对称”把BIC放出来又要在x和y两个方向同时、等效地放出来。设计难点就在这里。本文后面给出的建模思路就是围绕“什么样的破缺既能打破BIC的暗态又能保持x/y等价”来展开。2. 非偏振BIC结构的设计逻辑对称性与扰动怎么搭配2.1 从“过度对称”的暗态出发最省事的起点是周期金圆盘阵列。单元胞结构和电磁模式都是C∞旋转对称场的多极展开里包含偶极、四极、六极等成分。只有偶极项能与平面波直接耦合四极以上因辐射矩阵元为零而成为暗态BIC。COMSOL里给这种结构设置一套合理的默认参数就可以作为工作起点周期P600 nm金圆盘半径R160 nm盘厚度t50 nm底层二氧化硅衬底厚度取300 nm上方空气覆盖层高度取500 nm以上为什么要取这些数我通常按近红外波段来定。金盘局域等离激元共振大致在P/R比值约为3~4的时候出现在700~900 nmP600nm、R160nm时阵列瑞利异常约在600nm附近它与局域等离激元耦合能形成表面晶格共振。BIC往往就在这个共振附近产生具体位置靠特征频率扫描找。你不需要一上来就追求参数完美先把几何建出来、把模式找到再根据实际共振位置微调周期和盘径比空想参数靠谱得多。2.2 怎么扰动才能“偏振不敏感”打破圆盘对称性的方式很多但效果完全不同圆盘改成椭圆只打破了一组正交方向的对称性响应偏振依赖明显。圆盘上偏离圆心挖一个圆孔单孔偏移同样只在某一方向上破缺。圆盘上挖两个孔分别沿x和y方向偏离整体结构保持x/y互换对称是两个正交偏振都激发、且激发率相同的最直接方式。圆盘上挖一个沿45度方向偏移的椭圆孔同时把椭圆长轴放在与偏移方向垂直的方位也有机会做到x/y等效但需要对参数精细调整不建议新手一上来就试。我在类似项目中用得比较顺手的是“双偏心孔”方案。单元胞内金属圆盘半径R200 nm厚度70 nm盘上挖两个圆孔半径r均为45 nm第一个孔中心相对于盘心沿x正方向偏移d第二个孔沿y正方向偏移d。整体结构对交换x和y坐标保持不变而相对于x轴、y轴的对称性都被破坏了。这样原本在严格对称结构中完全暗态的BIC模式被同时“放行”给x偏振和y偏振入射光而且两边放行的强度相同。d从0开始增大到80~100 nm的过程就是Q值从趋于无穷回落到几百、几十的过程。2.3 用COMSOL参数化这一步要做的准备COMSOL里做参数扫描前先把全局参数建好。我习惯把所有几何尺寸用参数表示后续优化不重建模型。参数写法如下P600e-9 R200e-9 r45e-9 d0 t70e-9 h_sub300e-9 h_air600e-9建立几何时用“块Block”做空气背景和衬底“圆柱Cylinder”做圆盘两个“圆柱”做孔用“差集Difference”把圆盘上的孔挖掉。关键点是孔的偏移方向第一个孔的圆心坐标参数化表示为(d,0, t/2)第二个为(0,d, t/2)坐标原点放在盘心z方向为厚度方向。这样一个简单的参数化操作就把“破缺正交方向保持等价”这个设计条件固化在几何里了。下一步是检查几何对称性是否真的满足x/y等价。最简单的验证办法在COMSOL里把单元胞复制一份绕z轴旋转90度再和原胞对比若完全重合则x/y偏振一定等价。我一般通过“变换Transforms”里的旋转功能做一次直观检查比看代码可靠得多。这一步看似多花几分钟但它能提前暴露很多“明明应该对称却怎么也算不出对称结果”的问题。3. COMSOL建模的完整流程物理场、边界、网格与求解3.1 物理场接口与周期性条件怎么选模型用“电磁波频域”接口接口标识ewfd三维。对BIC研究分两个阶段。第一阶段找模式用特征频率求解。单元胞四周定义“周期性条件”类型选Floquet周期在Γ点kx0ky0下求解。此时顶部和底部不用PML直接留一定厚度的空气和衬底即可特征频率求解器会自动计算开放空间中的辐射衰减有辐射泄漏的模式虚部大BIC模式的虚部极小。第二阶段算响应用频域求解。把激励方式改成“周期性端口”在单元胞的一个面上定义入射端口另一个面定义出射端口。为了分别验证x和y偏振端口模式需要分别设置极化。COMSOL的周期端口在使用时需要和“边界对”配套端口类型选“衍射光栅”阶数按实际设置的布洛赫周期对应。需要注意的一点是在特征频率求解模式下周期性端口不一定能很好地配合所以常见做法是先用特征频率找模式再用频域加端口做透射谱验证。3.2 材料设置金属色散必须用实验数据金属的相对介电常数关系紧密。建议不要在COMSOL里直接用简化的Drude常数一路算到底因为金属膜的实际损耗、衬底附着效应会让结果偏移。更稳妥的做法建立插值函数epsilon_real(λ)和epsilon_imag(λ)数据来自公认的材料光学常数手册。或者在COMSOL材料库中找到金或银的“光学”材料但注意部分库材料属性覆盖波长范围有限超出范围会外插出奇怪结果。如果做的是近红外波段要用包含该波段实测值的色散数据。金在近红外区的典型行为是实部绝对值逐渐减小从约-12到-30左右变化虚部在10左右。如果你拿一个只覆盖可见光的材料库数据外插到1200 nm结果很可能偏离物理。在频域和特征频率两种求解里金属损耗的处理要分开看。特征频率求解时如果金属损耗打开BIC的虚部被欧姆损耗钳制很难看到“数值上趋于零”的BIC现象。为了分析辐射通道我通常先关掉损耗把金属色散里的虚部暂时设成0让虚部仅由辐射决定等验证完模式归属再打开损耗做最终结果。这是建模中的一个关键实用技巧能让你把“辐射Q”和“材料Q”两种物理机制分开来看。3.3 网格划分BIC仿真成败的关键网格这一块我要多说几句BIC数值结果对网格的敏感度比普通透射问题高一个量级。原因很简单虚部这个量本身是几乎为零的小数任何离散误差都会被放大成相对很大的虚假损耗。实际做法金属盘和衬底界面处用边界层网格第一层厚度1~2 nm层数5层增长率1.2。金属内部最大单元尺寸设为20 nm。单元胞整个计算域最大单元尺寸设在40 nm左右。空气区域可以放到80 nm但不要超过λ/10近红外800 nm时约80 nm。孔口、盘边缘是场集中处必须局部细化为10 nm以下。算完首先要做收敛性检查把最大单元尺寸分别设成40、30、20 nm对比BIC模式虚部。如果虚部随网格加密而显著变化说明网格还不够。通常要等到相邻两次网格对应的Q值相对变化小于2%再用来做最终数据处理。三维模型在这步计算量很大我一般先用粗网格定位模式、确定大致频率再针对单个目标模式做局部细化重算而不是全程细网格。3.4 特征频率求解器的设置细节特征频率求解器设置不当最常见的现象是“找不到想要的模式”。我一般这样设置“搜索的特征频率”采用“指定实数/复数频率附近”把初始值设在估计的共振波长对应的角频率并设置一个较小的搜索半径。一次计算模式数量给20~30个保证包含多个高阶暗模式。周期性条件的布洛赫波矢在Γ点设置为00用“Floquet周期”的“从k向量”自动配置。另外一个很容易忽略的点单元胞的背景空气层和衬底高度不能太薄要把倏逝场给够。如果空气高度只有200 nm高频泄漏模式可能被顶面边界影响虚部失真。我建议空气厚度取0.5λ到1λ衬底厚度取0.3λ以上。对于800 nm工作波长空气厚度500~800 nm、衬底厚度300 nm是安全的。4. 结果怎么提BIC定位、Q因子拟合与偏振验证4.1 特征频率结果里怎么认出BIC在d0的严格对称结构里特征频率列表会出现一批模式其中大部分是亮模式虚部在MHz到GHz量级Q值几百到几千。真正的BIC通常虚部最小小到10^3 Hz甚至10^2 Hz以下但同时你要看场分布确认它不是数值伪模式。我在看结果时常做两件事。第一件是打开特征频率对应的电场模分布BIC模式沿金属盘圆周有多极节点或者双孔结构的局域场围绕两个孔对称分布。第二件是看表面电荷密度分布反对称分布说明它不能和均匀平面波耦合。这种模式下电场增强集中在金属表面很薄的一层里这正是做灵敏检测的价值所在。需要提醒的是数值上d0时虚部不会严格等于0。由于离散化BIC的虚部会是一个很小的有限值。只要虚部比同一频率附近其他模式小2~3个数量级就可以认定是数值意义上的BIC。不要为了追求虚部绝对为零去调网格那样会浪费大量算力而且根本做不到。4.2 Q因子随扰动参数的变化规律d从0开始增大BIC变成准BIC虚部逐渐增大Q值下降。准BIC理论给出Q与无量纲扰动参数α的平方成反比αd/R。也就是说双对数坐标下Q-α关系是一条斜率-2的直线。扫描d从10 nm到100 nm每步5 nm用特征频率求解提取实部ω0、虚部γ算Qω0/(2γ)。然后把Q值随α的变化画在对数坐标里。如果斜率在-1.9到-2.1之间说明模式行为符合准BIC物理如果斜率明显偏离就说明可能有其他泄漏通道参与比如衬底高阶衍射或模式杂化需要重新审视模型。这个Q与α的规律也是文章里最有力的“证据图”比直接贴一个透射谱更有说服力。建议把特征频率结果导出到表格在图形界面或者外部工具里做拟合COMSOL自带的结果绘图功能也够用。4.3 非偏振响应的验证方法非偏振不是说“对任何偏振都不敏感”而是对x和y两个正交线偏振的响应一致。验证分几步取一个工作点比如d60 nm的结构。在频域求解中先把入射端口模式设为x偏振扫描650~1000 nm记录透射谱或吸收谱。再把入射端口模式设为y偏振做完全相同的频域扫描。对比两条透射曲线主共振峰位置应当重合半高宽应当一致谷底深度相对偏差建议控制在1%以内。若两条曲线高度重合就可以认定这个结构是非偏振响应的准BIC超表面。如果想更严格再做一个“斜入射角扫描”测试。虽然标题里说非偏振实际应用里通常还会关心入射角稳定性这里可以加一点角度扫描的结果但别把范围扯太远。除了透射谱还可以对比两种偏振下共振波长的近场增强。比如提取金属盘表面的最大电场增强两个偏振方向计算结果之差小于5%就进一步坐实了“偏振不敏感”的结论。5. 常见问题与实操排查5.1 “特征频率里根本没有虚部极小的模式”排查顺序看是不是周期条件没配对。Floquet周期条件要求单元胞相对两个面的网格和物理场完全匹配如果其中一对边界没配对计算域就变成“开口”结构所有模式都带大虚部。看几何是否有微小的无意破缺。比如孔偏移参数d虽然是0但布尔运算后的边界有细微gap或碎面也会破坏对称性。看搜索脉宽设置是否太窄。BIC模式可能不在你猜测的初始频率附近把初始频率范围放宽多算一批模式再筛选。5.2 虚部始终下不去像被什么东西托住了这种情况先区分是“物理损耗”还是“数值假损耗”。把金属损耗关掉再算一遍把色散虚部设成0看虚部是否明显降低。如果关了损耗还是托在10 MHz级以上那就是网格或边界问题。网格问题上首要怀疑金属孔口和盘面的解析度。此时把局部最大单元尺寸从10 nm降到5 nm重新算一次观察虚部变化幅度。若变化超过10%继续细化直到变小到可以忽略。边界问题上把背景空气厚度从500 nm加到800 nm看虚部是否下降。若显著下降说明原先是计算域截断引起的近似泄漏。把边界条件改成PML或者直接增大厚度都可以。5.3 两条偏振透射谱“不合拢”偏振不敏感的设计没做到位最常见的低级错误是几何只在一侧偏移。比如圆盘上只挖了一个偏心孔这种结构主响应其实只在一个偏振另一条透射谱显示不出同样的深谷。解决办法就是换成双侧偏心孔或使单元胞在x/y旋转90度后完全重合。还有一种隐蔽情况几何看似对称但材料属性各向异性或衬底的晶轴方向让x/y不再等价。COMSOL里如果材料定义为“各向异性”尤其需要注意折射率张量方向。仿真时统一用各向同性材料除非确有必要。5.4 参数扫描里模式“跳变”特征频率随d增大时和附近亮模式发生交叉/避让追踪会突然跳到另一条分支上。解决办法把一次模式数量提高取出全部特征频率后用手工判断。更可靠的是通过模式场重叠积分来追踪同一条分支。计算第n步与第n1步之间两个模式的场重叠积分乘积大的就是同一分支。如果只是画Q-d曲线可以限制在某个频率范围并固定搜索半径COMSOL会按接近参考值的顺序排列特征频率。5.5 模型太大、算得太慢三维等离激元BIC模型确实吃资源。我的建议是先算二维截面结构验证物理规律再扩展三维。参数扫描时粗网格多步扫找到趋势后再精算若干关键点。特征频率计算用稀疏直接求解器内存不足时改用多重网格预条件迭代求解器但要注意收敛性。近红外模型如果不需要太细节可以先关掉金属损耗用约十万自由度的网格看模式趋势最后对工作点做一次高精度细网格计算。提示如果你用的是双偏心孔方案d扫描范围别开太大。d超过孔半径之后两个孔可能在盘边缘重叠几何变形成非物理形状。我一般把d上限限制在R - 2r附近这样既能扫出Q值几个数量级的变化又不会惹出几何报错。坦白说我在这类模型上踩过最多坑的就是“网格”和“损耗”两者纠缠不清的那一步。BIC的虚部在理论上就是要等于零但数值上它永不会严格为零你必须用“相对其他模式的虚部小几个量级”这个判据而不是死磕绝对零值。另一个体会是非偏振设计不是“做一个圆对称结构再让它BIC被激发”——完美对称下BIC根本激不出来真正的物理在于“让破坏对称的那个扰动同时具有两个正交方向的等效性”。如果一开始就围绕这个原则选几何后面的建模和验证都会顺很多。最后给一个小技巧做d扫描之前先把d0的BIC模式在粗网格下固定好并在求解器中把“特征频率搜索”的参考值设置成该模式的频率附近。这样后面每一步扫描时求解器都会优先返回接近BIC分支的模式模式跳变的问题会少很多。希望这篇经验能帮你少走弯路。