做超表面或者光子晶体的同学最近大概率绕不开一个词BIC连续谱束缚态。我在给一个周期硅柱结构算能带的时候发现Q因子和远场偏振投影是验证BIC最直观的方式但COMSOL里这三件事想要串起来跑通中间其实有不少坑。这篇把我实际操作的整个链路写清楚包括能带计算、Q因子提取、远场偏振投影最后落到怎么把两个BIC在k空间合并这个具体任务上。1. 为什么说COMSOL是算BIC最顺手的工具1.1 BIC是什么一句话版本连续谱束缚态字面上看很奇怪——一个态被嵌在允许辐射的连续谱里却完全不往外辐射。类比说就是演唱会现场大家都站在开放区域里听歌声音四面八方传出去结果有一个人原地站着周围形成了一道看不见的隔音墙他的声音一点没漏出去。在光子学里这个隔音墙来自结构的对称性保护或泄漏通道之间的干涉相消。对称保护BIC最简单的例子一个上下对称的介质平板某个模式的本征磁场分布正好和辐射通道的对称性完全正交远场积分为零泄漏率就是零。另一种是非对称点off-(\Gamma)的偶然BIC它靠两个辐射通道的远场干涉相消实现。实验上我们一般看Q因子——共振峰的尖锐程度。BIC本身Q趋向无穷大但在真实结构中只能在某个k点附近看到Q疯子一样往上涨这就是准BIC。1.2 合并BICmerged BIC是什么为什么值得专门做单独一个BIC在k空间是一个点它很脆弱。实验上以有限数值孔径的光入射入射角稍有偏移k就离开BIC点了Q立刻掉下来。合并BIC的思路是通过调节结构几何参数让两个BIC在动量空间逐渐靠近、最后合并成一个点。合并之后Q因子在k点附近的下降变得极其平缓也就是说对入射角的容忍度大大提高。这个性质对实验极友好。你不用把光对准到纳米级别的角度误差也能测到很高的Q。我在仿真里复现过这个过程同一个结构不合并前(\Delta k0.01,\mu m^{-1})就掉一半Q合并后同样(\Delta k)下Q只下降不到1%。1.3 为什么要用COMSOL而不是FDTD或RCWA算能带和Q因子常用的有三条路RCWA严格耦合波分析算反射、透射谱很快做参数扫描很爽但物理场可视化弱BIC附近的相位和偏振信息不好提取。FDTD时域有限差分直接看时域衰减曲线拟合Q直观但周期性结构的Bloch边界条件在斜入射时比较麻烦而且计算资源消耗大。COMSOL特征频率法直接在频域中求复本征频率一步到位得到Q因子同时能拿到完整的近场分布、能带色散以及远场辐射方向图非常适合分析BIC的模式物理。我自己的经验是COMSOL的Floquet周期边界特征频率研究是辨识BIC最快的方式。你扫一条k路径看到某一支模式虚部断崖式归零那基本就是BIC了。2. 仿真初始化几何、材料、边界条件的配置2.1 结构建模一个晶胞就够以最简单的方形晶格介质柱阵列为例。晶格常数(a800,\text{nm})柱高(h400,\text{nm})柱半径(r200,\text{nm})。在COMSOL的几何里建立一个长为(a)、宽为(a)的2D矩形作为晶胞底面。在矩形的中心画一个圆形用这个圆从矩形中减去得到一个中间有孔的矩形区域——这个孔就是柱体的位置。用Extrude拉伸这个二维截面高度设为(h)得到柱体和背景介质区域。上下再加上空气层空气层高度建议至少(2\lambda)这样PML才不会因为离结构太近而扰动近场分布。以上是介质柱刻在薄膜上的常见结构。如果你做的是柱体立在衬底上那就把背景介质换成衬底材料就行原理一样。提示在工作平面上画几何时COMSOL默认的尺寸单位是m记得统一改成nm或(\mu m)后续设置材料折射率时才不会搞混量纲。2.2 材料参数与色散材料用最简单的无损耗介质即可背景用SiO2折射率(n1.45)柱体用Si折射率(n3.45)。如果研究的是可见光波段Si的色散不可忽略要使用COMSOL材料库中的Palik数据。但如果是近红外比如通信波段1550nm左右在单一目标频率附近用一个恒定折射率造成的误差很小。在BIC合并的研究中折射率的微小差异通常只造成BIC合并点的移动不会改变合并机制本身。先用常数折射率跑通流程再逐步加上色散是稳妥的做法。2.3 Floquet周期边界条件和布洛赫k扫描在晶胞的四个侧面分别设置Floquet周期边界条件。这是算能带的核心。两个互相平行边界成对设置周期矢量分别是((a,0))和((0,a))。COMSOL的Floquet边界条件里可以直接设置波矢(k)的分量。这个(k)就是Bloch波矢。当(k_xk_y0)时对应(\Gamma)点(k_x\pi/a, k_y0)时对应X点设置好之后可以通过参数化扫描来扫整个布里渊区。我习惯在模型的全局参数里定义a 800[nm] r 200[nm] h 400[nm] kx 0 ky 0然后在Floquet边界条件的波矢设置里直接填(k_x)和(k_y)。这样后续做扫描时只需要把(k_x)、(k_y)作为参数就行了。2.4 网格划分的实际策略网格是精度和速度之间的钥匙。周期性超表面的结构相对简单但要注意几个关键点结构内部的网格不能太稀柱体内部和柱体周围的折射率变化区域决定了场的正确分布至少要满足每个波长区域内6个二阶单元。这里由于介质折射率是3.45实际波长为(\lambda/n)所以网格要加密到空气中的1/3。空气区域使用较粗网格远离结构的区域场变化平缓可以用最粗的网格大大节省自由度。PML区域用扫掠网格PML内部通常使用扫掠网格Swept层数建议5层以上每层厚度递增。实际经验是使用物理场控制网格的常规Normal预设再手动在柱体区域加一个细化的尺寸子节点就足够让Q因子从1e3级别的精度提升到1e5级别了。3. 能带计算从特征频率到色散曲线3.1 特征频率研究的设置这是核心一步。物理场接口选择电磁波频域ewfd研究类型选择特征频率。在代码实现上COMSOL实际上会求解[ \nabla \times \mu^{-1} \nabla \times E - \omega^{2}\epsilon E 0 ]特征频率研究解出的是谐振频率的平方对应一系复特征值(\omega \omega_r i\omega_i)实部就是模式的振荡频率虚部代表衰减或增长。对于BIC附近的模式虚部绝对值会非常小甚至接近数值噪声水平。研究设置里需要指定所需特征频率个数比如6个让求解器一次性找出前6个最低频的模式。在周期性结构中特征频率的数目和你设置的周期边界条件严格对应——每个k点集合上一个固定数目的带。3.2 k路径扫描的两种做法能带图需要沿布里渊区的边线扫描(\Gamma-X-M-\Gamma)。第一种做法是在COMSOL里加参数化扫描研究节点扫(k_x)、(k_y)。但这种方式会生成一个巨大结果数据集非常占内存。我常用的做法是写一个扫参脚本支持LiveLink for MATLAB或直接用COMSOL的Java API。不过对于大多数场景直接在参数化扫描里扫描就够用了。参数设置原则是(\Gamma)点(k_x0, k_y0)X点(k_x\pi/a, k_y0)M点(k_x\pi/a, k_y\pi/a)间隔取100个k点。线路径太稀疏会漏掉BIC附近的窄特征太密则计算量成倍上升。100-150个k点是我实际用的平衡值一个波段扫描大概需要10-20分钟。3.3 识别BIC看虚部和线宽扫完k路径后处理里画能带图。直接把特征频率的实部作为纵轴k位置作为横轴可以看到几个带。这时BIC并不直接可见因为能带图上它是一条普通的色散曲线混在众多模式里。识别的关键在虚部。在结果里创建另一个图把特征频率虚部的绝对值或Q因子作为纵轴k位置作为横轴。BIC附近的模式会显示一个窄而尖锐的峰——Q因子突然飙升几个数量级。这就是准BIC的识别标志。我在长期实践里发现一个规律在能带图上两支模式在某个k点交叉交叉附近Q因子同时飙升通常意味着这里存在两个偶然BIC或者对称保护BIC与偶然BIC的耦合。这种情况就是合并BIC的前兆。3.4 导出能带数据画图COMSOL自带绘图工具可以直接出图但是想发表级的效果最好把数据导出来用Python或Origin处理。导出步骤在结果节点下创建一维绘图组。添加全局绘图选择特征频率数据集。手动勾选数据系列里的全部频率。右键导出格式选CSV。导出的CSV中包含特征频率实部、虚部以及对应扫描参数的k值再在Python里画出干净的能带图即可。4. Q因子计算的几道坎4.1 Q因子的公式与COMSOL实现Q因子的标准公式是(\omega_r/\Delta\omega)但对于本征模式计算更直接的方法是用复本征频率[ Q \frac{|\omega_r|}{2|\omega_i|} ]在COMSOL中特征频率求解结果是一个复数单位Hz因此可以直接在派生值→全局计算中设置表达式-real(freq)/(2*imag(freq))就算出Q因子了。注意这里的freq是COMSOL内置变量表示特征频率的复值。4.2 为什么BIC的Q因子发散但数值上却不是无穷大理论上BIC的(\omega_i 0)Q无穷大。但在实际仿真中你永远得不到真正的无穷大。原因有三层数值离散误差网格把连续的物理空间离散化了虚部的最小值被网格尺寸和单元阶次限制。你的网格越密虚部下限越小。PML反射即使PML设计得很好仍然有极微弱的反射这会造成泄漏模式的虚部偏大。求解器容差特征值求解的迭代收敛阈值也会给虚部一个底噪。所以如果你看到某个k点Q的计算值是1e8那基本已经到数值精度的极限了。别奢求它变成真正的无穷大关键看它相对周围模式是否高出几个数量级。4.3 网格和PML对Q因子计算精度的影响这是最容易踩坑的地方。我在算一个薄膜超表面时最初用默认PML3层比例1算出的Q因子在BIC附近只有1e3。后来仔细检查发现是PML参数设得不对导致泄漏通道被PML人为堵死或者放大了。正确的PML设置在结构上下方各加一块PML区域高度建议大于1倍波长。PML类型选笛卡尔吸收方向设为z方向。为什么要笛卡尔因为在周期平板结构中周期方向x、y的场分布已经由Floquet边界条件约束PML只需要在z方向吸收辐射使用笛卡尔坐标系下的各向异性PML就能准确做到。PML层数建议5-8层并且打开在PML中使用扫掠网格选项。PML到结构的距离至少保留一个波长的空气间隔。另一个影响Q精度的大因素是网格对称性。BIC附近模式往往具有较强的对称性如果网格不对称会人为引入对称性破缺导致虚部偏大Q被低估。因此我在划分网格时会使用贴片旋转法保证结构的旋转对称性在网格中也被保留。简单做法是在网格节点设置中打开对称检测选项或者对半格结构用镜像网格划分。4.4 用参数扫描追踪Q因子的峰值在识别BIC之后通常还需要看Q因子随某些结构参数的变化这就是合并BIC研究的关键步骤。比如扫描柱体半径r在全局参数里把r设为变量。添加参数化扫描节点扫r从180nm到220nm间隔5nm。在扫描中加入k点扫描比如在Γ点附近取30个k点。全局计算每个(r, k)组合下的Q因子。这里要注意参数化扫描产生的数据量很大。我的做法是在研究设置里勾选在每次扫描迭代后生成数据集让每次扫描结果独立存储然后在后处理里用派生值分组挑选特定参数的值避免数据集之间互相干扰。扫描完了之后把Q作为r和k的函数画出来你会看到Q的峰值在某个r处突然变得更尖锐、更窄——这正是两个BIC合并的signature。我自己做的时候观察到Q曲线从双峰两个BIC各自独立变成单峰但峰值更高BIC合并这个转变很直观。5. 远场偏振投影看看BIC往哪辐射5.1 为什么要看远场偏振BIC的拓扑指纹能带和Q因子能告诉你BIC的位置和Q值但不足以确认BIC的物理机制。远场偏振投影是一个重要的判据。对称保护BIC的远场偏振场在动量空间呈现涡旋结构——远场偏振矢量在k空间围着一个点转一圈会有(\pm 2\pi)的相位/偏振旋转这就是偏振涡旋。这个拓扑性质的实际意义在于如果远场偏振投影图上有清晰的偏振涡旋中心那几乎可以确定对应位置存在BIC而且可以通过涡旋的缠绕数来判断是何种BIC。5.2 COMSOL中设置远场域Far-Field Domain远场计算的思路是先求解近场然后把近场边界上的电场分布当作等效源再来一次瑞利-索末菲衍射积分得到无穷远处的辐射场。在COMSOL里比较简单的方式是在PML外面再加一块足够大的空气域。设置散射边界条件然后勾选远场域选项。在结果节点添加远场辐射模式绘图直接画出来。但这里有个局限性COMSOL内置的远场绘图通常是直接在给定角度(\theta, \phi)上显示远场强度(|E|^2)如果想要得到偏振投影比如把远场电场矢量映射到波矢空间的(k_x, k_y)平面内置功能就不够灵活了。我的做法是导出PML内边界的近场数据用外置脚本做远场积分。5.3 导出近场数据到Python做远场变换把PML内边界上的电磁场导出后用Python进行远场近场变换。核心公式为[ E_{ff}(\theta, \phi) \frac{ik}{2\pi} \int\limits_{S} \left[ E_t - \hat{r}(\hat{r}\cdot E_t) \right] e^{-ik \hat{r} \cdot \rho} dS ]其中(E_t)是边界面上的切向电场(\hat{r})是远场方向的单位矢量(\rho)是位置矢量。导出数据方法在结果里创建二维剪切图选PML内边界。显示Ex, Ey, Ez三个分量每个分量单独创建。用导出节点的数据功能保存为文本文件。导入Python按公式做积分得到(E_\theta)和(E_\phi)。需要注意导出的近场数据包含Floquet周期条件留下的相位因子(e^{ik\cdot r})在做远场积分时要把这个相位因子正确包含进去否则远场分布会错。5.4 绘制偏振投影图得到远场(E_\theta, E_\phi)之后把远场方向映射到波矢空间。在平面波垂直入射附近(\theta)和(k_x, k_y)的关系近似为[ k_x k_0 \sin\theta \cos\phi,\quad k_y k_0 \sin\theta \sin\phi ]在(k_x)和(k_y)网格上计算每一点的偏振椭圆轴比ellipiticity和长轴方向。偏振椭圆的长轴方向可以做为一组短线或箭头绘制在2D图中。BIC对应的就是一个偏振涡旋奇点箭头方向绕点旋转360度总缠绕数为(\pm1)。我这里给个Python画法的思路import numpy as np import matplotlib.pyplot as plt # ex, ey, ez: 近场导出数据(二维数组)这里直接假设已有 # theta, phi: 远场方向网格 # Etheta, Ephi: 远场电场分量 # 计算偏振长轴倾角 psi 0.5 * np.angle(Etheta / Ephi) # 画在kx-ky平面 kx np.sin(theta) * np.cos(phi) ky np.sin(theta) * np.sin(phi) plt.quiver(kx, ky, np.cos(psi), np.sin(psi), colork)这个方法做出来就是一张偏振矢量场分布图。把(|E|^2)作为背景颜色再叠加上偏振箭头BIC附近的偏振涡旋结构一目了然。这张图是论文里最具有说服力的BIC证据之一。它的好处是哪怕你的Q因子计算精度不够高比如某些杂散模式干扰远场偏振涡旋的拓扑性质也不容易被数值误差抹掉所以可以用来可靠地验证BIC。6. 合并BIC的几何参数优化与后处理6.1 参数扫描目标和BIC合并的判断合并BIC的核心调节参数通常是晶格常数a、柱半径r、柱高h三者之一。不同结构会不同需要根据能带上的两支BIC随参数的变化趋势决定。具体操作是固定一个参数比如r扫描另一个参数比如柱高h同时保持k扫描不变。观察两支BIC的位置在k空间是往一起靠还是互相远离。如果随着h增大两支BIC逐渐靠近最终在某个临界h处合并成一支那这个临界h就是合并BIC的工作点。在COMSOL里的实现是把两个参数同时放在参数化扫描中循环但要注意顺序先扫h再在内部循环扫k这样后处理时以k为横轴、h为颜色变量来画Q因子图合并的过程会非常直观。6.2 合并之后看什么Q因子在k空间的分布合并BIC是否成功不能只看能带上重合还要看Q因子的分布。两个独立的BIC合并前Q在k空间呈现两个独立的高峰峰与峰之间的Q相对较低。合并后两个峰合而为一并且峰宽变宽、峰顶Q变得极高。我用一个实际的例子来说明。一个晶格常数为760nm、半径180nm、高度380nm的结构两个准BIC分别位于k点的(-0.02)和(0.02)处Q峰值各为5e4。逐步增大柱高到420nm两个峰逐渐靠近最终在h410nm时完全合并成一个峰中心处Q达到9e7而且峰的半高宽从0.04增大到0.12。也就是说入射光的角度容忍度增加了三倍。这张图出来后基本就可以确认合并BIC的实现了。6.3 验证远场偏振投影在合并后的一致性合并BIC的远场偏振投影也有一个明显特征原本两个分离的偏振涡旋合并成一个单一涡旋拓扑电荷从二分之一对变成单一的(\pm1)。在远场图中你能清楚地看到两个涡旋在参数扫描中逐步靠近、最后融合。这时候我用导出的近场数据做了远场变换观察到的就是偏振涡旋的合并。这算是对合并BIC的交叉验证。这个验证不要省。如果只靠Q因子图有时数值误差也能造成类似的趋势但偏振涡旋的拓扑结构变化是特征性的不依赖数值精度。6.4 容易忽略的坑参数扫描时的简并模式在COMSOL特征频率计算中有很多简并模式频率相同或极接近的模式会造成特征值的排序不稳定。比如某个k点附近有一个模式虚部特别小准BIC另外还有一个模式虚部较大两个模式频率接近。求解器可能在某次扫描中把这两个模式顺序搞反导致Q因子曲线跳变。解决方法有几种增加特征频率求解数目比如从6个增加到10-12个让准BIC附近的所有模式都被找到然后再在后处理中按频率排序。用模式跟踪功能设置好参考模式后COMSOL会沿扫描参数追踪同一个模式。手动检查所有模式的近场分布确定哪个是你关心的目标模式然后只提取那个模式的指数。实话说模式跟踪选项在有些版本的COMSOL里并不稳定特别是在强烈简并的情况下。我更多使用的方法是在扫描前先对目标模式的近场分布做一次确认记住它的空间分布特征然后在后处理时根据场分布筛选模式而不是光看频率顺序。6.5 资源消耗与求解时间优化最后聊一下算力问题。一个800nm晶格、5nm精度的参数扫描k点150个每个k点解10个特征频率单个k点的求解时间大约10秒整体扫描一次需要25分钟左右。如果再加上参数扫描18组参数总耗时就是7个多小时。这是很现实的时间成本。优化技巧先用粗网格跑一遍k路径确定BIC的大致位置和参数方向。在BIC附近加密k点在其他区域用稀疏k点。关闭不需要的模式数目如果只关心最低频两支模式求解数目设为4就够。使用MUMPS直接求解器对于中等规模周期结构比SPOOLES更稳定速度也更快。参数化扫描时尽量使用单机多核并行COMSOL的扫描节点天然支持并行计算设置好并行核数能明显加快速度。提示在参数扫描时如果同时扫描k点和几何参数把几何参数设为外层循环、k点设为内层循环可以充分利用不同参数组合下的初始解加速COMSOL会用前一个解作为新参数的初始猜测整体时间能缩短30%。尾声我对这套流程的实际体会把这套流程走通之后最大的感受是能带计算、Q因子和远场偏振投影绝不是三个独立的任务。能带告诉你BIC在哪里Q因子告诉你BIC的泄漏有多低而远场偏振投影告诉你BIC为什么不会泄漏——三者的信息拼在一起才构成了对BIC的完整物理图像。如果你也打算在COMSOL里做BIC相关的研究我的建议是先把单个BIC的能带和Q因子算明白再往合并方向扩展。合并BIC虽然是热点但对数值计算的要求更高PML设置、网格收敛性、模式跟踪这三关缺一不可。最后分享一个写论文的小技巧远场偏振投影图不要用COMSOL默认的灰度图把偏振箭头叠加上去之后再把Q因子等高线叠在背景上三者叠加的图在审稿人眼里会显得信息量极大而且物理图像非常完整。我第一篇BIC论文里用了这张图审稿人几乎没有在BIC的验证上提任何问题。