简介在光子晶体研究中平带合并与连续态束缚态BIC是高效率光学器件设计的关键课题。这份COMSOL仿真资源包面向微纳光学领域的研究生、工程师及科研人员聚焦实现二维能带、三维能带、品质因子和远场偏振的复现计算可帮助读者掌握光子晶体仿真分析的完整流程。包内共10个文件包括6个txt技术笔记与解析文档、2张jpg结果图、1份html说明及1个doc文档整体仅383KB内容紧凑且易于检索。已有243人学习下载。资源不仅详细阐述了平带合并与BIC的原理还给出了COMSOL建模参数设置、能带结构与品质因子的计算结果并对远场偏振特性进行了深入分析配合可视化图片与分章节技术总结读者可以系统复现二维/三维能带计算理解BIC现象在光学器件设计中的意义。 前阵子一直在折腾一个光子晶体的课题——用COMSOL复现平带合并BIC连续谱束缚态的工作顺手把二维能带、三维能带、品质因子和远场偏振计算全走了一遍。这活儿看着只是“跑个仿真”实际做起来从建模、边界条件、特征频率扫描到后处理每一步都有不少讲究。这篇东西打算把这套流程完整拆开把我调过的参数、踩过的坑、以及那些“文档里找不到但实测很有用”的经验都写出来。适合正在做光子晶体、微纳光子器件或者刚接触BIC想快速上手COMSOL特征频率计算的研究生和工程师参考。先交代一下背景。BIC简单来说就是本征频率落在辐射连续谱里但不往外辐射的模式理想情况下Q因子趋向无穷大。光子晶体平板是研究BIC最常用的结构周期性调制让导模和辐射模之间产生耦合可以在某个波矢点形成对称保护BIC或者偶然BIC。平带合并BIC比单个BIC更进一步通过调控结构参数把多个BIC在动量空间内推到一起形成一段非常平坦的色散同时还能保持近BIC的高Q特性。这在慢光器件、非线性增强、窄带滤波和低阈值激光方面都很有前景。COMSOL在这类问题里的定位非常明确它不需要你手写麦克斯韦方程的求解器把几何、材料、周期边界和特征频率研究设置好它直接给出复本征频率。实部对应模式频率虚部对应辐射损耗Q因子就是两者比值。远场偏振则需要额外处理用远场域计算远场分量后推断偏振态。整套流程建模和物理图像很紧凑但坑也多我一个个说。1. 先说清楚BIC、平带合并和COMSOL的求解思路1.1 光子晶体平板里的BIC和平带合并光子晶体平板不同于完全二维的光子晶体它在面内是周期结构面外方向通常是z方向则是有限厚度光可以从上下表面辐射出去。正因如此平板里的导模会与自由空间中的辐射模发生耦合形成“泄漏模式”。如果某个模式在某个波矢点上由于对称性匹配不了辐射模或者多个模式之间的干涉让辐射分量完全抵消它就无法往外辐射这就是BIC。平带合并BIC的关键思想是把原本分散在不同波矢位置的多个BIC通过调节某个几何参数比如孔的半径、平板的厚度、占空比让它们在动量空间里靠近甚至重合。一旦BIC合并色散关系会变得非常平群速度很低同时模式的场约束能力很强这相当于在一个宽频谱范围内都拿到了高Q响应。这个“合并”不是碰运气它背后通常是多极子模式之间的干涉条件起了变化。1.2 COMSOL做这件事的计算路径COMSOL里计算这类问题最核心的是“特征频率”研究。你不需要设置入射光源直接求解波动方程的特征值问题得到的就是结构的本征模式。针对周期性结构面内方向使用Floquet周期边界条件通过设置Bloch波矢 $k_x, k_y$ 来扫描布里渊区面外方向要么用匹配层吸收要么在只算能带时用开放式边界近似。得到复数特征频率后[ Q \frac{\mathrm{Re}(f)}{2|\mathrm{Im}(f)|} ]注意COMSOL的特征频率结果表达方式在不同版本里略有差别有的直接给虚部带负号有的给的是正虚部需要根据结果确认一下符号约定。远场偏振计算则是另一回事。它需要在结构外一定距离处做一个积分面由近场数据外推出远场方向的复振幅进而得到两个正交偏振分量的相位关系。BIC的远场强度为零但它周围的偏振分布携带拓扑信息这也是判断BIC性质的重要依据。2. 几何建模与材料参数设置2.1 基础结构和参数选择我复现的是常见的方形晶格光子晶体平板介质材料用硅折射率3.48背景空气平板上是圆柱形空气孔。这类结构在现有文献里已经研究得很透适合拿来做对照验证。归一化参数取法如下晶格常数 (a 1\ \mu m)平板厚度 (t 0.55a)空气孔半径 (r) 作为调控参数基线值取 (r 0.25a)。实际建模时把a、t、r都定义为全局参数方便后续扫描r来观察BIC合并。单位上用微米COMSOL里几何尺寸和波长的单位一定要统一不然归一化频率和实频对应关系很容易出错。建模细节上空气孔贯穿整个平板上下各加一层空气域高度至少留到半个波长以上。如果只算能带不加PML空气域顶部用散射边界条件SBC就可以如果要算远场则需要把空气域顶部改成远场域并搭配完美匹配层PML。这两个方案对应的边界条件不一样建议分开建两个模型或用一个模型做两个研究。2.2 边界条件和网格设置面内方向必须用周期性条件。COMSOL的“周期性条件”节点里有Floquet周期选项需要指定两个周期矢量和Bloch波矢。对方形晶格周期矢量就是 ((a,0)) 和 ((0,a))Bloch波矢设成 (k_x, k_y) 两个参数后面扫描时赋值。网格经验是每波长至少6到8个网格单元。但特征频率计算和频域计算还不太一样你得先估算目标模式的归一化频率。假设目标模式在 (\omega a/2\pi c \approx 0.5) 左右那么介质内波长约为 (\lambda/n)据此设置最大网格尺寸。孔边界和圆角处别偷懒至少要两到三层边界层网格否则孔边缘的电场奇异会污染本征频率虚部Q因子计算直接失真。这里插一句最常被新人忽略的点几何建模不要用完整的三维大结构去算能带。COMSOL的Floquet边界条件只需要建一个原胞unit cell也就是一个晶格周期对应的最小单元算出来的就是无限大周期结构的结果。很多同学一开始按阵列建了几十个孔不但慢而且非周期边界会产生伪模式方向就跑偏了。3. 二维能带计算从单个模式到色散曲线3.1 特征频率研究与波矢扫描实操二维能带图指的是沿布里渊区某个高对称方向比如 (\Gamma)-X扫描 (k_x)画出模式频率随波矢变化的曲线。在COMSOL里操作如下第一步在“研究”里添加“特征频率”研究。求解器设置中要把“所需特征频率数”设成一个稳妥的值比如10到20个。搜索频率范围设置成目标区间比如0到1归一化频率或者对应真实频率范围 0 到 3e14 Hz。这一步不要贪多模式数设得太少会漏掉BIC设太多则求解变慢且容易出现无关模式。第二步添加参数扫描。在“研究”里对全局参数 (kx)以及对应的 (ky)沿 (\Gamma)-X 时 (ky0)做扫描。步距在 (0.01\pi/a) 到 (0.02\pi/a) 之间比较合适粗扫时可以用0.02找到感兴趣区间再用0.005加密。第三步求解完成后通过“全局计算”将特征频率的实部和虚部导出。COMSOL里有“特征频率”这个结果量可以直接取实部、虚部。导出的数据重新整理后用MATLAB或者Python画 (\omega a/2\pi c) 对 (kx \cdot a / 2\pi) 的散点图就是二维能带图。3.2 怎么从能带图里认出平带和BIC平带在能带图上非常醒目——一段频率几乎不随波矢变化的模式看起来像一条水平的线条。但要注意不是所有平带都是BIC。有些平带可能是模式在面内方向几乎没有色散但仍然有垂直方向的辐射泄漏。要确认一个模式是BIC最直接的办法是看它的Q因子如果在某个 (k) 点上Q因子突然飙升到 (10^6) 以上数值上通常只看有限值并且在频率上高于其他模式的辐射连续谱起始位置那它大概率就是BIC。实际操作里有个经验BIC通常伴随模式对称性与辐射模不匹配所以在COMSOL里可以直接画该模式的电场分布。如果是TE类BIC看 (E_z) 分量在孔中心附近是不是呈反对称或者满足某些特定节点结构与周围辐射模的对称性对比一眼就能看出“这个模式对外面的世界无感”。另外在扫描过程中BIC模式的本征频率虚部会趋向于0数值上可能出现跳跃或符号抖动的现象别慌这正是数值精度接近极限的表现加密网格并增大求解精度后虚部会稳定在极小值。3.3 后处理导出与绘图经验COMSOL自带的绘图功能适合快速看场分布但处理大量的扫描数据还是导出到外部脚本更方便。我个人习惯用“派生值 → 全局计算”选择“特征频率”的实部和虚部按参数扫描的每个 (k) 点分别计算导出成CSV。格式上每行对应一个模式和波矢后面接频率实部、虚部。画图时归一化公式建议统一用 (\omega a/2\pi c)这样与文献对比时不需要换算。如果是单元尺寸不是1微米记得换算频率时把晶格常数代进公式。这个看似简单的地方往往是复现文献数据对不上的首要原因我一开始就栽在这。4. 三维能带与品质因子计算4.1 三维能带的扫描设置与数据管理所谓三维能带字面上是频率作为 (k_x, k_y) 两个变量的函数得到的是一个曲面。它对应整个二维布里渊区内的色散关系而不只是某一条高对称线上的曲线。复现这类结果时关键在于扫描策略和数据量控制。理想的方案是对 (k_x) 和 (k_y) 做二维网格扫描常见做法是在 (k_x \in [0, \pi/a])(k_y \in [0, \pi/a]) 的不可约布里渊区内做均匀网格。不做太多模式的话比如只关心频率最低的5到6个带网格取 (20 \times 20) 时一次参数扫描就是400个特征频率求解数据量和耗时都非常可观。因此我建议先用14到16个点粗扫找出BIC所在的位置再在局部区域加密到 (20 \times 20) 甚至 (40 \times 40)。三维能带的可视化在COMSOL里可以直接用“三维绘图组”里画频率曲面。不过我个人偏好把数据全部导到Python里用matplotlib的plot_surface或者scatter颜色对应频率。这样控制坐标轴、标注高对称点都更方便。导数据时注意把 (k_x) 和 (k_y) 以及模式编号都保留下来否则后面做模式追踪时根本分不清哪条带是哪一个模式。4.2 模式追踪与BIC位置的确认三维能带最大的麻烦不是算而是“认”——不同波矢点上的特征频率排序会交换直接连接散点会画出很多跳变的折线。BIC模式在接近特殊波矢位置时虚部迅速变小与其他模式在实频上可能靠近甚至交叉如果不做模式追踪很容易把两条带在图上连错。一个实用的方法是用场分布做辅助识别。COMSOL的特征频率求解结果里可以存下每个模式的电场分布通过检查 (\Gamma) 点附近特定模式的场图记录它的对称性特征比如是否有偶极、四极特征再到其他 (k) 点去匹配相同对称性的模式。另一种方式是改用“按固定顺序的模式编号”配合小步距扫描在模式编号稳定时才计算更密的网格。实测下来粗扫判断大致分布细扫追踪模式效率最高。4.3 品质因子的提取与数值精度控制Q因子的提取本身不复杂公式前面已经给过了。但想要Q因子算得准远比提取过程更费心思。第一本征频率虚部是数值量它对网格疏密非常敏感。当真实Q达到 (10^5) 以上时虚部可能只有实部的十万分之一普通网格根本分辨不出来算出来的虚部会被网格离散误差淹没导致Q因子卡在 (10^3) 到 (10^4) 之间上不去。解决办法是对孔边缘和平板上下表面做局部细网格细化至少让最小网格尺寸达到目标模式介质内波长的1/20到1/30。同时把求解器的相对容差调低比如1e-6。第二PML的存在会影响虚部。如果PML离结构太近它会吸收模式本身的渐近场导致虚部偏大、Q值偏低。经验值是PML内边界到平板的距离至少留半波长到一倍波长的自由空间且PML厚度要大于内部波长的三分之一。这个距离并不是越远越好太远则计算域太大浪费资源太近则精度受损。第三Q因子曲线可以用参数扫描来看BIC合并过程。做法是固定一个接近BIC的波矢扫描空气孔半径 (r)记录模式Q值。合并BIC的特征在Q值曲线上会表现为在某个 (r) 值附近原先偏离的模式Q值集体向无穷大方向拉升而单带Q则对 (r) 的敏感度下降形成一段“平台”。这个趋势比单纯某一组能带图要直观得多。5. 远场偏振计算与BIC的拓扑特征5.1 远场域和PML的配置方法计算远场偏振需要在COMSOL里开启“远场域”。基本思路是在一个封闭边界内部求解近场在边界上把近场做等效源积分外推出无穷远处的远场分布。操作上在物理场设置里右键选择“远场域”指定积分边界和远场计算方向即可。对于光子晶体平板的半空间问题常见的做法是在模型上方加一段空气域空气域顶部设置为“远场”边界再在其上添加PML作为吸收层。为了避免PML反射和远场计算互相干扰PML应放在远场边界的外侧。如果只是想要某个方向的远场也可以在“远场域”节点的设置里限定计算方向比如从 (0^\circ) 到 (360^\circ)仰角方向固定为 (90^\circ)即垂直于平板表面的方向。这里要特别提醒远场偏振计算和能带计算对边界条件的要求相互矛盾——能带计算要周期边界远场计算要开放边界。虽然COMSOL可以在同一个模型中同时设置周期边界和远场域但两者的物理近似不同计算结果容易互相污染。我更推荐的做法是把模型拆成两个研究一个用周期边界算能带和Q因子另一个用较大空气域加PML算远场两者可以共享同一套几何参数。复现文献中的远场偏振图时务必在论文的方法部分确认作者使用的是哪种边界条件否则你算出来的偏振分布可能跟文献差很远的。5.2 提取远场偏振的两种方式远场偏振怎么描述最直接的是看远场电场的两个球面分量 (E_\theta) 和 (E_\phi)。两者的复振幅比决定了偏振椭圆的长短轴方向和旋向。实际操作时COMSOL可以在派生值里对远场电场求分量得到每个方向上的 (E_\theta) 和 (E_\phi) 的实部和虚部。拿到这组数据后偏振态就可以用Stokes参数来刻画[ S_0 |E_\theta|^2 |E_\phi|^2 ] [ S_1 |E_\theta|^2 - |E_\phi|^2 ] [ S_2 2 \mathrm{Re}(E_\theta E_\phi^) ] [ S_3 2 \mathrm{Im}(E_\theta E_\phi^) ]归一化后可以用 (S_1, S_2, S_3) 在邦加球上画点或者在 (k_x)-(k_y) 动量空间里画偏振椭圆。BIC附近的一个显著特征是包围BIC的路径上一圈走下来偏振角的变化总量为 (2\pi) 的整数倍这个环绕数就是拓扑荷。合并BIC的偏振结构中能看到多个拓扑荷相互靠近甚至抵消这是判断是否发生“合并”的最重要远场证据。5.3 计算偏振时的坐标系陷阱偏振计算里最隐蔽的坑是坐标系。COMSOL远场计算的默认投影是基于全局坐标系的球坐标当观察方向不在 (z) 轴上时(E_\theta) 和 (E_\phi) 的定义会随方向变化导致同一个模式在不同方向上看起来“偏振状态不同”。如果你只关心垂直方向上远场比如 (k_xk_y0) 附近这问题不大但如果想做整个远场半球映射建议自己写一小段脚本把COMSOL输出的笛卡尔远场分量(E_x, E_y, E_z)在正确的局部坐标系下投影到垂直于传播方向的平面内再计算偏振态。另外BIC本身的远场强度是零在BIC所在 (k) 点计算Stokes参数会出现 (S_0 \to 0) 而其余分量不稳定的情况。这是正常的绘图时通常把该点留空只画它周围一圈方向的偏振分布。不要试图直接算出BIC点的“偏振”那个量在物理上没有意义。6. 常见坑与排查清单每个项目跑下来我基本都踩过这些坑列出来方便自查症状可能原因排查与解决特征频率求解不收敛或模式数不够求解器搜索频率范围过窄、模式数设置太小扩大搜索区间模式数设到目标数量的1.5倍以上Q因子虚部突然变负或剧烈跳变特征频率编号切换、简并模式混在一起做模式追踪输出电场分布比对对称性能带曲线在某个区域断截波矢步距过大错过了模式交叉点在断截位置局部加密扫描步距远场偏振图严重不对称周期边界或PML设置不对称网格不对称检查周期边界两条边是否严格对应网格是否对称细化网格细化后Q因子反而降低PML离结构太近吸收了渐近场增加PML到结构的间距适当降低PML吸收性能还有一个非常容易忽略的点三维模型里用特征频率研究时默认的求解器是SPOOLES或MUMPS。BIC模式的虚部极小数值上接近零这可能导致特征值求解中实部和虚部精度差异过大。如果发现虚部总在 (10^{-6}) 附近抖动试试把求解器换成PARDISO并打开“舍入误差保护”选项。这个操作我在几个版本上都验证过虽然不太明白底层原理但对结果稳定性有明显改善。网格这一步再强调一次。BIC计算属于“高频高精度”问题一个均匀的粗网格能给出正确的实频但虚部往往不可靠。建议做一次收敛性测试固定波矢把最大网格尺寸从 (\lambda/4) 逐步缩小到 (\lambda/16)看Q因子是否趋于一个稳定值。如果Q值始终随着网格细化而上升说明还没有达到收敛需要继续细化。这个测试很费时间但比跑完一遍才发现Q值全部偏小要高效得多。最后关于COMSOL版本兼容性我能给的经验是物理场设置在不同版本之间差异不大但“远场域”节点和“特征频率”后处理量的名称偶尔会有变化。早期版本里远场电场分量叫ewfd.Efarx等新版则统一了变量命名。不确定时直接用菜单里的“派生值”选物理量不要手敲变量名。整套流程走下来我最深的感受是“计算BIC”这件事本身门槛不高难的是“证明你算出来的是BIC”。能带图、Q因子、远场偏振三样缺一不可它们分别回答了“模式在哪里”、“损耗有多小”、“为什么辐射被抑制”。做复现的时候我也建议你把每一步的数值判据都跟文献对照一次尤其是Q因子的数量级和远场偏振的拓扑荷方向。凡是能对上的后面再调参数做新结构心里就有底气得多。本文还有配套的精品资源点击获取