最近在做一个关于一维光子晶体边界态的项目，最卡人的地方不是建模型，也不是仿真计算，而是算完能带之后别人问了一句话：“你这个结构的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边界干净，尤其当研究斜入射时会出错。我用的参数给个参考：

当然你完全可以根据自己的材料体系改参数，但只要保证d1+d2=a，模型就是周期性的。

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_x=0)处求解一次，然后把上一次计算得到的特征频率作为下一次的初始猜测。在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_{j+1}}\rangle ]

然后累加相位：

[ \theta = -\mathrm{Im}\left[\log\left(\prod_{j=0}^{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)

这里的本质是计算(u=E 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' % kxList：k点列表，单位 rad/um % ax：晶格常数 um % 输出： % theta：Zak相位，单位 rad Nk = numel(kxList); theta = 0; for j = 1:Nk-1 data1 = readmatrix(sprintf(filePattern, j)); data2 = readmatrix(sprintf(filePattern, j+1)); % 假设每列是 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(j+1)*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_{j+1}}\rangle ]

然后累加这个矩阵的乘积，Zak相位从乘积矩阵的行列式或者展开后的相位里取。对应Matlab里，需要导出每个k点下的多个模式场，还要识别哪些模式属于同一个简并空间。这个工作量会大不少，但如果你要做的是普通无简并一维结构，一开始可以不碰。

4. 数据后处理里最容易翻车的五个坑

4.1 坑一：分不清导出的是总场还是周期场

这个坑我在2.3节提过，但值得单独拿出来说，因为它的后果非常隐蔽。你算出来的Zak相位会整体偏差一个随k线性变化的相位项，更难受的是这个偏差不稳定，Wayne会造成结果看起来有非零值，但实际是伪相。我在第一次计算时就吃过这个亏，结构反射相位验证死活对不上，后来用场周期性的检查才发现Comsol导出的是总场。

避免办法是在Matlab里把所有场都乘上(e^{-ikx})，然后检查它在x=0和x=a处是否相同。如果相同，说明你现在手头的是周期场；如果不同，继续调整相位约定。不同版本甚至同一版本不同模块都可能不一样，所以这一步不要省。

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的时候，相位提取更敏感。

判断收敛的办法很粗暴：算一次(N=40)，再算一次(N=80)，看结果差是否小于0.01 rad。对于需要精确判断0或(\pi)的系统，这个容差往往比相位本身对边界态的影响更值得关注。

4.4 坑四：模式排序在k扫描时发生跳变

Comsol特征频率求解器在参数扫描时，会对每个k值列出所有特征频率，但模式编号不一定按能带顺序。某个模式在k=0是第3条带，到k=0.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相位当作一个常规工具去用了。