一维光子晶体的zak相位计算 内含comsol文件和matlab程序 注意这个是重复别人文章的结果方法是论文中所提到的今天咱们来唠唠一维光子晶体Zak相位的计算实操。这玩意儿听起来挺玄乎其实就是个描述拓扑特性的数学量但算起来可真是让人头秃。先扔个重点用COMSOL建模Matlab后处理完整复现文献里的结果。文末附文件自取先别急着关网页看完保证你少踩三个坑。先上COMSOL建模核心设置。在波动光学模块里定义周期结构的时候记得把晶格常数a设为变量后面扫参方便。材料交替堆叠的介电常数要严格对应文献参数这里有个骚操作用矩形函数定义空间依赖的ε分布比手动画几何省事十倍。关键代码长这样% 结构参数 a 1e-6; % 晶格常数 d 0.4*a; % 高介电层厚度 epsilon1 13; % 高介电材料 epsilon2 1; % 低介电材料 % 空间坐标函数 epsilon (x) epsilon1*(mod(x,a)d) epsilon2*(mod(x,a)d);建模时边界条件最容易翻车。上下边界用完美电导体PEC左右用周期性边界条件。重点来了端口激励要设置成布洛赫边界条件kx从-π/a扫到π/a这步错后面全白给。网格划分建议用极端细化特别是介电突变处别心疼算力精度第一。一维光子晶体的zak相位计算 内含comsol文件和matlab程序 注意这个是重复别人文章的结果方法是论文中所提到的数据导出后上Matlab处理。Zak相位计算的核心是积分布洛赫波函数的导数但数值计算要转成离散求和。这里有个神仙操作用unwrap函数处理相位跳变比手动判断省心。关键代码段phi angle(eigenvectors); % 取本征矢量相位 phi_unwrap unwrap(phi,[],2); % 解相位包裹 dphi diff(phi_unwrap,1,2); % 差分求导 Zak_phase sum(dphi,2)/(2*pi); % 积分归一化注意本征矢量排序问题COMSOL输出的模式顺序可能乱序必须按频率排序后再处理。建议用sortrows搭配eigenvalues别问我怎么知道的都是血泪教训。最后验证环节当结构参数对称时Zak相位应该是0或π。如果算出0.5π这种妖孽值先检查本征矢量的归一化是否正确八成是边界条件里的kx范围设反了。用文献中的参考图对比时注意频段缩放比例特别是高频段容易出现模式交叉导致相位跳变。文件包里已经调好参数的.mph和.m文件把晶格常数改成你的实际值就能跑。遇到计算结果抽风时重点检查1. 周期性边界是否真周期 2. 端口激励的kx扫描方向 3. 介电分布函数的模运算是否正确。祝各位算力全开相位不乱