最近在整理手头一个跟光学超构表面相关的课题需要把一维光子晶体能带里的拓扑不变量——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 方向的波矢分量kxy 方向给 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); % 这里只演示 solnum1实际需要循环所有模式 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必须是周期函数所以在x0和xa处的值应当一致误差在数值容忍范围内。如果两端的幅值一样、相位差接近零说明符号选对了如果出现接近2k a的相位差那就把指数符号反过来再试。这一步我当时花了很长时间才想明白。原因是 Comsol 的文档对 Floquet 因子的写法写在了不起眼的边界条件说明里很容易略过。所以这个经验必须记下来不要盲信任何教程永远用周期函数这个物理约束去检验自己的符号约定。3.3 能带追踪模式顺序交错是最大干扰每次特征值求解返回的模式在solnum里的排列顺序本应该按照特征频率从小到大排。但问题在于有限元求解器返回的顺序在正常情况下是按频率排序的可是当你把不同 k 点的结果放到一起时某个能带在 k 点 1 是第 5 个求解模式到了 k 点 2 可能变成第 6 个。如果直接按模式序号去连接能带Wilson loop 会算出一堆乱七八糟的内积。解决办法是做重叠矩阵最大匹配。对相邻两个 k 点的本征模式集合{u_m(k_j)}和{u_n(k_{j1})}先计算所有模式对之间的重叠积分矩阵S_{mn} ⟨u_m(k_j) | u_n(k_{j1})⟩然后对每一行取绝对值最大的列索引作为下一 k 点应该对应的模式序号。这个过程等价于人眼追踪能带走向但用矩阵运算自动完成。当能带远离简并时非常稳定在带边附近接近简并时会误判这时候需要提高 k 点采样密度。具体到实现我会先按能量从低到高排序一次然后逐段做最大匹配。注意能带在布里渊区边界处可能存在级数交换比如第 1 带和第 2 带在kx π/a附近若发生反交叉追错就会导致后续相位全错。处理原则是宁可多增加 k 点也不要让相邻 k 点频率差太大。我一般把Nk设为 60 到 120具体看带边复杂程度。3.4 Wilson loopZak相位的离散化计算一旦模式排序搞定了Zak 相位就可以通过相邻 k 点之间的内积连乘得到这就是 Wilson loop 的离散形式W ∏_{j1}^{Nk-1} ⟨u(k_j) | u(k_{j1})⟩对单带情况W是一个复数Zak 相位取θ -Im(log(W)) -angle(W)为什么可以这样算因为 Berry 联络的积分可以写成相邻态重叠积分的连乘这是平行移动思想在数值上的直接实现。当采样点足够密相邻重叠积分趋于 1连乘的对数虚部就趋于连续积分的结果。一个关键点是最后首尾要不要连接。我们扫描的区间是从-π/a到π/a这是整个布里渊区但由于k -π/a和k π/a实际上是同一个物理点相差一个倒格矢所以严格来说应当把最后一个点和第一个点也做一个重叠积分构成一个闭合环。不过实际计算中两个端点对应的模式可能因为解算器给出的规范不同在布洛赫函数上差一个相位闭合积分可以吸收这个相位差取对数得到的是模2π的相位。这正是拓扑不变量需要的。因此我的代码会额外把jNk和j1连起来算不能漏。下面是一个简化版的 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 层界面处开始一组完整的 AB 周期所有结果都干净了。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 函数还原是否正确再做批量扫描。这个前置检查能帮你省下大半天排错时间。之后把整套流程封装成函数以后换结构、换材料一行参数改动就能重新算一遍效率会高很多。