企业官网建设流程全解析

计算一维光子晶体的Zak相位，最近不少做拓扑光子学的朋友都在问我同一个问题：Comsol里能算出漂亮的能带图，可怎么把它变成具有物理意义的拓扑不变量？本文就围绕这个项目，完整记录我从Comsol建模到Matlab分析Zak相位的全过程，包括原理、参数设置、代码实现和踩坑记录。这套流程适合正在入门拓扑光子学、或者需要复现一维光子晶体边界态研究的人参考，看完之后你可以直接迁移到自己的体系里。

1. 项目思路拆解：为什么选Comsol加Matlab这条路线

1.1 算Zak相位到底在算什么

一维光子晶体本质上就是折射率沿某个方向周期调制，形成布拉格散射和能带带隙。它和拓扑扯上关系，是因为在2014年前后，研究者发现两个带隙结构相同的一维光子晶体拼接在一起时，界面处是否出现带隙内的局域态，取决于两条能带在带隙边界的Zak相位差。如果带隙两侧的Zak相位差为π，界面处就一定会出现拓扑边界态；如果相位差为0，则没有。

Zak相位是Berry相位在一维周期性体系里的特殊形式，数学定义是布洛赫波函数中的周期部分在布里渊区内沿闭合路径积分得到的相位。它的核心价值在于：虽然能带本身只告诉我们频率和波矢的关系，但Zak相位额外携带了波函数在动量空间演化过程中的几何信息，这才是判断拓扑性质的关键。

对于一维光子晶体，Zak相位在结构具有中心对称性时表现为量子化的0或π，这也是为什么那么多工作都围绕着一维体系的拓扑边界态展开。把Zak相位算出来，你就能解释实验里观察到的反射相位异常、边界态存在性、以及透射谱中的共振峰来源。

1.2 Comsol和Matlab的分工逻辑

Comsol在光子晶体计算里的优势非常明显：几何随意、材料参数随便设、边界条件丰富，网格局部加密也方便。但它有一个死穴——它默认给你能带频率和模式场，却不直接给拓扑不变量。你无法在Comsol的界面里找到一个叫“Zak相位”的输出选项

Matlab则正好补上这一环。把Comsol算出的每个波矢处的本征模式场提取出来，在Matlab里做归一化、做相邻k点内积、累加相位，就能得到Wilson loop，进而推出Zak相位。这个过程本质上是一个数值线性代数问题，Matlab做这种事轻车熟路。

很多人可能会问，能不能直接用传输矩阵法或者平面波展开法，在Matlab里把能带和Zak相位一次性搞定？当然可以，那种方法编写出来之后运行也很快，但可扩展性差，换个复杂几何就得重写。Comsol建模的优势在于：后期你如果想要加入增益损耗、各向异性材料、非线性介质、或者把模型推广到二维三维，代码框架不用大改，只要调整Comsol里的物理场配置就行。这也是我把Comsol和Matlab结合起来作为一个完整项目来做的原因。

2. Zak相位的原理与数值计算要点

2.1 从Berry相位到Zak相位的完整定义

在固体物理里，一条能带对应的布洛赫函数可以写成psi(k,r) = e^{ikr} u(k,r)，其中u(k,r)是周期函数，满足u(k,r+a)=u(k,r)。Zak相位定义为：

θ_n = i ∮_BZ ⟨u_n(k)|∂_k u_n(k)⟩ dk

这里的积分区间是整个一维布里渊区，从-k0到k0再闭合回来，对于一维体系通常取[-π/a, π/a]。|u_n(k)⟩是第n条能带在波矢k处的布洛赫周期函数。

这个积分形式上和Berry相位完全一致，只是维度变成了一维。一维的特殊性在于，布里渊区是一个闭合环，波函数在k空间绕一圈之后回到原点，但由于周期规范的存在，会多出一个相位因子，这个累积相位正是Zak相位。注意Zak相位是规范依赖的——原点平移会导致Zak相位发生改变。实际操作中，只有当结构本身具有中心对称时，Zak相位才会被限制为0或π两个离散值，这也是做一维拓扑光子学的工作都选择中心对称结构的原因。

2.2 Wilson loop：数值计算Zak相位的标准做法

解析计算Zak相位需要知道u(k,r)的完整表达式，这对数值计算不友好。数值上最常用的方法是Wilson loop方法，原理是把积分离散化，变成相邻k点之间内积的连乘。

把布里渊区分成N个点，k_m = -π/a + m·Δk，其中m=0, 1, ..., N-1，Δk = 2π/(aN)。离散形式的Zak相位可以写为：

θ_n = -Im Σ_m log⟨u_n(k_m)|u_n(k_{m+1})⟩

这里|u_n(k_{m+1})⟩是从k_m出发的下一个波矢处的模式，整个连乘再取log和虚部，就能得到相位。为什么这个式子成立？因为当Δk很小时，⟨u_n(k_m)|u_n(k_{m+1})⟩约等于1 + Δk⟨u_n|∂_k u_n⟩，取对数后虚部正好对应Berry联络的累积。

实际计算中需要特别注意两点。第一，相邻k点之间的内积模应当非常接近1；如果内积模显著小于1，说明网格不够密，离散化误差太大，需要加密k点采样。第二，每个k点的模式场在导出时往往带有任意相位，这会导致内积结果不稳定，因此每次内积前需要先归一化，或者用一个固定的参考相位把模式场校正到同一规范下。

2.3 边界闭合条件：最容易出错的一步

前面提到布里渊区是一个闭合环，但实际采样时，你不可能同时把-π/a和π/a作为独立的采样点各算一遍，因为这是同一个物理点。Wilson loop的正确做法是只采样N个点，比如从k_0=-π/a到k_{N-1}=-π/a + (N-1)Δk，然后利用周期规范把最后一个点与第一个点连接起来。

周期规范的表达式是：

u_n(k+G) = e^{-iGr} u_n(k)

在一维情况下，G = 2π/a，所以当你处理边界项、计算⟨u(k_{N-1})|u(k_N)⟩时，k_N对应的模式场不是独立计算出来的，而是用起点处k_0的模式场乘以e^{-i2πx/a}得到的。这个边界闭合项如果不加，积分的路径就是开的，算出来的相位会差一个边界贡献，数值上可能既不是0也不是π，导致错误结论。

这是我实测下来整个项目里最容易出问题的地方。很多教程在代码里没有明确处理这一项，结果算出来的Zak相位看起来乱跳，但很难找出原因。建议你在写代码时单独输出每一步的相位增量和累计相位，检查是否存在从接近+π跳变到-π的情况，这种跳变往往说明k点采样过粗或边界项处理有误。

2.4 符号约定与结果检验方法

Zak相位的符号取决于两个地方：一是Berry联络的积分方向，二是Wilson loop公式里取log前的正负号。不同文献的约定可能不同，如果你发现自己的结果和文献上的相位符号相反，不要急着改代码，先统一约定再判断。

我建议在做完整计算之前，用一个结构简单、结果已知的模型来验证代码。比如用二分之一占空比的一维光子晶体，每个周期内两种介质各占一半宽度，通过改变原点位置观察Zak相位的变化，或者直接和传输矩阵法得到的反射相位对比。理论上，当结构从中心对称变为非中心对称时，Zak相位会从量子化的0或π变成非量子化的任意值。拿这个性质来验证代码是最直接的。

3. Comsol建模实操：参数、几何、边界条件与扫描设置

3.1 建模前想清楚你算的是哪种偏振

一维光子晶体的模式分为TE和TM偏振。在Comsol二维模型中，不同偏振对应不同的场分量。我的项目里选择电场沿面外方向，也就是电场只有z分量E_z，磁场在面内。对于这种偏振，波动方程变成标量形式，处理起来最简洁。

如果你算的是面内偏振（即电场在x-y平面内），那需要同时求解E_x和E_y两个分量，方程组更大，Floquet边界条件的设置也更繁琐。对于一维周期性结构来说，算标量形式足够说明问题。实际课题研究中如果涉及斜入射或者偏振变换，再去考虑面内偏振不迟。

3.2 几何建模与材料参数设置

我用的是二维模型，几何是在x方向取一个晶格常数a的长度，y方向取一个很小的矩形高度h。比如a = 1微米，h = 0.05微米。两个介质层并排放在一个周期内，材料采用无损介质，折射率分别为n1和n2。

这里给出我实际用的一组参数：

之所以选n1=3.45和n2=1.0，是因为折射率对比度足够大，带隙会比较宽，后续在带隙边界处讨论Zak相位时特征更明显。如果你用二氧化硅和空气组合，折射率对比只有1.45左右，带隙窄，计算时需要更密的能带采样才能分辨带隙边界，对数值精度要求高。

材料参数中要注意：Comsol中相对介电常数是ε_r = n^2，如果材料有损耗，还需要设置电导率或复介电常数。算Zak相位时建议先关掉损耗，因为损耗会模糊能带边界，也会让本征模式场的相位提取变得困难。

3.3 Floquet周期边界条件与波矢扫描

这是Comsol建模最核心的一步。选择电磁波、频域接口，二维模型。左边界和右边界设置成Floquet周期边界条件，K矢量设为(kx, 0)。上下边界设置成周期性边界条件，这样y方向无限延伸，等价于一维体系。

在全局参数里定义一个kx变量，然后在研究设置中的辅助扫描里让它变化，扫描范围取[-π/a, π/a]，步长按需要的k点数量确定。我通常先扫41个点做验证，加密到81个点得到最终数据。

需要特别注意的是，Comsol的特征频率求解器中，k矢量通过Floquet周期边界条件输入，而kx是作为参数扫描的。这种情况下，每个kx值都会触发一次特征值求解，得到的每个特征频率都对应一条能带上的点。Comsol会把所有特征值结果按求解器默认方式排列，但这不保证是按能带顺序排列的，后期需要用频率排序或重叠积分来整理数据。

3.4 特征频率研究设置与网格无关性测试

研究类型选择特征频率，要设定想要的模态数。我的项目中扫描前10条能带，所以特征频率搜索基数设为10或12，留一些余量。物理场中，电磁波频域接口的特征值方程本质上是关于频率ω的特征值问题，解出来的特征频率经过后处理可以直接画能带。

网格建议先用较粗的网格跑通流程，再逐步细化。因为你的主要目标不光是能带频率，还要提取模式场，网格密度直接影响Wilson loop计算内积时场矢量的离散表示。我在这个项目中测试过最大网格尺寸从0.1a细化到0.02a，前几条能带的频率变化很小，但Zak相位的计算结果在粗网格下会出现明显偏差。最终我采用最大网格尺寸0.03a，既保证频率精度，又不会让导出文件太大。

一定要关闭自适应网格细化。特征频率求解器中如果有自适应网格，每个kx扫描点的网格会变化，导致相邻k点模式场所在的节点位置不一致，后续Matlab算内积时会出现严重问题。这个坑我踩过，如果你在数据后处理时发现内积模远小于1，先检查是不是网格随参数扫描变化导致的。

3.5 导出模式场数据的具体操作

计算完成后，我们要导出每个kx处的本征模式场。在Comsol的导出数据选项中，选择数据集为特征频率解的某一索引，表达式填E_z的实部和虚部，再额外导出x、y坐标。

这里的关键问题是：Comsol导出的是物理场E_z，也就是布洛赫函数ψ_k(r)本身，而不是周期函数u_k(r)。由于ψ_k(r) = e^{ikx}u_k(r)，因此导出后需要在Matlab里乘以e^{-ikx}来得到u_k。另一种做法是在Comsol的派生值中直接定义表达式E_zexp(-ikx*x)，导出的就是周期函数部分。推荐后者，省去在Matlab里处理的麻烦。

每个kx点导出一个单独的文件，文件名包含kx索引，例如data_0001.txt，data_0002.txt。导出时选择所有特征频率索引，这样每个文件都包含多条带的数据。实际文件里会包含大量行，每条带在一个时间步索引中给出，需要按特征值索引分拆。建议导出前在Comsol中先按特征频率排序，或者用脚本生成一个包含特征值顺序信息的表，方便Matlab读取时对应。

4. Matlab读取与预处理：从导出文件到可用的模式向量

4.1 文件目录规划与批量读取

Comsol导出数据的组织方式直接影响后续代码的复杂度。我的做法是建一个data目录，里面存放所有kx点的导出文件，命名格式是kx_01.txt这样。每个文件的列顺序设为：x坐标、y坐标、Re(E_z)、Im(E_z)、Re(E_z_bloch)、Im(E_z_bloch)，其中E_z_bloch就是前面说的剥离布洛赫指数后的场分量。

Matlab里读取很简单，不依赖任何额外工具箱：

function [x, y, Ez, u] = loadFieldData(fileName) data = readmatrix(fileName); x = data(:, 1); y = data(:, 2); Ez = data(:, 3) + 1i*data(:, 4); u = data(:, 5) + 1i*data(:, 6); end

读取后需要检查每个文件的节点数是否一致。由于我关闭了自适应网格，整个扫描过程的网格没有变化，所以所有文件的节点数应该完全相同。如果节点数不一致，后面做内积时直接按元素相乘就会出错。万一因为某些原因网格变了，解决办法是先在Matlab里用griddata插值到统一网格上，但这个方法会引入误差，能不用就不用。

4.2 模式场的归一化处理

Wilson loop公式中使用的态矢量应当是归一化的。Comsol导出的模式场通常有自己的归一化方式，但不同kx点之间的模长尺度可能不同，因此必须在Matlab里对每个模式重新归一化。

对于一个复场向量u，归一化的代码是：

u = u / norm(u);

norm函数默认计算L2范数，即sqrt(sum(abs(u).^2))。这个归一化操作对后续内积计算至关重要，因为如果不归一化，每个内积都会带入一个未知的比例因子，最后取log时这些比例因子不会抵消干净，会污染相位累积。

还有一个细节：如果模式场在空间不同位置上的幅度分布差异极大，比如局域在某一层中，那么简单的L2归一化可能对远场区域的小幅度噪声过于敏感。实际中我发现，网格质量足够好时这个问题不会太明显；如果确实出现异常，可以只取结构内部区域的节点做内积，但不能完全丢弃周期边界附近的点。

4.3 节点一致性与内积模检验

在开始计算Zak相位之前，我强烈建议先做一步诊断：计算相邻k点模式间的内积模|⟨u(k_m)|u(k_{m+1})⟩|。这个量在理想情况下应当非常接近1，偏差通常小于千分之一。如果发现某个k点对之间的内积模为0.5甚至更小，说明模式排序错乱或网格有问题。这时直接算Zak相位肯定出错，先把诊断问题解决再继续。

内积模诊断的Matlab代码：

for m = 1:Nk-1 inner = dot(u{m}, u{m+1}); fprintf('k index %d to %d, |<u|u>| = %.6f\n', m, m+1, abs(inner)); end

当你遇到内积模明显偏离1的情况，多数原因是模式排序错乱，就是第m个k点的第n条带和第m+1个k点的第n条带实际上不是物理上的同一条能带。这种问题在带隙较窄、能带交叉区域经常出现，需要回到Comsol里查看能带回线，或者改用频率排序之外的场重叠积分方法重新匹配能带顺序。

5. Matlab计算Zak相位的核心代码与结果分析

5.1 单条能带的Zak相位计算主程序

假设你已经把所有kx点的模式场存在一个cell数组u_list中，数组长度是Nk，每个元素是一个Nnode×1的复向量。下面这段代码实现了完整的Wilson loop计算：

clear; clc; % 参数定义 a = 1e-6; % 晶格常数 Nk = 41; % 布里渊区采样点数 G = 2*pi/a; % 倒格子基矢 kx = linspace(-pi/a, pi/a, Nk+1); kx = kx(1:end-1); % 去掉最后一个重复点，闭环由周期规范处理 % 这里需要先加载所有模式场数据到 u_list % u_list{m} 是第m个k点处某一条能带的模式场向量 ZakSum = 0; % 相邻点内积累积 for m = 1:Nk-1 u1 = u_list{m}; u2 = u_list{m+1}; u1 = u1 / norm(u1); u2 = u2 / norm(u2); inner = dot(u1, u2); if abs(inner) < 1e-8 error('内积模接近0，模式匹配失败'); end ZakSum = ZakSum + log(inner / abs(inner)); end % 边界闭合项：最后一个点与第一个点的连接 u_first = u_list{1}; % k = -pi/a u_last = u_list{Nk}; % k = -pi/a + (Nk-1)*Δk u_first = u_first / norm(u_first); u_last = u_last / norm(u_last); % 这里假设导出时已经是剥离了布洛赫指数的u场， % 所以闭合项需要乘上e^{-iGx} x_coord = x_list{1}; % x坐标，从导出文件中读取 phaseFactor = exp(-1i * G * x_coord); inner_boundary = dot(u_last, u_first .* phaseFactor); ZakSum = ZakSum + log(inner_boundary / abs(inner_boundary)); % Zak相位（弧度） Zak = -imag(ZakSum); % 归一化到 [0, pi) Zak = mod(Zak, pi); fprintf('Zak phase = %.4f rad\n', Zak);

这段代码里我特意用log形式而不是直接连乘后取arg，原因是log形式可以逐步观察相位累积过程，更容易定位问题出在哪两个k点之间。

5.2 多条能带同时计算的循环结构

如果要计算前N条能带，需要在外层加一个能带索引循环。每个能带使用各自对应的模式场数据。由于Comsol导出的特征频率结果可能没有按能带排列，Matlab侧需要先做一个能带分拣。

一种简单有效的分拣方法是：读取每个kx点的全部特征频率，按频率大小排序后，依次对应能带1、能带2等等。但在能带交叉区域，这种方法会出错。更可靠的方法是使用模式匹配：从k_0出发，用第n条带在k_m的模式与k_{m+1}的所有模式计算重叠积分，选择重叠积分模最大的模式作为同一条带的延续。这个思路实现起来也不复杂：

for band = 1:Nbands % 存储该条带所有k点模式 u_band = cell(Nk, 1); for m = 1:Nk if m == 1 % 第一个k点，按频率排序取第band条 u_band{1} = mode_data{1}{band}; else % 后续k点，找与上一k点同带模式重叠最大的模式 maxOverlap = -1; bestIndex = 1; for n = 1:Nmodes overlap = abs(dot(u_band{m-1}, mode_data{m}{n})); if overlap > maxOverlap maxOverlap = overlap; bestIndex = n; end end u_band{m} = mode_data{m}{bestIndex}; end end % 对u_band运行Wilson loop计算 zakPhase = computeZakPhase(u_band, x_coord, G); end

这种基于模式重叠的能带追踪方法在带隙较宽、模式区分明显的体系中非常稳健。如果体系出现近简并，两组模式的重叠积分都很接近1，就需要进入简并子空间做Wilson loop，这是一个稍微复杂的升级版，本文不展开。

5.3 结果判定：怎么看Zak相位算对了

算出来的Zak相位应当在0或π附近，因为中心对称结构的Zak相位是量子化的。我这里给出一个实际计算例子。用之前说的参数结构，晶格常数1μm，高低折射率层厚度分别为0.7a和0.3a，折射率3.45和1.0，取41个k点。前四条能带的计算结果如下：

与文献对照，这个二分之一结构在折射率对比足够大的情况下，Zak相位按带隙排列为0, π, π, 0，符合预期。如果换成非中心对称结构，例如改变第二层介质的位置使原点的选取不再对称，Zak相位会明显偏离0和π，这也是一个很灵的敏感性测试。

5.4 采样点数对结果的影响

k点采样数量对Zak相位的影响比很多人想象的更大。我用同一模型测试了Nk从11到121变化时的结果。Nk=11时，计算误差比较大，有些能带的Zak相位偏差达到0.1 rad级别；Nk=41时，基本稳定在0或π附近，偏差小于0.01 rad；再加密到121时，偏差进一步减小到0.001 rad量级。

原因是Wilson loop中相邻k点之间的相位增量|Δθ|不能超过π，否则log函数的虚部会混叠。k点越密，相邻内积相位差越小，累积越准确。建议至少取41个k点作为起步，若要追求高精度或处理近简并能带，建议取81~121个点。求解时间和数据量会相应增加，但对单个一维模型而言完全在可接受范围内。

6. 常见问题与排查技巧实录

6.1 能带模式排序错乱导致内积跳变

这是我自己做这个项目时遇到最多的一个问题。现象是：诊断内积模时，大部分k点对的内积模在0.999以上，但某个点特别低，只有0.3甚至0.1。原因几乎都是Comsol特征频率求解器在每个kx扫描点返回的模式顺序不是按能带排列的，频率接近的两条带交叉时，求解器返回顺序会交换。

解决办法就是在Matlab里用模式匹配法做能带追踪，而不是直接按频率排序取第n条模式。另外还有一个技巧：在Comsol中设置辅助扫描时，勾选“使用前一个解作为初始猜测”，可以让相邻kx点的模式顺序更加一致，减少后期处理工作量。

6.2 边界闭合项缺失导致相位不量子化

如果Zak相位计算结果不是0或π，但它们分布在某个中间值附近，比如0.3π，首先不要怀疑物理模型，先检查边界项。边界项是整个Wilson loop里唯一体现布里渊区周期性的地方，漏掉它，算出来的是一个开路径积分，物理上没有任何意义。我在代码里特意把边界闭合项单独列出来，就是希望读者看清楚这一项的作用。

还有一个容易搞混的地方：周期规范算符到底是乘以e^{-iGx}还是e^{+iGx}。这取决于你对正k方向的定义和Comsol中Floquet周期条件的k矢量符号设置。我建议你在小范围内用解析模型验证一次符号。最简单的验证对象是自由空间（均匀介质），此时Zak相位应当为0，如果符号写反，会得到一个与路径长度相关的非零值。

6.3 内积模总是略小于1：网格和节点不一致

如果诊断显示所有内积模都在0.95左右，虽然能算出Zak相位但总觉得不放心，通常是两个原因。一是网格不够细，离散化误差大；二是上下边界或内部界面处的场在相邻k点之间发生了微小变化，而这种变化是因为Floquet边界条件的数值实现不是完全一致的。

解决方法：把网格加密一个量级再跑一次，如果内积模明显上升，就是网格问题。如果加密后内积模还是0.96，就要检查是不是上下边界条件选得不对。我在之前的模型中上下边界用周期性边界条件，换来的是y方向的均匀场，内积模轻松到0.999以上。换成PEC边界会引入y方向的横向模式变化，内积模就会下降。

6.4 Zak相位符号和文献相反

出现这种情况，大概率不是计算错误，而是约定不同。有些文献定义Zak相位时积分方向是从0到G，有些是从-G/2到G/2，符号差一个负号。还有的文献在Wilson loop里取log之后用+imag而不是-imag，结果也是反号。我在文章里给出的公式和代码采用的是最常见的约定：θ = -Im Σ log(...)，方向沿k正方向。如果你的情况特殊，统一改掉符号即可，物理结论不受影响。

6.5 特征频率中出现不想要的杂散模式

算特征频率时，设置的搜索基数越大，越容易混入一些和物理问题无关的模式。比如上下边界条件造成的横向高次模、求解器数值产生的非物理模式。这些模式在能带图上表现为一些偏离主能带的光滑曲线之外的零散点。

判断是否为杂散模式有一个简单方法：在Comsol后处理中查看该模式场的空间分布。物理模式应该主要集中在介质结构内，场分布沿x方向周期变化均匀；杂散模式往往在场分布上有明显的横向振荡，或者能量集中在边界上。处理时把这些模式从导出列表中排除，或者在Matlab里根据模式场的傅里叶成分做一个过滤。

6.6 常见问题速查表

7. 实际操作中最后想说的话

这套流程跑通以后，我的最大感受是：算Zak相位比算能带要敏感得多。能带数据差一点还能看出趋势，Zak相位有一处弄错就直接不量子化了，反而逼着我把整个数值链路从头到尾查了一遍。这是好事，因为通过这个排查过程，我对Comsol场导出、Floquet边界条件的数值实现、以及Wilson loop的离散形式都有了更扎实的理解。

如果你刚开始做类似项目，我建议先不要急着上复杂结构，就用最简单的两个介质层交替的一维光子晶体，把整个流程跑通，再用一个已知结果来验证你的代码，确认无误后再推广到多层结构、渐变结构、或者带损耗的体系。另外，Comsol的LiveLink for Matlab可以省去文件导入手动操作的麻烦，算是锦上添花的工具；没有的话，纯文本导出加readmatrix也完全够用。

最后分享一个调试技巧：在Matlab的循环里每隔几个k点打印一次当前累计相位和步进相位，画出来看看。如果累计相位在某个位置出现了接近π的跳变，那就是内积相位越过分支切线了，需要增加k点采样密度或者调整数据处理方式。这种逐步观察的习惯，比最后只看一个最终数字有用得多。