做光子晶体仿真的人，这几年绕不开一个词：BIC（Bound states in the continuum，连续谱束缚态）。尤其是Merging BIC这个概念火了之后，做激光、超表面、非线性器件的人都在往这个方向靠。我这次用COMSOL跑了一套光子晶体平板的仿真流程，目标很清晰：通过结构调整把普通BIC变成合并BIC，同时把三维能带、Q因子、远场偏振这几件事一次性算明白。这篇文章就按这套流程的先后顺序拆开讲，模型怎么建、能带怎么扫、Q怎么算才靠谱、远场偏振图怎么解出来，以及中间踩过哪些坑，全部摊开说。适合正在做光子晶体、超表面、拓扑光子学研究的朋友，尤其是有COMSOL基础但被BIC类计算折磨过的人。

先说清楚这套仿真的物理背景和计算目标，再往后谈操作才不糊涂。

1. 先想明白这几个物理图像

1.1 光子晶体和能带

光子晶体说白了就是介电常数周期排布的微纳结构。电子在周期势场里有能带，光子在周期介电结构里同样有能带，这个类比是整类器件的理论基础。COMSOL算光子晶体能带，本质就是在周期边界条件下解麦克斯韦方程组，得到布洛赫波矢k对应的本征频率ω(k)。跟电子能带的算法思路几乎一样，只不过这里求解的是电磁场而非薛定谔方程。

一个典型的介质平板光子晶体，会同时存在导波模式和辐射模式。导波模式局域在结构内部，随着距离指数衰减；辐射模式则与外部自由空间耦合，能把能量带到远场。在动量空间里，光子晶体的结构周期性会把能带折叠进第一布里渊区，于是本来在光锥外的导波模式也可能被折叠到光锥以内，跟连续的辐射模式共存于同一个动量位置。这种导波与辐射连续态在同一个（k, ω）区域里共存的局面，就是BIC出现的先决条件。很多人一开始不搞明白这个，上来就扫能带，看到一条带就以为是BIC，其实大概率只是普通导波。

1.2 BIC到底是什么

BIC这个名词看起来很高深，物理图像其实很直观。它就是一个包含在连续辐射谱内部的束缚态，本征模式本身的频率落在辐射连续谱范围内，理论上却不向外辐射，Q因子趋向无穷大。听起来矛盾，但没有违背物理规则，因为它和所有辐射通道之间的耦合恰好被干涉相消了。

光子晶体平板里最常见的BIC是Γ点处的对称性保护BIC。在Γ点，模式具有反演对称性，而外部辐射场具有另一种对称性，两者对称性不匹配，辐射通道就被强制关闭。另一类叫偶然BIC（accidental BIC），出现在k不等于零的位置，靠参数调谐让两个辐射通道的泄漏相互抵消。这两种BIC的本质都是远场辐射通道的耦合为零，但在动量空间里的位置和鲁棒性差别很大。

有一个容易混淆的点：BIC不是缺陷态，不是某个局部结构产生的局域态，而是整个周期性结构的延伸态，它的能量完全局限于平板内部，但模式在空间上依然是布洛赫波扩展的。这一点在做COMSOL仿真时很重要，因为边界条件、网格尺寸、PML处理稍有不对，BIC就会被数值噪声“打开”辐射通道，Q因子变成有限值，看起来就像个高Q共振，而不是真正的BIC。

1.3 Merging BIC的价值

单个BIC在理论上是完美的，但实际加工出来的结构总有尺寸误差、侧壁倾斜、材料吸收。一旦结构对称性被破坏，对称性保护BIC立刻就变成有限Q模式。这导致BIC在实际器件中非常脆弱。Merging BIC的思路，就是通过调节几何参数，把动量空间不同位置的全新BIC（某种意义上是“可移动”的BIC）移动到同一位置，合并成一个增强了拓扑特性的BIC。

为什么合并以后的BIC更鲁棒？文献里通常解释为：合并后的BIC在动量空间里对应一个更高阶的偏振涡旋，结构微扰对它的影响是二次幂关系，而普通BIC是一次幂关系。反应在宏观结果上，就是你制造误差变大，Q因子依然能维持在十万量级甚至更高。因此Merging BIC在单模激光、高灵敏传感、窄带滤波方面特别有价值。

2. 在COMSOL里搭一个可信的光子晶体平板模型

2.1 几何参数与材料设定

实际建模时我采用一个典型的介质柱光子晶体平板：硅柱阵列放在一个薄的介质平板上，或者直接在硅平板上刻空气孔，两者都可以，仿真设置几乎相同。下面以一个通用结构为例：晶格常数a = 800 nm，背景平板是硅，折射率n = 3.48，平板厚度h = 500 nm，圆柱半径r从200 nm到280 nm范围做参数扫描。

几何建模尽量用参数驱动，不要直接在几何里填死数字。因为后面做Merging BIC调控，必然要扫描某个结构参数，比如半径、孔偏移量、椭圆率。把所有尺寸都设成全局参数，扫参数时只需改一个变量，模型重建自动完成。我习惯在参数表里一次性定义a（晶格常数）、r（柱半径）、h（厚度）、g（PML与结构距离）等，并且用r/a这种无量纲组合做扫描，这样换材料换波段时所有几何自动缩放，不用重做模型。

材料参数要特别注意色散问题。对半导体材料，仿真波段较窄时可以忽略色散，固定折射率即可。如果要做可见光到近红外的宽谱能带，建议加入材料色散数据，否则二维能带和Q因子趋势都会失真。常见做法是在组件内定义材料函数，允许分段插值。

这里还有一个重要经验：算BIC相关问题，材料损耗先设为0或很小的值。如果一开始就加入材料吸收，特征频率虚部会被材料吸收主导，辐射通道贡献的信号完全被淹没，Q因子分析根本没法做。分开两步走：先做无耗仿真看清辐射Q，再单独估计材料Q对总Q的影响。

2.2 Floquet周期边界与开放条件

单元胞的四条侧边设置Floquet周期性边界条件（在COMSOL里叫Periodic Condition或者Floquet Periodicity），并给定布洛赫波矢k = (kx, ky)。这里有个细节：Floquet边界要求两个相对边界上的场满足相位关系，也就是u_dst = u_src * exp(i * k · r)。其中k的取值直接决定你算的是哪个布洛赫模式。

很多新手在这里犯的典型错误是：把周期性边界条件设为零相位（即kx=0、ky=0），以为这就是Γ点。实际上设置波矢分量时，需要在周期边界条件属性里填写kx和ky的具体表达式。我直接把kx和ky定义为全局参数，在扫描研究特征频率时用“辅助扫描”扫这两个参数，每扫一组(kx, ky)就能得到一个本征频率列表。

光子晶体平板是开放结构，上下表面不能设成完美电导体或完美磁导体，否则辐射通道被人工关闭，算出来的模式全是“伪BIC”。正确做法是上下边界用散射边界条件或完美匹配层。做Q因子提取时，强烈建议用PML而不是简单的散射边界条件。本质上，只有PML才能可靠地吸收辐射场，从而在复数频率中准确提取出辐射损耗虚部。

PML的设置也讲究：厚度至少要覆盖一个波长，与结构表面留空约半个波长的间隙，PML内部网格要各向异性拉伸。COMSOL自带PML坐标缩放参数，通常依赖内置设置就行，但间隙层和厚度这两个值别省。我默认取h_gap = 0.5 * λ0，h_pml = 1.5 * λ0，精度和计算量都比较平衡。

2.3 网格剖分：决定仿真实力的地方

对于三维光子晶体计算，网格往往是大坑。BIC类模式对网格非常敏感，因为你需要精确解析场分布在周期性单元内的细节，特别是结构边缘附近的电场增强区域。我用的是物理场控制网格基础上手动加密：结构内部与平板区域的网格最大尺寸设定为a/8；结构锐利边角处用角部细化；PML区域用映射网格拉伸。

计算Q因子时网格不足的直接表现是：本应趋于无限的Q因子，算出来却只有几百、几千，而且随着网格细化还会系统性地缓慢上升。这是典型的未收敛症状。碰到这种情况，先固定其它所有条件，把网格尺寸从a/16、a/20、a/24各自算一遍，看Q值变化趋势。如果Q与网格尺寸基本无关，说明结果收敛了；如果Q一直在涨，你的结果不能直接引用。

经验值上，三维全波仿真中，BIC附近模式的网格密度通常是普通能带模式的2~4倍。一个900 nm尺度晶胞，网格数量大概在几十万自由度的规模，单次特征频率求解一两分钟能完成。这也是为什么后续做全布里渊区扫描时，计算量会非常巨大，所以网格质量和扫描步长之间要做一个妥协。

3. 三维能带扫描：从一条线到一张面

3.1 高对称路径扫描（一维能带）

先不要直接上二维扫描。我的习惯是，先沿布里渊区高对称路径做一维能带扫描，把所有模式的行为摸清楚，再开展全布里渊区的二维扫描。高对称路径在正方晶格中是Γ(0,0) – X(0,π/a) – M(π/a,π/a) – Γ(0,0)。这个路径上的模式身份比较清晰，可以快速确认目标频率附近有哪些导波带、有哪些交叉点、有几个BIC。

操作方式：研究步骤添加特征频率，把Floquet边界中的kx、ky用路径公式参数化。比如从Γ到X，令kx = π/a * s，ky = 0，s从0到1扫描。这里的s用辅助扫描的索引参数，步长取20~30步就够看到趋势。

关键点是，一维能带扫描时要把“特征频率搜索区域”设置得窄一点。如果求解器搜索范围太宽，比如从0到几个THz，它会返回一大堆高次模式，杂模非常多，无法区分哪些是物理模式哪些是数值模式。我通常先把目标频段估计好（用带结构的平均折射率公式估算），再设置搜索中心在目标频率附近，频率跨度设为带宽的1~2倍。

一维能带图出来后，判断BIC的线索有三个：模式在Γ点出现退化；能带与连续辐射态的色散交叉；模式的本征场在平板内部高度局域，外部场几乎没有泄漏。通常配合Q因子计算一起看。

3.2 全布里渊区扫描与三维能带面

一维能带确认物理图像后，再做二维扫描，得到的就是三维能带面：ω = ω(kx, ky)。标题里的“三维能带研究”说的就是这个——不是三维实空间模型（模型本来就是三维的），而是把动量空间两个维度都展开，看频率面在(kx, ky)空间中的分布。

COMSOL的辅助扫描支持嵌套循环，把kx定义成-pi/a到pi/a，ky同样定义成-pi/a到pi/a，步长建议取20×20到30×30。注意这里的扫描步长直接决定了能带面的分辨率。Merging BIC分析中，两个BIC在临界参数附近会靠得很近，分辨率不够就看不出合并过程。基础筛选用25×25就够，参数精调阶段再用40×40。

这种扫描的计算量非常大。25×25就是625个特征频率求解，每个求解又有好几十个本征模，总体规模和一次中等的三维电磁场参数扫描相当，单机要算几小时到十几个小时。建议先把几何网格适当放宽到a/12得到一个粗结果，锁定模式行为后再在关注的动量区域局部加密扫描。这个局部扫动量空间的技巧特别省时间，尤其在参数临界点附近追踪BIC位置时。

算完后把所有本征频率按(kx, ky)顺序导出，在外部脚本里重新排成二维频率矩阵。COMSOL自带的绘图也可以直接画三维曲面，但数据量太大时交互会卡。我通常把数据导出为表格后交给数据处理脚本，绘制ω(kx,ky)曲面和等频率图，等频率图在判断能带分离、连续谱交叉时更直观。

3.3 模态追踪，防止模式跳变

这是整个流程中最容易出冷汗的环节。特征频率求解器在每个k点返回多个特征值，顺序是按频率大小排的，所以模式会“跳线”——同一支物理模式在不同k点可能是第5个、第6个、第7个特征值，直接按索引画能带会得到一堆断线和跳跃。这个问题在Merging BIC研究中尤其致命，因为BIC往往处在能带交叉点附近，模式高度简并。

解决思路有几个层面。第一，在同一条能带曲线上做扫描时，用“前一个k点的解作为当前k点的初始值”。COMSOL特征频率求解器允许启用“使用上一个解”作为初始猜测，这样能大大提高模态追踪成功率。第二，提取每个特征值的参与因子和电场能量占比，数值假模参与因子明显偏低，可以直接剔除。第三个方法，凡是平面上有交叉危险的地方，用更小的k步长加密，手动分辨模式轨迹。

我踩过一次很深的坑：二维扫描结果看起来能带图非常漂亮，但Merging BIC的“合并”现象始终不明显。排查才发现，模式追踪在某个参数突变点上发生了跳线，我追踪的是一支残余的普通导波，不是BIC。从那以后我所有参数筛选都先跑一维路径验证，再跑二维。

4. Q因子计算的完整链路

4.1 为什么用复数特征频率

Q因子的定义是谐振腔存储能量与每个周期损耗能量的比值，数学上可以写成ω0/Δω。在COMSOL的特征频率研究里，求解器直接给出复数的本征频率，实部就是共振频率，虚部代表损耗率。对无源无增益系统，虚部为负值，Q因子公式是：

Q = Re(ω) / (2 * |Im(ω)|)

没有PML的问题在于，开放式结构的辐射损耗没有物理机制被吸收，求解器把表面看作反射边界，虚部被严重低估，Q因子虚高。所以前面才强调PML不是可选项，而是计算Q因子的必需项。PML把辐射波吸收掉，等效于开放空间，虚部才包含真实的辐射损耗。

对比一下常见的宽频法（S参数扫描再Fano拟合）和特征频率法：宽频法需要先扫频响，再用三参数或五参数Fano公式拟合出线宽，步骤多、拟合质量对初始值敏感；特征频率法一次求解直接得到复数频率，信息量更大，特别是可以同时拿到本征模场分布。表格对比如下：

实际项目里两种方法经常配合：先用特征频率法锁定目标模式和Q因子，再用宽频扫描验证透射谱里对应的共振线宽。

4.2 计算环节的细节

PML设置好之后，固定k点，运行特征频率研究。求解器设置有几个关键参数要手动调整：

• 特征频率数：不要贪多，默认会求6~10个。BIC通常落在目标频段，特征频率数量设置成10~14比较保险，太多反而浪费求解时间。
• 搜索中心频率：填目标共振频率附近的预估值，CTRL键逐个尝试，直到目标模式出现。这样解析出来的模式是最接近物理真实的弱辐射本征模。
• 相对容差：默认1e-6在小虚部BIC模式上不够用，改成1e-8或更严格，尤其当Im(ω)比Re(ω)小四五个量级时。
• 特征频率的研究扩展项：启用“计算参与因子”，协助甄别假模。

求解完成后，在派生值里提取实部频率f_real和虚部频率f_imag，Q = f_real / (2 * abs(f_imag))。这个计算可以直接做一个全局表达式，这样在参数扫描结果里每个参数组合都能自动算出Q。我习惯在结果里同时输出f_real、f_imag和Q三列，方便后续判断模式是否进入BIC状态（虚部趋于零、Q趋于大值）。

4.3 BIC的Q因子：要收敛不要迷信

BIC理论上Q是无穷大的，但COMSOL里你永远不会得到无穷大。原因很简单：网格把连续空间离散化了，辐射通道无法精确关闭；PML的吸收也不是理想吸收；求解器迭代收敛的数值噪声也会贡献到虚部。于是你看到的BIC模式的Q往往是几万、几十万甚至更高，但不会严格是无穷大。

碰到这种“高但不无穷”的结果，要做收敛性分析，不能说看到一个大Q就宣布找到了BIC。方法是固定k点，将网格最大尺寸从a/12细化到a/24，看Q值变化。如果Q随加密几乎不变，说明数值误差已经足够小，可以认为该模式在物理上是BIC。通常我们以Q大于1e6且随网格加密升高不超过10%作为一个工程判定标准。

还必须把辐射Q和材料吸收Q分开。如果材料折射率有虚部（吸收），那么总Q是辐射Q和材料Q的并联关系：1/Q_total = 1/Q_rad + 1/Q_abs。在实际仿真里加入材料损耗后，虚部几乎全被材料吸收主导，你看到的总Q会远小于纯辐射Q。这也是为什么我前面强调，正式研究中要把材料损耗单独处理后合并分析。

5. 调控Merging BIC与远场偏振验证

5.1 结构调整：让两个BIC走向合并

Merging BIC的物理实现，本质是通过改变某个几何参数，让动量空间中的多个基态BIC位置发生移动，最终在临界参数处汇聚。比较常见的方案是改变介质柱半径、高度、椭圆度或横向偏移量，把原本不在Γ点的偶然BIC推向Γ点，跟对称性保护BIC合并。

我以改变圆柱半径r为例：固定a和h，把r从200 nm到280 nm逐次扫描。每个r值都做二维能带扫描，记录非Γ点BIC在(kx, ky)空间里的位置。典型的行为是：r较小时，两个在kx轴上（或对角线方向）对称分布的BIC离Γ点较远；随着r增大，两个BIC逐渐向Γ靠拢；在某个临界半径rc处，两个BIC移到Γ点并与原本的对称性保护BIC合并，这时候Γ点附近的能带色散关系会发生明显变化：频率面变得平坦，同时Q因子出现一个极大尖峰。

这里最花时间的是定位两组BIC的准确位置。我通常用粗扫描看清楚大致坐标，再用局部加密扫描精确定位。如果你观察到BIC位置随参数单调变化，说明追踪方向是对的。如果有突然的跳变，先回一维能带图检查是否出现了模式身份交换。这种身份交换在结构参数扫描中极其常见，是Merging BIC模拟做得怀疑人生的主要来源之一。

5.2 远场偏振涡旋：BIC的“身份证”

怎么确认你看到的BIC真的是BIC，而不是一个带边上恰好漏到外面的模式？最可靠的方法之一是看远场偏振分布。BIC在动量空间的远场偏振分布会形成一个涡旋（vortex），涡旋中心对应BIC位置，它周围的偏振方向绕中心旋转一圈，累计旋转角度等于2π的整数倍，这个整数就是拓扑荷。

实际操作中，我通常用两种方式在COMSOL里展开远场偏振分析。第一种是用COMSOL自带的远场计算：在模型外缘设置远场域，然后计算本征模在该域的场值，输出随角度变化的远场方向图和偏振分量。第二种更灵活的是把本征模在单元胞上表面的电场分布导出，在外部脚本里做二维傅里叶变换，从而得到平面波展开系数。每个平面波分量对应一个远场辐射方向，其矢量的两个分量Ex、Ey就构成该方向的线偏振（或椭圆偏振）态。

在(kx, ky)平面里扫一片区域，把每个k点的远场偏振方向用短线箭头画出来，你就得到一张“偏振地图”。普通辐射态的偏振箭头方向是均匀的或缓慢变化的，而BIC附近的偏振箭头会绕着一个中心点旋转，形成一个清晰的偏振涡旋。这个涡旋中心就是BIC所在的动量位置。顺便提一句，拓扑荷一般取±1，合并后的BIC可以呈现更高阶的拓扑荷，所以远场偏振图还能定量告诉你合并的完成度。

5.3 合并判定与参数临界点

判定Merging BIC是否成功的标准，通常看三点：

• 第一，二维能带面上，Γ点附近的频率面是否在临界参数处出现上下分离或凹槽，形成孤立点；这不是开口率很大的能带，而是模式频率与连续谱交叉点处的行为突变。
• 第二，Q因子在临界参数处出现尖峰，偏离临界值后Q迅速下降到百级或千级，这个峰越尖说明合并越干净。
• 第三，远场偏振图上，原本分开的涡旋在临界参数处合并在同一个中心，拓扑荷变成叠加值。

三个判据里至少满足前两条，才能下“Merging BIC实现”的结论。我实际做下来，最稳定的是Q因子判据：因为能带交叉点附近的模式本来就难追踪，偏振涡旋的数值噪声也偏大，而Q因子的局部极大值非常明显，不容易误判。当然这只是作为工程判定，物理上严谨的结论还是要以磁场分布、偏振拓扑为主。

为了准确锁定临界半径rc，我在rc附近用极细的步长做局部扫描，比如r从240 nm到255 nm，步长0.5 nm。这个量级的扫描会让总计算时间成倍上升，但整个过程是值得的，因为Merging BIC的临界行为常常只出现在非常窄的参数窗内。扫得太粗，很容易直接把临界点跳过去。

6. 踩坑实录：常见问题速查

6.1 特征频率假模泛滥

做周期结构计算时，特征频率求解器经常会返回一堆能量集中在PML区域或结构内部奇异点上的假模。判据很简单：查看参与因子，物理模式的参与因子通常接近1，假模偏低；同时看电场能量分布，假模往往大量集中在网格奇异处或边界附近。

处理方法：加参与因子筛选；把搜索频段限制在目标频率附近；检查边界条件是否设置成周期性。还有一招很有效：建模时把PML区域材料损耗单独设一个微小的值，假模在PML里的损耗特别大，虚部会明显偏高，一眼能区分开。

6.2 Q因子虚部为负或者忽正忽负

特征频率模式虚部为负表示损耗，为正表示增益，但数值噪声会让本应为负的虚部在某些模式下翻成正号，尤其在BIC模式虚部极小、只有1e-7量级甚至更低的时候。这种时候Q值会从正无穷跳到负无穷，数据完全没法用。

解决办法是回到收敛性分析，加密网格并提高求解器容差。如果加密之后虚部的正负号依然不稳定，说明这个位置其实没有严格意义上的BIC，只是一个高Q准BIC模式。这种情况下就不要强行提取Q了，改报“Q > 某个阈值”更诚实。我在一篇数据分析里，Γ点附近的虚部在不同网格尺寸下从1e-8变化到1e-6，我最终报告的是“Q因子下限为1e6”，而不是写一个具体数字。

6.3 远场偏振图噪声大

远场偏振分析遇到的最大问题是取样窗口不够宽或k空间步长过粗，导致傅里叶变换的旁瓣渗入，偏振箭头方向显得杂乱。这时把导出场分布的范围扩大两个单元胞，并对窗函数做窗化处理，旁瓣会显著减弱。另外偏振图的分辨率取决于你扫描的k点密度：要分辨涡旋中心区域，至少要在中心周围扫描15×15个k点以上，否则拓扑荷判读不准。

还有一个容易忽略的点：远场偏振图的涡旋中心可能在两个k网格点之间，而不是恰好在求解点上。这时需要用插值或局部加密扫描找到准确的涡旋中心，否则你提取的拓扑荷可能不完整。

回头看这套流程，从物理图像、建模、能带扫描、Q因子到远场偏振，每一环都是相互嵌套的。单纯把COMSOL模型跑通不算什么，难的是在计算噪声里判断哪些结果是物理的，哪些只是数值伪影。我自己的体会是，做Merging BIC这类题目，至少60%的时间要花在验证计算的收敛性和模态追踪上，真正建模求解的时间反而只占一小部分。另外，参数化设计不只是为了调出一个好看的结果，更是为了理解整个系统如何响应结构变化——这一点对后续做实验和器件优化比一个临界参数值有用得多。希望这份实操记录能帮你绕开我踩过的坑。