手撕BIC：从能带仿真到拓扑电荷计算
nature文章中的BIC能带仿真计算 ，包括能带计算Q因子计算，拓扑电荷计算，包括询问等，所见即所得。
搞过光子晶体的老铁都知道，BIC（连续体束缚态）这玩意儿就像二维材料里的幽灵——明明应该辐射损耗却稳如老狗。今天咱们直接上手撕代码，用Python+MPB（MIT Photonic Bands）复现Nature级别的BIC仿真。
第一步：搭建光子晶体平板
先整个三角晶格光子晶体结构，介质柱半径r=0.3a（a是晶格常数），介电常数ε=12。MPB脚本里核心参数这样写：
`python
geometry = [mp.Cylinder(radius=0.3, material=mp.Medium(epsilon=12))]
geometry_lattice = mp.Lattice(size=mp.Vector3(1,1,0.3)) # 平板厚度
k_points = mp.interpolate(9, [mp.Vector3(0.5,0), mp.Vector3(0,0), mp.Vector3(0.5,0.5)])
`
这里geometrylattice的z方向厚度控制着平板模式，kpoints沿着Γ-M-K高对称点扫频。跑完能带后重点关注TM模式在Γ点附近的反交叉现象——这就是BIC的藏身之处。
Q因子暴力计算法
找到BIC频率后，直接在时域用FDTD算Q值：
`python
sim = mp.Simulation(resolution=20, cell_size=mp.Vector3(5,5,3),
sources=[mp.ContinuousSource(frequency=freq_BIC)],
boundary_layers=[mp.PML(1.0)])
sim.run(until=200) # 先跑稳态
decay_data = []
def save_decay(sim):
decaydata.append(sim.getfield_energy())
nature文章中的BIC能带仿真计算 ，包括能带计算Q因子计算，拓扑电荷计算，包括询问等，所见即所得。
sim.run(mp.aftertime(20, mp.atevery(0.1, save_decay)), until=500)
`
抓取场能量衰减曲线做指数拟合，Q=ωτ/2（τ是衰减时间常数）。注意要等谐振稳定后再开始采样，否则会吃到瞬态响应的噪声。
拓扑电荷的相位手术
计算拓扑电荷得玩转Berry相位，这里有个取巧的办法——直接在动量空间画圈：
`python
def calccharge(kxlist, ky_list):
phases = []
for kx in kx_list:
for ky in ky_list:
# 计算该k点的本征模式相位分布
phasemap = getphase_map(kx, ky)
phases.append(np.angle(phase_map))
# 计算绕数
delta_theta = np.diff(phases, axis=0)
winding = np.sum(np.abs(delta_theta) > np.pi)
return winding // (2*np.pi)
`
核心是提取TE/TM模式的电场相位分布，沿着闭合路径积分相位变化。当绕数不为零时，恭喜你逮到了非平庸拓扑电荷！
避坑指南
- MPB计算二维结构时记得设置
num_bands足够多（至少8个），防止模式遗漏 - Q因子超过1e6时改用频域法，时域衰减法会吃到数值误差
- 拓扑电荷对网格分辨率极其敏感，建议用自适应k点采样
跑完这三板斧，你的BIC基本就能在Nature子刊混个Supplementary Fig了。不过要真发顶刊，记得在结构不对称性上多做文章——毕竟完美的BIC就像理想气体，只存在于教科书里。