简介围绕 COMSOL 仿真光子晶体平带合并与连续谱束缚态BIC的研究复现资料面向需要开展二维/三维能带、品质因子及远场偏振计算的光学方向研究人员和高年级学生。内容覆盖从模型建立、参数设置到结果解释的完整思路配套技术博客文章、HTML 说明文档、TXT 笔记及 JPG 结果图共 10 个文件压缩包约 383KB便于按需查阅。已有 243 人学习下载。重点给出光子晶体平带合并的复现方法并结合能带图、品质因子与远场偏振计算结果帮助读者理解 BIC 的形成条件及评估器件应用潜力适合作为课题入门、仿真复现和教学参考资料。 上个月我接到一个复现任务用COMSOL计算光子晶体平带合并BIC要求把二维能带、三维能带、品质因子和远场偏振全部算齐。当时我觉得这活儿并不难无非是搭一个平板光子晶体单元胞跑特征频率扫描再把Q值和偏振后处理一下。真正在COMSOL里过了一遍才发现能带计算从扫k这个环节就开始出幺蛾子后处理庞杂模式时更是容易南辕北辙。我在这里把整个复现过程整理出来给同样准备用COMSOL做BIC和能带计算的人做个参考。文章会按物理图景、建模策略、二维/三维能带实操、Q因子提取、远场偏振计算和踩坑记录来展开所有参数和工作流都来自我实际跑过的算例可以直接抄作业。1. 为什么平带合并BIC值得在COMSOL里折腾1.1 BIC、平带和“合并”到底指什么连续谱束缚态Bound state in the continuumBIC这个名字容易吓到人实际物理图景并不复杂。想象一个水池里有一个点源在振荡正常情况下水波会向四面八方传播能量不断流失。但有一种特殊情况由于波的对称性、拓扑性质或干涉相消某个模式虽然在频率上落在连续辐射谱范围内却始终无法把能量“交”给往外传的波能量就一直锁在结构里。这就是BIC辐射损耗严格为零对应品质因子Q趋于无穷大。平带则是指光子晶体能带中某一段频率对波矢k的变化极不敏感群速度v_g dω/dk接近零。平带意味着光在结构里走得极其慢光与物质的相互作用时间被拉长对非线性增强和传感都很有利。把这两个概念放在一起的“平带合并BIC”是近几年高Q光子晶体器件里很热的一个方向。说的是在一个光子晶体板里通过调节某个几何参数比如双孔半径差、孔位移原本在Γ点附近分开的多个BIC或者BIC与准BIC发生融合合并成一个高阶BIC。合并之后Q因子随动量偏离的衰减律会变慢从|Δk|⁻²变成|Δk|⁻⁴甚至更慢相当于在更大的动量范围内都能保持高Q。这对实际器件非常关键因为真实加工中总会有动量展宽单一的BIC很容易被破坏合并后的BIC则皮实得多。1.2 为什么这类问题适合用COMSOL做我周围不少同行做这类结构习惯用FDTD或者RCWA。FDTD做能带和Q值虽然方便但复本征值提取不够直接精细分辨高Q模式时需要极长的仿真时间和很小的Courant步长RCWA对垂直厚度分层处理很强但要做模式场分布、辐射偏振这类后处理会绕不少弯路。COMSOL的优势在于把建模、能带扫描、本征模式识别和远场后处理放在同一个模型里。特征频率研究直接给出复数本征值Q值计算几乎是免费的Floquet周期边界条件对任意k点都能处理需要的时候还能在边界上做近远场变换提取远场偏振。虽然FEM在三维光子晶体全结构模拟上计算量偏大但利用周期性单元胞加PML计算量是可控的。如果你恰好在做超表面、光子晶体微腔或者BIC激光器件这套流程可以直接复用。2. 建模前必须敲定的几何、材料与求解策略2.1 单元胞、材料和几何参数选择先别急着开新模型把结构参数先钉死后面会省大量时间。我用的是一块悬浮的氮化硅Si₃N₄光子晶体板方形晶格晶格常数a 500 nm平板厚度t 220 nm空气孔半径r₀ 0.25a。这个结构在可见光和近红外波段有丰富的导模材料损耗低加工也很常见。材料设置上如果只看本征频率直接用折射率定义相对介电常数εᵣ n²即可。Si₃N₄在近红外波段n ≈ 2.0玻璃衬底n ≈ 1.45空气n 1。我不建议一开始就把衬底加进去衬底会引入额外的模式杂化和泄漏路径BIC的判断会变浑浊。先把自由悬浮膜体系跑通再按需加入衬底做对比这是最稳妥的路径。几何上只需画一个单元胞不要画整块晶体。平板区域在x、y方向各取一个周期z方向从板表面向外各留出至少一个a的空气层再在空气层外侧放PML。这样既保证了周期性边界条件可以正确设置又给了辐射模式一段“飞出”的空间。2.2 Floquet周期边界与k矢量设置COMSOL里周期边界条件选“周期条件”子类型选“Floquet周期”。它和普通周期边界不同允许边界两侧之间存在一个相位因子用户需要指定布洛赫波矢的分量k_x、k_y。这个波矢分量一定要设成全局参数而不是直接写数字否则后面做k路径扫描时没法参数化驱动。我在模型中定义了两个参数kxg和kyg单位是rad/μm。在Floquet周期边界设置里把波矢分量分别指向这两个参数就行。注意周期边界要同时施加在x方向和y方向的对应边界对z方向不要加周期边界这里后面要放开放边界或PML。2.3 特征频率研究的求解器配置研究类型选“特征频率”这一步是整篇文章的地基。特征频率范围设在150 THz到450 THz大约覆盖近红外到可见光几个导模频段。搜索基准值设在你感兴趣的带边频率附近比如250 THz这比让COMSOL盲扫整个范围高效得多。特征频率数量建议一次算30到50个。不要担心算多因为能带扫描时会有大量高阶泄漏模混进来只有把候选模式都算出来才能在能带图上做筛选。关键是求解器迭代求解器在网格大、模式密集时容易漏根我直接改用PARDISO直接求解器虽然在内存上压力大一点但鲁棒性提升明显。网格这块多说一句光子晶体板的场在板内振荡平板厚度方向至少需要3层网格x、y方向的最大网格尺寸不要超过λᵉᶠᶠ/6其中λᵉᶠᶠ 2π/(k₀n)。我用板厚方向5层、面内最大尺寸60 nm的三角形网格全局自由度在十几万量级单k点求解时间约1到3分钟属于可接受范围。3. 二维能带计算从Γ到M的完整操作链路3.1 k路径扫描的实现思路二维能带的意义在于快速确定模式频率范围、识别带隙和平带位置为三维仿真缩小搜索范围。倒空间路径选Γ → X → M → Γ这三个高对称点就够了Γ (0, 0)X (π/a, 0)M (π/a, π/a)。实际操作时我引入一个路径参数s用条件表达式来定义k_x和k_y。假设s从0到3s ∈ [0,1]Γ → X kxg s * π/a kyg 0 s ∈ [1,2]X → M kxg π/a kyg (s-1) * π/a s ∈ [2,3]M → Γ kxg (3-s) * π/a kyg (3-s) * π/a在COMSOL“全局定义”里把kxg和kyg写成关于s的解析函数然后在研究里对s做参数扫描。这里的坑是扫描步长第一遍先用每段10到20个点跑粗扫找到能带密集的区域后再对对应段加密到每段50个点以上。如果全程都用很细的步长三维仿真时求解时间会成倍上涨没必要。3.2 从一大堆本征模里揪出BIC模式特征频率研究跑完之后第一直觉是直接画频率对s的散点图。但真正跑过的人都知道这个图往往是一团乱麻因为里面混了大量泄漏模和数值伪模。我的筛选顺序是这样的先在结果里按本征频率虚部排序。COMSOL返回的是复频率格式是f i·f_imag对于随时间衰减的泄漏模虚部通常是负的绝对值代表了辐射损耗的大小。普通泄漏模的虚部可能在GHz量级而BIC附近模式的虚部会小好几个数量级甚至接近数值噪声水平。然后再看模式场分布。把电场模画在xz或yz截面上BIC模式的场被牢牢局域在板内向外衰减很快泄漏模则会在空气层PML区域留下明显的传播条纹。两者结合起来判断基本不会认错。有一个经验如果某个模式在k点扫描中虚部忽大忽小那多半是数值问题不是物理BIC。真正的BIC即使在不同网格下虚部有波动量级也始终远低于普通泄漏模。3.3 二维能带里平带出现的判据平带在能带图上非常好认一段模式的频率几乎不随s变化看起来就像一横排密集的点。此时群速度∂ω/∂k趋近于零。做复现时要注意平带往往是“看”出来的但容易误判。比如网格不够密时原本缓变的能带会被离散点搞得像阶梯状看起来像平带。我建议对疑似平带区间加密k点单独扫描一次如果加密后能带斜率稳定在一个很小的值才确认是平带。另一个辅助手段是直接看模式场平带模式在广阔的面内区域有近乎均匀的场分布粒子数密度高局域性强。平带和BIC的关联在能带图上往往表现为平带在Γ点附近出现一个极值点且该点模式的Q值远高于同一频段其他模式。此时就值得把参数变量设回三维模型去算精确的Q和远场偏振了。4. 三维能带计算厚度方向带来的模式本质变化4.1 三维单元胞建模与周期边界二维模型在物理上对应无限厚的介质柱或空气柱光沿着面内传播z方向场分布不变化。真实光子晶体板是有厚度t 220 nm的光在z方向有驻波结构所以三维建模是必须的。把二维几何拉伸成三维时只需在z方向设定板厚板的上方和下方各加一层空气和PML。x、y方向依然用Floquet周期边界z方向用PML吸收层。三维模型的网格压力比二维大得多但因为有周期性和薄板结构网格可以做得比较巧妙。平板区域z方向5层网格面内使用扫描式网格让周期边界两侧节点对齐这会显著降低周期边界带来的数值误差。空气层和PML区域用较粗网格即可总自由度控制在百万以内。4.2 如何区分板内导波模式与泄漏模式三维能带里会混入更多“杂兵”因为z方向自由度增加了。区分方法还是老两样虚部和场分布但第三个维度让场分布判据更直观了。看yz或xz截面如果模式场在PML边界处发生明显畸变或者有出射波特征可以判定为泄漏模式。BIC或导波模式在PML区域场强几乎为零。需要提醒一点PML本身就是人为设置的“数值悬崖”如果模式场延伸到PML边缘说明空气层高度不够应该把空气层从a增加到1.5a到2a再试。还有一个从文献里沿用的判据比较同一频段不同模式的Q值分布。BIC在能带图上是一个孤立的高Q点周围的模式Q都很低在Q对k的曲线上会形成一个尖锐的尖峰。如果Q值尖峰被宽化可能是网格太粗或PML吸收效果不理想。4.3 用参数微调实现BIC合并这是整个复现里最有意思的部分。我的做法是把单孔晶格改成2 × 1超胞里面有两个半径分别为r₁和r₂的空气孔。当r₁ r₂时体系具有C₂对称性Γ点存在对称性保护的BIC当r₂逐渐变化Δr r₂ - r₁从0增大体系的对称性被破坏原来的BIC会移出Γ点另一个准BIC或原本靠近连续谱的模式会向Γ点移动。扫描Δr从0到0.05a每步0.005a在Γ点和Γ点附近几个小k点各做一次特征频率计算。你会看到两条模式频率在某个Δr处发生接近和交换这正是BIC合并的前兆。合并点的典型特征是能带色散从线性变成二次甚至四次接触同时Q值在Γ点附近的平台区间变宽。这个参数扫描比较耗时我强烈建议先用粗k扫描定位到候选Δr区间再在该区间内加密。加密时不仅要扫Δr还要把k点步长同时加密否则观察不到Q平台的展宽效应。5. 品质因子与远场偏振从复数频率到远场涡旋5.1 Q值提取复数本征值的正确解读Q值不需要额外的仿真步骤它直接写在特征频率的虚部里。对每个本征模式COMSOL给出复本征值f_c f_r i·f_i计算公式为Q f_r / (2|f_i|)很多初学者会在这里犯一个错误直接把COMSOL显示的“本征频率”当成实数频率然后单独去另外一个模型里算Q。其实本征频率虚部一直在结果里在“派生值”或“全部特征值”列表里就能看到。但有一个前提Q值可信的前提是模型里有正确的辐射通道。如果不在z方向放置PML泄漏模的能量无法从计算区域流出虚部会被严重低估Q值虚高。PML的位置和厚度我建议设为板面到PML内边界距离1.5aPML厚度0.5a到1a。太薄的PML会反射辐射导致虚部出现振荡。另一个经验是Q值的网格收敛行为。真正的BIC在网格加密时Q值单调上升且不收敛到某个有限值理论上发散如果Q值随网格加密收敛到某个常数那说明它只是高Q准BIC不是严格BIC。做复现时可以把网格从粗到细做三档看Q值走向。5.2 远场偏振计算的二次仿真流程远场偏振不能在特征频率研究里直接算COMSOL的远场节点是在频域分析中配合电磁波、频域接口使用的。我把这个过程拆成两步一次特征频率研究锁定目标模式的频率和k点后再建一个单独的频域研究同一个模型物理场还是电磁波、频域激励给一个弱电流源点偶极子放在板上表面场强最大的位置。频率设定为目标模式的实部频率放宽扫描范围后可以做单频点计算节省时间。启用远场计算时在“电磁波、频域”接口下添加“远场”节点选择模型最外侧的闭合边界作为远场计算边界。这个边界必须包住所有散射体和非均匀结构且最好离结构至少一个波长。计算完成后在派生值里请求远场电场E_θ和E_φ的复数值相位基准点选为结构中心。5.3 偏振奇点与BIC位置的对齐验证得到远场E_θ、E_φ之后可以先算Stokes参数来刻画偏振态S0 |Eθ|² |Eφ|² S1 |Eθ|² - |Eφ|² S2 2 · Re(Eθ · Eφ*) S3 2 · Im(Eθ · Eφ*)用S1、S2可以确定线性偏振方向角用S3/S0确定圆偏振度。BIC的拓扑性质会在远场偏振矢量场上留下明显的印记在动量空间这里的观察角θ、φ近似对应面内波矢方向中偏振矢量场会出现一个涡旋奇点奇点附近偏振方向绕一圈转2π这就是V点。单BIC情形下远场偏振图上只有一个孤立的奇点。当BIC合并发生时两个奇点会先靠近再在某个临界参数上合并成一个复合奇点或者退化为一条闭合的暗线偏振椭圆率变为零的轨迹。观察这个变化过程最直观的方式是在极坐标图上用箭头叠加椭圆率填充色。角度分辨率要够细我在θ方向取了0.25度步长φ方向取180个点才能把奇点看清楚。偏振计算最容易翻车的地方是相位基准。BIC模式的全局相位本身没有绝对意义但远场相位差E_θ相对E_φ的相位是物理的。只要两次计算用同一个相位基准最终偏振图就有意义。6. 复现过程中踩过的坑与可以直接照抄的工作流6.1 四个容易忽略但会浪费两天的隐性Bug第一能带“跳带”问题。特征频率求解器不会按照能带顺序排列本征值所以在画能带图时如果直接把所有模式按频率排序画出来交叉点附近会出现一团乱麻。正确做法是借助模式场分布的重叠积分或对称性分类来连接能带。我在Γ → X路径上输出了每个模式的x/y/z偶极子辐射投影再手动分类能带图立刻清爽很多。第二PML与周期边界的交界处容易出现伪模式。有段时间我在k空间远角处看到一条异常的平坦模式以为是物理上的平带后来发现是PML网格不够平滑导致边界处反射并在PML里形成驻波。解决方法是PML区域用映射网格并让PML内边界紧贴空气层网格不产生网格错位。第三特征频率的虚部符号。COMSOL采用e^{jωt}约定时泄漏模的虚部为负但也有版本或设置下虚部为正。判断泄漏方向不靠符号靠虚部的绝对值是否远大于数值噪声。如果发现某个模式的虚部恰好是零附近的小量不要惊喜先加密网格看看它是否稳定。第四参数扫描时数据集排序会变。在COMSOL里做完参数扫描后“一维绘图组”选择“所有特征值”时曲线的颜色和样式会乱跳。解决办法是在计算后把“全部特征值”数据导出成表格文件用外部脚本按s值和频率排序后再画图或者在主绘图里按s值着色来区分扫描点。6.2 复现工作流总览最后按照我实际跑通的经验给出一个可以复制的分阶段工作流阶段目标关键设置1. 2D粗扫确定模式频段与平带位置粗k步长30个特征频率开Floquet周期边界2. 2D加密扫描识别BIC候选模式加密k点观察虚部数量级筛选局域场模式3. 3D单k验证确认BIC存在与Q量级加PML固定目标k点对照虚部收敛性4. 参数扫描定位BIC合并参数扫描Δr观察Γ点附近能带与Q峰值变化5. 频域远场计算提取远场偏振与涡旋弱电流源激励远场边界距结构1波长6. 后处理画出能带、Q曲线、偏振图Stoke参数计算奇点位置追踪这个流程里唯一需要反复调试的是第3和第4步之间的衔接。3D全参数扫描计算量很大我建议先用第3步的3D模型在几个关键Δr值上验证趋势确认物理规律后再批量扫描。我在最终跑完整个参数区间时用了约60个k点×15个Δr值总计算时间大约一个周末属于可接受范围。这套流程跑顺了之后换个材料体系比如Si、TiO₂或者换成三角晶格、椭圆孔只需要改几何和材料参数绝大部分求解器设置和后处理逻辑都可以复用。量子点增益、非线性增强或传感方向的应用也都可以从这套能带和Q值框架往下延伸。本文还有配套的精品资源点击获取