comsol BIC本征态计算，支持comsol直接出图。 2019PRL。 此为通用算法

comsol BIC本征态计算，支持comsol直接出图。 2019PRL。 此为通用算法

COMSOL这货玩本征态计算是真能整活。前阵子搞光子晶体的时候发现它的BIC（连续谱束缚态）求解功能有点东西，直接调参数就能出图。2019年PRL那篇论文里提到的通用算法，其实在软件里已经内置了现成的套路，咱不用自己造轮子。

comsol BIC本征态计算，支持comsol直接出图。 2019PRL。 此为通用算法

先整点几何模型试试水。在模型向导选电磁波->频域，然后随手画个二维光子晶体结构。记得在全局参数里塞几个变量，比如晶格常数a=500nm，圆柱半径r=0.3*a这种，后面改参数不用重新建模贼方便。代码片段大概长这样：

% 全局参数定义 a = 500e-9; r = 0.3*a; k0 = 2*pi/(1550e-9); % 工作波长

边界条件要特别注意，完美匹配层(PML)必须安排上。但BIC的特性决定了它对外辐射得是零，所以得在周期性边界里玩花样。有个骚操作是设置Floquet周期边界时把波矢kx设为扫描变量，这样能自动抓取Gamma点的模式：

% Floquet边界设置 floquet = mphphysic('floquet', 'emw', 'Floquet1'); floquet.set('k', {'0' '0' 'kz'}); % 这里kz用来扫参数空间

求解器配置才是重头戏。本征频率计算选频域分解法，记得勾选"Shift"选项。这个偏移量要设在工作频率附近，比如1550nm对应约194THz，不然算到地老天荒也抓不到模式。代码配置里藏着个彩蛋——把特征值缩放模式改成"Auto"，软件会自动处理发散问题。

% 本征求解器参数 study = mphstudy('Eigenfrequency'); study.set('shift', 2*pi*194e12); % 频移量 study.set('numbmodes', 5); % 算前五个模式

后处理阶段要看电场分布里的玄机。真·BIC模式的Q值会爆表，辐射损耗几乎为零。在截面图上用箭头图叠加等高线，能清晰看到电场局域在结构内部。有个冷技巧是导出模式场后做傅里叶变换，看看k空间有没有落在光锥外，这才是验证BIC的终极姿势。

% 模式验证脚本 Ex = mphinterp(model, 'emw.Ex', 'dataset', 'dset1'); [Px, Py] = fft2shift(Ex); % 自定义傅里叶变换函数 contourf(kx_grid, ky_grid, abs(Px)); % 绘制k空间分布

最后说个坑：网格剖分别太糙！BIC对结构对称性敏感得很，建议在关键区域用极端细化的三角形网格。算完记得检查模式能量是否主要分布在介质区域，别让PML把结果带偏了。用这个套路，基本能复现大部分文献里的BIC结果，不用跪求原作者代码的感觉真香。