一维光栅拓扑BIC的COMSOL模拟与单向辐射设计

我最近在搭一个光子晶体超表面的单向辐射模型时，撞上了那个经典现象：扫参数的过程中，某个模式的特征频率虚部突然跌到接近零，Q值在图上像坐了火箭一样往上冲。这个现象就是连续谱束缚态（Bound States in the Continuum，简称BIC），而把它放到一维光栅里研究，配合COMSOL模拟，恰好是理解拓扑BICs和设计单向辐射器件最顺手的路径。这篇文章想用做仿真的人之间最直接的交流方式，把物理图像、COMSOL建模细节、Q值提取、单向辐射量化这几个环节全部拆开讲清楚。内容适合正在做超表面、周期光栅、导模共振器件的朋友，也适合刚接触COMSOL波动光学的同学当进阶案例参考。

1. 先把物理图像捋顺：一维光栅里的拓扑BIC和单向辐射

1.1 一维光栅为什么是研究BIC最顺手的模型

先回答一个最基本的问题：为什么研究拓扑BICs要选一维光栅，而不是直接上三维超表面？

一维光栅在y方向平移不变，这意味着整个电磁问题可以用2D模型完整描述。2D模型在参数扫描里的优势太明显了。BIC在本质上不是一个单点现象，它是在参数空间中出现的拓扑奇点，你需要同时扫描占空比、光栅高度、布洛赫波矢kx等多个变量才能把它的行为看清。3D超表面模型即使能跑，一轮扫描也要几天甚至一周，而2D模型几小时就能出一批数据，迭代速度完全不在一个量级。

物理上，周期性光栅支持布洛赫模式。模式的色散关系落在真空光锥之外时，会被全内反射限域在波导层内；落在光锥之内时，原则上可以耦合到外部平面波并向外辐射，这就是导模共振（Guided-Mode Resonance）的基本条件。BIC最反直觉的地方在于：它明确落在光锥里，有开放的辐射通道，却依然不向外发射能量。一维光栅中多个衍射级次和多条导模色散相互交织，正好是观察“辐射通道被干涉锁死”的最简单结构。

对想用COMSOL练手的初学者，一维光栅还有一个额外的好处：偏振分离非常干净。TE偏振对应电场只有Ey分量，TM偏振对应磁场只有Hy分量，在2D模型里两者完全解耦，可以分别建立模型。你不会遇到三维超表面里那种模式混叠导致“算出来一个模式却不知道它是什么”的尴尬。

1.2 对称保护BIC、Friedrich-Wintgen BIC和拓扑BIC的区别

这个话题容易越讲越玄，我尽量用做仿真的语言把它落到可操作层面。文献里常出现的BIC其实可以分成三类，实际复现时你会发现它们经常耦合在一起。

对称保护BIC最直观：当光栅结构上下/左右对称、布洛赫波矢在Γ点（kx=0）时，某个本征模式的场分布与所有开放辐射级次的对称性失配，耦合积分严格为零，模式因此被“锁死”在结构内部。在COMSOL里算一次特征频率你就能看到：对称结构中，某个模式的虚部在kx=0处突然趋近于零，同时Q值出现一个尖锐的峰。这种BIC不需要精细调参，但只能出现在高对称点，实际器件用起来不够灵活。

Friedrich-Wintgen型BIC靠的是两条模式之间的干涉相消。一维光栅中相邻导模的色散在参数变化时会发生反交叉，反交叉点附近两条模式同时耦合同一个开放辐射通道。如果参数合适，两条模式的辐射路径刚好相位相反，总辐射被完全抵消。因为这类BIC不要求kx=0，它可以出现在布里渊区内部，需要你自己在参数空间里把它找出来。从仿真角度讲，这就是在log10(Q)二维图里找亮斑或亮线的过程。

拓扑BIC的说法来自远场偏振矢量的拓扑结构。把布洛赫波矢和某个结构参数（比如占空比）当成参数空间，每个模式对应的远场偏振可以表示成一个矢量场。BIC所在的位置远场强度为零，偏振方向失去定义，相当于矢量场里的一个涡旋奇点，带有整数拓扑荷±1。这种奇点有一个重要性质：一般参数扰动只能移动它的位置，不能凭“微调”把它消除，除非两个拓扑荷相反的奇点碰到一起发生湮灭。在Q值的具体表现上，拓扑BIC附近满足 Q ∝ 1/Δk²，即偏离BIC的波矢越大，Q按平方反比下降。这条规律可以直接用来验证你在COMSOL里看到的Q发散到底带不带拓扑保护性质。

1.3 单向辐射的本质：把完全锁死变成受控泄漏

BIC最直接的用途是提供超高Q值，但单向辐射需要的是一个有方向性的发射通道。这个逻辑我一开始也绕了很久：BIC完全不辐射，怎么可能做单向辐射？

关键是把BIC变成准BIC。通过轻微破坏结构的对称性，辐射通道被打开，但因为耦合效率仍然很低，模式可以保持很高的Q值，同时泄漏能量被引导到指定的辐射方向。对称结构的上/下辐射一定是相等的，想要真正的方向性必须同时打破上下镜像对称性，让向上的泄漏路径和向下的泄漏路径不等价。

工程上实现单向辐射有两条主流路线：一是做几何非对称光栅，让上下界面形态不一样，比如矩形齿改成梯形或锯齿形；二是在光栅的一侧加反射结构，把泄漏到某一侧的波重新引回另一侧，实现能量单向汇聚。后者往往能把单向性推到极致，但结构复杂度上去之后，仿真和工艺的麻烦都会增加。

对做器件的人来说，这一步是整个设计的甜蜜点：Q决定光谱线宽，单向性决定能量利用率，而在准BIC框架里这两个指标是互相独立的。靠轻微破缺调整单向性，不会把Q压得不可接受，这正是拓扑BICs在单侧激光器、集成量子光源、定向传感器件里频繁出现的原因。

2. COMSOL模拟方案设计：从物理场到边界条件的取舍

2.1 特征频率法用来找BIC，频域法用来测响应，两条腿都要走

COMSOL里没有专门“算BIC”的按钮。通常的路线是先用“电磁波，频域（ewfd）”接口下的特征频率研究找出模式的复频率，再用频域研究验证谱线形状。这两种方法定位不同，不能互相替代。

特征频率法直接给出复频率的实部wr和虚部wi，Q = |wr| / (2|wi|)。在BIC附近虚部趋近零，Q值飙升，效率非常高。频域法则模拟真实实验条件：入射平面波、扫频率，在准BIC附近观察反射/透射光谱的Fano线型，用线宽算Q = λ0/Δλ。问题在于高Q对应的线宽只有皮米量级，比如Q=10^5、中心波长1μm时，线宽大约0.01nm，频域扫描步长要取得很小，计算量很恐怖。我的个人习惯是先特征频率法扫参数空间锁定BIC的位置，再用频域法只在目标波长附近做窄带确认，两者互相印证。

偏振处理也要注意。2D模型中，TE偏振意味着只有电场面外分量Ey有意义，TM偏振则对应磁场面外分量Hy。在ewfd接口中默认求解全部三个场分量，但建模前最好明确约束偏振，比如直接设置电场分量为Ey方向。否则特征频率求解器会混入大量无用的面内模式，后续筛选成本很高。我试过不约束偏振直接算，出来的前几十个特征频率里真正关心的目标模式被埋在大量噪声里，极其痛苦。

2.2 周期边界与PML的布置是成败关键

建模时x方向取一个周期单元，左右边界设置Floquet周期性边界条件，指定布洛赫波矢kx。上下边界怎么处理，取决于你要算的东西。

找BIC时，上下外侧要放PML吸收层，让模式向外辐射的能量被吸收，特征频率的虚部才会包含真实的辐射阻尼项。做反射/透射谱时，上下边界改用周期端口，COMSOL可以自动处理各个衍射级次的透射反射，直接输出S参数。

这里有一个我踩过好几次的坑：PML区域本身也必须参与Floquet周期性边界条件，而不是只在物理区加周期边界。如果不给PML设置周期条件，在kx≠0或者斜入射时，PML的角点会产生强反射，严重污染虚部，Q值的物理趋势会被噪声盖掉。解决方法是把左右两条边界上的PML区域一并选中，和物理域一起设置周期边界条件。

PML的厚度和放置距离也有讲究。PML太靠近结构会耦合到倏逝近场，一般物理域上下各留至少0.5λ到1λ的空气层，PML厚度取0.5λ到1λ。2D模型里PML坐标拉伸类型要选择与模型边界方向对应，常见的是z方向拉伸。方向选错的话，轻则吸收效率低，重则直接报错。每次新建模型我都会反复确认这个坐标对应关系。

2.3 材料与几何参数的初始取值：先让BIC出现再谈工程化

材料选择上，我强烈建议用低损耗甚至无损介质起步。氮化硅（折射率约2.0）或二氧化钛（约2.4）都是很好的选择，衬底用二氧化硅（约1.45）。不要一开始就用硅或带复折射率的金属，材料损耗会给Q设一个上限，拓扑奇点的“虚部趋零”特征会被吸收淹没，你根本看不到预想中的现象。

一组能快速出结果的初始参数：周期Λ=600nm，光栅高度h=180到220nm，占空比f从0.3到0.7扫描，光栅齿用矩形截面。在近红外约900nm到1200nm范围内，这组参数支持一阶导模，TE偏振下可以看到明显的模式反交叉，非常适合复现Friedrich-Wintgen型BIC。

折射率对比度也要仔细选。如果包层和衬底折射率非常接近，模式上下泄漏近似对称，单向性的设计空间会完全不同。我的建议是先做“光栅在顶层、空气包层、SiO2衬底”的标准配置，把物理机制研究清楚之后，再切换成全对称平板结构去验证拓扑荷特性。工程化的事情放后面，先把现象复现出来。

3. 实操流程：从几何建模到Q值提取的完整套路

3.1 几何与网格：高Q模式最挑剔的就是网格

几何建模本身没有太多技术含量：一个矩形光栅齿、一个衬底矩形、上下空气层、PML区域。但我强烈建议把所有尺寸用全局参数管理，比如Lambda、h、w、d_sub、d_air、d_pml。这样后续参数化扫描占空比时，只需要改w这一个参数，省去大量重复建模操作。

网格是BIC仿真里的重头戏，我在这上面浪费过大量时间。高Q模式的场主要集中在波导层内，所以光栅区域和近场缓冲层必须优先加密。我的经验性规则是：波导层内最大网格尺寸设为λ_min / (8 * n_eff)。比如波长1000nm、有效折射率2.0时，网格大约60nm。空气区域可以放宽到λ_min/6。初次扫描阶段可以先放宽到λ/4，锁定候选BIC位置后再把网格加密一倍做收敛性验证。

需要提醒的是，不要在BIC位置单独加密网格。正确做法是记录同一参数下网格尺寸减半前后的Q值变化：普通准BIC模式的Q会有一定波动，但拓扑BIC附近的Q会随着网格加密持续上升，说明Q发散是物理真实的。反过来，如果加密后Q下降并趋于某个平台，那很可能只是一个普通的高Q模式，或者Q已经被边界反射、PML截断限制住了。

3.2 特征频率求解器配置与复频率的符号陷阱

特征频率求解建议使用“在预定频率附近搜索”的选项，而不是全局搜索。以目标波长1μm为例，角频率约1.88e15 rad/s，可以把搜索区间设为1.7e15到2.1e15 rad/s。如果选择“所有特征频率”，PML伪模式和空气模式会被一起扫出来，筛选过程非常费神。

COMSOL默认的时间因子是exp(-iωt)，复数特征频率写作ω = wr - i*wi，一般wi>0表示衰减。但不同版本、不同接口之间符号约定可能有差异，算完Q值后建议用已知参数手动验证一下：把一个文献里的模式参数代进去，看Q数量级是否合理，然后再继续后续工作。否则你有可能得到虚部为正的“增长模式”，实际上只是符号约定颠倒，白白排查几天。

另外不要把特征频率的实部直接当成真实频率画色散。特征频率的实部在损耗较大时会偏离频域结果，一般需要Q>100之后两者才比较接近。初学阶段建议在一个固定的kx下同时算特征频率和频域窄扫，对比中心波长是否重合。

3.3 参数化扫描、logQ图与Q~1/Δk²拟合

拓扑BIC在参数空间里的识别方法非常直接：扫描kx和占空比f，计算每个参数点的Q值，然后画log10(Q)的二维图。BIC在图上会表现为一条或几点明亮的尖峰，Q突然冲到10^6以上，离开这个点马上掉回10^3到10^4。如果两个拓扑荷相反的奇点都在参数区里，甚至能看到脊线断裂、涡旋对合并的图像特征。

单个kx方向的发散规律要用拟合验证。在BIC所在的kx0附近取一组点kx = kx0 + Δk，计算Q值。如果是拓扑保护BIC，Q与Δk的关系应该满足Q ∝ 1/Δk²，在log-log坐标下斜率就是-2。这一步能把“算出来了一个高Q”和“这是一个拓扑BIC”区分开，是论文里很有说服力的证据。

手动一个个跑参数点显然不现实。参数化扫描研究功能可以直接批处理，但涉及多参数和后续拟合时，我建议把COMSOL的全局参数接口和MATLAB或Python脚本结合起来，用外部脚本控制批处理、导出复特征频率并做后处理。实测下来，一个kx取40步、f取30步的扫描，2D模型在普通工作站上大约几小时能跑完，足够判断BIC位置和拓扑性质。我自己后来甚至把整个扫描-拟合-画图流程写成了脚本，一键出logQ图，效率提升非常明显。

4. 把单向辐射作为设计目标：量化与优化

4.1 用坡印廷矢量积分量化上下辐射比

单向辐射需要一个明确的量化指标。我用的方法是在准BIC模式下用平面波从顶部入射，然后在距离结构大约一个波长的上下两条水平积分线上，分别求坡印廷矢量的法向分量积分，定义单向比为：

uni = P_top / (P_top + |P_bottom|)

其中P_top是向上辐射功率，P_bottom是向下辐射功率。2D模型对应单位宽度功率，单位W/m。对称结构的uni天然等于0.5，典型的非对称结构可以做到0.85到0.95。对于更严格的前后比指标，后续可以再结合远场积分，但前期优化用这个近场功率比值已经非常稳定，并且与远场趋势一致。

这里有个要注意的坑：功率积分线必须离结构足够远，同时不能太远导致PML吸收掉太多能量。我一般把积分线放在距离结构3λ/4的位置，此时气场基本进入远场区域，已经可以反映实际辐射功率。如果积分线放在离结构太近的位置，倏逝场会污染结果，算出来的“高单向性”到实验里完全对不上。

4.2 三类破缺对称性的方案与仿真取舍

想要拉开上下辐射比，结构上常见的方案有三类，各有各的适用场景。

第一类是几何非对称光栅，把矩形齿改成梯形或锯齿形，让上下界面的散射效率不一致。优点是只改几何，不引入额外材料和工艺，在COMSOL里做参数扫描非常方便。缺点是实际加工中侧壁倾角不好精确控制，仿真中的最优角度在工艺上可能实现不了。

第二类是底部加反射镜，在光栅下方加介质或金属镜层，把向下泄漏的能量反射回上方，然后与顶部直接泄漏干涉相长，把能量都赶到同一侧。这类方案需要给反射镜使用真实金属材料（比如金、银），计算量上去了，而且金属趋肤深度决定了光栅到镜面的缓冲层厚度必须仔细设计。

第三类是双层光栅或布拉格反射端，用第二层小周期光栅相当于反射器，把某一侧的泄漏禁掉。这个方案的优势是反射带宽可控，缺点是层间对准公差在仿真中体现不出来，设计时要留较大的工艺裕度。

为了便于对比，我把三种路线的关键差异列在下面：

我自己偏好的顺序是：先用几何非对称方案扫参数，确认方向性可调之后，再在需要极限单向性时换成光栅加金属镜的方案。在COMSOL里从前者改到后者只是新增几个矩形和边界条件，但要注意金属层的网格分辨率，以及金属与PML之间的距离，避免出现非物理反射。

4.3 远场方向图：从近场模式走向实验可测的量

单向性最终要在探测器上看到，所以远场方向图必须算。COMSOL中求远场需要先定义远场计算边界，通常选在PML内侧的截断边界上，然后使用ewfd接口的远场计算功能，输出角分辨辐射强度。

一维光栅的辐射主要发生在x-z平面，方向图可以画成极坐标角分布曲线。准BIC模式在不对称参数下远场往往只有一个主瓣，位置和强度随几何参数改变。如果主瓣指向斜上方，并且上/下方向图明显不对称，说明单向辐射设计已经生效。实验上更常用的上/下半球总能量比，可以直接用前面4.1节的功率积分替代，两者趋势一致，不必重复计算。

需要注意的是，2D模型远场对应的是线源辐射，不能直接套用三维远场的功率公式。如果需要与实验绝对数值对比，要么把2D结果按单位长度辐射功率处理，要么扩展到3D模型重新计算。在前期参数优化阶段，相对比较已经完全够用，不必纠结绝对辐射强度的换算。

5. 常见报错、数值陷阱与排查速查表

5.1 特征频率求解器混入的PML伪模式怎么筛

用PML算特征频率经常冒出一批虚部很大的伪模式，它们的特点是场全部聚集在PML层内部，完全是数值产物。我的筛选方法是：逐个查看每个模式在物理区的电场模分布，物理模式的场会集中在光栅波导层里，伪模式则趴在PML层或者空气层角落，一眼就能分辨。

更省力的验证方式是用参数对比：保持其他条件不变，把PML厚度从0.5λ改成1λ，观察目标模式的复频率是否显著变化。物理模式的辐射阻尼应该基本稳定，伪模式则对PML参数非常敏感，复频率会大幅漂移。这样一轮就能筛掉九成以上可疑模式。

5.2 网格导致的“伪BIC”判别

网格太粗时，数值耗散太强，高Q模式根本显现不出来；网格只要恰到好处，Q可能冲到特别大，但继续加密之后Q反而剧烈跳变。遇到这种情况不要急着下结论，先回粗网格看趋势，再做收敛性曲线。

真正的BIC在细化网格时有明确特征：Q持续上升，频率收敛，虚部单调趋近于零。如果Q到了某个平台不再变化，说明边界反射或PML截断成了新的泄漏通道；如果Q上下抖动，大概率是网格没有解析出角落奇异点附近的场分布。我个人的处理经验是在光栅矩形齿的四个角做局部细化，那些位置是电场奇点最集中的地方。之前我犯过一个错误：只依赖全局加密，网格从5万加到50万，Q还是在抖，后来对四个角单独细化之后才稳定下来。

5.3 单位、偏振与周期性条件：三个不起眼却致命的细节

2D模型默认z方向厚度是1m，所有功率、损耗值都按单位厚度归一化。计算比值没问题，但如果需要做端口功率和PML吸收功率的绝对校验，必须把这个厚度系数考虑进去，否则能量守恒怎么都对不上。

偏振问题再强调一次：TE模式下只有Ey有物理意义，如果算出来的模式里Ex、Ez分量占比不小，说明结构或边界条件把偏振混杂了，或者端口激励类型设错了。一维光栅中TE和TM的BIC条件完全不同，kx位置、色散形状都不一样，不要指望一个模型同时覆盖两种偏振。

Floquet周期边界条件里需要填写的波矢是kx = k0 * sin(θ)。这个θ是入射角还是衍射角，取决于你的端口设置。我有一次复现时把θ定义反了，整个色散图左右颠倒，Q发散点全部跑到负kx侧，排查了很久才反应过来。建议每次建模时把入射角和kx的换算公式写进模型的描述文本里，免得隔几天自己都搞混。

5.4 一套稳妥的复现checklist与个人建议

按下面这个顺序走，我个人踩坑最少：

1. 先跑对称结构，在Γ点找对称保护BIC，确认模型整体设置正确；
2. 再扫描kx找到模式反交叉，定位Friedrich-Wintgen型BIC；
3. 画log10(Q)在(f, kx)二维参数空间的分布图，锁定拓扑BIC位置；
4. 在BIC附近做Δk收敛拟合，验证斜率是否为-2；
5. 之后才引入非对称几何，扫描单向辐射比；
6. 最后加损耗和真实材料折射率，评估实际可达到的Q值，确认样品能够测到对应的窄线宽。

最后想分享一点个人体会：BIC仿真最大的风险不是算不出高Q，而是算出一个数值上极其漂亮、但对工艺偏差极其敏感的奇点。所以仿真阶段务必多做容差分析，看看Q值和单向性如何随侧壁角、厚度、占空比扰动而变化。一维光栅的优势就是参数少、扫描快，特别适合做这类鲁棒性分析，这也是我强烈建议新手从一维结构切入拓扑BICs研究的原因。别一上来就冲三维超表面，那种模型的调参周期太长，容易把物理直觉磨没，也很难体会到参数空间里那个奇点从出现到消失的全部过程。