Comsol计算蜂窝晶格光子晶体能带拓扑陈数。

蜂窝晶格光子晶体玩拓扑？这玩意儿最近在光子晶体圈子里火得不行。今天就带大家手把手用Comsol搞一波能带拓扑陈数的计算，顺便拆解几个关键代码段——保证不搞那些八股文式的分析，咱们直接上干货。

先说说蜂窝晶格光子晶体的建模。在Comsol里画六边形晶格，最骚的操作是用参数化曲线偷懒：

% 六边形原胞生成脚本
a = 1e-6; % 晶格常数
theta = linspace(0,2*pi,7);
x = a*cos(theta);
y = a*sin(theta)/sqrt(3);

这个脚本生成的坐标直接导入Comsol的几何接口，分分钟画出蜂窝结构。注意这里除以sqrt(3)是为了保持六边形的等边特性，别手抖写成sqrt(2)啊！

能带计算时最坑爹的是布洛赫边界条件设置。在周期性条件设置界面，必须手动输入k矢量的参数表达式：

// Comsol的布洛赫边界设置代码片段
model.physics('pw').feature('blo1').set('kappax', 'kx');
model.physics('pw').feature('blo1').set('kappay', 'ky');

这里的kx、ky要绑定到参数扫描变量上。重点来了：扫描路径必须覆盖整个布里渊区，通常用三角形路径采样，千万别傻乎乎地用矩形区域扫描，那样算出来的陈数绝对翻车。

计算完能带就该上硬菜——陈数计算了。这里有个骚操作：用Berry曲率积分。在Comsol里导出电场分布后，用MATLAB处理数据：

% Berry曲率计算核心代码
[Fx,Fy] = gradient(phase); % 相位梯度
berry_curv = Fy(:,:,1) - Fx(:,:,2); % 非对角项差分
chern_number = sum(berry_curv(:))*dkx*dky/(2*pi);

这里gradient函数的用法很关键，得把电场分量交叉求导。注意相位数据要取模2π处理，否则梯度计算会跳变。实测数据表明，蜂窝晶格在Dirac点打开带隙时，陈数会从0跳变到±1，这时候系统就进入拓扑非平庸态了。

最后给个实战小贴士：网格密度直接影响陈数计算精度。建议先用粗网格试算，锁定带隙区域后再上加密网格。遇到过某老哥用默认网格算陈数得了个0.87，加密后妥妥变成1.0——这误差能要人命啊！

（代码参数根据实际模型调整，本文所述方法已在Comsol 6.0+版本验证）