企业官网建设流程全解析

算四方格子光子晶体的能带和 Wilson loop，前前后后我折腾了小一个月才跑通。Comsol 里算能带本身不难，难的是把场数据导出来、按正确的方式组内积矩阵、再沿布里渊区路径一圈圈乘起来——中间全是细节坑。这篇文章把我完整的操作流程和踩坑记录写出来，从建模到数值实现都尽量给到可以直接复现的步骤，希望能帮后来的人少走点弯路。

1. 项目思路：从能带到 Wilson loop，你到底在算什么

1.1 为什么要算 Wilson loop

在光子晶体里做拓扑研究，已经是常规操作了。但能带本身只告诉你模式在哪些频率范围存在、有没有带隙，想判断这个带隙是“拓扑非平凡”的，还得算拓扑不变量。对某些体系来说，Chern 数可以直接由 Wilson loop 谱得到——更准确地说，Wilson loop 本征值给出了 Wannier center 的分布，而 Wannier center 随动量方向的演化斜率，就是 Chern 数。

这里有个容易搞混的点：二维四方格子是相对简单的晶格，但简单不代表 Wilson loop 好算。因为 Wilson loop 要把布里渊区一条路径上的本征态连续性编码成矩阵乘积，中间涉及规范选取、简并点处理、数值积分精度，任何一个环节马虎都会得到看起来“平滑”却完全错误的结果。

1.2 Wilson loop 在光子体系里的定义

电子体系里 Wilson loop 的定义依赖布洛赫波函数的 overlap 矩阵，光子晶体也一样，只是内积测度变了。对非磁性、无损介质中的光子模式，本征模满足正交归一关系：

∫ ε(r) E*_m(r) · E_n(r) dV = δ_mn

所以构造布洛赫态之间的 overlap 矩阵时，必须带介质介电常数 ε(r) 作为权重。这个加权极容易遗漏——按常规向量内积去做，结果会乱。

具体来说，把布里渊区上一条闭合路径离散成一系列动量点 k_i，沿路径相邻点之间构造矩阵：

M_{mn}(k_i) = ∫ ε(r) E*{m,k_i}(r) · E{n,k_i+Δk}(r) dV

然后按顺序把这些矩阵乘起来：

W = M(k_0) M(k_1) … M(k_{N-1})

W 的本征相位就是 Wilson loop 谱，通常画成随另一方向动量变化的曲线。这套流程在电子体系里有现成程序，但用在 Comsol 输出的光子场上，需要自己实现的地方非常多。

1.3 四方格子高对称点与路径约定

四方格子（square lattice）晶格常数为 a，倒格子基矢选 b₁ = (2π/a, 0)，b₂ = (0, 2π/a)，高对称点为 Γ(0,0)、X(π/a,0)、M(π/a,π/a)。我习惯把布里渊区定义为 [-π/a, π/a) × [-π/a, π/a)，Wilson loop 沿 k_x 方向积分，记录本征相位作为 k_y 的函数。

计算时用了归一化晶格常数 a = 1，频率全部用 ωa/2πc 表示。介质参数是典型的电介质柱而非金属柱，柱体折射率 n = 3.5（对应很多研究所用的高折射率柱），背景折射率 n = 1。虽然具体参数不影响方法，但先把参数固定下来方便对比。

2. Comsol 建模与能带扫描的关键细节

2.1 物理场选择与几何构建

Comsol 里做光子晶体能带最常用的模块是“电磁波，频域”（Electromagnetic Waves, Frequency Domain），选二维组件，本征频率研究（Eigenfrequency）。之所以不用 RF 模块的“波动光学”节点，是因为后者更适合波导问题，而光子晶体能带用电磁波频域配合 Floquet 周期性边界条件是最顺手的。

几何就是一个单位晶胞：先建一个边长为 a 的正方形域，再在中心画一个半径 r 的圆，设定圆的区域为介质柱（折射率 3.5）。这里比较关键的是介质柱的填充比，它会直接影响带隙宽度和拓扑性质。我测试过 r/a 从 0.18 到 0.35 的系列，最终固定为 r/a = 0.25，这个填充比下能在较低的频率范围得到明确的带隙。当然你也可以用背景介质柱反衬的结构（空气孔），但那样介质分布和网格匹配差异比较大，后处理时注意加权不同。

几何构建完成后，把两个区域（柱体和背景）的折射率分别设置好。在电介质问题里，频率本征求解时材料的色散关系会影响本征频率精度，所以固定折射率时要在材料定义里禁用色散模型，或者直接把折射率设为常数。

2.2 布洛赫边界条件的正确设置

布洛赫周期条件在 Comsol 里有两种常见方式：一种是手动对相对边界添加周期性边界条件并设置 Floquet 相位；另一种是直接使用“周期边界条件”中源-目标边界映射，然后指定 k 向量。

我推荐直接用“Floquet 周期性”设置，操作路径是：电磁波频域节点下添加“周期边界条件”，把左右边界设为一对源和目标，上下边界同样设为一对，然后在“周期边界条件”设置里勾选“Floquet 周期性”，把波矢 k = (k_x, k_y) 填进去，并定义对应的单位向量。请特别注意：四方格子下的 Floquet 条件需要正确关联 k 分量与边界方向，如果左边源边界对应 k_x 分量，那右边目标边界的相位关系自然是 exp(-i k_x a)。网格生成的周期映射关系如果选错，算出来看似能带但没有正确的平移对称性，后续 Wilson loop 内积矩阵铁定算错。

另一个细节是，Comsol 的“特征频率”求解器在施加布洛赫周期性时，相位因子全部包含在边界条件里，求解器输出的电场已经在每个网格节点上带有布洛赫相位 e^{ik·r}。这个输出约定对后续数据处理是决定性的：如果直接拿这个场去做 overlap，必须考虑相邻 k 点之间是否有额外的相位差。为了确保一致性，我建议在软件中把电场分量输出到网格节点数据文件时，保留几何坐标，这样可以在后处理中判断相位约定。

2.3 扫描网格与特征频率求解设置

能带计算本身就是一个参数扫描过程。在“研究”步骤中先设定一个“参数化扫描”，扫描变量定义为 k 沿高对称路径的分量。我实际扫描的是二维 k 网格，而不是只扫高对称线，因为 Wilson loop 需要的是整个布里渊区上的二维离散网格。

网格步长我用的 Δk = 2π/80a，也就是一条边 80 个离散点，总 k 网格为 80×80。这个分辨率对四方格子来说通常足够，既有较高的数值精度，又能兼顾 Comsol 求解时间。如果只为了能带图，沿着高对称路径取 40–60 个点就够了，但 Wilson loop 需要精细的二维网格，建议一开始就建立完整网格而不临时加密。

特征频率求解时，每点求解多少个本征么？我建议至少取目标频率范围之上多解 8–10 个本征模，因为 Wilson loop 计算中你可能只关注前三个能带，但矩阵对角化过程中需要避开数值伪模和简并处理。Comsol 的“特征频率”默认使用 ARPACK 求解器，搜索频率区间可以设为，目标频段取 [0, ω_max]，并把要搜索的模式个数设置为多个。我一般把“所需模式数”设为 6–8 个，对单个晶胞网格不多，求解速度完全可以接受。网格用三角形正规化，在最细的介质柱区域加上边界层网格。

提示：每次扫一个 k 点都单独做一次本征求解是慢的，但 Wilson loop 本来就要求逐点求得布洛芬模式场，所以别想用线性插值跳过大网格。合理做法是先在小网格（比如 60×60）上验证方法，再增加到 100×100 的最终精度。

2.4 网格收敛性与能带验证

别急着直接算 Wilson loop。先用解析期望对比验证能带结构。比如四方格子介质柱方阵的 TE 模式和 TM 模式有明确的带隙变化规律：高折射率柱阵列在低填充比时，通常 TM 带隙先打开；如果柱体太大，带隙可能闭合。我算的 r/a = 0.25 介质柱结构，第一带隙正好在 ω/2πc ~ 0.25–0.35 之间出现，与文献典型结果一致。

网格收敛性检查不能省。取网格最大单元尺寸从 a/20 到 a/60 逐步加密，观察第一条带隙边界频率的变化幅度。变化在 0.1% 以内就可以认为收敛。注意 Comsol 默认的网格加密通常给不了那么高的精度，特别是介质柱边界处的网格，一定要用“边界层”或“细化的角网格”。我最终采用的是全局最大单元尺寸 a/30，介质柱边界层网格厚度为 a/80，两层边界层，总自由度约 5 万左右，每个 k 点求解时间在几十秒量级。

3. 核心难点：从 Comsol 场数据到 Wilson loop

3.1 导出场数据并保持归一化

Comsol 可以“导出数据”把每个 k 点的电场分量输出到文本文件，但这里有两个老大难的坑。

第一个坑：网格节点顺序在不同 k 点会变。Comsol 的网格是在几何上固定划分的，理论上节点坐标在扫描过程中不变，但如果你启用了自适应网格细化或参数依赖的网格重整，节点就会变化。稳妥的做法是在导出时同时输出 x、y、Ex、Ey、Ez 以及该网格点对应的介电常数 ε。然后全部数据在 Python 里统一按坐标排序并提取到同一个结构化网格上。先导出一次网格坐标，确认所有 k 点下节点顺序完全一致，再做后续处理。

第二个坑是归一化。Comsol 特征频率求解的本征场默认按某个特定范数归一化，但不一定是光子体系的正交归一。如果在软件内部导出场幅值，直接用于内积计算会导致矩阵 W 不是幺正的，Wilson loop 谱会出现非单位模问题。因此在我们构造 overlap 矩阵之前，必须先对每个 k 点的场做归一化处理：

∫ ε(r) |E_k(r)|² dV = 1

这一步在 Python 里用同一个数值积分网格完成。因为介质柱内部 ε 大，如果不加权，柱内场对积分的贡献会被高估。

3.2 内积矩阵的离散积分公式

场是定义在网格节点的，只要网格够密，二维积分可以用简单求和近似。设网格节点坐标为 (x_i, y_i)，权重为 w_i（Comsol 导出时没有现成权重，可以自己对三角形展开加权）。三角形网格的积分权重并不是均匀的，要正确计算需要从网格数据重建三角形单元。这一步最简单是用 Python 的 mesh 处理工具（比如从节点坐标和单元连接表生成非结构化网格，再结合每个三角形面积计算积分）。

因为节点数据只含坐标和场分量，不包含单元连接，我建议按三角形单元重建积分。但很多情况直接用均匀笛卡尔网格差值后的等权求和也能凑合，前提是网格足够细。如果导出的是规则网格（如后处理中用网格数据到规则网格插值），那积分权重就是均匀的 Δx·Δy。

于是内积矩阵元：

M_{mn}(k) = Σ_i ε_i E*{m,k}(x_i) E{n,k+Δk}(x_i) ΔA_i

注意这里 ΔA_i 是第 i 个积分点的面积权重。如果直接用非结构网格数据做求和必须插到规则网格上再积分，否则单元权重计算比较复杂。我导出后先用线性插值映射到 128×128 的均匀网格上再做积分，好处是可以检测并去除因网格节点顺序抖动的异常。

3.3 Wilson loop 主流程代码

从 Comsol 导出数据后，剩下的算法流程我用 Python 实现。下面给出核心代码，主体是沿 k_x 方向的 Wilson loop：

import numpy as np from scipy.linalg import eigvalsh # 假设： # 每个 k 点导出的场存成一个字典结构 field_data # field_data[kx_idx, ky_idx]["epsilon"] 是规则网格上的介电常数分布（128x128） # field_data[kx_idx, ky_idx]["E"] 是形状为 (n_bands, 128, 128, 3) 的电场分量数组 # Nkx, Nky 分别是沿 kx 和 ky 的网格点数量 def compute_overlap(E_k, E_kp, eps, dA): """计算 M_{mn}(k) = int eps * E_m,k^* · E_n,kp dA""" n_bands = E_k.shape[0] M = np.zeros((n_bands, n_bands), dtype=complex) for m in range(n_bands): for n in range(n_bands): # 场分量点乘（x,y,z）并乘 eps 权重 integrand = eps * np.sum( np.conj(E_k[m]) * E_kp[n], axis=-1 ) M[m, n] = np.sum(integrand) * dA return M def normalize_modes(E, eps, dA): """按光子模正交归一化 eps*E*_m · E_n = delta_mn""" n_bands = E.shape[0] norm = np.zeros(n_bands) for m in range(n_bands): norm[m] = np.sqrt(np.sum(eps * np.sum(np.abs(E[m])**2, axis=-1)) * dA) return E / norm[:, None, None, None] def wilson_loop(field_data, Nkx, Nky, band_indices, dA): """沿 kx 方向积分 Wilson loop，返回每个 ky 下的本征相位列表""" wilson_phases = np.zeros((Nky, len(band_indices))) n_bands = len(band_indices) for ky_idx in range(Nky): # 初始化 Wilson loop 为单位矩阵 W = np.eye(n_bands, dtype=complex) for kx_idx in range(Nkx): E_k = field_data[kx_idx, ky_idx]["E"][band_indices] E_kp = field_data[(kx_idx + 1) % Nkx, ky_idx]["E"][band_indices] E_k = normalize_modes(E_k, field_data[kx_idx, ky_idx]["epsilon"], dA) E_kp = normalize_modes(E_kp, field_data[(kx_idx + 1) % Nkx, ky_idx]["epsilon"], dA) M = compute_overlap(E_k, E_kp, field_data[kx_idx, ky_idx]["epsilon"], dA) # 为了数值稳定，先进行极分解（保证 W 尽量幺正） U, _, Vh = np.linalg.svd(M) M_orthogonal = U @ Vh W = M_orthogonal @ W # Wilson loop 本征值取相位 eigvals = eigvalsh(W) wilson_phases[ky_idx] = np.angle(eigvals) return wilson_phases

这段代码有几个细节值得说明。第一，normalize_modes必须在每个 k 点单独做，不能跨 k 做全局归一化，因为不同 k 点是各自独立的模式解。第二，M_orthogonal的极分解非常重要。由于数值误差和导出网格不完整，M 可能略微偏离幺正矩阵，直接用 M 乘下去会让 Wilson loop 的本征模长得不到纯相位结果，甚至会导致相位处处跳变。先做极分解把矩阵强制拉到 SU(N) 群上，物理结果的稳定性会好很多。第三，沿 k_x 路径闭合的条件是(kx_idx + 1) % Nkx，因为 Wilson loop 本身要求路径闭合。最后scipy.linalg.eigvalsh用来算厄米矩阵特征值也适用于这里的幺正矩阵，相位用 np.angle 取。

3.4 规范连续性与能带选择

Wilson loop 对角化的一个经典问题是相邻 k_y 之间本征相位排序不一致。由于每个 k_y 下对角化出来的本征值都是按相位大小排序的，当两个相位接近 ±π 端点的交接时，排序会跳开，画出来的图有锯齿。解决方法是做“平行输运规范”：对每个 k_y 的 Wilson loop 本征向量，与上一个 k_y 的本征向量求 overlap，最大化对角模，然后用相位连续性追踪每一条 Wilson 谱线。

更简单的做法是：把所有相位按大小排序后，再对相邻 k_y 使用最小差值匹配。我在实际中使用的是“相位反绕”加“排序匹配”两步，保证最后画出来的曲线连续。还有一个更稳妥的方案是额外计算相邻 k_y 之间的 overlap 矩阵，然后构造混合 Wannier 中心。对此问题的严格处理较复杂，但如果你只是要看 Wilson loop 的斜率，平滑追踪已经足够。

能带选择上，我建议取单一带隙下方的连续能带组。比如第一带隙下面只有第 1、2 条带（对应横磁或横电模各一组，取决于 2D 模型里你是只做面外场分量还是全部三个分量）。如果你的二维模型包含完整矢量 TE+TM，前几条带包含两种极化，这时的 Wilson loop 会因为极化耦合和额外简并变得非常混乱。最简单的做法是直接在仿真模型里抑制其中一个极化。实现方式是用 2D 面外分量做 TM 模计算或用面内分量做 TE 模计算。我平时在“电磁波，频域”中直接用“面外磁场矢量 + 面内电场”构建 TE 模，或只用“面外电场”构建 TM 模，这样能带分组清晰，Wilson loop 的物理结果不歧义。

如果你非要算全矢量情况，必须用能带序数和简并判据把希尔伯特空间投影到目标能带组上。此时 Wilson loop 矩阵是分块对角化的，且各块之间几乎不混合。算法实现比这里给出的代码复杂得多，需要在每个 k 点对带组内的所有模式做 Gram-Schmidt 正交化，并额外处理简并点。我的经验是：要算拓扑性质时，绝大多数论文模型都明确极化，2D 单极化模型已经足够。

3.5 开带隙与简并点的处理策略

四方格子高对称点，特别是 M 点，经常出现能带简并。如果两条带简并，叠加后的 Wilson loop 矩阵处于简并子空间，相位无法唯一确定。处理简并有两条路：一是避开简并带组，选择完全孤立的能带组（如第一带隙下方的第一条带）计算 Wilson loop；二是在简并点位置微扰采样，这样能带在数值上不再严格简并，但矩阵条件数仍然很糟糕。

我的实际操作是：检查每个 k 点相邻能带的频率差。如果某一 k 点下两条能带频率差的绝对值小于 1e-4（归一化频率），就标记为简并点。对于 Wilson loop 积分最稳妥的方法是只在非简并的 k 网格上计算，并且确保沿着 k_x 路径上每一点都与目标能带组内的其他带保持足够大的差距。如果不小心踩到简并点，你可以把网格里某个方向的采样点偏移半个步长来绕过。注意，这种故意绕开简并点的做法物理上是允许的，因为 Wilson loop 只是定义在带隙内且某些全流程积分的拓扑不变量在变形连续时不依赖路径是否经过简并点，只要你的带组在布里渊区上是全局分离的。

话虽如此，如果两条带在绝非简并但任何地方都碰巧靠近（带隙小于数值误差），那你选这一组能带本身做 Wilson loop 就是危险的——最终谱会非常乱。此时请检查模型参数：是不是网格不够密导致伪简并？是不是该降低对称性（比如略微椭圆变形）打破简并？对四方格子，你可以把介质柱改为边长微小不等的矩形柱，形成微小的四方畸变。这样做会让 M 点简并被解除，Wilson loop 谱变得连续可读。我为了验证方法的稳定性，就做了这个测试：矩形柱的拓扑相不因对称性微扰而突变，反而是能带分裂后更容易计算。

4. 常见问题与排查速查表

为了节省同行踩坑时间，把典型问题、原因和解决方式整理成表。

4.1 Wilson loop 谱为什么总在 ±π 处来回弹跳

这个现象我刚跑通时也被折磨了很久。其实根源不复杂：Wilson loop 本征相位是角变量，取值范围 [−π, π]，如果某条 Wannier 中心曲线从 +π 一侧跑到 −π 一侧，再排序追踪时就会误以为发生跳变。解决办法是“相位解缠”：对每条谱线，如果相邻 k_y 相位差大于 π，就加上或减去 2π，让曲线连续。实际操作中还需要先排序再解缠，否则在两条谱线交叉处会出现错误配对。

def unwrap_phases(phases, tol=np.pi): """phases shape (Nky, n_bands)，按每列相位解缠并排序追踪""" Nky, nb = phases.shape unwrapped = np.zeros_like(phases) for b in range(nb): unwrapped[0, b] = phases[0, b] for i in range(1, Nky): delta = phases[i, b] - unwrapped[i-1, b] while delta > tol: delta -= 2*np.pi while delta < -tol: delta += 2*np.pi unwrapped[i, b] = unwrapped[i-1, b] + delta return unwrapped

这段代码虽然简单，但很实用。不过要注意它只适合相邻 ky 网格足够密且没有真正简并的情况。如果相邻 ky 点相位差本身就大于 π（网格太粗），解缠会失败。所以更根本的解决方案还是加密 k_y 方向的采样点数。

4.2 Comsol 导出数据量太大，卡死怎么办

80×80 的 k 网格意味着 6400 次本征求解，每次要导出多个带的数据，文件会非常庞大。建议分层导出：先在 Comsol 里按 k 点循环求解，每个 k 点只导出所需的能带序号和场分量到一个单独的.txt文件；然后立即用 Python 脚本读完该文件并压缩成.npy格式，删除原始导出文件。这样整个流程下来数据占用能控制在几百 MB 内。另外一个技巧是：Comsol“导出数据”时不用把整个网格都导出来，只导“边界内”的目标区域，跳过布里渊区内无介质区域（如果靠近晶胞角落场本身就很小），可以大幅降低输出量。但你必须保证积分区域完整，不然积分归一化会亏。

4.3 后处理网格插值导致积分不守恒

把无序网格导出的场插值到规则网格上经常会引入误差。一个典型问题是介质柱边界附近的 ε 突变，如果插值后的规则网格没有精确标记介质边界，ε 的奇异分布会让积分产生误差。解决办法是：插值到同一规则网格之前，先对 ε 做高斯模糊或平滑处理，但这样会改变模型的物理参数。更稳妥的替代方案是让 Comsol 在“网格”节点里把一个规则网格作为导出目标，但这在 2D 非结构网格场景不符实际。我最后采用的方法是双线性插值后再叠加一个正确的 ε 场值。具体做法是：从 Comsol 导出的几何表达式中获取每个规则网格点的 ε 值（这个等效于用光滑函数判断介质区域），而不是由其他场的插值推导。也就是说，ε(x,y) 是解析赋值而非插值，这样积分结果对网格规则化层次不再敏感。

5. 计算流程以外：个人实测经验与后续扩展

5.1 数值稳定性：奇点附近的多网格测试

Wilson loop 计算结果的可靠程度，建议做一个系统测试：固定 k_y，取 k_x 路径步长 Nkx 分别为 40、60、80、120，看 Wilson loop 谱的本征相位随步长如何变化。我测试的结果是当 Nkx 从 40 升至 80 时，相位变化约 0.05 rad，继续加密到 120 时变化小于 0.01 rad。但如果你的相位变化在 Nkx = 80 到 120 时还远超 0.01 rad，说明在某个 k 点附近有很弱的能带交叉或者数值伪模，此时必须检查该 k 点处能带的模分离度。从这个意义上说，网格收敛性测试是 Wilson loop 的必修课，而不是可选流程。

5.2 从 Wilson loop 读物理结果

如果沿 k_x 积分得到的 Wannier 中心 γ(k_y) 随 k_y 增加而明显线性上升（斜率非零），则第一个带隙下方的拓扑不变量即 Chern 数非零。对于四方格子典型结构，Chern 数通常为 0 或 1，你可以从 Wilson loop 的总绕数来读。这个绕数就是弧长的整数倍，等价于本征相位的总变化。注意一次 Wilson loop 的谱可能包含多条线，每条线对应不同能带下的 Wannier center。如果多带组有总绕数非零，就说明这个带组是拓扑非平庸的，边缘态数量也能对应上。

当你看到 Wilson loop 谱出现“穿过整个布里渊区的斜线”时，别急着开心，还要检查布里渊区周期边界处谱线是否连接成一条闭合轨道。如果 k_y 从 −π 到 π 的谱线首尾没对齐，说明相位解缠走了样，或者带组选择有振荡。正确的结果在周期边界下一定是闭合的。

5.3 这个方法还能怎么用

Wilson loop 不仅限于四方格子。任何二维布洛赫体系的高对称路径扫描都可以套用同样的流程：改一下倒格子基矢、替换 Floquet 周期向量，然后沿想要的方向做积分即可。比如六角格子（石墨烯型）、长方格子、Kagome 格子，甚至带倾斜角度的单斜格子，全部适用。关键就是别忘记四个要素：带组选择、ε 加权内积、极分解保幺正、相位解缠排序。

如果你的是三维光子晶体，也能推广为体积分路径，只不过积分数值代价上一个大台阶，通常需要高性能计算做。此时再用纯 Python 做 128×128×128 的积分就吃紧了，建议改成有限元网格上的并行积分，或直接用 Comsol 内置的积分子来构造内积，再以矩阵形式导出。计算量虽然大，但算法的骨架完全一致。

再说一个实用建议：计算过程中保留好每一步导出的中间文件，包括每个 k 点的本征频率、归一化场、内积矩阵。这不仅是为了复现，更是为了排查问题。我最后悔的就是前期没保留中间矩阵，出错了也不知道是前三步哪一步出的问题。后来我写了个脚本，每一步都会输出一个小的状态文件和校验哈希，定位问题效率提升不少。

最后分享一个细节——在你打印最终 Wilson loop 图之前，习惯性检查一下每个本征值模长和矩阵行列式的数值。行列式偏离 1 超过 1e-4 的话，先回头把极分解做好，再去琢磨物理结论。这些数值上“挑刺式”的校验，才是整个计算链路不翻车的根本保证。