nature文章中的BIC能带仿真计算 包括能带计算Q因子计算拓扑电荷计算包括询问等所见即所得。搞过光子晶体的老铁都知道BIC连续体束缚态这玩意儿就像二维材料里的幽灵——明明应该辐射损耗却稳如老狗。今天咱们直接上手撕代码用PythonMPBMIT Photonic Bands复现Nature级别的BIC仿真。第一步搭建光子晶体平板先整个三角晶格光子晶体结构介质柱半径r0.3aa是晶格常数介电常数ε12。MPB脚本里核心参数这样写pythongeometry [mp.Cylinder(radius0.3, materialmp.Medium(epsilon12))]geometry_lattice mp.Lattice(sizemp.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值pythonsim mp.Simulation(resolution20, cell_sizemp.Vector3(5,5,3),sources[mp.ContinuousSource(frequencyfreq_BIC)],boundary_layers[mp.PML(1.0)])sim.run(until200) # 先跑稳态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)), until500)抓取场能量衰减曲线做指数拟合Qωτ/2τ是衰减时间常数。注意要等谐振稳定后再开始采样否则会吃到瞬态响应的噪声。拓扑电荷的相位手术计算拓扑电荷得玩转Berry相位这里有个取巧的办法——直接在动量空间画圈pythondef 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, axis0)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就像理想气体只存在于教科书里。