最近在做一个关于一维光子晶体边界态的项目最卡人的地方不是建模型也不是仿真计算而是算完能带之后别人问了一句话“你这个结构的Zak相位是多少”很多教程讲到能带图就停了可拓扑性质恰恰藏在波函数的相位信息里Zak相位就是解决这个问题的一把钥匙。这篇内容我会结合Comsol文件和Matlab程序把一维光子晶体里Zak相位的计算流程完整走一遍从原理讲到数据后处理再讲到结果验证适合那些已经能算出一维光子晶体能带、但还不清楚拓扑不变量怎么从数值解里提取出来的朋友。1. 先搞明白Zak相位到底在算啥从Berry相位到一维布洛赫系统的拓扑不变量1.1 一个看似简单但很容易绕晕的定义Zak相位是Berry相位在一维周期系统里的特殊形式。对于一个能带指数为n的一维布洛赫态周期部分是(u_{n,k}(r))在倒空间沿着整个布里渊区走一圈Berry联络的积分就是Zak相位[ \theta_n i\int_{-\pi/a}^{\pi/a}\langle u_{n,k}|\partial_k u_{n,k}\rangle dk ]其中内积是在一个原胞内对空间坐标积分。这里的(u_{n,k})不是完整波函数(\psi_{n,k}u_{n,k}e^{ikx})而是周期部分。为什么必须是周期部分因为在做k微分的时候完整波函数里的平面波因子(e^{ikx})会产生额外的贡献这个东西不是规范不变的而Zak相位作为物理量必须规范不变。很多人在这一步就开始犯错后面我会专门讲怎么从Comsol导出数据里把周期部分剥出来。1.2 能带图不包含的信息才是拓扑问题的核心一维光子晶体的能带图(E(k))只给了频率与波矢的关系里面没有波函数的相位信息。两个结构可能能带结构非常像但边界态性质完全不同。最经典的类比是紧束缚模型里的SSH链二聚化结构会给出两类拓扑Zak相位为0或为(\pi)。光子晶体虽然是连续介质不是严格的紧束缚系统Zak相位依然可以反映能带的拓扑性质而且在一维结构中还能联系到表面阻抗、反射相位、边界态的出现条件。从另一种角度看一维光子晶体的Zak相位相当于电介质极化的一维形式跟“表面电荷”或“表面态”有对应关系。这也解释了为什么很多做多层膜、一维光栅的人现在都开始算Zak相位因为它能直接判断两个光子晶体拼接起来之后带隙里会不会出现拓扑边界态。1.3 为什么不能用完整的场直接算Zak相位数值计算里经常有人图省事直接用Comsol导出的总电场(E(r))去套Berry联络公式。这样算出来的相位并不是Zak相位而是一个换上带规范依赖的混合量。原因很简单[ \langle\psi_k|\partial_k\psi_k\rangle \langle u_k|\partial_k u_k\rangle i\langle u_k|x|u_k\rangle ]第二项在周期结构里没有明确的物理意义而且它依赖于整体上的(e^{ikx})相位约定。如果你在Comsol里设置的Floquet周期条件相位约定不同这一项就会变化。所以所有Wilson loop算法都必须先确保你手里的是周期函数(u_k(r))。2. Comsol建一维光子晶体模型几何、材料和Floquet边界的正确打开方式2.1 用什么模型表示“一维结构”“一维光子晶体”指的是介电常数只在一个方向上周期变化其他方向均匀。以多层膜为例x方向是周期方向y方向是均匀延展方向。实际仿真不需要做成三维用二维模型就可以。原因是加入z方向的面外传播并用模式分析求解时同一个色散关系可以通过选取合适的极化在二维模型里表达出来。我在Comsol里一般用二维模型x方向取一个晶格常数(a)y方向取一个合适长度然后x方向加Floquet周期边界y方向也加周期边界但把Bloch波矢设成0等效于在y方向没有相位变化。这样做的思路是让y方向“看起来像个无限均匀介质”避免因为边界硬截断出现虚假的导波模式。当然y方向边界选PEC或PMC也能出结果但不如Floquet边界干净尤其当研究斜入射时会出错。我用的参数给个参考参数数值说明晶格常数 a1 μm周期方向的基本单元介质1厚度 d10.3 μm高折射率层介质2厚度 d20.7 μm低折射率层折射率 n12.5高折射率层折射率 n21.5低折射率层求解频率范围0.1–1.5 THz覆盖前几条能带当然你完全可以根据自己的材料体系改参数但只要保证d1d2a模型就是周期性的。2.2 几何与材料设置的小细节画几何时不要真的画几十个周期只画一个晶胞就够了Floquet边界条件已经包含了周期复制的信息。我做的是在矩形区域中间画一条分界线左侧材料1右侧材料2这样材料界面是竖直的。网格用映射网格x方向加密y方向均匀就好因为一维结构在y方向没有变化。材料设置里可以直接用折射率定义也可以定义介电常数。Comsol RF模块或波动光学模块都支持从折射率自动计算介电常数。要注意频率单位在求解特征频率时Comsol输出的频率单位是Hz但光子晶体能带图通常用归一化频率(a/\lambda)所以要自己在后处理里转换。2.3 Floquet周期边界的相位约定必须第一次就搞清楚选中x方向的一对相对边界在“周期性条件”里选择“Floquet周期性”。Comsol会让你填周期矢量通常是一个二维矢量。因为我们要扫的是x方向的波矢(k_x)所以周期矢量可以写成[ (k_x, 0) ]然后把这个(k_x)设成模型参数后面用参数扫描来扫。这里最关键的坑是Comsol的不同版本对“周期场”和“总场”的约定不一样。有些版本里因变量本身就是周期部分解出来的电场(E)就是(u)有些版本里因变量是总场(E\psi)。判断方法很简单查看解在x方向的相位分布。如果(E)沿x方向呈线性相位变化并且在周期边界处不满足(E(0)E(a))那就是总场如果内部相位基本不变只有在材料界面处突变边界处满足相等关系那就是周期场。这个检查决定了你后面要不要额外除以(e^{ikx})我在第4节会专门讲。2.4 特征频率扫描的一个实用技巧用参数扫描直接扫描(k_x)时会发现特征频率求解不稳定尤其在高频段。比较好的做法是先在(k_x0)处求解一次然后把上一次计算得到的特征频率作为下一次的初始猜测。在Comsol里可以在“特征频率”研究设置中打开“特征频率搜索”填入上一次的某个频率附近的值或者使用“辅助扫描”并选择“延续”。还有一个经验扫描步长尽量小一点。我一般在一个布里渊区内取41个点即从(-\pi/a)到(\pi/a)分成40段。如果结构折射率对比大有些带比较陡41个点可能不够至少要到81个点。Wilson loop对(k)的分辨率很敏感后面数值上会有体现。3. Matlab计算Zak相位Wilson loop程序的核心思路与关键代码3.1 Wilson loop在离散网格上怎么走通Zak相位计算现在主流做法是Wilson loop。把布里渊区分成(N)个小段令(k_j -\pi/a j\Delta k)(\Delta k 2\pi/(aN))。对每个能带单独定义重叠量[ M_j \langle u_{k_j} | u_{k_{j1}}\rangle ]然后累加相位[ \theta -\mathrm{Im}\left[\log\left(\prod_{j0}^{N-2} \frac{M_j}{|M_j|}\right)\right] ]相当于把连续积分离散成有限差分再通过(\log)提取虚部。这里关键是内积要完全按照原胞内的体积分计算不是随便对网格节点求加权平均。很多人栽在“直接求和”上因为有限元网格节点密度不均匀简单求和会把网格密的地方权重放大。所以我建议先用Comsol把场导出到规则网格再在Matlab里用二维梯形积分做内积。这样既绕开了有限元积分的复杂实现又保证空间积分的权重是正确的。3.2 从Comsol导出特定k点场的完整流程在Comsol的“结果”节点里使用“导出”功能把某一维k点下的解导出成文本文件。为了后续方便文件名里带上k索引比如k_001.txt、k_002.txt。导出的变量至少包含x坐标y坐标电场分量的实部电场分量的虚部如果你在Comsol里已经确认因变量是总场那么还得额外定义导出变量E_periodic_re Re(E) * cos(kx*x) Im(E) * sin(kx*x) E_periodic_im Im(E) * cos(kx*x) - Re(E) * sin(kx*x)这里的本质是计算(uE e^{-ik_x x})。注意符号约定要和Comsol的Floquet条件里的相位项保持一致。导出时可以在“数据集”里选择“网格”然后设置“数据点”为“规则网格”或者直接在Matlab里做插值。我个人更倾向于Comsol直接导出到规则网格这样后面省事很多。规则网格范围就是原胞也就是一个周期单元。3.3 Matlab脚本的主体框架下面是一段可运行的思路不是完整封装但主干逻辑足够你改写function theta zak_phase_via_wilson(Npoint, filePattern, kxList, ax) % 输入 % Npoint每个场文件里规则网格点数比如 [Nx, Ny] % filePattern文件路径模板如 field_%03d.txt % kxListk点列表单位 rad/um % ax晶格常数 um % 输出 % thetaZak相位单位 rad Nk numel(kxList); theta 0; for j 1:Nk-1 data1 readmatrix(sprintf(filePattern, j)); data2 readmatrix(sprintf(filePattern, j1)); % 假设每列是 x,y,E_re,E_im x1 data1(:,1); y1 data1(:,2); E1 data1(:,3) 1i*data1(:,4); E2 data2(:,3) 1i*data2(:,4); % 如果Comsol导出的是总场需要转到周期部分 % E1 E1 .* exp(-1i*kxList(j)*x1); % E2 E2 .* exp(-1i*kxList(j1)*x1); % 规则网格化可以用griddata或者scatteredInterpolant [Xg, Yg] ndgrid(linspace(0, ax, Npoint(1)), ... linspace(min(y1), max(y1), Npoint(2))); F1 scatteredInterpolant(x1, y1, real(E1), linear, none); G1 scatteredInterpolant(x1, y1, imag(E1), linear, none); H1 F1(Xg,Yg) 1i*G1(Xg,Yg); F2 scatteredInterpolant(x1, y1, real(E2), linear, none); G2 scatteredInterpolant(x1, y1, imag(E2), linear, none); H2 F2(Xg,Yg) 1i*G2(Xg,Yg); % 二维梯形积分计算重叠 dy mean(diff(Yg(:,1))); dx mean(diff(Xg(1,:))); ovl sum(conj(H1(:)) .* H2(:)) * dx * dy; % 归一化并累加相位 ovl ovl / abs(ovl); theta theta - imag(log(ovl)); end theta mod(theta, 2*pi); if theta pi theta theta - 2*pi; end end这个脚本有几个假设第一每个k点的网格坐标一致所以可以直接用第一个文件的坐标第二导出时已经确保每个文件覆盖同一个原胞第三没有处理能带简并。实际使用中这三个假设都需要验证尤其第三条。3.4 如果简并能带出现了怎么办一维光子晶体在布里渊区边界有时会出现两支模式简并或者由于结构对称性在某些k点出现交叉。这时候单独的Wilson loop按“每支能带分别算”不再严格有效应该处理简并子空间。做法是把每支模式的场并成列向量即对第j个k点构一个矩阵[ W_{mn}^{(j)} \langle u_{m,k_j}|u_{n,k_{j1}}\rangle ]然后累加这个矩阵的乘积Zak相位从乘积矩阵的行列式或者展开后的相位里取。对应Matlab里需要导出每个k点下的多个模式场还要识别哪些模式属于同一个简并空间。这个工作量会大不少但如果你要做的是普通无简并一维结构一开始可以不碰。4. 数据后处理里最容易翻车的五个坑4.1 坑一分不清导出的是总场还是周期场这个坑我在2.3节提过但值得单独拿出来说因为它的后果非常隐蔽。你算出来的Zak相位会整体偏差一个随k线性变化的相位项更难受的是这个偏差不稳定Wayne会造成结果看起来有非零值但实际是伪相。我在第一次计算时就吃过这个亏结构反射相位验证死活对不上后来用场周期性的检查才发现Comsol导出的是总场。避免办法是在Matlab里把所有场都乘上(e^{-ikx})然后检查它在x0和xa处是否相同。如果相同说明你现在手头的是周期场如果不同继续调整相位约定。不同版本甚至同一版本不同模块都可能不一样所以这一步不要省。4.2 坑二Wilson loop里的“内积”没有按物理积分来内积的定义是在原胞上的二维积分。如果直接对有限元网格节点做sum(conj(E1).*E2)很容易因为网格不均匀而出现偏差。我做过一个对照测试同样的场用直接节点求和算出的Zak相位和用规则网格梯形积分算出的结果差了0.2 rad左右这足以让一个原本(\pi)的边界态相位判断错掉。为了避免这个问题Comsol里面的网格明明很密却还是建议你在Matlab里做一次规则化。线性插值虽然会引入一些数值平滑但对前几条能带的影响远小于积分权重错误。4.3 坑三k点采样太粗Wilson loop相位没有收敛Zak相位是路径积分离散化误差和采样密度直接相关。如果你用20个k点算出来是1.3π50个k点变成0.95π那肯定是采样问题。我一般至少用40个点起步再多取一倍验证收敛性。不是每个结构都收敛那么快特别是带比较平缓、相邻k点场重叠接近1的时候相位提取更敏感。判断收敛的办法很粗暴算一次(N40)再算一次(N80)看结果差是否小于0.01 rad。对于需要精确判断0或(\pi)的系统这个容差往往比相位本身对边界态的影响更值得关注。4.4 坑四模式排序在k扫描时发生跳变Comsol特征频率求解器在参数扫描时会对每个k值列出所有特征频率但模式编号不一定按能带顺序。某个模式在k0是第3条带到k0.2π/a时可能跳到第4条带因为你没有跟踪同一条能带的演化。如果你直接按模式编号取场Wilson loop会把两条不同能带的波函数混在一起结果毫无意义。解决方法是做“模式连续追踪”用上一k点的场和当前k点每个模式的场做内积重叠模最大者判定为同一条能带延续。这个可以在Matlab里自动完成核心代码类似[~, idx] max(abs(overlap_matrix));但如果遇到简并点这个贪心判断会失效需要用简并子空间的方法。4.5 坑五符号约定不统一导致Zak相位反号不同文献对Zak相位的定义有时相差一个负号。比如有的定义是(\theta i\int\langle u|\partial_k u\rangle dk)有的则取负号。Comsol的Floquet条件相位项符号也会影响导出场的定义。所以在你开始写论文或者对比文献之前先做一个已知结构的验证比如SSH模型对应的一维介质结构在标准参数下Zak相位应该从0跳到(\pi)。如果你的验证结果明显反了先把整体符号调过来不要继续往下摸索。符号检查还有一个技巧看重叠矩阵的幅角在带隙中的连续性。物理Zak相位在能带内部是连续变化的如果计算出来的相位出现无意义的跳变先怀疑符号约定或模式混叠。5. 用反射相位和边界态验证计算结果5.1 为什么必须做交叉验证Zak相位是个积分量你算出一个数很容易但怎么知道对不对最好的办法不是信程序而是用另一个独立物理量去交叉验证。对一维光子晶体反射相位和拓扑边界态是两个最常用的验证窗口。如果你有传输矩阵法的解析结果可以直接把Wilson loop算出的Zak相位和反射相位在带隙内的行为做对比。当表面截断位置变化时带隙内的反射相位会和Zak相位发生关联尤其是半无限晶体表面为高折射率层和低折射率层时往往对应不同的拓扑相位标记。这个关系能非常敏感地暴露数值计算错误。5.2 构建两个拓扑不同的结构在Comsol里找边界态更直观的验证方式是找边界态。选两个一维光子晶体A和B它们的带隙有重叠但Zak相位在某个共同带隙里不同。把A和B沿x方向拼接在一起构成一个超原胞x方向两端还是用周期边界条件然后做特征频率扫描。如果A和B的Zak相位不匹配理论上会在界面附近出现局域在带隙中的模式频率落在A和B的共同带隙内。具体在Comsol里操作时模型里画两个区域左侧是一定周期的A结构右侧是一定周期的B结构中间界面贴合。整个超原胞的长度可以是两个结构周期的最小公倍数或者直接取足够长的组合然后加周期边界。网格不需要太密但界面附近需要加密。我一般会取左右各4到5个周期否则超原胞带会被有限尺寸效应污染边界态不够干净。找到边界态后检查场分布应该集中在界面上且频率在带隙内。这一步如果匹配不上八成是Zak相位计算的符号或模式追踪出了问题回到第4节逐项排查。5.3 一个快速验证的小例子假设你已经算好了单个原胞的能带。取一个结构参数计算前三条带的Zak相位。然后改变表面层的厚度重新算同一能带范围的Zak相位。比较两个结构的相位值再拼接这两个结构观察边界态。一个常见的规律是如果两个结构在公共带隙内的Zak相位相差(\pi)边界态会出现在带隙中央附近如果相位相同边界态会消失或移出带隙。这个规律不是普适的绝对判据但用来做自洽检查非常有用。你在自己的结构里如果发现“相位差(\pi)但没有边界态”或“相位差0却有边界态”那基本可以判断有一处计算流程出问题了。6. 我对Zak相位数值计算的一点习惯这套流程跑通之后我现在已经形成固定工作流Comsol扫完能带Matlab跑Wilson loop再用反射相位或边界态做交叉验证。整个链路里最消耗时间的往往不是仿真本身而是排查那些微妙的相位约定和模式排序问题。所以我建议你第一次做的时候不要急着追求高精度先用20个k点跑通再加密到40个点对比一次。顺手把每个阶段的结果存下来这样后面检查问题时能快速定位到底是场导出错了、内积错了还是模式追踪错了。最后分享一个实用小技巧在Comsol里做特征频率扫描时打开“使用上一解作为初始猜测”这个选项能够大幅减少模式跳变的情况。配合Matlab里的重叠追踪能带连续性基本不会出大问题。等这一步稳了你就可以比较放心地把Zak相位当作一个常规工具去用了。