这几年只要做光子晶体几乎绕不开拓扑光子学。前一阵我接手一个验证项目要在 Comsol 里模拟光子晶体谷霍尔效应并且把谷陈数算出来。这个任务光看标题好像不算难真正做起来才发现能带图只是第一步后面陈数怎么算、拓扑相变怎么判断、边界态怎么验证每一步都有不少隐形的坑。这篇博文我把整条流程完整复盘一遍从模型选型、Comsol 参数设置到 Wilson loop 提取谷陈数再到最终的边界态验证全是我实际跑通过的做法希望能给正在做类似仿真的人省点时间。1. 光子晶体谷霍尔效应到底在算什么1.1 谷自由度与 K/K 点从电子体系到光子体系谷霍尔效应这个名字是从电子体系借过来的。石墨烯里布里渊区角落有两个不等价的狄拉克点通常叫 K 和 K电子在 K 附近的动量和能量关系跟 K 附近几乎完全一样只是手性相反。这两个点就叫“谷”相当于给电子多了一个离散的自由度。后来人们发现如果打破某种对称性两个谷会表现出不同的贝里相位电子在谷间会产生方向选择性输运这就是谷霍尔效应的大致图景。光子晶体里完全可以做同类的事。把介电常数周期排布成蜂窝状或者三角晶格横磁TM模式下电场只有 z 分量能带图里同样会出现 K/K 这类谷点也具备线性的色散关系。既然方程组在数学上同构电子体系里的谷自由度、谷陈数这些概念就能搬到光子体系里。所以做这个项目第一件事不是急着开 Comsol而是要想清楚我们准备用哪一对谷点做文章以及用什么手段把谷简并解开。这里要强调一个容易混淆的地方很多文章里说的“谷”不是能量本征态的“谷”而是动量空间里的特殊高对称点。光子晶体能带图里K 点处有三重旋转对称和镜像对称如果体系保持完整的空间反演对称狄拉克点会稳定存在。只有把反演对称破坏掉原本简并的狄拉克锥才会劈裂出带隙形成真正可用的“谷带隙”。这个带隙的大小和位置直接决定了后面所有拓扑性质是否能被观察到。1.2 反演对称破缺把“质量项”加进光子晶体理论框架里K 点附近的有效哈密顿量可以写成一个二维狄拉克哈密顿量反演对称对应的就是哈密顿量里的质量项 m。m 0 时狄拉克锥不打开m ≠ 0 时带隙劈开符号不同对应两种不同的拓扑相。在真实结构里最常用的做法是让元胞内两个相同尺寸的散射体变成一大一小比如三角晶格里两个圆孔半径不同或者让六边形单元胞里的两个柱体高度不同。我在实际仿真里用的是一套比较经典的结构三角晶格每个元胞里有两个空气孔其中一个孔半径是 R1另一个是 R2。当 R1 R2 时体系具备反演对称性K/K 点保持狄拉克简并当 R1 ≠ R2 时反演对称性被破坏能带在 K 和 K 打开带隙。这个模型的优点是几何构造简单、Comsol 建模容易、参数变化时能带演化很直观非常适合用来复现谷霍尔效应。这里有一个选择逻辑为什么不用最简单的正方晶格因为正方晶格没有六重/三重旋转对称K 谷和 K 谷的概念不自然而且能带简并的拓扑保护不够稳健。三角晶格天然有三重旋转对称K 和 K 在布里渊区对称位置时间反演把它们互相关联这套对称性正是谷物理的基础。后面计算谷陈数时对称性会直接体现在数值结果里。1.3 谷陈数为什么是 ±1/2以及算它有什么用对于连续二维狄拉克模型K 谷附近的贝里相位是 π对应的谷陈数是 ±1/2。这里的 1/2 不是错误也不是分辨率不够它本身就是单谷贡献。整个体系如果保留时间反演对称两个谷的陈数加起来必须为零所以一个谷是 1/2另一个就是 -1/2。看上去很“碎”但它足以决定边界态是否存在在谷畴壁两侧如果谷陈数从 1/2 变成 -1/2体边对应关系就要求界面处存在手性的边界模式这可以用超胞模型或者有限大结构直接验证。算谷陈数的主要意义在于它比单纯看能带图更可靠。能带打开带隙只能说明存在带隙不能说明带隙是否具有拓扑非平凡性。两个参数不同但都打开带隙的结构可能一个是普通绝缘体一个是谷拓扑绝缘体单靠“哪些模在带隙里”看不出来。陈数把波函数的相位信息浓缩成一个整数或半整数告诉我这个带隙到底“拓扑”在哪。这也是为什么这个项目不能止步于能带图而一定要把陈数算出来。2. Comsol 建模仿真与能带计算的四个关键步骤2.1 结构怎么选能稳定出结果的模型参数集我用的参数集如下晶格常数 a 1 μm空气孔背景介质取硅相对介电常数 ε_r 12.25。元胞内两个空气孔分布在元胞中心两侧初始半径 R1 0.30a、R2 0.20a注意二者所在位置的局部坐标保持镜像对称即可。后面扫参时固定 R1改变 R2就可以观察能带演化。COMSOL 建模时建议直接在“全局定义”里把参数全部写清楚后面所有几何尺寸、扫描变量都引用这些参数。这样做的好处是参数化扫描时只需要改 R2 一个量不用手动重建几何。Comsol 版本差异会影响具体菜单名称但思路一致先建一个二维元胞几何在元胞边界上加 Floquet 周期边界条件然后扫波矢得到本征频率。网格部分二维模型自由度不大但要注意三角形空气孔或圆形孔边界处网格要细分一点。我一般用“自由三角形网格”最大单元尺寸设为 a/30 左右孔边界处加一个“尺寸”节点最大单元尺寸设为最小孔半径的 1/5。这个网格密度跑能带图完全够算边界态图也不会慢到无法接受。如果开 4 核并行一个扫描点大概几十秒到几分钟整条能带路径 30 个点大概一小时上下属于可以接受的量级。2.2 Floquet 周期边界条件布洛赫波矢的设置细节光子晶体能带计算的核心是布洛赫定理周期性体系中电场满足 E(r a) E(r) e^{ik·a}。Comsol 里对应的功能是“Floquet 周期边界条件”使用时要给出一组波矢分量例如 kx kx0ky ky0然后计算对应布洛赫波矢下的本征模式。实际操作中我会在“全局定义”里定义两个扫描变量 sx、sy分别作为 kx 的系数和 ky 的系数然后让 Floquet 周期条件里的波矢等于 kx sx * 2pi/a、ky sy * 2pi/a。这里容易犯的错误是把波矢直接定义成“无量纲的 0 到 1”却忘记乘 2π/a导致算出来的频率和文献对不上。我建议统一用物理波矢单位 rad/μm这样带隙频率单位是 THz和实验也容易对应。设置周期边界时还有一个细节Comsol 的 Floquet 条件需要把相对的边界匹配起来源边界和目标边界的选择会影响相位方向。如果不小心把一对边连反了结果会出现本征频率全变复数或者能带形状不对。所以建议每次建模后先扫一个高对称点比如 Γ 点验证一下前几个本征频率是否可靠。Γ 点布洛赫波矢为零周期条件退化成普通周期边界最容易排查问题。2.3 特征值研究与参数化扫描把能带“扫”出来能带计算本质上是一个特征值问题给定波矢 k求满足布洛赫边界条件的本征频率 f。Comsol 里我建议用“特征值研究”物理场选“电磁波频域”研究类型选“特征值”。这里的关键是特征值搜索基准如果研究的频率范围在 0.3 THz 到 0.8 THz搜索基准就填一个中间值比如 0.5 THz然后设置“所需特征值数量”我一般取 6 或 8保证在目标频率附近能拿到足够多的能带分支。参数化扫描设在“研究”里的“参数化扫描”节点扫描参数就是定义好的 sx、sy或者直接扫描路径参数 s。为了让能带图方便画我通常会直接定义一条动量空间路径比如沿着 Γ-M-K-Γ 走把路径长度参数化成一个参数 s 从 0 到 1在“全局定义”里用分段函数把 sx、sy 写成关于 s 的表达式。之后研究里只需要扫描 s一次就把整条高对称路径扫完后处理也方便。需要注意的一点是扫描参数如果跨越高对称点比如经过 Kk 点会经过简并点附近特征值求解器偶会报“缺少特征值”或者特征频率丢失。这不是模型错了而是扫描点太靠近简并点或者特征值数量不够。我的处理办法是把路径上 K 点附近加密比如 K 点前后各多放 3~5 个点同时把特征值数量从 6 改成 10基本能解决。2.4 能带图后处理不要只靠 Comsol 画图Comsol 自带的后处理可以画能带图但它把每个参数点单独作为一组解直接画出来的图是散点而且不同分支不会自动连接看起来像一堆乱点。我用的是导出数据到本地再用 Python 画的流程。在“结果”里新建“一维绘图组”用“全局”绘图把纵轴设为特征频率值比如使用变量 f 或 “ewfd.freq”横轴用参数 s。如果你只是快速预览可以选所有扫描解一次性显示所有本征频率点。对于正式论文或报告建议在“派生值”里选择“全局计算”把每个参数点的所有特征频率导出成表格再用 Python 或 Matlab 重画成连续能带图。绘制能带时可以把频率归一化成 a/λ也就是纵轴用 a*f/c这样与文献对比时不用管具体晶格常数。我对比过不同 a 值下的结果归一化后几乎完全重合这也可以用来验证建模有没有单位错误。3. 从能带图到谷陈数Wilson loop 与相位提取3.1 为什么用 Wilson loop 而不是直接积分贝里曲率教科书中陈数被定义为贝里曲率在动量空间某个闭合流形上的积分。对二维体系把整个第一布里渊区当成一个环面陈数就是贝里曲率在环面上的积分除以 2π。理论上可以直接在每个 k 点算贝里曲率再积分但数值上非常麻烦因为贝里曲率的表达式包含波函数对 k 的导数而 Comsol 给出的本征场是离散点在固定网格上的值求导会引入很大误差。Wilson loop 是更稳的做法。它在动量空间的一个闭合路径上把相邻 k 点的波函数内积乘起来得到一个幺正矩阵的本征相位谱这些相位谱的绕数winding number就等于陈数。对应到谷体系我们可以绕着 K 或 K 画一个小圆环圆环上的贝里相位除了 2π就能得到谷陈数 ±1/2。这个方法不需要显式求导只需要本征函数在离散 k 点上的值因此特别适合和 Comsol 这类有限元工具配合。绕圈半径的选取也要注意。理论模型里有效哈密顿量在 K 点附近是线性的所以圈半径越小越接近理想的 ±1/2但数值上圈太小波函数差异太小内积矩阵接近单位阵相位噪声会变大。我的经验是绕圈路径半径取布里渊区尺寸的 5%~10%比如在倒空间中 r 0.05 * (2π/a) 左右再配合足够密集的采样点结果会很稳定。3.2 从 Comsol 导出本征场与构建 Wilson loop 的流程既然有了绕圈方案后面就是数据提取。我以 K 谷为例给出完整流程。第一步定义一条闭合路径。在倒空间中以 K 点为中心取 N 个点把角度等分比如 N 36每个 k 点就是一个扫描参数记为 path_i。第二步在 Comsol 里把这 N 个 k 点按顺序扫一遍每个点保留目标能带分支的本征场。这里需要你对能带图已经比较清楚知道带隙上沿或下沿哪条分支是我们关心的。我建议在带隙上方的第一支导带和下方第一支价带分别做 Wilson loop对应不同的谷陈数。第三步在结果里导出该本征模的电场 z 分量复数分布。导出时最好固定一套网格点的坐标这样不同 k 点之间的内积才有意义。Comsol 数据导出支持“网格”导出选“所有域”导出变量选“实部(Ez)”和“虚部(Ez)”得到每个网格点的坐标和复数值。这里要注意导出的是包含布洛赫相位因子的总场还是周期部分取决于 Comsol 内部定义。我的做法是在表达式里手动写成周期部分 u Ez / exp(i*(kxx kyy))再导出 u 的实部和虚部确保得到的是布洛赫函数的周期包络。第四步用 Python 读取这些文件。每个 k 点对应一个文件里面是网格点坐标 (x, y, re_u, im_u)。将复数场 u_k(x, y) 存入数组。然后计算相邻 k 点之间的内积矩阵M_{mn} ∫ u_m^*(r) u_n(r) dr积分通过在网格点上求和近似注意要乘上每个网格点对应的面积权重对于非均匀网格建议用非均匀权重不然结果会偏向网格密集的区域。最简单的办法是把场差值到均匀网格上再积分精度够而且省事。第五步把所有相邻内积乘成一个矩阵乘积然后求这个矩阵乘积的行列式的幅角除以 2π 就得到该闭合路径对应的陈数。用代码表示大概是import numpy as np # U_list: list of complex fields for each k, shape (N_points, N_grid) U [u0, u1, u2, ..., uN_minus_1] M [] for i in range(N): M_i U[(i1) % N].conj().T U[i] # 注意顺序 M.append(M_i) # Wilon loop matrix W np.eye(N_bands) for mat in M: W mat W # 取行列式相位 phase np.angle(np.linalg.det(W)) chern phase / (2 * np.pi)这段代码我没有把归一化写进去实际使用中每个场都要归一化内积矩阵也要做正交化处理可以用 Gram-Schmidt 或者极分解不然幅角会受模长影响。另外路径闭合后最后一个点和第一个点是同一个 k也就是 k_N k_0这一步才是“闭合”的关键千万别漏掉。3.3 三个必须处理干净的数值细节第一个细节是相位对齐。Comsol 每次独立求解时特征模式的整体相位是完全随机的同一模式在不同 k 点可能差了一个任意相位 e^{iθ}。Wilson loop 本身对这种“规范相位”其实是协变的不会影响行列式的总幅角但如果波函数没有归一化或者内积矩阵不正交随机相位会放大数值误差。所以我在每个 k 点都会先对波函数做归一化然后用“平行输运”思路把相邻 k 点的相位对齐也就是让内积 u_k|u_{k1} 的幅角尽可能接近 0 或连续变化。这个步骤不是理论必须但对数值稳定性帮助很大。第二个细节是模式排序。Comsol 输出的本征模不是按我们关心的能带分支排好的同一 k 点附近可能出现频率很近的多个模式如果挑选分支时不仔细Wilson loop 会混入其他带算出来的相位谱就是乱的。我通常用频率排序先初步判断是哪一支再通过场分布的空间对称性做二次确认。带隙上下两支的场对称性差异在 K 谷附近很明显一支像偶极子一支像四极子。第三个细节是导出的网格点顺序。Comsol 导出的网格点坐标不是按照空间顺序排列的不同 k 点文件里网格点顺序也不保证相同。这在做内积时是致命问题。我建议在 Python 里先按 (x, y) 排序再构造网格点列表如果你用“均匀网格重采样”这个问题会自动解决。之前我就是没在意顺序算出来的陈数一会儿 0.2 一会儿 0.8整整排查了两天才发现是网格点错位这个坑相当隐蔽。4. 谷拓扑相变的仿真验证参数扫描、边界态与单向传输4.1 参数扫描看带隙闭合-重开的过程谷陈数不是从带隙打开就自动固定的。把 R2 从 0.1a 一直变到 0.4a保持 R1 0.3a能带图会经历一个明显的过程带隙先缩小在 R2 R1 处闭合此时体系恢复反演对称K 点重新变成狄拉克点然后 R2 继续增大带隙再次打开。看似对称实际上两次打开的能带顺序是翻转的谷陈数也差一个符号。具体表现如下R2/R1K点带隙状态导带底与价带顶的空间对称性谷陈数K谷0.5打开导带底为一种涡旋手性1/2或 -1/2取决于定义1.0闭合狄拉克点简并01.5再次打开两种手性互换-1/2或 1/2如果你观察到的“闭合再打开”过程中能带在闭合点之前一直打开着说明可能没有经过对称点比如 R1 和 R2 的差异没有扫到零或者周期边界条件设置有问题导致反演对称一直没有恢复。这个检查非常重要因为“带隙闭合再打开”本身就是体陈数变化的直接证据。4.2 能带反转的直接判据只扫带隙宽度还不够最好从模式的对称性上给出“反转”的直接证据。在 K 谷附近带隙上下两支本征模在元胞内的场分布都有明显的旋转对称性可以用角动量或者镜像对称性来区分。常见的做法是看两个模式在元胞中心附近的电场相位分布一个模式表现为顺时针旋转的涡旋另一个表现为逆时针旋转或者两个模式的镜像对称性相反。R2 大于 R1 和小于 R1 时上下两支的对称性会互换这就是“能带反转”。我用一个更直接的方法验证画某个高对称点的 Ez 场分布计算出它关于镜像面的对称投影或者直接右键“派生值—体积分”积分 Ez 在元胞左右两半的强度差。如果带隙闭合重开后这个量的符号发生翻转说明拓扑相变确实发生了。这个信号比单纯看频率差更容易定量对比。4.3 边界态与单向传输验证谷陈数判据的最终手段算完陈数最好再用边界态做一次交叉验证。方法是在 Comsol 里做一个“谷畴壁”结构左边用 R2 R1 的元胞右边用 R2 R1 的元胞中间形成一条直的界面左右两侧整体仍保持周期性。用频域求解在界面附近放置一个偶极子点源频率设在体带隙内。如果拓扑分析正确会看到场只能沿着界面单向传播反向几乎没有能量。实际操作时我不建议一上来直接建无限大超胞而是先建一个有限大矩形区域左右各 6~8 个周期上下面加“散射边界条件”或 PML中间放点源。然后在带隙频率范围内做频域扫描找到边界态的传输峰。边界态在透射谱上会表现为一个明显的尖峰对应手性边界模式。为了确认单向性可以把点源放在界面中间左右两端各放一个探针比较两端能量。由于边界态是手性的朝一个方向传播的场应该远大于另一个方向。建有限区域时网格和内存仍可控边长在 10a 量级时自由度大概几十万普通工作站都能跑。这一步千万记得用和能带计算相同的材料参数几何否则边界态频率会和体带隙对不上一个是真空/介质界面的表面波一个是谷边界态很难区分。5. 常见问题与排查技巧实录5.1 能带扫出来像乱麻怎么办这是最常见的问题。很多人扫完参数以频率为纵轴把每个参数点上的所有特征频率画出来结果一堆曲线交叉错乱根本看不出带隙。原因通常有三个一是特征值数量太少某些分支在扫描过程中漏掉导致相邻参数点的频率点不能连续连接二是扫描步长过大尤其经过高对称点附近时分支变化很快步长太大看起来就像跳变三是模式顺序在 k 点之间发生了交换你没有手动识别直接连线时就会交叉。我的做法是先用比较粗的扫描找到大致带隙范围然后在带隙附近加密扫描步长并把特征值数量提高到 10 以上。如果某些频率明显“断档”单独取那个参数点打开“特征值求解器”的日志检查求解器是否提示“未收敛”或“特征值缺失”。这一步俗称“在频带图上做标记哪条线缺了就补哪个点”看起来土但很有效。5.2 带隙打不开或者对称性破坏了却还是简并如果 R1 和 R2 已经不同但 K 点带隙始终为零原因很可能是周期边界条件的对称性根本没被破坏。比如两个孔的相对位置选得不对虽然半径不同但元胞里存在一个平移加旋转的对称操作能把 R1 孔变成 R2 孔这相当于体系仍有某种反演对称。这种情况尤其容易出现在“中心对称但两个孔离得不够远”的几何里。另外还有一个隐蔽因素Floquet 周期边界的波矢定义错误。K 点的坐标应该是 (2π/3a, 2π/√3a) 附近如果设置成了 Γ 点或 M 点即便反演破坏带隙也不会出现在预期的位置。建议先画出布里渊区标出 Γ、M、K 三个点然后在扫描路径里显式地把这几个点包含进去不要凭感觉设 k 范围。5.3 内存不足和求解变慢怎么办二维模型内存压力通常不大但如果做带隙扫描时选择的域非常大网格数会快速增加。常见的优化手段有三个一是把最小网格尺寸从 a/30 放宽到 a/20先验证趋势再加密做最终结果二是把特征值数量减少到 6只在必要的范围里保留更多分支三是关闭“自动重新划分网格”当几何形状不变只是扫描 k 时网格完全不需要重画这个设置在 Comsol 里默认会重建网格有时特别浪费时间。如果用了“有限大结构加 PML”的验证模型PML 区域的网格不能太粗否则会产生虚假反射。我会在 PML 区域使用“扫描”网格生成器沿厚度方向划分 5~8 层面内网格和结构内部的边界网格尺寸一致这样既保证吸收效果又不至于让自由度爆炸。5.4 Wilson loop 结果不稳定怎么查Wilson loop 跑出来陈数是 0.3、0.8 这种中间值而不是接近 ±1/2优先检查几件事第一个是导出场时是否除以布洛赫相位因子。Comsol 解的因变量是总场如果你直接用总场做内积等于在每个 k 点额外引入了一个空间相位 e^{ik·r}相邻点内积里会出现一个很大的相位差把 Wilson loop 的相位完全扰乱。第二个是路径采样点是否太少。如果绕 K 点一圈只取 12 个点相位差离散化误差很大我建议至少 24~36 个点。第三个是模式选择是否发生跳变。如果绕圈过程中需要跨过一条能带你选的是“同一序号”的模式但在某个 k 点附近的频率变成了另一条分支Wilson loop 就会混入完全不同的态。这种情况可以通过画出每支模式的频率随角度变化曲线来排查理想状态下应该是光滑连续的曲线。另外每次导出后可以用一个 sanity check在两个完全相同的 k 点分别做两次导出计算两次场的内积结果显示归一化后的模长应该接近 1相位应该接近 0。如果做不到说明数据读入或网格对齐有问题后面算出来的任何陈数都不可信。5.5 若干“经验型”建议最后分享几个项目过程中积累的实用习惯。第一所有关键参数都用全局变量文件名里也带参数名比如R2_020.txt方便回溯。第二每次改动几何参数后第一步先看 Γ 点频率再扫 K 点附近一小段能带不要直接全路径扫描可以快速判断设置有没有错。第三如果同一台机器上要同时跑大量扫描建议把 Comsol 的“参数化扫描”拆分成多个独立任务比如 10 个一组避免单次任务超过内存限制导致整个研究中断。还有一点算陈数的过程虽然看起来独立于 Comsol但我强烈建议把能带图和 Wilson loop 结果放在同一个文档里对比不要分开记。只有当带隙闭合重开、模式对称性反转、谷陈数符号翻转三者全部对上才能比较确定地说这个体系发生了谷拓扑相变。单看其中任何一个都有可能被数值假象骗过去。如果让我重新做这个项目我会用不超过半天的时间先把一个已知文献里的基准结构完整复现一遍确认能带图和谷陈数都对得上再开始改参数、加边界态。这样后面所有新结果都有参照不会在错误模型上越跑越远。这一步多花的时间通常能省下后面好几天瞎折腾的时间。
