简介：这份资料围绕光子晶体能带计算与陈数分析，面向使用Comsol Multiphysics和MATLAB开展光子晶体研究的初学者与科研人员，重点解决从三维建模、周期边界设置到能带结构提取的全流程问题。压缩包共5个文件，包含2份PDF文档（原理说明与文件说明）、1个mph仿真模型、1个m文件（陈数计算脚本）及1个txt数据文件，整体大小约20MB，文件搭配完整，可边看边练。该资源已有1043人学习下载，关注度较高。通过对照PDF中的步骤与脚本，读者可掌握光子晶体能带计算、布里渊区分析、能带图绘制以及陈数提取的实现方法，并能借助LiveLink接口完成Comsol与MATLAB的数据交互，为后续设计光子晶体光纤、激光器等器件提供理论与算法支持。
1. 陈数计算：为什么你的 COMSOL 能带图还差一个整数
做光子晶体、声子晶体或电磁超材料仿真的人，大概率都有过这种经历：COMSOL 里跑完能带，带隙位置对得上文献，却说不上这条带“有没有拓扑数”。陈数（Chern number）就是答案，它给每条能带贴一个整数标签，决定边界态的数量和存在性。“陈数的数怎么写”落到 COMSOL 里，其实是一个工程问题：从 Floquet 周期条件解出的复本征场出发，在布里渊区上离散积分，再把结果稳定地取整到 0、±1、±2。这篇内容按我从几何建模到拿到陈数、再用边界态反推验证的完整路径来写，适合正在做拓扑光子学、声学拓扑和电磁超材料的人；新手可以照着复现一遍，熟手可以直接跳到第 5 章对照踩坑清单。
2. 陈数的物理框架：从 Berry 曲率到那个整数
2.1 Berry 相位与 Chern number 的几何含义
陈数不是从能带图上“数”出来的，它是在动量空间积分出来的。考虑一个二维周期结构，布洛赫态写成 \(\psi_{n,k}(r) = e^{ik \cdot r} u_{n,k}(r)\)，其中 \(u_{n,k}\) 是满足周期边界条件的格点周期函数。随着 \(k\) 在布里渊区里走一圈，\(u_{n,k}\) 会累积一个几何相位，这就是 Berry 相位。
Berry 相位可以写成沿闭合路径的积分：
\[ \gamma_n = \oint_{\mathcal{C}} i \langle u_{n,k} | \nabla_k u_{n,k} \rangle \cdot dk \]
对二维布里渊区，定义 Berry 曲率 \(\Omega_n(k_x, k_y)\) 为 Berry 联络的旋度。单条带的陈数就是 Berry 曲率在整个布里渊区上的积分：
\[ C_n = \frac{1}{2\pi} \int_{\mathrm{BZ}} \Omega_n(k_x, k_y) , dk_x dk_y \]
这个积分结果必须是整数，物理上它描述布洛赫态在动量空间“缠绕”的次数。如果算出来是 1.997 或者 -1.003，不要怀疑物理，先怀疑自己的积分域、波矢单位或模式追踪。
很多人第一次听到“陈数”这个名，误以为它像特征值一样是某个方程的解。实际上它更像绕数：给你一条闭合曲线绕某点转了几圈，你数出来几就是几。只是这里“绕”的对象不是曲线，而是量子态在布里渊区上的相位变化。
2.2 三种数值计算路径，COMSOL 该选哪种
算陈数有三大类主流数值路径，网上教程多半只讲其中一种，但实际工程里要根据你的 COMSOL 模型结构选。
第一种是“格点 Berry 曲率积分”，也就是 Fukui-Hatsugai-Suzuki（FHS）方法。做法是把布里渊区切成 \(N \times N\) 的网格，在每个小方格上计算相邻格点布洛赫态的重叠矩阵，再把每个方格的“曲率相位”累加。这个方法的优点是实现简单，只需要本征场本身，不需要对场求导，COMSOL 导出的复场数据直接能用。缺点是网格必须覆盖整个布里渊区，而且每条带单独算，带交叉附近容易出错。
第二种是 Wilson loop 方法。沿 \(k_x\) 方向做一条直线路径，固定若干 \(k_y\) 值，把路径上每一小步的连接矩阵乘起来，得到一个幺正矩阵，再对这个矩阵求对数特征值。对数特征值的斜率就是陈数。它天然适合多带一起算，能处理带交叉，是现代拓扑计算的主流做法。
第三种是实空间方法，把周期结构换成有限大超胞，直接对角化实空间 Hamiltonian，再用实空间纠缠谱或位置算符的本征值分布来提取陈数。好处是对无序和边界开放系统更稳，坏处是 COMSOL 里建大超胞的计算量迅速爆炸，而且超胞的边界处理很容易掩盖拓扑信号。
对 COMSOL 用户，我一般建议首选 FHS：研究类型选“本征频率”，扫参扫出每个 \(k\) 点的本征场，导出后在 MATLAB 或 Python 里做积分。跑通一条带之后，再用 Wilson loop 做交叉验证。实空间方法只在涉及无序或缺陷时才考虑，且通常需要换工具。
2.3 格点上的陈数离散公式与采样密度
FHS 方法的核心是把 Berry 曲率的积分换成格点上的规范场累积。在二维 \(k\) 网格上定义两个连接矩阵：
\[ U_x(k_x, k_y) = \frac{\langle u(k_x + \Delta k, k_y) | u(k_x, k_y) \rangle}{|\langle u(k_x + \Delta k, k_y) | u(k_x, k_y) \rangle|} \]
\[ U_y(k_x, k_y) = \frac{\langle u(k_x, k_y + \Delta k) | u(k_x, k_y) \rangle}{|\langle u(k_x, k_y + \Delta k) | u(k_x, k_y) \rangle|} \]
每个小方格的格点曲率定义为：
\[ F_{xy}(k) = \ln \left[ U_x(k) , U_y(k + \Delta k_x) , U_x(k + \Delta k_y)^{-1} , U_y(k)^{-1} \right] \]
将所有方格的 \(F_{xy}\) 求和，再除以 \(2\pi i\)，取实部，就得到该带的陈数。这里log取复对数主值，所以每个方格的相位必须落在 \((-\pi, \pi]\) 内。如果网格太粗，相位跳变超过 \(2\pi\)，复对数的割线就会造成回绕数丢失。
网格密度怎么选？我做过测试：对简单的三角晶格光子晶体，\(N=20\) 时陈数已经能稳定在 ±1；对带隙窄、色散平的系统，\(N=60\sim80\) 才收敛。COMSOL 的代价是每次扫参都要解一遍特征值问题，\(N=60\) 意味着 3600 次特征值求解，即使每次只有 1 秒也要一个小时。所以不要一上来就开大网格，先用 \(N=12\) 跑通全流程，再加密验证。
提示：COMSOL 扫参是按顺序逐个求解的，参数个数直接影响总耗时。先想清楚布里渊区路径要不要全部覆盖，再决定网格规模。
3. 在 COMSOL 里跑出能带与本征场：模型建立与导出
3.1 物理场接口、研究类型与单胞几何怎么选
COMSOL 里做陈数计算，物理场接口的选择没有唯一正确答案，但选错会直接导致后续导出字段不可用。我做光子晶体时常用“电磁波，频域”接口；做声学拓扑时用“压力声学，频域”；如果研究的是抽象模型或自定义偏微分方程，就用“弱形偏微分方程（Weak Form PDE）”接口。
关键不是选哪个物理场，而是确认两点：第一，接口支持复值本征场；第二，接口内置“周期性条件”里能选“Floquet 周期性”。以电磁波频域接口为例，因变量是 \(E\) 场的三个分量，本征求解后每个分量都是复数，完全够用来算内积。弱形式 PDE 接口则需要自己定义因变量的实部虚部分量组合成复场，比如把因变量设为u1和u2，构造 \(u = u1 + i u2\)，重叠积分的表达式要写u1*u1c + u2*u2c再加交叉项。
几何上只建一个单胞，注意单胞的边界要严格周期对应。对六边形蜂窝晶格，COMSOL 里直接画六边形单胞再裁剪即可；对矩形单胞，把相对边标记为源和目标即可。网格用“周期性网格”功能，保证源边界和目标边界的网格节点一一对应。这一点非常关键，如果两侧节点不对应，后续导出场做内积时每个 \(k\) 点的矩阵维度都会不一致，根本没法积。
下表是我常用的模型参数：
| 项目 | 常见取值 | 说明 |
|---|---|---|
| 物理场接口 | 电磁波/压力声学/弱形式 PDE | 必须支持复本征场 |
| 研究类型 | 特征频率（Eigenfrequency） | 特征值虚部对应衰减，取实部画能带 |
| 单胞尺寸 | \(a = 1\) 归一化或按物理单位 | 会影响导出波矢的换算 |
| 网格 | 周期性网格，最大单元尺寸 \(\lambda_{\min}/8\) | 高带必须加密，否则特征值漂移 |
| 特征值搜索区域 | 从 0 开始到目标频率之上 | 多设几条带，给模式追踪留余量 |
3.2 Floquet 周期条件的相位约定与参数扫描设置
Floquet 周期条件是整个计算里最容易翻车的部分。COMSOL 的“周期性条件”节点下选“Floquet 周期性”，然后指定波矢分量变量名，常见写法是flq1.kx和flq1.ky。如果你用的是弱形式 PDE 接口，周期性条件里也是同一套设置。
默认的相位约定是：目标边界上的场等于源边界上的场乘以一个复指数，指数中的符号在 COMSOL 帮助文档里写得很清楚，不同版本符号可能相反。我第一次算蜂窝晶格光子晶体时，能带走势和文献完全相反，就是因为在目标边界上少乘了一个负号。
我的建议是：建立模型后先不急着扫整个布里渊区，只扫一个高对称路径，比如 Γ-M-K-Γ，画出能带和文献对比。若能带位置、简并性都对得上，再开二维扫参。别跳过这一步，等到陈数算完才发现符号反了，所有网格都要重跑。
参数扫描有两种方式：辅助扫描和普通扫描。计算陈数必须用辅助扫描（Auxiliary sweep），因为普通扫描对每个参数组合重新初始化求解器，可能会打乱特征值顺序；辅助扫描则保留前一参数点的解作为初始猜测，模式连续性更好。
波矢范围要换算准确。COMSOL 里 Floquet 波矢可以填归一化值 0 到 1，也可以填物理波矢 \(-\pi/a\) 到 \(\pi/a\)。我习惯填归一化值，导出后在外部分脚本里乘上 \(2\pi/a\)，这样网格索引不容易出错。但要时刻记住：积分域是真实的布里渊区面积，如果单位搞错，陈数会被乘上一个 \((2\pi/a)^2\) 之类的因子，最后取整会得到一个荒谬的整数。
3.3 把本征场导出成可计算的矩阵
COMSOL 本身不直接输出陈数，所以必须把本征场导出到 MATLAB 或 Python 里。三种常用做法：
第一种是 LiveLink for MATLAB，直接通过 COMSOL API 把model.result().dataset()里的场数据拉回 MATLAB 工作区。好处是不需要重复读文件，坏处是需要 LiveLink 许可证，且 MATLAB 版本和 COMSOL 版本要匹配。
第二种是“导出-数据”节点，把每个 \(k\) 点的本征场写成 CSV 或文本文件。在导出设置里选“数据”作为数据集源，表达式填实部和虚部分量，文件名用参数表达式区分 \(k_x\)、\(k_y\) 和带序号。这个过程很机械，但文件数量是 \(N^2 \times \mathrm{bands}\)，管理起来比较痛苦。
第三种是网格插值法，利用 COMSOL 的“网格数据”数据集，把网格节点编号固定，再在每个 \(k\) 点导出节点上的场值。我只推荐这种做法给有大量内存和处理批量文件经验的用户。
无论哪种方式，导出的场必须包含网格节点坐标，因为后续内积是逐点加权的。更稳妥的做法是直接在导出表达式里让 COMSOL 先算好某个目标分量（比如 \(E_z\)），再在外部脚本里把同节点的实部和虚部合成复数。节点顺序一致性是底线，否则重叠矩阵完全错乱。
4. 把 COMSOL 结果变成陈数：MATLAB 脚本与参数解析
4.1 数据读入与布里渊区网格表示
假设你从 COMSOL 导出了psi_kx_ky_mode.csv文件，每行包含x, y, real(Ez), imag(Ez)，文件名里的kx、ky是 1 到 \(N\) 的网格索引。先把数据读入 MATLAB，整理成四维数组：psi(nx, ny, mode, np)，其中np是单胞内的网格节点数。
% 参数设置 N = 40; % kx, ky 格点数 Nmodes = 6; % 每条能带取 6 个特征模，后面按频率排序 a = 1e-6; % 单胞尺寸，单位 m（记得和 COMSOL 几何一致） k = linspace(-pi/a, pi/a, N); % 物理波矢，COMSOL 里归一化值要乘 2*pi/a % 读取一个 k 点、一条模式的数据文件（示例文件名格式） % 数据结构：每行 = x, y, real(Ez), imag(Ez) for ix = 1:N for iy = 1:N for m = 1:Nmodes fname = sprintf('psi_%02d_%02d_%02d.csv', ix, iy, m); data = readmatrix(fname); x(:, ix, iy, m) = data(:, 1); y(:, ix, iy, m) = data(:, 2); psi(:, ix, iy, m) = data(:, 3) + 1i * data(:, 4); end end end这个读入过程的关键是固定x和y的节点顺序。COMSOL 导出时默认按网格节点编号排序，只要单胞网格不重新划分，每个文件的行序是一致的。数据读入后要做一遍归一化：对每个 \(k\) 点和每条模式，把 \(\psi\) 除以自身模长，消除求解器对整体幅度的任意缩放。
k向量的构造要特别注意：如果 COMSOL 里 Floquet 波矢用的是归一化值 0 到 1，那么实际物理波矢必须乘 \(2\pi/a\)。很多 FHS 脚本算出的陈数是准确的 2 或 3，就是因为这里少乘或重复乘了因子。
4.2 Fukui-Hatsugai-Suzuki 积分主循环
FHS 积分是全文的核心，代码反而很短，难的是把索引方向搞对。下面这段是我常用的单带版本，适用于目标带在频段上与其他带相隔较开的情形：
% FHS 陈数计算：遍历 k 网格上的每个小方格 C = zeros(1, Nmodes); for n = 1:Nmodes % 对每一条带单独计算 total_phase = 0; for ix = 1:N ix_p = mod(ix, N) + 1; % 沿 kx 的周期邻居 for iy = 1:N iy_p = mod(iy, N) + 1; % 沿 ky 的周期邻居 % 四个连接矩阵 Ux1 = overlap(psi(:, ix, iy, n), psi(:, ix_p, iy, n)); Uy2 = overlap(psi(:, ix_p, iy, n), psi(:, ix_p, iy_p, n)); Ux3 = overlap(psi(:, ix_p, iy_p, n), psi(:, ix, iy_p, n)); Uy4 = overlap(psi(:, ix, iy_p, n), psi(:, ix, iy, n)); % 格点曲率 = 复对数的圈乘 F = log(Ux1 * Uy2 * Ux3 * Uy4); total_phase = total_phase + F; end end % 陈数 = 总相位 / (2*pi*i)，取最接近的整数 C(n) = round(real(total_phase / (2 * pi * 1i))); end function ol = overlap(a, b) ol = dot(a, b) / (norm(a) * norm(b)); end代码逻辑讲解：Ux1是当前 \(k\) 点沿 \(k_x\) 方向平移一格后与当前态的重叠，Uy2是从平移后的点沿 \(k_y\) 方向再平移一格的重叠，Ux3是沿 \(k_y\) 方向平移后的点沿 \(k_x\) 负方向返回的重叠，Uy4是沿 \(k_x\) 方向平移后的点沿 \(k_y\) 负方向返回的重叠。四个连接矩阵乘起来，相位围绕了一个闭合的小回路，这个回路对应的几何相位就是该格点的 Berry 曲率贡献。
这里的overlap函数做了归一化，作用是把每个 \(k\) 点任意施加的整体相位消掉。FHS 格式的一个明显好处是“规范不变”：只要重叠矩阵的模长归一化，整体相位偏移会在闭合圈乘中自动对消。这也是你不必逐个 \(k\) 点做相位对齐、脚本也能跑出稳定整数的原因。
但有一个前提：相邻 \(k\) 格点的态必须是同一条带。如果你发现某个带在布里渊区某处和其他带发生交叉或简并，那么模式排序一旦换位，这个闭合回路的相位就断了。此时用单带 FHS 就是白费功夫，必须切换到 Wilson loop。
4.3 Wilson loop 校验与“数是什么”的最终判定
FHS 出结果后，我从不直接信，一定会用 Wilson loop 再验一遍。Wilson loop 的做法是在固定 \(k_y\) 的行上，把 \(k_x\) 方向所有相邻格点的连接矩阵连乘起来，得到一个幺正矩阵，对角化后取相位特征值。每行 \(k_y\) 对应一组相位，这些相位关于 \(k_y\) 的缠绕次数就是陈数。
% Wilson loop：对固定 ky 行做 kx 方向累积 % 这里假设只考察单条带（带间有 gap 的情况） WilsonPhase = zeros(N, 1); for iy = 1:N Mprod = 1; % 累积极矩阵 for ix = 1:N ix_p = mod(ix, N) + 1; M = overlap(psi(:, ix, iy, n), psi(:, ix_p, iy, n)); Mprod = Mprod * M; end WilsonPhase(iy) = log(Mprod); % Mprod 的模接近 1，相位即几何相位 end % 陈数 = 相位绕数，也就是对 iy 方向差分后求和除以 2*pi C_wilson = sum(diff(wrapToPi(cumsum([real(WilsonPhase); real(WilsonPhase(1))])))) / (2 * pi);严格来说，Wilson loop 在多带时应该是矩阵连乘，而不是标量连乘。上面的标量版本只适合单带且与邻近带完全分离的情况，多带系统的正确做法是把连接矩阵换成 \(N_{\mathrm{occ}} \times N_{\mathrm{occ}}\) 的重叠矩阵块。初次验证可以做简化版，但论文级别的计算不要用这个版本。
Wilson loop 的好处是把“积分整数”转化成“相位曲线的绕数”，肉眼也能判断。你画出的 Wilson 相位谱如果是一条从 \(-\pi\) 到 \(\pi\) 斜穿过的直线，斜率就是陈数；如果是一条水平线，陈数就是 0。FHS 和 Wilson 的结果互相对得上，再去填访谈结果；对不上，首选怀疑模式排序。
5. 陈数计算常见踩坑与排查：从能带翻车到取整翻车
5.1 简并点上的 Berry 曲率爆炸与 NaN
现象：FHS 积分结果不是整数，而是 NaN 或一个随网格密度漂移的浮点数。检查单个格点的F值，发现某些 \(k\) 点上的曲率相位异常大。
原因：这些点通常对应能带的简并点，也就是两条带的色散在某个 \(k\) 上重合。简并处本征态不唯一，任何微小的数值误差都会让重叠矩阵的方向突变，曲率在数学上发散，数值上表现为相位折叠。
解决：先画能带图，找出简并点在 \(k\) 空间的位置，然后让网格偏移半个步长，即 \(k_x = (i - 0.5) \Delta k\)，避开简并点所在的格点。另一个办法是改用 Wilson loop，因为 Wilson loop 是在一条线上累积连接矩阵，单个简并点不会直接毁掉一整行。
5.2 本征态随机相位导致积分漂移
现象：同一个模型，只改了 COMSOL 求解器的容差，陈数就从 1 变成了 0。
原因：COMSOL 特征值求解器在每一个 \(k\) 点求解时，会为本征向量引入一个任意的全局相位。这个相位在不同 \(k\) 点之间没有任何联系，FHS 的闭合圈乘虽然能消掉整体相位，但前提是重叠矩阵里的两个态确实来自连续的能带。求解容差变松后，特征向量的方向被污染，重叠矩阵的模偏离 1，闭圈相位也随之失真。
解决：把求解容差收紧，特征值求解器的相对容差从默认 \(10^{-5}\) 调到 \(10^{-8}\) 或更低。另一个实用技巧是检查重叠矩阵的模长：如果 \(|U_x|\) 明显小于 1（比如 0.97），说明相邻 \(k\) 点的态已经不对齐，要么加密网格，要么降低容差。
5.3 COMSOL 参数扫描中模式跳变
现象：扫参后导出的第 1 条带在第 40 个 \(k\) 点上突然变成另一条带，能带图出现“断带”。
原因：COMSOL 默认按特征频率从小到大排序。两条本征频率曲线接近或交叉时，同一序号的解会在交叉点前后发生交换，也就是模式跳变。普通扫描比辅助扫描更容易跳变，因为普通扫描每个参数点都重新初始化。
解决：改用辅助扫描，并在研究设置里开启“特征值模式排序”或手动指定搜索区域，让求解器始终搜索目标频率附近的模式。导出后还可以做一次基于频率连续性的自动重排：对每个 \(k\) 点，把本征频率与上一个点最接近的模式重新编号。
5.4 网格不对称导致能带错位
现象：单胞在几何上是对称的，但扫 \(k_y\) 方向时能带图明显左右不对称，陈数算出来是 2，理论预期却是 -1。
原因：源边界和目标边界的网格节点没有严格一一对应，或者单胞本身在网格划分时被各向异性处理。Floquet 周期条件要求源和目标边界在离散层面完全相同，如果一侧加密另一侧没加密，重叠积分内部的一致性就被破坏了。
解决：在网格设置里选择“周期性网格”，确保相对边界使用相同的网格序列。更常用的做法是先用映射网格手动控制相对边的边界单元数，再生成整体网格。检查方法也简单：导出后比较同一节点编号在两份文件里的坐标，若坐标偏差超过机器精度，网格就有问题。
5.5 k 空间单位与导出坐标不一致
现象：脚本里 G 点（Γ 点）不在网格中心，积分域算出来总面积多了一倍；陈数结果乘以 2 才等于理论值。
原因：COMSOL 里 Floquet 波矢变量可以是归一化值 0 到 1，也可以是 \(-\pi/a\) 到 \(\pi/a\)，不同物理场接口默认不同。如果你在模型里设置了归一化波矢，又在脚本里用了 \(-\pi/a\) 到 \(\pi/a\) 的网格，积分就是重叠覆盖，陈数被放大两倍。
解决：在模型设置里统一约定——推荐在 COMSOL 参数里把kx定义为归一化坐标，脚本里再用kx_phys = kx * 2*pi/a做换算。同时检查布里渊区面积是否为 \((2\pi/a)^2\)：把格点数乘单元面积，与理论面积对比，对不上就先修这个再谈陈数。
6. 验证陈数：用体边对应看边界态是否对得上
陈数算出整数只是第一步，真正让人信服的验证是体边对应：如果某条带的陈数是 \(C\)，那么把周期结构沿一个方向截断后，在两个边界上出现的拓扑边界态数目至少为 \(|C|\)。这个验证在 COMSOL 里做起来不复杂，还能反过来确认第 4 章脚本没有把符号搞反。
做法是把单胞沿 \(x\) 方向复制 \(N_x = 12\) 个，\(y\) 方向保持单胞宽度，并保留 \(y\) 方向的 Floquet 周期条件。研究类型仍然是特征频率，扫 \(k_y\) 从 \(-\pi/a\) 到 \(\pi/a\)。由于 \(x\) 方向被截断，原本的体带会在带隙中出现边界态模式，色散曲线横贯带隙。下面的表格是我常用的验证设置：
| 项目 | 取值 | 作用 |
|---|---|---|
| 超胞尺寸 | \(N_x = 12\)，\(N_y = 1\) | 足够消除两边界之间的耦合 |
| 边界条件 | \(x\) 方向默认自由，\(y\) 方向 Floquet | 模拟无限长边界 |
| 网格 | 边界处加密一层边界层网格 | 避免边界态场集中在边界上时分辨率不足 |
| 扫描参数 | \(k_y\)，32 个点 | 画出连续边界态色散 |
如果陈数计算正确，带隙内的边界态不会与体带交叉，而是从一条体能带的边缘连续延伸到另一条体能带的边缘。边界态的传播方向还和陈数符号相关：符号相反时，边界态的群速度方向在同一个 \(k_y\) 处是相反的。
我自己在第一版计算时，FHS 给出来是 \(C = 2\)，超胞里却数出 4 条边界态，怎么都对不上。排查到最后发现是 \(y\) 方向的 Floquet 相位符号填反，导致边界态来自两个不同带隙的叠加。这个经历给我的教训是：陈数永远不是一个孤立的数，必须和边界态数量一起看；两套验证对不上时，先怀疑符号约定，再怀疑模式排序。希望帮到你。
本文还有配套的精品资源，点击获取
