做光子晶体仿真做了快一年最让我头疼的并不是网格怎么剖而是在COMSOL里把一个本该很局域的缺陷模算清楚之后还要判断它到底会不会漏光。最近一个月我把精力集中在二维光子晶体缺陷模BIC的数值探索上从建模、扫描、到Q因子分析踩了不少坑也把很多原本理解得模糊的概念重新捋了一遍。如果你也在做COMSOL光子晶体计算尤其是对缺陷模、连续谱束缚态Bound States in the ContinuumBIC以及高Q谐振器有兴趣这篇文章就是我基于实际操作的完整记录不是理论讲义更偏怎么一步步算出来。1. 为什么缺陷模BIC值得折腾光子晶体的核心看家本领是光子带隙PBG周期结构里某些频段的电磁波不能传播。传统的光子晶体功能器件比如波导、微腔大都利用带隙来约束光。引入缺陷后缺陷态会落在带隙中场被局域在缺陷附近理论上不会向背景泄漏这是大家都很熟的质量因子做得很好的路径。但BIC完全是另外一回事。它落在连续谱里也就是它的频率和可传播的连续扩展态重合按直觉必然向外辐射可对称性或干涉效应让这些态完全不耦合到辐射通道Q值理论上发散。用光子晶体结构去实现BIC比用简单平板腔难但也更能找到一些让人惊讶的模场分布。缺陷模BIC还有一个额外好处它是局域的携带的信息往往比扩展态更可控。从应用角度BIC可以把Q因子做到极高所以近些年很多研究都在往这个方向靠。低阈值激光器、高灵敏传感、非线性谐波增强、手性光发射只要需要强场约束和窄线宽BIC都常被提及。普通的侧向光栅腔Q做到几千几万就需要精细加工了BIC结构却在理论上可以直接冲上十的六次方、七次方以上这种数量级的差异让实验组和仿真组都很难忽视它。对我的工作来说这个问题的计算目标非常明确在COMSOL里构造一个带缺陷的二维光子晶体超胞找到缺陷模然后判断这个缺陷模在什么情况下变成BIC定量给出它的Q因子。这里要提前给新手打个预防针BIC是在严格意义上存在于某个参数点的理想情形而我们数值模拟能得到的是Q因子尖锐峰网格剖分越细、超胞越大峰越高。把数值峰和严格BIC区分开是后面所有后处理工作的出发点。2. 动手前先把物理图像和计算目标拆干净2.1 缺陷模带隙里的局域态还是带内的寄生虫我最初对缺陷模的理解就是禁带中间单独的那一条平坦能带。这个理解只对了一半。引入点缺陷后如果缺陷模频率落在完美晶体的带隙范围内确实与所有布洛赫模式都不共振场指数衰减是严格局域态。但有些缺陷模的频率并不在带隙里而是落在允许带中这时它可以和连续扩展态共存并不自动衰减。这种允许带内的局域态通常是有泄漏的也正因如此它才有机会在特定条件下变成BIC。搜索文献的时候被一些概念绕晕是因为不少人把缺陷模和带隙缺陷模混着用。真实物理上带隙缺陷模的Q由结构损耗和工艺误差决定而BIC则是Q→∞的极限数值上只是在参数空间里一个孤立的奇异点。缺陷模BIC的研究动机恰恰是让这个寄生在连续谱里的局域模变成可用的高Q态而不是躲避到禁带里求安稳。2.2 BIC与高Q模的本质区别高Q模是工程上可以实现的BIC是理想数学极限。为了后面计算时不误判我给自己整理了三类常见BIC的画像。第一类是对称性保护BIC最容易在高对称点上出现比如Γ点的偶极子模和辐射通道奇偶性不匹配耦合积分严格为零。它不需要精细调参换个结构只要对称性不变BIC就一直在。第二类是Friedrich-Wintgen型BIC靠两个模式在同一个连续谱通道中干涉相消即使结构不对称也能产生。它通常需要调节某个几何参数使一个模式的远场与另一个模式发生完全相消所以出现的是参数空间里的孤立点。第三类是Fabry-Perot型BIC多与扩展态在边界处多次反射形成的束缚有关在光子晶体平板里也比较常见。缺陷模BIC常常与超胞的结构对称性绑定在一起计算上最便宜也最适合作为第一次的复现对象。判断是否BIC有三个层次第一场分布局域在缺陷附近是局限于单个缺陷的局域态第二向外辐射的功率为零或理论上严格为零第三Q因子在k空间表现出发散尖峰。单看任何一个指标都会误判必须三管齐下。2.3 在COMSOL里怎么观察这些东西COMSOL能做两件事一是本征值计算二是频域响应。对BIC来说本征值计算更直接——寻找复频率本征值虚部对应泄漏损耗Q对应比值。光子晶体带结构本质上是一系列k点下的本征频率。缺陷模BIC就是在某个k点附近那条缺陷色散线扫过连续谱时突然出现一个复频率虚部急剧减小点。所以我的完整流程是先算完美晶体的带结构再改超胞引入缺陷扫描标准k路径并从结果中筛出缺陷模频带随后加入PML层提取复本征值计算每个k点的Q值最后分析远场确认泄漏为零。下面每一节对应一个环节。3. COMSOL建模模块、几何、边界条件一次到位3.1 物理接口和研究类型别再纠结选哪个COMSOL中做光子晶体带计算常用RF模块或波动光学模块的电磁波频域接口。二维情况下如果电场沿z方向是TM模磁场沿z方向是TE模。带隙位置、模式顺序和偏振关系很大不要选错。推荐在模型向导里选模型空间二维物理场选RF模块→电磁波频域或波动光学模块→电磁波频域研究选特征频率。这里有个细节接口默认的因变量是电场分量Ex、Ey、Ez和磁矢势对2D问题若你用TM偏振只需保留Ez这一个分量运算量小很多。在组件中把不必要的分量关掉或者在波动方程电节点中设置面外维度即可。在设置中指定折射率或介电常数。频率单位无所谓但我后处理时习惯统一转成归一化频率u f·a/c这样横向对比文献更方便也便于跨结构比较。单位制的混乱会造成后面所有实验对比出错建议从建模第一分钟就把坐标系、尺寸单位、频率单位写清楚。3.2 参数和几何用归一化尺度才省心我的做法是先定义全局参数参数含义参考值a晶格常数520 nmr_hole空气孔半径0.26*an_mat介质折射率3.48n_air空气折射率1.0kx, ky归一化波矢分量0 到 0.5几何其实非常简单一个边长为a的正方形域中心挖掉一个半径为r_hole的圆。如果多孔结构可以用COMSOL的阵列功能铺开n×n个圆孔再做一个外边界为n倍a的大正方形最后用差集得到平板。把圆边界单独命名后面网格要加密。孔板的核心性能由r/a这个比例决定r/a太小带隙会消失太大孔之间材料太薄计算不稳定常见范围是0.2~0.4。材料与域设置平板区域设成ε n_mat²孔内设成空气ε 1。如果是柱状结构介质柱在空气中反过来设置即可。不管哪种边界条件都一样唯一差别是带隙的位置和模式分布。材料参数尽量用常数折射率先不用色散等模式找对之后再考虑材料吸收带来的Q变化。3.3 周期边界从元胞到超胞核心是Floquet波矢完美晶体计算用单位元胞边界条件用Floquet周期边界。在COMSOL的电磁波接口下添加周期性条件边界选择类型Floquet周期然后设定波矢k的实部分量kx k_bx * π/aky k_by * π/a。单位这样定义是为了方便扫描第一布里渊区的边界。参数扫描的路径我一般取Γ(0,0) → X(0.5,0) → M(0.5,0.5) → Γ(0,0)把参数k_bx、k_by按线性插值打点。一个常见做法是用扫掠参数名定义变量k_path从0到3然后给k_bx、k_by赋上分段线性表达式。COMSOL参数扫描支持参数化扫描但不方便在一个列表里直接写分支表达式我会在参数定义处提前写好这两个条件表达式再用参数扫描去扫k_path导出数据后就能直接画出带图。超胞计算时周期边界的元胞尺寸从a变成了Na倒格子也随之缩小。如果还把k定义成原胞时的π/a单位必须小心缩放。比如5×5超胞的Γ点对应的其实是一组折叠后的模式缺陷模可能出现在任意一个折叠带上。更方便的做法是默认扫描小的k范围比如kx从-0.1到0.1单位π/a看缺陷模在Γ点附近的色散和Q变化。3.4 网格划分影响结果的第一道坎光子晶体仿真最容易翻车的就是网格。圆孔表面如果网格太疏等效折射率就错了带隙位置和缺陷模频率都会漂。我一般用三角形网格最大单元尺寸设为λ/(5√ε_eff)其中λ c/fε_eff取背景材料折射率的平方做初步估计即可圆边上再加一个边界层网格至少2~3层每层厚度约λ/(80√ε)圆孔内部的空气区域网格也要保留不要让孔内部完全没有体单元。网格的最大尺寸不宜只设一个全局固定值COMSOL会根据波长自动适应只要在网格节点里设置多个尺寸等级全局最大单元尺寸、圆边边界层、以及孔内单独最大单元尺寸。如果算完发现特征频率对网格密度还很敏感那说明网格还不够细加大密度重新算。同样的模型网格稀疏时算出的Q数值可以差一个数量级以上这是仿真新手最不容易察觉的地方。4. 从完美晶体到缺陷晶体完整计算路径4.1 先扫出完美晶体的带隙图完成单胞建模后将研究设定为特征频率特征值范围先给一个足够宽的窗口比如归一化频率u的范围在0.1到0.8对应的角频率ω(rad/s) 2πc·u/a。注意COMSOL特征频率的单位可能是rad/s或Hz后处理换算公式固定为u f·a/c ω·a/(2πc)。设置参数扫描k_path从0到3求各点前若干阶本征值。求解器采用特征值求解器线性求解器选MUMPS因为网格较大时SPOOLES会慢。完成后在结果中生成一维绘图组横坐标k_path纵坐标为第一个到第N个本征频率每个本征频率一组数据。把最多12条带画出来找出带隙范围。跑完一版典型参数a 520 nmr/a 0.26n 3.48后TM模带隙大致出现在u 0.2~0.4区间的某个位置具体由r/a决定。这一步的目的主要有两个一是确认带隙的频段方便后面设置特征值搜索范围二是存下来作为之后判断缺陷模是否进入连续谱的参照系。如果后面算出来的缺陷模频率落在带隙内那它就是传统局域模如果落在带隙外却又局域在缺陷中心就有BIC潜质。4.2 超胞构造与缺陷引入确定带隙后用5×5或者7×7的超胞替换原单胞。这里不建议直接用大超胞而不用周期性条件否则四周边界上会出现假态。超胞加上同样的Floquet周期边界物理含义是让超胞本身成为新的晶格基元无限重复下去从而模拟缺陷阵列。缺陷的引入方式有很多把中心孔半径改成0完全删除、改成0.1a缩小、改成0.4a增大或者把中心孔的位置稍微偏移。圆孔半径变化最简单因为改动参数后COMSOL自动重新剖网格。偏移中心孔会破坏对称性可能把对称性保护的BIC直接破坏掉所以第一次尝试最好只改半径保持对称性不变。超胞大小选择很讲究。理想状态下缺陷模指数衰减5×5足够时缺陷间的耦合很小如果算出来模式场在超胞边界处没有衰减到噪声水平就是超胞不够大要换7×7甚至9×9。超胞越大特征值求解的未知量越多内存也跟着涨。我的机器是16G内存算7×7加PML已经有点喘5×5则非常流畅所以先用5×5扫盲区、再用7×7验证关键Q值是合理省钱的方案。4.3 缺陷模怎么从结果里筛出来把超胞模型的参数扫描从原来的k_path简化成只扫Γ点附近的一小段比如kx 0、0.01、0.02直到0.15单位是π/a。然后求特征频率。带缺陷的超胞本征频率列表很大不能直接看。筛选缺陷模有三个实用特征第一色散曲线几乎是平的不随k明显变化这是因为缺陷局域态几乎不受周期边界波矢影响第二场图能量集中在中心区域周围按指数规律衰减第三出现在带隙内时它上下没有别的扩展带非常显眼出现在允许带内时需要与扩展模交叉区分——就是看场分布扩展模是均匀铺开的缺陷模是集中在中心的。后处理中我会用电场模的二维表面图扫一遍各本征模按场局域程度排序。将局域程度定义为中心区域比如以圆心为中心、半径为a的圆内的电场能量占总能量的比重。这个量可以自己做积分后处理步骤是在结果节点中定义一个积分算子分子是中心圆域内的能量积分分母是全域能量积分。如果一个模的局域化占比超过70%大概率是缺陷模接下来才谈得上BIC。5. 判断BICQ、远场与对称性三级论证5.1 复本征值与PML的正确姿势判断BIC最直接的数值证据是Q因子发散。对无吸收材料的结构若边界完全周期且超胞足够大BIC本征频率的虚部应为零普通泄漏模则会被数值吸收写成虚部较大的复值。关键问题在于普通特征值求解默认假设无损耗这时频率是实数的。要让向外辐射的泄漏变成可量化的损耗需要人为引入吸收边界。在COMSOL中主流做法是在超胞外面加一层完美匹配层PML。几何上做一个外框区域把整个超胞包起来外框的厚度通常取半个工作波长到两个波长PML内部各向异性材料参数由COMSOL自动填充。添加PML的具体操作组件中将PML区域设为完美匹配层域节点指定厚度和吸收方向。把PML外部边界设为散射边界条件或者默认的低反射边界。重新求解特征频率从复特征值得到Q Re(ω) / (2|Im(ω)|)。注意一个反直觉之处PML不能太薄否则对掠射角分量吸收不足Q会被低估也不能太厚太厚会有数值反射。我一般用厚度0.5~1.5λ做几次Q收敛测试取Q对该厚度不敏感的区间作为最终设置。5.2 Q因子在k空间的扫描峰就是印记选定缺陷模后扫描kx从0到0.1或0.2每隔0.005采样算出每个点的Q。正常泄漏模的Q随k平滑变化而在BIC点会出现窄而尖的峰。要判断这个峰是不是真的发散最简单的办法是看收敛性加密网格后如果峰顶Q进一步提高且峰的位置逐渐稳定说明趋近BIC若峰顶不再变化说明只是一个普通高Q泄漏模不是BIC。这里还要注意超胞的折叠效应。超胞的倒格子比原胞小了N倍因此Γ点附近很小的k取值就可能已经跨越原胞布里渊区的若干个高对称点。我在5×5超胞上看到的现象就是模式数量暴涨原因就是折叠需要依托原胞带图去对照别把折叠模当缺陷模。把原胞带图和超胞模式频率叠在一起画能很快认出哪些是折叠带哪些是真正的缺陷模。5.3 远场分布和对称性辅助判定复频率提供Q值场分布提供物理图像。BIC可以不向外辐射但近场仍然可以有一个很好的局域模。在COMSOL后处理中我可以画出远场辐射功率分布做法是添加远场节点或者手动在PML外侧边界做边界能流积分。远场图里如果出现完全为零的辐射角而近场又局域这就和BIC的物理图像对上了。在Γ点分析对称性最重要。假设超胞和缺陷本身关于x轴和y轴都对称缺陷模如果是某种奇对称模式而相同频率处所有可辐射扩展模都是偶对称那么二者耦合积分为零就是对称性保护BIC。COMSOL里可以分别选择x0和y0平面附近做模式分解观察电场分量在镜像操作下的符号变化给模式贴上偶/奇标签。如果缺陷模标签与同一频率处扩展模的标签正交就判BIC。我实际试的一组结构是Si孔板a 520 nmr/a 0.265×5超胞中心孔半径缩小为0.1a。在这个结构里Γ点附近存在一个清晰的缺陷模kx 0时Q可以到10的6次方以上kx偏离0.02π/aQ掉到几千。这种急剧下降的曲线配合远场几乎没有辐射基本可以认定为对称性保护BIC。6. 实际计算中踩过的坑做这种计算踩坑是必然的我列几个对我影响最大的都是真金白银的时间换来的。坑一特征值搜索范围太窄缺陷模被遗漏。超胞之后模式数成倍增长一开始我只搜前10个本征值缺陷模根本没出现。实际上缺陷模可能在频带中位于第20~30阶。解决办法是先把搜索范围扩大搜前50~100个本征值或者先用带隙位置做预估再用搜索接近频率的方式求解。坑二k参数归一化搞错。单位原胞的k值定义和超胞不一样。我把超胞的kx也写成k_bx*π/a但这个π/a是原胞单位在5×5超胞中应该除以Na否则色散曲线完全乱掉。我的经验是始终在参数注释里写明以原胞基矢长度归一还是以超胞基矢长度归一别靠脑记。坑三PML与周期边界打架。PML要吸收泄漏而周期边界要求场在超胞两侧相等。PML区域不能放在周期边界之间否则周期条件被PML截断。我一开始把PML放在超胞内部结果求解出大量虚假模式。正确做法是周期边界加在超胞四周超胞外再包一层PMLPML最外边界用默认低反射边界。坑四网格各向异性带来虚假高Q。如果网格在缺陷区特别密、在扩展区特别粗算出的Q会因数值泄漏而偏低反过来如果边界层网格把本应辐射掉的场也包住了Q会被虚假抬高。我试过同一个模型把圆周边界层加厚Q夸张地从几千跳到十万但那只是网格反射造成的假象。统一判据是PML厚度和网格密度同时翻倍若Q仍然稳定收敛数据才可信。坑五纯2D模型对实验的局限。二维光子晶体忽略了垂直方向的泄漏BIC在真实平板中可能因为垂直辐射通道而不再严格存在。二维模型是研究机理和教学的好工具但如果你要对标实验至少得换成3D平板模型并扫描面外辐射Q。文章里可以写二维结论但心里要清楚它的边界。坑六特征值求解器内存溢出的处理。7×7以上超胞加PML时特征值求解放大内存可能不够我会把线性求解器改成MUMPS并将特征值求解器的重启动向量数从默认值调大。再不行把直接求解改成迭代求解或者降低待求特征值数量。批量扫描时尽量单次求解再逐个取结果避免一次把所有k点都塞进求解器。坑七复频率变量提取顺序。COMSOL返回的复特征频率有时实部虚部顺序会把后处理搞晕。要严格用全局计算提取real(freq)和imag(freq)再做Q real/(2*abs(imag))。别在表达式里写错符号否则Q值变成负的都不知道。坑现象对策搜索范围太窄缺陷模不出现扩大搜索特征值数量或用频带预估k单位混淆色散乱、模式错位注释清楚归一化基准统一用原胞单位PML位置错大量虚假模式PML在超胞外部周期边界保持完整网格反射Q被虚假抬高同时加密网格和PML厚度做收敛测试2D局限纯2D Q偏大对标实验时换3D模型内存不足求解崩溃改MUMPS减少特征值数量或缩小超胞复频率提取乱Q为负或异常用real/imag单独提取再计算7. 最后一点实操经验这篇文章从动机到建模、计算、判断基本走完了一个缺陷模BIC探索的闭环。以我的个人经验做这种仿真最重要的不是熟悉COMSOL菜单而是每一步都知道自己在验证什么。先算完美晶体带图再引入缺陷找局域模然后加PML提取Q最后用对称性和远场交叉验证这套固定流程跑熟了无论换成哪种晶格、哪种缺陷形式都能快速上手。如果你只是想复现一个BIC案例我建议从5×5超胞、Si孔板、TM模开始先把Γ点的对称性保护BIC试出来。等看到Q峰出现再慢慢引到连续谱中的Friedrich-Wintgen型BIC。到时候你会发现COMSOL里能做的玩法其实非常多难度主要在如何设置一个既算得动又有物理意义的模型。我再分享一个习惯每一个参数批量计算完成后第一时间把频率、Q、局域占比三个量导出成txt或csv之后画图和处理都用它不反复重新打开模型。这些看似琐碎的做法在连续几百组参数扫描时能救命。祝你的BIC探索之路少踩坑、多出好图。