BIC单向辐射损耗计算 COMSOL光子晶体超表面计算
光子晶体超表面的BIC特性最近在光学圈子里特别火，但真要自己上手算单向辐射损耗的时候，总感觉参数调得脑壳疼。今天咱们直接拿COMSOL开刀，手把手教你如何在仿真里揪出那个神秘的辐射损耗值。
先搞个最简单的二维光子晶体结构练手——正方形晶格介质柱阵列。建模时记住这组黄金参数：晶格常数a=500nm，介质柱高度h=220nm（氮化硅材料），基底用二氧化硅。重点来了，想要触发BIC，得先让结构对称性玩点花样。比如把介质柱从正圆改成椭圆，x轴半径给120nm，y轴半径调成100nm，这种对称破缺操作直接关系到辐射损耗的计算精度。
边界条件设置最容易翻车。上下边界用完美匹配层（PMPL），左右用Floquet周期边界。特别注意端口设置要开启衍射阶数分析，COMSOL里的"wave vector port"选项记得勾选"compute diffractive orders"，这里藏着辐射损耗的关键数据。跑完频域计算后，在结果里翻出散射参数的S11曲线，Q值计算公式直接拍脸上：
lambda0 = 1550e-9; % 中心波长 FWHM = abs(lambda_peak(2)-lambda_peak(1)); % 半高宽 Q_factor = lambda0/FWHM;但这时候你可能会发现Q值高得离谱，八成是网格划分太糙。试试在介质柱边缘加自定义网格，边界层数给到3层，最大单元尺寸压到λ/8。这步做完至少能让计算结果从玄学走向科学。
想验证是不是真BIC？上参数扫描！把介质柱椭圆度从0到0.3扫一遍，观察Q值随不对称度变化的曲线。这里祭出COMSOL的批处理大法：
for (double eta=0; eta<=0.3; eta+=0.05){ model.param().set("eta", eta); model.study("std1").run(); double Q = model.result().numerical().getReal(); exportData(Q); }当曲线出现明显拐点时，恭喜你找到BIC的藏身之处。这时候的远场辐射图应该像个害羞的刺猬——能量死死锁在平面内，基本没有向上泄漏的趋势。
最后给个实用小技巧：计算损耗时别死磕频域，试试时域仿真看能量衰减曲线。虽然计算量暴增，但能直观看到电磁场能量如何随时间泄漏。用这个数据拟合指数衰减曲线，辐射损耗系数直接现原形：
from scipy.optimize import curve_fit def decay_func(t, gamma): return np.exp(-gamma*t) popt, pcov = curve_fit(decay_func, time_points, energy_data) print(f"辐射损耗系数: {popt[0]:.3e} 1/s")这个方法特别适合验证那些处在BIC临界状态的结构，毕竟眼见为实嘛。调参调到怀疑人生的时候，不妨泡杯咖啡让时域仿真慢慢跑，说不定就等来个柳暗花明呢？
