COMSOL仿真指南：准BIC增强复合波导光栅的古斯汉森位移

上个月帮课题组跑了一个复合波导光栅的电磁仿真，目标就一句话：在1550 nm通信波段，用准BIC把古斯汉森位移做大。模型本身不算复杂，但真正把Q因子、反射相位和横向位移这几个量串起来，中间有不少容易翻车的地方。这篇把完整的仿真思路和关键设置在COMSOL里的落地方式记下来，给后面做导模共振、BIC或者位移传感器的同学一个可以直接参考的流程。内容会涉及物理背景、几何参数、特征频率分析、频域后处理、参数扫描和几个我实际踩过的坑，不需要你提前懂BIC理论，跟着一步步做就行。

1. 为什么是“准BIC+复合波导+光栅”这个组合

1.1 从BIC到准BIC：先搞懂这个“幽灵模式”

连续谱中的束缚态，也就是BIC，这个词听起来很玄，但物理图像其实很直接。一个波导模式本来可以通过光栅的倒格矢耦合到自由空间辐射连续谱，产生损耗；但如果这个模式与辐射波的对称性恰好“不匹配”，比如电场分布关于光栅中心对称轴是奇函数，而正入射平面波在这个轴上是偶函数，那么两者耦合积分为零，模式就不会向外辐射。它把自己锁在结构里，既不吸收也不散射，形式上是一个没有线宽、Q值为无穷大的暗态。这就是理想BIC。实际器件里，只要稍微打破这个对称性，比如把光栅占空比从0.5改成0.6，模式就能和辐射通道耦合，获得一个虽然有限但通常很大的Q值，成为准BIC。准BIC的本质是“泄露得很少的高Q导模共振”，但它比普通导模共振更可控，因为你可以通过对称性破缺量来连续调节辐射损耗。这个特性在传感、滤波和非线性增强里都很有用，我的仿真里主要用它来放大相位响应。

1.2 GH位移放大机制：相位梯度是关键

古斯汉森位移是全反射时反射光束相对于几何反射点沿界面发生的横向偏移。它不是什么新现象，一百年前就有人研究了。稳态相位理论给出的表达式是：

Δ = - (λ / 2π) · (dφ / dθ)

其中φ是反射系数相位，θ是入射角。这说明位移的大小直接取决于反射相位对入射角的梯度。普通介质界面的全反射，相位变化比较平缓，位移通常只有波长量级；但在共振附近，反射相位会在极窄的角度范围内发生急剧跳变，dφ/dθ 可以比本底大几个数量级，位移也就被放大了。那为什么非得用准BIC？因为准BIC的共振线宽可以做到比普通导模共振窄得多，Q值高意味着相位变化更陡峭。理论上一旦进入BIC，Q发散、线宽趋于零，相位梯度趋于无穷，但实际不可能无限放大：准BIC总要有一定的辐射尾巴才能和入射光耦合，否则它就是个纯暗态，你根本看不到反射峰。所以实际设计要找一个平衡点，让Q值足够高但又没有被网格和角度步长“抹平”。

1.3 复合波导比单层波导多了什么

单层波导光栅也能做准BIC，但复合波导有额外的设计自由度。所谓复合波导，一般指两个高折射率波导层中间夹一层低折射率间隔层，或者光栅层本身与另一平面波导层相距很近。两层波导各自支持模式，当两个模式的传播常数接近时会发生杂化，形成对称和反对称的super mode，色散曲线上出现反交叉。这个反交叉点的位置可以通过层厚、间隔层折射率来调节，也就等于给准BIC的共振波长和角度增加了一个“调谐旋钮”。我的经验是：如果你想在某个固定通信波长下工作，复合波导能让你把模式钉在目标波长附近。后期做参数容差分析时，这一点特别重要，因为工艺上光栅周期的误差很难避免，但层厚的调节相对容易。

2. 几何建模与参数化：从物理结构到COMSOL模型

2.1 一个可复现的结构参数起点

仿真不可能从零开始拍脑袋。你的第一个任务是用平板波导本征方程估算有效折射率，再用光栅耦合条件确定周期量级。我这次用的是“空气 - 部分刻蚀Si3N4光栅层 - Si3N4波导层 - SiO2间隔层 - Si3N4第二波导层 - SiO2衬底”这样一个复合波导结构。因为光栅和波导同材料，通过一次刻蚀就能同时形成，实际加工相对友好。参考参数如下：

这些数不是随便写的。周期1000 nm配合Si3N4波导的有效折射率约1.7，可以让导模共振发生在近正入射附近，方便用端口和PBC处理。如果你用的材料是TiO2或Si，折射率不同，周期要相应回落或增加。建模前用解析平板波导算一下TE0模的neff，再代进相位匹配条件 k0·sinθ + m·2π/Λ = k0·neff，基本能把周期确定在正负50 nm以内。

2.2 COMSOL中的几何搭建要点

我用的是COMSOL 6.4的电磁波频域接口，二维模型，x方向取一个周期长度，y方向从下到上依次是：PML、SiO2衬底、第二Si3N4层、SiO2间隔层、第一Si3N4层、Si3N4光栅层、空气层、PML。左右边界设置Floquet周期条件，特别注意周期条件里的波矢是 kx = k0·sinθ，而不仅仅是默认的0。空气层高度至少留2个波长，不然倏逝波会被PML边界影响；PML厚度取一个波长左右就行，我一般设1500 nm，比例系数用1。

光栅部分用两个矩形画“齿”和“槽”，然后通过差集得到光栅截面。占空比直接作为参数驱动，方便后面扫描。这里有个几何细节：光栅齿两侧必须严格平行于y轴，齿顶要平整，不要有微小斜面。因为准BIC对对称性极其敏感，COMSOL几何容差里如果出现零点几纳米的误差，虽然不会根本改变物理，但会干扰你对Q值的判断。

2.3 物理场与边界条件怎么配

• 偏振：二维模型里，TE偏振对应电场沿面外方向，也就是z分量，选择“面内磁场”还是“面内电场”要看清接口。我要算的是TE模式的GH位移，所以设置电场Ez分量。
• 入射：用周期性端口作为激励，端口类型选“周期性”，模式指定为TE平面波。出射端口放在衬底下方，用来同时读取透射。端口会自动给出S参数，S11就是反射系数。
• 上下边界：必须用PML，不能用完美磁导体或完美电导体代替，否则反射波会在边界来回振荡，相位根本不准。
• 网格：光栅齿和波导层用映射网格控制，最大单元尺寸设到 λ/12 左右；空气区域用自由三角形；PML区域独立划分。第一次粗扫可以放宽到 λ/8，锁定共振峰后再加密。

建模阶段最容易犯的错误是把端口放得离结构太近。端口在COMSOL里会默认把参考面设在端口边界位置，如果你后面要用S11的相位算位移，必须把端口参考面到光栅表面的传播相位减掉。我后来干脆把端口面直接放在光栅表面上方约50 nm的“近场”位置，配合足够高的空气层，省去了一堆手动扣除相位的麻烦。

3. 特征频率研究：从模场分布里把准BIC认出来

3.1 为什么不直接扫反射谱，先做模式分析

你当然可以用频域扫描直接找反射峰，但只看反射峰你无法判断它到底是导模共振还是准BIC，也无法获得Q值信息。我的做法是先用特征频率研究，把目标模式的复本征频率算出来：实部对应共振频率，虚部对应辐射损耗。Q = Re(f) / (2·|Im(f)|)。理想BIC的虚部理论上为零，数值上会给出一个很大的Q；准BIC的虚部不为零但很小。这一步的价值是让你确定准确的共振波长，之后再回到频域扫描，你会知道峰应该出现在哪，网格和角度扫描范围该怎么设置。

3.2 特征频率研究里的边界条件设置

在电磁波频域接口下添加“特征频率”研究，注意特征频率求解时不加端口激励，所以要把端口条件去掉，左/右Floquet周期条件保留。周期波矢同样设为 kx = k0·sinθ，而不是扫描频率后自动变化。上下边界用散射边界条件或者PML都可以；如果想快速摸清模式分布，我建议先用散射边界跑一遍，等确认模式后再换PML验证Q值。特征频率求解范围设置在 150 THz 到 200 THz 附近，对应1550 nm附近的频率是约193 THz，给一点带宽余量。求解器可能会找到一堆背景模式，不用慌，看模场图逐个排除。

3.3 怎么判断“这就是我要找的准BIC”

判据有三个。第一，模场能量绝大部分局域在波导层内，空气和衬底里的泄露尾巴很小。第二，电场分布关于光栅中心线呈反对称或具有“四极子”特征，和正入射平面波的对称性不匹配。第三，在对称占空比0.5附近，本征频率虚部趋于零；当你把占空比改成0.55，虚部会增大但仍在可接受范围。我用的是一个复合波导模式，模场在两个Si3N4层里都有分布，中间SiO2间隔层里的场较弱但并非零。这种杂化模的对称性比单层波导更复杂，反而更容易找到满足BIC条件的模式。实际操作时，我会在“结果”里画Ez的实部分布，再用“全局计算”输出 Q = 0.5 * real(f) / abs(imag(f))，把多个本征模式的Q列成表，一目了然。

这里有一个容易踩的坑：特征频率求解器默认用实数搜索还是复数搜索，不同版本有差异。一定要确认输出的本征频率带有虚部，虚部会以“+ i * 1.2e-4 [THz]”的形式写在结果里。如果你的结果全是纯实数，那说明求解设置把虚部丢了，需要调整特征频率研究的“搜索频率周围”或者手动指定复数搜索范围。Q值计算错误会直接影响后续参数扫描的判断。

4. 频域反射谱与相位梯度：古斯汉森位移的后处理计算

4.1 反射相位怎么取才准

确定准BIC频率后，回到频域研究，扫描入射角。比如在共振角附近从 -1° 扫到 +1°，步长先取0.01°，锁定峰后再用0.002°加密。计算完成后，在“派生值”里选择端口S参数，可以得到S11的实部和虚部。反射系数振幅 |r| = sqrt(Re² + Im²)，相位 φ = atan2(Im, Re)。这里最大的坑是相位包裹：atan2给出的相位在 -π 到 π 之间，共振处相位变化可能跨过 ±π，如果你直接对包裹后的相位求导，会在跳变点产生巨大的虚假尖峰。解决办法是用后处理里的“解包裹”功能，或者自己写一段循环把相位相邻差值超过π的地方加上或减去2π。我在COMSOL里用“全局计算”导出数据，再在外部用Python做unwrap和数值微分，这样更灵活。

4.2 GH位移公式和数值微分实现

使用稳态相位公式：

Δ = - (λ0 / (2π)) · (dφ / dθ)

注意角度必须是弧度制。数值微分采用中心差分：

dφ/dθ ≈ [φ(θ_{i+1}) - φ(θ_{i-1})] / (θ_{i+1} - θ_{i-1})

如果角度步长是0.002°，换算成弧度约3.5e-5 rad，相位差可能是几十度，差分结果会比较稳定。我实测在共振角附近 dφ/dθ 可以达到数千rad/rad，代入公式后位移到几十微米到几百微米量级，比普通界面的波长级位移大了三四个数量级。可以画两条曲线：一条是反射率随角度的变化，一条是GH位移随角度的变化。你会发现位移峰值并不正好在反射率峰值处，而是偏向相位梯度最大的位置，也就是反射率峰的侧翼。如果发现位移峰和反射峰位置完全重合，一般是相位处理出了问题，回头检查unwrap。

4.3 高斯光束直接模拟：交叉验证

相位梯度法虽然快，但它是基于平面波稳态相位理论的近似。准BIC共振线宽极窄，入射的高斯光束如果角谱太宽，共振响应会产生复杂畸变，实际位移可能偏离公式预测值。严谨一点的做法是额外建一个模型：入射端口不再用理想平面波，而是用背景场设置为高斯光束，在距离结构表面一定高度处观察反射场分布，取反射光束的质心位置减去几何反射位置。这个验证不需要做全参数扫描，只选两三个关键参数点看一下就行。COMSOL里设置高斯光束背景场不算麻烦，但网格需要更密，因为要分辨光束空间分布。两种方法如果偏差在20%以内，说明相位梯度法可用；如果偏差大，需要缩小入射角谱宽度，或者检查模型边界是否引入了额外反射。

4.4 用MATLAB或Python批量导出

手动点界面只能算单点。我在做参数扫描时习惯用COMSOL LiveLink for MATLAB写循环：修改占空比、重新求解、提取S11、计算相位和位移，把结果存成CSV。实际上如果你熟悉COMSOL 6.4的Java或Python API，在外部脚本里批量控制模型会更方便，整个参数空间几十个点可以无人值守跑完。这个流程不复杂，核心就是这个循环：更新参数、求解、取复S参数、后处理。别在COMSOL界面里手动点几十次，容易点在错误的数据集上，后处理代码写脚本也能保证每次计算口径一致。

5. 参数扫描与结果解读：让准BIC和GH位移真正“对上”

5.1 扫描占空比：从完美BIC到准BIC

占空比0.5对应结构关于z=0面镜像对称，这个点通常存在对称性保护的BIC。在特征频率研究里，你会看到Q值非常大，但频域反射率几乎看不到峰，因为正入射平面波无法激发它。把占空比从0.5调到0.52、0.55、0.6，Q值会逐步下降，反射峰逐步明显，相位梯度先增大后减小。为什么不是Q越大位移越大？因为位移来源于实质的反射过程，你需要足够的耦合把光“送进”准BIC再“放出来”；Q无穷大时光根本不进去，位移反而趋于零。所以存在一个最优Q区间。我这次扫描占空比0.5到0.7，Q从10^6附近降到10^3附近，GH位移在占空比0.57附近达到峰值，约80 μm。你可以用这个规律反推实验容差：如果工艺误差让占空比偏了0.03，位移掉多少，一看曲线就知道。

5.2 间隔层厚度与模式反交叉

复合波导比单层波导多出来的那个“旋钮”，就是两层Si3N4之间的SiO2间隔层。扫描间隔层厚度从200 nm到500 nm，你会发现两个模式的共振频率出现反交叉：原本随厚度单调变化的两个峰，在某个厚度附近互相排斥、交换模式特征。反交叉点附近，模式杂化最强，准BIC的电场分布会从“偏上波导”变成“上下都有”，这对GH位移的影响非常明显。实际操作时，可以固定占空比，把间隔层厚度作为第二个扫描参数，做二维扫描。结果可以整理成一张二维色图：横轴是间隔层厚度，纵轴是波长或入射角，颜色是GH位移。你会看到一条明亮的“位移脊线”，沿着它选取工作点，就是兼顾制造容差与增强效果的位置。

5.3 反射谱中的Fano线型

准BIC共振在反射谱上通常不是对称的洛伦兹峰，而是Fano线型，因为连续背景通道和离散准BIC通道之间的干涉。做参数扫描时不要只看峰值，要关注线型的不对称性。GH位移的峰值通常出现在Fano线型的快速上升沿或下降沿，具体在哪一侧取决于准BIC的耦合相位。我这次算的是反射率在共振角附近先跌到接近零再冲高，位移最大值出现在反射率上升沿偏右侧。在写论文或工程报告时，最好把反射率、反射相位、GH位移三条曲线放在同一张图里，能直观看出它们的关联，也方便你判断数据的自洽性。

6. 实际跑下来最容易翻车的几个地方

6.1 相位参考面和端口的距离

第一个坑我在前面提过：端口参考面不在结构表面时，S11会携带一段传播相位，且随角度变化。这段相位是随角度平滑变化的“背景”，会混入你要求的相位梯度。如果空气层高度有10个波长，背景相位梯度可能足以淹没共振信号。要么像我把端口贴近结构表面，要么在提取S11后人为减去 k0·cosθ·h 对应的相位路径。

6.2 共振峰太窄，网格不够密会“抹掉”准BIC

准BIC的反射峰线宽可能只有0.1 nm甚至更窄，如果网格尺寸按空气波长λ/10划分，共振会被人为展宽和降低，计算出的GH位移偏低。我的经验是：先用粗网格找到峰位，然后只对光栅和波导层加密，重点加密两个材料界面附近；再对比λ/12、λ/16、λ/20三种网格密度下的位移，如果峰值变化小于5%，就认为收敛了。网格加密后内存占用上升很快，二维模型还好，如果做三维模型，建议利用对称性减少计算域，或者只做单个入射面。

6.3 角度扫描步长太粗，相位跳变被漏掉

如果角度步长大于共振角宽度的一半，你可能根本采不到相位跳变的中间点，差分计算结果会严重偏小甚至出现锯齿。先扫0.01°确定共振角，再在共振角附近用0.002°加密，这是最稳妥的做法。更高效的方法是用COMSOL的自适应网格和参数化扫描配合，但手动观察曲线趋势仍然不可少。

6.4 不要只看正入射，倾斜入射可能更优

很多准BIC研究喜欢强调正入射，因为实验光路简单。但古斯汉森位移对角度敏感，倾斜入射时模式耦合条件变化，相位梯度可能更大。我的建议是扫描入射角从0°到2°，同时观察模式演化。倾斜入射还会引入额外的模式对称性破缺，使Q值对角度变得更敏感，这既是麻烦也是机会：如果你能精确控制角度，等于多了一个调位移的旋钮。

6.5 后处理脚本要统一数据集

手动操作时很容易选错数据集，比如用了“特征频率”研究的解去算S参数，自然什么也得不到。写脚本时强制指定解的标签和数据集名称，把S11提取、相位unwrap、位移计算封装成同一个函数，每次参数扫描都调用，避免人为错误。这一点虽然听起来很基础，但我见过不止一个模型因为数据集选错而“算出”几十米量级的离谱位移。

做这类仿真，我个人的体会是：不要在建模初期追求复杂结构，先用对称占空比0.5把特征频率和模场跑通，再逐步打开非对称参数；每改一个几何参数都要重新确认Q值、反射谱和位移三者的一致性。等你真正理解了准BIC的“暗-亮”转换过程，后面换材料、换波段、换器件结构都只是参数问题。最后一个建议：所有关键结果都做一次网格收敛性验证，并在报告里写上使用的网格尺寸和角度步长，这对论文审稿和后续复现都特别有帮助。