1. 蜂窝晶格为什么值得算狄拉克锥与拓扑光子学的切入点做拓扑光子学研究的人迟早都会碰到蜂窝晶格光子晶体。这个体系几乎就是为演示“如何从能带计算走向拓扑不变量”而生的。它的晶格结构和石墨烯一模一样由两套三角子格穿插构成所以天然具备狄拉克锥而只要破坏了空间反演对称性狄拉克点处就会打开一个带隙进而可以定义陈数、研究边缘态。很多初学者以为能带图算出来就万事大吉了实际上从Comsol模型到Matlab里的陈数结果中间隔着不少容易踩坑的环节。我在这篇文章里会把这套流程完整地拆开讲蜂窝晶格单胞怎么建、布洛赫周期边界怎么设、k路径怎么扫、本征场怎么导出来最后又怎么用Matlab脚本从实空间场分布出发算出Berry曲率和陈数。内容适合有有限元基础、想做拓扑光子学数值验证的研究生也适合想快速搭一套“能带计算拓扑不变量”工作流的工程师。里面所有参数和脚本逻辑都是我在实际项目中反复验证过的方案。1.1 从能带折叠到狄拉克锥蜂窝结构的关键是两套格子蜂窝晶格在实空间里的形象是一张六角网但计算时并不推荐直接画一个六边形单胞而是把它看成两套三角子格。设晶格常数为a两个基矢取为a1 (√3/2·a, 1/2·a)a2 (√3/2·a, −1/2·a)A子格放在原点B子格放在最近邻位置比如(0, a/√3)。两套子格在几何上完全等价时系统的空间反演对称性没有被破坏倒易空间第一布里渊区的K和K′点会出现两重简并能带形成线性交叉也就是狄拉克锥。对应到光子晶体里如果A、B位置放的介质柱半径相同、介电常数相同那么TE或TM模式的能带在K点附近一定会呈锥形。这个结构在Comsol中非常好验证把A、B柱半径设成相等扫出一条能带曲线你会在带隙关闭的状态下看到两条带恰好接触。这里有一个很容易忽略的前提——几何对称性必须被数值网格忠实地保留。如果网格剖分破坏了C3旋转对称性K点的简并度会被人为地打开狄拉克锥就会变成一条假带隙拓扑结论自然全部作废。所以后面讲网格剖分的时候我专门花了一节来说这个事。1.2 对称性破缺决定带隙能否打开和陈数直接相关一旦A、B两套格子的参数不再相同比如让A柱半径rA大于B柱半径rB空间反演对称性就被破坏。K点的两重简并上升原先线性交叉的两条能带被分开形成一段带隙。这段带隙的重要性在于它给了我们一个干净的“频率窗口”可以定义某条占据带通常取带隙下方的一组带的拓扑性质。这里需要特别说明光子晶体中“陈数”的定义方式。光子本身是玻色子不存在电子那样费米面和占据态但我们可以把布洛赫模式当成一个参数空间中的向量场对某个频带的态进行投影构造出等价于量子Hilbert空间的结构。具体做法是只看带隙以下那条带或者那组互相简并的带把它当成“占据带”在布里渊区内对Berry曲率积分得到的整数就是陈数。这个整数直接和体边对应关系挂钩带隙内存在单向传输的边缘态边缘态的数量和方向就由陈数决定。所以整个项目的逻辑链其实很简洁蜂窝晶格给出狄拉克锥对称性破缺打开带隙带隙的非平庸性质由陈数描述。下一步Comsol负责给出能带和本征场Matlab负责把场的相位信息提炼成拓扑判据。1.3 这个体系的“陈数”到底在描述什么很多人刚接触陈数时总被抽象公式吓到。我习惯用一个类比来解释把布里渊区想象成一个封闭曲面每个k点处都有一个本征向量也就是布洛赫函数的周期部分。这些向量会随着k变化而指向不同方向。陈数描述的就是这个向量场在封闭曲面上缠绕的“圈数”。就像莫比乌斯带和普通纸带的扭转数量差一个整数没法通过连续形变改变陈数也是一个无法连续变化的整数拓扑不变量。对于蜂窝晶格光子晶体陈数通常不止一个值。如果介质柱半径rA略大于rB某几条频带可能得到1或−1的陈数而rA小于rB时符号反转。符号的意义直接映射到边缘态的传播方向——正陈数的边界上边缘态沿某一方向单向传播负陈数则反过来。这个特性就是拓扑保护的根源即使路径上有弯折、有缺陷只要带隙没有被填满单向传播就不会被背散射破坏。我在后面会用这个性质和边缘态能带做交叉验证。2. Comsol建模与MPH工作流单胞怎么画、边界怎么设才能不出伪模2.1 先澄清“MPH”这个说法在Comsol生态里的实际含义很多人以为“MPH算法”是什么高深的学术名词实际上在Comsol的语境中.mph就是模型文件的扩展名。一份.mph文件包含了整个项目的一切几何、物理接口、网格、研究步骤、结果和后处理。所谓“含MPH算法”合理的解读就是依托Comsol的.mph模型文件来组织整套计算流程的方法。我个人喜欢把“MPH”拆成三个动作Modeling建模与参数化、Processing求解与数据提取、Harvesting后处理与拓扑判据计算。这样工程文件的结构会很清晰模型目录下放几何和物理研究目录下放参数扫描结果目录下放导出数据Matlab脚本单独放一个文件夹。如果你后续要换体系、换晶格类型只需要改参数和几何研究步骤几乎不用动。2.2 几何搭建和参数化用晶格矢量而不是普通直角坐标在Comsol中创建一个二维模型物理接口选择“电磁波频域”Electromagnetic Waves, Frequency Domain。几何上我们建一个平行四边形单胞而不是正六边形单胞。原因是平行四边形单胞的两条边正好对应a1、a2两个平移矢量Floquet周期边界条件的两个相对边就沿这两个方向施加周期关系一目了然。具体参数可以这样设置全局参数a 800 nm作为晶格常数A柱半径rA 0.2aB柱半径rB 0.1a对称时令两者相等背景介电常数eps_bg 1柱介电常数eps_rod 12这是硅在近红外的典型值平行四边形单胞的两个相邻边向量就是a1、a2绘制时将左下角定在原点。介质柱在单胞内的位置需要严格对应子格坐标。我推荐用“几何”里的“圆”来添加两个圆域圆心分别放在(0,0)和(0, a/√3)。注意不要用普通坐标要用参数表达这样以后改a时整个模型自动更新。材料域分配上背景域设为空气圆形域设为高介电材料。为了后面网格收敛性验证我还会单独加一个参数mesh_size来控制网格最大单元尺寸。还有一个细节是单位问题。Comsol的几何尺寸一般带单位你用nm画图时特征频率结果会是Hz级别的高频数看起来很不直观。我习惯把几何单位设成μm微米频率出来后换算到常用的归一化频率f a / c也就是在Matlab后处理里统一处理这样能带图的纵轴才具有普适性。2.3 布洛赫边界条件的关键周期性结对与网格一致性这一步是整个Comsol建模中最容易出问题的地方。在物理接口中添加“周期条件”类型选择“Floquet周期”然后为两对相对边界分别指定波矢分量。Comsol中通常需要输入两个方向的波数kx、ky我建议先在全局参数中定义kx k1 * (2pi/(sqrt(3)a)) k2 * (2pi/(sqrt(3)a))ky k1 * (2pi/a) - k2 * (2pi/a)这里的k1、k2是倒格矢约化坐标后续k路径扫描会以它们为参数。这里有一个到k空间直角坐标的换算必须和你的基矢定义严格对应。最稳妥的做法是用约化坐标(k1, k2)作为扫描变量kx、ky都通过参数表达式实时计算这样最后导出的数据和Matlab脚本里的k网格也是一一对应的。周期性边界条件还有一个坑相对边界上的网格必须一一配对否则Comsol会在边界插值上引入误差产生大量伪模。你需要进入网格步骤打开“周期性结对”Periodic Pair功能让左右边界、上下边界的网格节点精确对应。我自己的经验是如果忘了配对特征频率会多出很多能带曲线看起来就像一团毛线完全没有规律而一旦配对正确哪怕网格密度一般低频段的能带也会非常清爽。另外在Floquet边界条件下Comsol内部的场变量自动带有布洛赫相位exp(ik·r)的约定。这一点对后面的Matlab相位修正特别重要——我建议在建模完成后先做一个自检算一个均匀背景没有介质柱的平行平板模式把解析色散关系和Comsol结果对比确认k空间相位约定到底差一个正负号。这个自检两小时就能做完能省掉后面调试陈数的三天时间。2.4 求解器设置特征频率研究、搜索基准、所需模态数研究类型选择“特征频率”Eigenfrequency。在研究设置里需要指定“所需特征频率数”我一般取8到12个覆盖带隙附近主要频段即可。搜索基准值也要设置如果目标带隙在归一化频率0.3附近可以根据公式f (0.3*c/a)换算成Hz填入。这个基准值不需要非常精确但能极大提升特征值求解器的收敛速度也能避免求解器抓住一堆高频模式。参数扫描要在“研究扩展”里使用“辅助扫描”Auxiliary sweep而不是简单地在研究里嵌套扫描。辅助扫描可以保证每个参数组合都执行完整的特征值求解步骤而且结果数据会用一个额外的维度标记每个扫描点方便后面一维绘图。我们把s设为扫描变量范围0到3步长0.02左右然后用参数表达式把s映射到k1、k2上。求解器设置里还有一个容易忽略的点特征值求解器默认使用SPOOLES或MUMPS对于二维小模型都够用。如果模式数多、网格细建议切换成MUMPS并打开“稀疏直接求解器”的“行预排序”选项速度通常能改善不少。对蜂窝晶格这个体系模型规模其实不大Intel八代以上处理器一般一两分钟就能扫完一个k点整条路径扫描在半小时到一小时之间。3. 能带路径扫描与数据导出把k空间“走”成一条连续的曲线3.1 波矢路径和参数映射表Γ-K-M-Γ的换算蜂窝晶格的倒空间高对称路径通常取Γ→K→M→Γ路径在倒格矢坐标下可以用一个统一的路径参数s来定义。s从0到3每段对应一个高对称段。我建议把k1、k2写成s的分段表达式s范围k1k20到1s/3s/31到21/3 (s−1)/61/3 − (s−1)/32到31/2 − (s−2)/20比如s0代表Γ点s1代表K点s2代表M点s3回到Γ点。把这些表达式直接写进Comsol的全局参数定义里k1、k2就变成s的函数了。此时也可以定义kx、ky的表达式如前文所述。需要注意的是这里的s的物理含义是无量纲路径比例而不是k的大小所以能带图的横轴直接用s或者累计路径长度都行。我个人习惯在Matlab里把s映射成实际倒空间距离这样横轴单位是1/a看起来更专业。3.2 在Comsol中设置辅助扫描的完整步骤我在Comsol界面的实际操作流程如下在“全局参数”中定义a、rA、rB、eps等几何和材料参数再定义s、k1、k2、kx、ky。在“研究1”中选择特征频率研究在“研究设置”里填所需特征频率数。展开“研究扩展”勾选“辅助扫描”在扫描参数里选择s填入范围0到3、步长0.02。确保“在所有参数组合上求解”被勾选。运行研究Comsol会自动对每个s做一次特征值求解。这里有个经验数值步长0.02对应的路径段上有大约150个k点已经足够画出一条光滑的能带曲线。如果你后面要做精确的Wilson loop建议路径扫描的s步长再细化到0.005不过这会让总求解时间拉长一倍以上前期调模型时保持0.02就好。3.3 导出频率和场数据的推荐方式数据导出有两种主流路线取决于你手头有没有LiveLink for MATLAB。如果你有LiveLink那最理想的做法是直接在Matlab里通过COMSOL Java API控制模型扫描s后把频率和本征场一次性读入工作区所有数据都不落盘。这个方案最适合反复迭代缺点是前期需要写一段比较长的调用脚本。我自己的项目用的是这个方案后面给Matlab脚本时也会给出对应的数据接口设计。如果你没有LiveLink也不用慌Comsol GUI导出的数据足够完成陈数计算。需要导出两个东西频率数据通过“派生值 → 全局计算”得到每个s下的所有特征频率把表格导出为文本文件每行是“s, f1, f2, …”。场数据在“数据集”中选择某个参数点的主特征模式然后“导出 → 解”中选择你关心的域勾选坐标列以及Ez的实部、虚部导出为CSV文件。文件格式大致是x, y, Re(Ez), Im(Ez)0.000, 0.000, 2.34e-1, -1.22e-3...要特别提醒的是Comsol导出本征场时实部和虚部通常是分离的而且归一化到能量或某种绝对尺度。这意味着不同k点导出的场幅度不一定可比但对Wilson loop重叠积分来说每个k点内部的归一化已经足够——重叠积分的值会在适当地归一化后变为酉矩阵元素。我建议在Matlab脚本中显式地对每条能带的场向量做一次单位化避免幅度不匹配导致的重叠矩阵退化。3.4 网格收敛性校验一个系数就能判断结果可不可信能带计算和有限元的所有计算一样必须先确认结果不受网格影响。蜂窝晶格光子晶体有一个天然的探针对称状态下K点的狄拉克频率。它必须和网格密度无关地收敛到某个值。我的做法是分别用最大单元尺寸为a/10、a/20、a/30三套网格各算一次对称状态下的K点频率和带隙关闭程度。如果两次加密之间狄拉克频率变化小于0.5%就认为网格足够。另一个更苛刻的验证是rA≠rB时带隙宽度的收敛。网格太粗通常会高估带隙因为数值色散会人为加大模式分裂网格加细之后带隙会单调下降到一个平台。另外圆柱边界处的网格需要单独加密。圆边界如果只用粗网格相当于用一个多边形去近似圆弧这等于在几何层面额外引入了对称性破缺。我一般会在“网格”里添加一个“边界层”沿圆柱边界设置5到8层边界层网格厚度因子设为0.2左右。这一步对狄拉克锥频率的精度影响非常大值得每次建模都加上。4. 陈数计算的Matlab脚本详解从本征场到Berry曲率积分4.1 核心公式和数值算法为什么用plaquette Wilson loop陈数的定义是Berry曲率在布里渊区上的积分C (1/2π) ∫ F(kx, ky) d²kBerry曲率F由Berry联络A i⟨u|∇_k u⟩的旋度给出。数值实现时我们不会去显式求导而是使用Wilson loop方法把k空间划成很多小方格plaquette每个小格子的四个角点各自对应一个本征场u(k)。然后构造相邻k点之间的重叠矩阵M(k, k′) ⟨u(k)|u(k′)⟩由四个角点围成闭环的Wilson loop相位就是该小格上的Berry曲率U_total M(k1→k2) M(k2→k3) M(k3→k4) M(k4→k1)F ≈ arg(det(U_total))把所有小格子的相位加起来除以2π就得到了陈数。这个算法的好处是回避了本征态相位规范选择的问题因为每条能带的全局相位在共轭内积中抵消了。4.2 Matlab读取Comsol数据的格式约定我建议建立一个统一的数据目录结构comsol_output/ klist.txt % 每一行: kx, ky, s freq_list.txt % 每一行: s, f1, f2, f3... fields/ k001.txt % 每个k点的场文件 k002.txt ...场文件的格式为x, y, Re(Ez_band1), Im(Ez_band1), Re(Ez_band2), Im(Ez_band2), ... 0.0, 0.0, 2.3e-2, -1.1e-3, ...Matlab读取脚本里我建议固定按行读入并reshape成网格矩阵然后用实际网格坐标做数值积分。为了让Wilson loop计算省事我推荐在Comsol导出场时使用统一的矩形网格坐标这样不同k点之间的场可以直接配对。如果拿到的场分布是三角形网格上的散点就要先在Matlab里用griddata插值到统一网格。这一步我强烈建议在Comsol导出时就做好在导出节点中多花几分钟指定网格坐标比之后在Matlab里插值省事得多也避免插值误差积累进Wilson loop。4.3 Wilson loop矩阵构造和相位追踪下面是核心的Matlab代码片段。我先给出主循环部分然后单独说明两个重要细节。% 载入k列表和场数据 [kx, ky] load_klist(comsol_output/klist.txt); fields load_fields(comsol_output/fields/, length(kx)); bands 1:2; % 占据带索引例如带1~2 Nx 30; Ny 30; % k网格数实际使用需根据数据覆盖调整 C 0; for i 1:Nx for j 1:Ny % 四个角点的索引 idxA sub2ind([Nx, Ny], i, j); idxB sub2ind([Nx, Ny], i1, j); idxC sub2ind([Nx, Ny], i1, j1); idxD sub2ind([Nx, Ny], i, j1); % 构造各边重叠矩阵 Uab wloop_matrix(fields(idxA), fields(idxB), bands); Ubc wloop_matrix(fields(idxB), fields(idxC), bands); Ucd wloop_matrix(fields(idxC), fields(idxD), bands); Uda wloop_matrix(fields(idxD), fields(idxA), bands); % 闭合环路相位 Uloop Uab * Ubc * Ucd * Uda; F angle(det(Uloop)); C C F; end end C C / (2*pi);wloop_matrix函数的实现如下function U wloop_matrix(f1, f2, bands) nb length(bands); M zeros(nb, nb); for n 1:nb for m 1:nb u1 f1.u{bands(n)}; u2 f2.u{bands(m)}; M(n, m) sum(conj(u1(:)) .* u2(:)); end end % M是方阵直接作为Wilson线算符 % 必要时可以做逆矩阵修正M / (M*M)^0.5 U M; end如果计算中curly的M矩阵不等于酉矩阵需要做极分解修正。我自己的经验是当场数据归一化且平面波近似时M矩阵都会很接近酉矩阵但严格起见可以加上奇异值分解修正[Uq, ~, Vq] svd(M); U Uq * Vq;这一步能有效消除数值噪声特别是网格较粗时。4.4 最终陈数计算与结果验证直接算出来的C未必是漂亮整数。比如你可能得到0.934或者1.06这是k网格离散误差的正常表现加密网格会逼近整数。我自己通常画一条“网格加密曲线”把k网格从10×10逐步增加到50×50观察C怎么收敛。如果C稳定收敛到1或者−1就可以放心如果C在某个网格密度下反复跳变大概率是能带交叉或者相位修正出了问题。另外一个重要的验证手段是和已知对称性对照在rArB的极限带隙关闭陈数失去定义在rA和rB的关系相反时陈数应当变号。用Matlab跑一遍参数扫描画出“rA/rB − C”曲线如果看到清晰的平台区比如比值大于1的地带一直保持C1小于1的地带C−1这个结果就和物理图像完全吻合。4.5 脚本中两个容易写错的细节第一个细节是布洛赫相位修正。Comsol的Floquet边界条件解出来的物理场是完整的布洛赫函数E(r) exp(ik·r)u(r)但Wilson loop中需要的是周期部分u(r)。所以从Comsol导出Ez后必须做一次换算u E·exp(−ik·r)。这里的正负号以上面提到的自检结果为准。如果在Matlab里把这个符号搞反重叠矩阵会多出一个随位移指数变化的相位因子直接导致Wilson loop相位不闭合陈数算出个非整数。第二个细节是能带排序。参数扫描时Comsol输出的特征频率是按大小排列的但相邻k点之间的模式身份可能互换尤其在高对称点附近。做Wilson loop时必须确保“第n条带”在相邻k点上是同一条物理能带否则重叠矩阵的行列会错位。解决的办法就是前面提到的用上一个k点的模式与当前k点所有模式做重叠选择重叠最大的索引作为当前能带顺序。我通常在load_fields之后就完成按顺序重排再送进Wilson loop避免在循环里反复处理。5. 踩坑实录假带隙、边界模式混叠和陈数符号问题5.1 狄拉克锥没有打开先查网格对称和边界条件如果你设置了rArB理论预期K点应该是狄拉克锥可是算出来的能带在K点附近有明显分裂这就是一个典型警告信号几何或者网格对称性被破坏了。第一个可查的项目是单胞的几何坐标确认A柱圆心在原点、B柱圆心在准确的(0, a/√3)不要有小数点误差。第二个可查的项目是网格剖分。我遇到过最隐蔽的一次问题出在周期性结对边界网格虽然节点数一致但两对边界的“源/目标”配对方向反了导致周期性条件施加的是反对称关系狄拉克锥被强行打开成一条假带隙。排查方法很简单——把网格隐藏掉单独看能带如果带隙依然存在就关掉周期条件改做镜像模式或者查看边界两侧的场确认相位关系。这种问题纯靠看频率数值很难发现关键时刻还是要看场分布。5.2 布里渊区路径端点不连续模式互换与能带穿越能带图上最让人头疼的现象是某条能带在跨过高对称点后“跳”到了另一条颜色对应的曲线上。这通常不是物理问题而是特征频率排序造成的。在K点附近群速度为零模式简并度高微小数值误差会让求解器把模式顺序打乱。处理思路有两条。第一条是把s步长加密比如从0.02改成0.005让每个模式在相邻两点间的频率差足够小排序更稳定。第二条是在后处理时按重叠积分重排方法在4.5节里已经说了。如果项目只关心带隙位置和陈数其实能带图横轴上的这种跳跃不影响最终结果但如果你要画漂亮的能带图给论文用就必须把重排逻辑加进去。5.3 陈数算出0.3或−0.7这不是非整数是网格和相位问题当你第一次跑通整个脚本算出来的陈数往往是0.31、−0.72这种不尴不尬的值。有时还会出现正负号漂移。出现这个问题的原因几乎总是以下三个之一布洛赫相位修正的正负号错了k网格太粗Berry曲率集中在K/K′尖点附近没有被充分采样场数据没有统一网格导致重叠积分误差大。我的建议是先花半天时间用紧束缚模型验证Matlab脚本本身没有Bug。写一个2×2或者4×4的紧束缚哈密顿量把本征向量按解析公式算好直接喂给同一个Wilson loop函数如果算出来结果与解析陈数一致那脚本就没问题剩下的锅全甩给Comsol导出和相位修正。这个方法帮我节省了大量排查时间。5.4 与其他结果交叉验证边缘态和群速度方向即使陈数算出了漂亮的整数我也建议做一步额外的验证计算超胞的边缘态能带。把单胞改成条带结构沿一个方向有限另一个方向仍用周期边界在同一个Comsol模型里再算一次特征频率。如果体态带隙内确实出现了跨过整个带隙的边缘态色散曲线而且边缘态在同一边界的两个方向的群速度方向相反对应不同手性那陈数结果基本可以盖章确认。这一步还有个额外好处可以直接观察边缘态场分布来判断拓扑保护是否成立。比如在带隙频率范围内激励一个源看波是否能绕过缺陷转角继续单向传播。我在实际项目中就是靠这个最终验证了设计的拓扑波导整个过程比单纯相信一个整数踏实得多。最后分享一点个人的体会能带计算和陈数计算虽然是两套东西但它们的错误模式高度耦合——能带有一点点不干净陈数就会显现成莫名其妙的数值。所以做这个项目时不要急着冲到最后一步先把“对称状态下狄拉克锥正确关闭”这个最基本的结果反复算稳、看稳再谈拓扑不变量。把这套流程走通之后以后再换三角晶格、Kagome晶格或者其他复杂体系框架完全可以直接复用只是换参数、换波矢路径而已。
