四方格子光子晶体能带与Wilson loop计算实践指南

1. 项目概述：四方格子光子晶体能带与Wilson loop计算

四方格子光子晶体是光子晶体研究中的经典模型结构，其周期性介电常数分布形成的能带结构对光场调控具有重要意义。Wilson loop作为拓扑光子学中的重要工具，能够有效表征光子晶体的拓扑性质。在COMSOL Multiphysics中完成这一系列计算，需要跨越电磁场仿真、能带计算和后处理分析三个关键环节。

我最初接触这个课题时，发现现有文献大多只给出理论公式或最终结果，对COMSOL具体操作细节和常见问题避而不谈。经过多次尝试和调试，总结出一套可靠的工作流程。本文将重点分享从模型建立到Wilson loop计算的完整过程，特别是那些容易出错的参数设置和数据处理技巧。

2. 模型建立与参数设置

2.1 四方格子光子晶体基本结构

四方格子光子晶体由介质柱在空气中周期性排列构成，典型参数包括：

• 晶格常数a = 1 μm（归一化单位）
• 介质柱半径r = 0.2a
• 相对介电常数ε = 12（模拟硅材料）
• 空气区域ε = 1

在COMSOL中建立模型时，几何构建需注意：

1. 使用"周期阵列"功能而非手动复制基本单元
2. 设置足够大的外围空气区域（至少3a）以减少边界效应
3. 明确区分材料边界，避免网格生成时的几何混淆

关键提示：介质柱边缘的网格密度直接影响计算精度，建议设置至少5层边界层网格

2.2 物理场和边界条件配置

选择"电磁波，频域"物理场接口，关键设置包括：

% 对应COMSOL中的材料参数表达式 epsilon = (x^2+y^2<=r^2)*12 + (x^2+y^2>r^2)*1; mu = 1; % 非磁性材料

边界条件配置要点：

1. 使用Floquet周期边界条件处理晶格周期性
2. 完美匹配层(PML)厚度设为1/2工作波长
3. 端口激励设置为禁用（仅特征频率研究）

3. 能带计算关键技术

3.1 布里渊区路径选取

四方晶格的不可约布里渊区路径为Γ→X→M→Γ：

k_path = [0,0; 0.5,0; 0.5,0.5; 0,0]; % 标准化k点坐标

在COMSOL中实现时：

1. 创建参数化扫描研究
2. 设置波矢量k为扫描参数
3. 每个k点计算6-8个模式以确保完整性

3.2 求解器配置技巧

特征频率研究的关键参数：

• 搜索频率范围：0.2-0.8 c/a（覆盖典型光子带隙）
• 搜索方法：shift-invert
• 网格数：至少10,000个自由度
• 收敛容差：1e-6

常见问题处理：

• 出现虚假模式 → 检查材料定义和边界条件
• 模式交叉 → 启用模式跟踪功能
• 收敛困难 → 调整初始猜测频率

4. Wilson loop计算实现

4.1 能带数据导出与处理

计算完成后，需要导出所有k点的本征模式和场分布：

1. 导出为.mat格式保持数据结构
2. 在MATLAB中重组为(k,E,ψ)三维数组
3. 对能带进行排序和编号

典型数据处理代码框架：

load('band_data.mat'); num_bands = size(E,2); num_k = size(E,1); % 能带排序 for k_idx = 2:num_k [~,order] = pdist2(E(k_idx-1,:)', E(k_idx,:)', 'euclidean', 'Smallest',1); E(k_idx,:) = E(k_idx,order); psi(:,:,k_idx) = psi(:,:,k_idx)(:,order); end

4.2 Wilson loop算法实现

Wilson loop计算核心步骤：

1. 沿k路径离散化采样
2. 计算相邻k点间的重叠矩阵：

   M = psi(:,:,k)' * psi(:,:,k+1);

3. 累积乘积得到Wilson loop算符：

   W = eye(num_bands); for k = 1:num_k-1 [U,S,V] = svd(M); W = W * U*V'; end

4. 计算本征相位得到Wannier中心：

   theta = angle(eig(W));

注意事项：相位缠绕(phase wrapping)问题需特殊处理，建议使用unwrap函数

5. 常见问题与解决方案

5.1 能带计算不收敛

可能原因及对策：

1. 网格太粗糙 → 加密网格特别是介质边界处
2. PML设置不当 → 调整PML层数和拉伸参数
3. 初始猜测不准 → 先用较大范围扫描再局部细化

5.2 Wilson loop相位跳变

典型现象：相邻k点相位差超过π 解决方法：

1. 增加k点采样密度
2. 实施相位连续性校正：

   for n = 2:length(theta) while theta(n)-theta(n-1) > pi theta(n) = theta(n) - 2*pi; end while theta(n)-theta(n-1) < -pi theta(n) = theta(n) + 2*pi; end end

5.3 计算资源优化

大型模型加速技巧：

1. 使用对称性减少计算域
2. 分布式计算参数扫描
3. 适当降低收敛精度要求
4. 采用渐进式网格加密策略

6. 结果分析与可视化

6.1 能带结构绘制

标准可视化流程：

figure; hold on; for band = 1:num_bands plot(k_points, E(:,band), 'LineWidth',1.5); end xlabel('Wave vector'); ylabel('Frequency (c/a)'); set(gca,'XTick',[0,0.5,1,1.5],'XTickLabel',{'Γ','X','M','Γ'});

6.2 Wilson loop结果展示

拓扑不变量计算示例：

figure; plot(k_points(1:end-1), theta/(2*pi), 'o-'); xlabel('Wave vector'); ylabel('Wannier center (2π)'); title('Wilson loop spectrum');

典型分析要点：

1. 寻找跨越整个布里渊区的Wannier中心
2. 计算Chern数：θ(2π)-θ(0)
3. 识别拓扑边缘态的存在

在完成整套计算流程后，我发现介质柱边缘的网格处理对结果精度影响最大。通过对比测试，当边界层网格达到7层时，Wilson loop相位的波动可以控制在0.05π以内。另一个实用技巧是在MATLAB中使用parfor并行处理不同k点的数据，能使整体计算时间缩短40%左右。