最近在整理手头一个跟光学超构表面相关的课题，需要把一维光子晶体能带里的拓扑不变量——Zak 相位——用数值方法算出来。坦白讲，这个量在拓扑光子学文章里出现频率很高，但真正落到计算上，比教材里那行积分公式要折腾得多。我最后采用的方案是 Comsol 和 Matlab 联合：Comsol 负责建模、扫 Bloch 波矢、求解本征模场，Matlab 负责导数据、追踪能带、算 Wilson loop 并提取相位。整套流程跑通之后效果很稳，今天把过程和关键坑都记录下来，给有同样需求的小伙伴一条能直接照着走的路。
1. 一维光子晶体的Zak相位：为什么值得专门写一套流程
1.1 Zak相位到底在描述什么
对沿 x 方向周期排列的一维光子晶体，布洛赫定理告诉我们本征场可以写成E_k(x) = u_k(x) exp(i k x)，其中u_k(x)是周期函数，k是波矢。Zak 相位就是给定能带在布里渊区上积累的 Berry 相位，写成公式就是：
θ_n = i ∮ ⟨u_{n,k} | ∂_k u_{n,k}⟩ dk
积分路径取整个布里渊区，也就是从-π/a到π/a。这里的a是周期。如果系统具有时间反演对称性，Zak 相位只能取两个值：0或π。这个离散化的性质非常重要，因为它直接和边界处是否出现局域表面态挂钩。两块相邻光子晶体如果带隙重叠、但对应能带的 Zak 相位不匹配，界面上会出现拓扑界面态；匹配则没有。这在设计拓扑光子学器件时是实打实的判据，不是可有可无的理论点缀。
很多入门教程会告诉你 Zak 相位能在纯解析模型里算出来，比如二元交替层状介质，确实可以手推。但实际操作中，往往会碰到各种变形结构——渐变层、缺陷腔、非平面界面，甚至各向异性材料。这时候基于简单解析解的代码就不够用了，需要一个能处理任意几何的电磁场求解器，这就是 Comsol 存在的意义。
1.2 纯写代码能算，但为什么要拉上 Comsol
纯用传输矩阵法或者平面波展开法也能得到一维光子晶体的能带和 Bloch 函数，Zak 相位自然也能算。问题是这些方法每一步都要自己维护色散关系、边界条件、模式归一化，一旦几何结构偏离"平板/真空/平板"的理想模型，开发成本会迅速膨胀。Comsol 的优势在于，你只需要把周期单元画出来，设置 Floquet 周期条件，它就能返回本征频率和完整的空间场分布。Matlab 再接手做拓扑不变量计算，正好互补。
不过这种分工有一个隐含门槛：Comsol 默认输出的电场包含 Floquet 相位因子，直接拿去算内积并不对，必须先还原出周期函数u_k(x)，并且在整个 k 扫描过程中保证同一能带不被搞混。这两个问题如果不处理好，算出来的 Zak 相位就会在0和π之间乱跳，甚至出现中间值。下面我的记录重点就是解决这两件事。
1.3 一个可复现的参考模型参数
为了后面讨论方便，先给出一组我实际用过的参数。结构选用最常见的 AB 二元交替介质：
- 材料 A：空气，折射率
n_A = 1.0，厚度d_A = 600 nm - 材料 B：硅，折射率
n_B = 3.48，厚度d_B = 400 nm - 周期常数
a = d_A + d_B = 1000 nm
计算采用电场垂直于模拟平面的偏振，对应 TE-like 模式，只需要关心电场在 z 方向的分量Ez。频带范围大致在 100 THz 到 350 THz，对应归一化频率a/λ大约在 0.33 到 1.17，足够看到前几条能带和带隙。
这套参数本身没有特殊的物理意义，只是因为结构简单、周期尺度适合在近红外波段做有限元仿真，同时对比传输矩阵结果也方便。你完全可以根据自己的波段改厚度和折射率，后处理流程不用动。
2. Comsol端建模：决定Zak相位成败的关键设置
2.1 用二维模型表示一维周期延伸结构
要模拟一维光子晶体，严格来说需要做一个沿 x 方向无限周期、y 和 z 方向均匀的模型。在 Comsol 中我建议直接用 2D 组件，画一个宽度等于周期常数a、高度可以任意取的矩形作为单元胞。为什么可以这样？因为对于 2D 电磁波模型，第三个方向默认是无限延伸的，我们只需要在高度方向加一个周期条件，让 y 方向也"假装"无限，就能消除有限高度带来的波导截止效应。
具体设置时，几何里创建一个矩形，宽度设为a = 1 μm，高度设为h = 0.1 μm。这个高度不要取得太大，否则 y 方向可能出现高阶横向模式干扰特征值搜索；也不要太小，太小会让有限元网格产生病态。之后在 x 方向的两条边设置 Floquet 周期条件，在 y 方向的两条边也设置 Floquet 周期条件，但波矢的 y 分量固定为 0。
物理场使用电磁波、频域接口。在 2D 中默认有面内分量和面外分量两种模式选择，这里选面外电场矢量，也就是只有Ez分量。这样求出来的本征模式是 TE-like，电磁场分布沿 y、z 均匀，只在 x 方向呈布洛赫振荡，恰好符合一维光子晶体的假设。
2.2 Floquet 周期边界条件的波矢怎么填
Comsol 的周期条件里有几个选项，一定要选 Floquet 周期条件（有的版本写成 Bloch-Floquet），而不是默认的"周期性"或者"连续"。在 Floquet 设置里需要指定空间相位因子。对一维结构，只需要填 x 方向的波矢分量kx，y 方向给 0。
这里的单位是rad/m，不是归一化波矢。我第一次建模型时直接把π/a填进去了，结果能量明显不对，翻回来才发现少了1/a的量纲。建议在全局参数里先定义一个归一化扫描量kk，取值范围[-1, 1]，然后在周期条件里填：
kx = kk * pi / a
这样扫描关系一目了然，也方便后面 Matlab 端统一计算。
2.3 特征频率研究的搜索范围和模式数
研究步骤选"特征频率"，这是求本征模式的正规路子。里面有几个参数需要认真处理：
- 特征频率搜索基准点：我会先做一次单点求解，比如
kk = 0.3，看目标频带落在哪个频率附近，然后把搜索基准点设在那里。 - 待搜索的特征数：至少要大于你想研究的能带数。因为有限元特征值求解会返回一堆模式，除了物理上成立的布洛赫模式，还可能出现角点奇异导致的伪模。我通常设成目标能带数的 1.5 到 2 倍。如果搜索范围太宽，伪模数量会爆炸，反而拖慢求解。
- 搜索范围：可以用频率区间限制，也可以直接搜基准点附近的若干个模式。建议先扫一两个 k 点看能带的大致频率范围，再决定区间，否则容易漏带。
还有一个非常重要的细节：在整个 k 扫描期间，网格千万不能变。如果每次修改参数后 Comsol 因为某些原因重新划分网格，本征场的插值位置就会漂移，后面算 Wilson loop 时相邻 k 点的重叠积分会引入不必要的噪声。因此在批处理设置里，要把网格划分排除在研究序列之外，只对参数kk做扫描。
2.4 网格无关性检查
一维光子晶体往往有高折射率对比，比如硅和空气，折射率差达到 3.48 倍，电场在界面上变化非常陡。这种情况下网格要保证每个介质层内部至少有 20 到 30 个单元，尤其界面附近要加密。
我的检查方法很简单：先在较粗网格下把 Zak 相位算一遍，再加密网格计算一遍，看结果是否稳定在0或π。如果两次结果一致，说明网格已经收敛；如果相位值飘在0.3π、0.7π这类中间值，基本是网格不够细或者 k 点不够密。网格数量不是越多越好，我测试下来每层 50 个单元加 10 层边界加密，就已经足够稳定，再多只是增加计算时间。
3. 从Comsol到Matlab：本征场到Zak相位的完整链路
3.1 用LiveLink把数据导到Matlab的正确姿势
Comsol 和 Matlab 联合有几种方式，最便捷的是 LiveLink for MATLAB。启动后在 Matlab 命令窗口输入mphstart，然后打开模型：
model = mphopen('pc_zak.mph');之后每个 k 点循环里，修改参数、运行研究、提取场数据：
Nk = 60; klist = linspace(-pi/model.param.get('a'), pi/model.param.get('a'), Nk); for j = 1:Nk model.param.set('kk', klist(j) * model.param.get('a') / pi); model.study('std1').run(); % 提取 x 方向均匀采样点上的 Ez 分量 xq = linspace(0, model.param.get('a'), 501); yq = zeros(size(xq)); Ez = mphinterp(model, 'ewfd.Ez', 'coord', [xq; yq], 'dataset', 'dset1', 'solnum', 1); % 这里只演示 solnum=1，实际需要循环所有模式 end有几个坑需要提前知道。第一，mphinterp返回的是一组按求解模式索引排列的数据，特征值研究里每个特征频率对应一个solnum，所以至少要循环所有目标模式次数。第二，comsol 的特征向量存在任意整体相位，也就是说每个模式的全局相位因子是随机的，这不会影响后面 Wilson loop 的相位结果，但会影响你能看到的场图。第三，如果模型里存在多个数据集，建议在 Comsol 端就固定一个默认数据集，否则 Matlab 端容易取错。
版本兼容性也是个容易忽略的问题。COMSOL 6.x 和 Matlab 各版本之间的联调接口并不总是开箱即用，建议在开始之前查一下当前 Comsol 版本官方支持的 Matlab 版本列表。我遇到过一次mphstart后连接不上的情况，换了匹配版本后一切正常。
3.2 还原Bloch函数以及符号约定
Comsol 求解出的Ez是完整布洛赫波场，其中已经包含了 Floquet 因子。为了得到周期函数u_k(x)，需要把该因子除掉。我这里采用的约定是全场形式为：
E_k(x) = u_k(x) exp(-i k x)
所以周期函数为：
u_k(x) = E_k(x) exp(+i k x)
注意这个符号不是绝对的。不同版本的 Comsol，甚至不同物理场接口，Floquet 条件里相位因子的定义可能差一个符号。如果搞反了，后面的 Wilson loop 会直接得到错误相位，而且不会自动修复。
怎么确认符号对不对？很简单：算完u_k(x)后，检查它在单元胞两端是否相等。因为u_k必须是周期函数，所以在x=0和x=a处的值应当一致，误差在数值容忍范围内。如果两端的幅值一样、相位差接近零，说明符号选对了；如果出现接近2k a的相位差，那就把指数符号反过来再试。
这一步我当时花了很长时间才想明白。原因是 Comsol 的文档对 Floquet 因子的写法写在了不起眼的边界条件说明里，很容易略过。所以这个经验必须记下来：不要盲信任何教程，永远用"周期函数"这个物理约束去检验自己的符号约定。
3.3 能带追踪：模式顺序交错是最大干扰
每次特征值求解返回的模式在solnum里的排列顺序，本应该按照特征频率从小到大排。但问题在于，有限元求解器返回的顺序在正常情况下是按频率排序的，可是当你把不同 k 点的结果放到一起时，某个能带在 k 点 1 是第 5 个求解模式，到了 k 点 2 可能变成第 6 个。如果直接按模式序号去连接能带，Wilson loop 会算出一堆乱七八糟的内积。
解决办法是做"重叠矩阵最大匹配"。对相邻两个 k 点的本征模式集合{u_m(k_j)}和{u_n(k_{j+1})}，先计算所有模式对之间的重叠积分矩阵：
S_{mn} = ⟨u_m(k_j) | u_n(k_{j+1})⟩
然后对每一行取绝对值最大的列索引作为下一 k 点应该对应的模式序号。这个过程等价于人眼追踪能带走向，但用矩阵运算自动完成。当能带远离简并时非常稳定，在带边附近接近简并时会误判，这时候需要提高 k 点采样密度。
具体到实现，我会先按能量从低到高排序一次，然后逐段做最大匹配。注意能带在布里渊区边界处可能存在级数交换，比如第 1 带和第 2 带在kx = π/a附近若发生反交叉，追错就会导致后续相位全错。处理原则是：宁可多增加 k 点，也不要让相邻 k 点频率差太大。我一般把Nk设为 60 到 120，具体看带边复杂程度。
3.4 Wilson loop：Zak相位的离散化计算
一旦模式排序搞定了，Zak 相位就可以通过相邻 k 点之间的内积连乘得到，这就是 Wilson loop 的离散形式：
W = ∏_{j=1}^{Nk-1} ⟨u(k_j) | u(k_{j+1})⟩
对单带情况，W是一个复数，Zak 相位取：
θ = -Im(log(W)) = -angle(W)
为什么可以这样算？因为 Berry 联络的积分可以写成相邻态重叠积分的连乘，这是平行移动思想在数值上的直接实现。当采样点足够密，相邻重叠积分趋于 1，连乘的对数虚部就趋于连续积分的结果。
一个关键点是最后首尾要不要连接。我们扫描的区间是从-π/a到π/a，这是整个布里渊区，但由于k = -π/a和k = π/a实际上是同一个物理点（相差一个倒格矢），所以严格来说应当把最后一个点和第一个点也做一个重叠积分，构成一个闭合环。不过实际计算中，两个端点对应的模式可能因为解算器给出的规范不同，在布洛赫函数上差一个相位，闭合积分可以吸收这个相位差，取对数得到的是模2π的相位。这正是拓扑不变量需要的。因此我的代码会额外把j=Nk和j=1连起来算，不能漏。
下面是一个简化版的 Matlab 片段，展示单能带 Wilson loop 的核心思路：
% u_all{n}(k_index, x_index) 已经存储了归一化的布洛赫函数 % 假设能带追踪已经完成，band_index 是当前要算的能带 W = 1; for j = 1:Nk jp = mod(j, Nk) + 1; % 环闭合 u1 = u_all{band_index}(j, :); u2 = u_all{band_index}(jp, :); ov = sum(conj(u1) .* u2); % 等距采样用简单求和近似积分 W = W * ov; end theta = -angle(W);实际使用中，等距采样点的内积应该带有积分权重，如果采样点足够多，均匀间距下权重因子相同，可归一化后直接求和，精度足够。更多高阶做法是每段都用梯形法积分，但 Zak 相位对采样密度的依赖不强，密集取样后结果基本一致。
3.5 归一化和相位展开问题
u_k计算出来以后先要做归一化，否则重叠积分的幅度不是 1，连乘之后会偏离纯相位。标准做法是：
u = u / sqrt(sum(abs(u).^2));这是每个 k 点、每个模式都要做的。做完后相邻重叠积分的模应该非常接近 1，如果偏离 1 太多，说明网格不够密或采样点不足。
相位展开也是个小坑。-angle(W)返回的值域在[-π, π]，由于我们期望结果在0或π附近，这个范围已经够用。但如果后续要计算带隙边界处伴随的反射相位，可能需要连续展开相位曲线，那时候建议用unwrap函数处理相邻频率点的相位值。
4. 联调中踩过的坑：从错误结果到稳定复现
4.1 Floquet因子符号导致的π误差
前面提过符号约定问题，这里再说一下它造成的典型症状。我第一次跑完整流程时算出来的 Zak 相位不是0/π，而是0.7π和1.3π这种奇怪值。检查了很久，最后发现是布洛赫函数还原时用错了指数符号。把exp(+i k x)换成exp(-i k x)后，所有结果立刻变成清晰的0/π。这个教训让我意识到，凡涉及相位类计算，最先排查的一定是约定符号，而不是网格或算法精度。
4.2 模式排序错误造成的虚假拓扑相变
另一个高频问题是模式排序错误。表现是：稍微改一下 k 点数量或者网格密度，Zak 相位就从0跳到π，似乎在某个参数处发生了拓扑相变。其实这不是物理现象，是模式追踪链条断裂了。解决办法是使用重叠矩阵排序，并且在排序后把能带图重画一遍，肉眼确认没有断裂。我会在算 Zak 相位前，先把相邻 k 点的模式顺序连线画出来，确认所有能带连续，再进入相位计算。这一步成本很低，收益却巨大。
4.3 单元胞边界切割位置是否影响结果
理论上，Zak 相位作为一个体态拓扑不变量，不取决于你选择单元胞的哪个位置作为边界。比如你可以把单元胞从 A/B 中间切开，也可以从某个介质层的中心切开，算出来的闭合 Wilson loop 应该相同。但数值上如果你抽样范围没有完整覆盖一个周期，或者坐标取点不准确，结果就会飘。我之前因为几何建模时把 A 层放在左边、B 层放在右边，但坐标原点恰好落在 B 层中间，导致u_k(x)在端点上不严格相等，虽然 Wilson loop 闭合后相位误差很小，但反应到反射相位验证时就对不上。后来我把模型统一调整为从 A 层和 B 层界面处开始一组完整的 A+B 周期，所有结果都干净了。
4.4 特征频率搜索范围里混入伪模
特征频率研究返回的模式，除了物理能带，还常常包含少数非物理数值模式。典型特征是在某些点出现局部场尖峰，或者场分布完全不符合布洛赫模式的空间轮廓。这些伪模混进 mode sorting 后，会严重干扰重叠矩阵的最大匹配。解决办法有两个方向：
- 计算之前限制频率搜索范围，避免把太高频率的伪模也搜进来。
- 计算之后根据场分布过滤掉那些空间变化异常剧烈的模式。比如一个模式的电场主要集中在一个窄条带上，而不是全周期分布，基本可以判定是伪模。
我通常会在 Comsol 端先画出几个候选模式的场图，大致确认哪些模式看起来"像物理模式"，哪些是伪模，然后再跑全扫描。千万不要闷头跑全批量，最后拿到一堆数据再清理，那会非常痛苦。
4.5 批量扫描的效率优化
用 LiveLink 循环跑 60 个 k 点、每个点求 8 个模式，如果每次都启动一个完整的特征值求解，耗时非常可观。我的优化建议：
- 每次求解前不要重建几何和网格，固定网格后只更新
kx参数。 - 使用上一步的解作为当前步的初始猜测，可以显著加速带边附近的求解。
- 将研究设置为只更新参数并求解，不要执行"初始化研究"之类的操作。
- 如果条件允许，直接把 COMSOL 批处理放到服务器上跑，输出结果文件再给 Matlab 做后处理。
实测下来，60 个 k 点、每个点 8 个模式的扫描，优化后耗时能降低一半以上。但要注意，使用上一个解做初始猜测时，如果某个 k 点模式变化太大，可能迭代到局部解，所以结果要在能带图上校验一下。
5. 验证和扩展：怎么确定算出来的Zak相位可信
5.1 和传输矩阵法对表
Zak 相位算完之后，最直接的验证是用独立方法对照。传输矩阵法是二维交替层状介质最经典的全解析算法。流程大致是：
把每层介质写成 2×2 转移矩阵，周期结构的单胞总矩阵为M_total，能带色散由cos(K a) = (M_total(1,1) + M_total(2,2))/2给出。这里K是布洛赫波矢。把K反解出来，画成能带曲线，和 Comsol 的特征频率曲线重叠，两者误差通常在 1% 以内。
TMM 同样可以给出 Bloch 函数的数值形式，进而用同样的 Wilson loop 流程算出 Zak 相位。我拿自己的 TMM 代码和 Comsol/Matlab 流程对照，前四条能带的 Zak 相位完全一致。这种对照虽然不能百分之百保证你 Comsol 模型没有对称性破坏之类的问题，但至少能确认后处理链路没问题。
5.2 用反射相位判别法做物理交叉验证
一个更偏物理的验证方法是用半无限光子晶体的反射相位。反射相位与 Zak 相位存在明确对应关系：某个带隙的反射相位在带隙两端的变化方式，取决于产生这个带隙的两条能带各自 Zak 相位的相对关系。
具体表现为，从无穷大结构侧面正入射一个平面波，计算反射系数相位随频率的变化轨迹。当频率扫过某个带隙时，如果反射相位连续光滑，说明带隙两端 Zak 相位匹配；如果反射相位在带隙中心附近发生 π 量级的跳变，说明两端 Zak 相位不同。用这个方法可以非常直观地判断"带隙是否拓扑非平庸"。
我当时在 Comsol 里额外建了一个截断结构模型，右端是有限周期光子晶体，左边是空气，然后频域求解反射系数。反射相位轨迹和 Zak 相位判断一致，这一步做完我才觉得结果可信。
5.3 超胞模型做表面态验证
还有一个御三家级别的验证方法：直接找表面态。从拓扑角度，截断光子晶体的带隙中若存在表面态，通常对应界面两侧的某个 Zak 相位不匹配。做法是把光子晶体切成有限长度，比如 9 个周期，左右都是空气，然后做一个超胞能带计算。特征值结果里，如果带隙内部出现平直的本征模式，说明这些频率就是被局域在表面的态。这个验证会在视觉上非常直观，也方便你后续做实际器件设计。
代价是超胞模型的计算规模比单胞大不少，网格数成倍增加。建议先在二维简单模型上验证，确认表面态存在后，再进入三维或者更复杂的结构。
5.4 后续扩展思路
这套 Comsol + Matlab 流程并不只限于二元平板型一维光子晶体。换一个角度思考，它还能用于：
- 含缺陷层的一维光子晶体，研究缺陷模和表面态耦合。
- 两个不同周期的一维光子晶体拼接，构造拓扑界面态和拓扑波导。
- 结构中加入非线性介质，在频率域扫描里追加 Kerr 项，研究非线性对 Zak 相位的影响。
- 把一维模型推广到二维六角晶格光子晶体，Wilson loop 从标量变成矩阵形式，这时 Matlab 后处理的灵活性优势会更加明显。
每一条扩展路径都需要在 Comsol 端调整模型、在 Matlab 端调整后处理脚本，但核心逻辑不会变：确保 Bloch 函数还原正确，确保能带追踪正确，Wilson loop 连乘之后查相位。只要这三根柱子立住，整个框架就是可移植的。
最后再说一个实际体会。整个流程里最花时间的不是 Comsol 建模，也不是 Wilson loop 公式，而是模式追踪和符号约定这类"看起来很小"的问题。它们不像网格加密那样可以靠算力硬顶，必须靠逻辑判断和经验去识别。所以建议第一次跑时，先不要急着追求完整结果，花十几分钟把单个 k 点的本征场导出来，用"周期函数"约束检查 Bloch 函数还原是否正确，再做批量扫描。这个前置检查能帮你省下大半天排错时间。之后把整套流程封装成函数，以后换结构、换材料，一行参数改动就能重新算一遍，效率会高很多。
