COMSOL光子晶体谷霍尔效应仿真：从单胞到单向传输全流程

做光子晶体的同学，今年基本绕不开一个词：光子晶体谷霍尔效应。简单说，就是利用晶格对称性破缺，在动量空间的K和K谷处打开带隙，然后在两个拓扑性质不同的光子晶体边界上，激发出一种由谷自由度锁定传播方向的光学边界态。我最近用COMSOL完整跑了一遍从单胞、超胞到谷单向传输的全链路仿真，从能带计算到实空间电场分布，一步步把现象做出来了。整个过程踩了不少坑，也攒下很多设置细节，写出来供同样想复现谷霍尔效应的朋友参考。这篇文章适合刚开始接触拓扑光子学、想用COMSOL做能带计算的读者，也适合已经算过能带、想继续做超胞和传输仿真的人。

1. 光子晶体谷霍尔效应的核心逻辑与仿真路线

1.1 从电子体系的谷自由度到光子晶体

谷自由度的概念最初来自电子体系。在石墨烯这类蜂窝状晶格中，能带在动量空间第一布里渊区的一对角顶点上形成不等价的能量极值，习惯上叫K和K'谷。电子在不同谷附近的行为不完全一样，这种额外自由度就被称为谷自由度。后来研究进一步发现，如果打破空间反演对称性，比如让晶格中两个子格位置处的原子半径不同，原来的狄拉克锥会打开带隙，两个谷里会出现符号相反的贝利曲率，这就是谷霍尔效应的起源。

光子晶体是由不同介电材料周期性排列构成的人工结构，电子能带理论里的许多概念，包括能带、狄拉克锥、拓扑相，都可以平移过来。区别在于，光子晶体中不需要真实的自旋和电子，只需要调控光频段的艾克斯韦场模式。只要让光子晶体的单胞同样具有蜂窝类晶格结构，再通过几何尺寸不对称打破反演对称性，就能得到K和K'两个谷，并在两谷附近获得方向锁定的边界态。

为什么大家都愿意用光子晶体做谷霍尔效应？一个很实际的原因是实验门槛低。微波波段加工大晶格常数样品非常容易，近红外波段也有成熟的介质柱或空气孔工艺。再加上电磁场分布可以直接用探头测、可以直接看仿真云图，观察拓扑边界态比电子体系直观很多。另一个原因是光子晶体几何灵活，单胞形貌、半径比、晶格常数都可调，拓扑相变往往只依赖某个几何参数是否穿过临界点，非常方便参数扫描。

1.2 为什么用COMSOL做这套仿真

很多做能带计算的人习惯用平面波展开法、紧束缚法或者现成的光子晶体计算程序。这些方法在简单周期性结构上确实快，但一旦涉及超胞、异质界面、有限器件和场分布可视化，修改几何和边界条件会非常绕。COMSOL的优势是几何建模自由度高，能直接从单胞切换到超胞再到有限结构，不需要换工具；电磁波频域接口天然支持周期边界、PML和特征频率扫描，后续还能做参数化优化。

缺点也很明显：特征频率扫描慢，模式排序容易错，网格设置不当会得到假的带隙。所以我的建议是先用数学模型把预期频率估算出来，再用COMSOL做确认。这篇文章里的所有截图和数值，都是基于一个完全虚构的模拟项目X，重点关注可复现的设置流程，而不是具体某个实验数据。

1.3 从单胞到传输仿真的整体路线

我采用的标准流程分三步：

• 第一步，单胞能带计算。通过改变两个不等价格点位置处的介质柱半径，打破空间反演对称性，在K和K'点打开带隙，确认谷态频率区间。
• 第二步，超胞色散计算。将两种半径参数相反的单胞拼接成含有界面的超胞，沿界面方向扫描波矢，寻找能隙内部局域在界面上的模式，也就是谷边界态。
• 第三步，有限器件传输仿真。在由两种拓扑相构成的大尺寸结构中放一个点源，观察界面上的电磁波是否单向传播，并定量评估单向传输比。

这给流程每一步都有独立验证指标：单胞看能带是否打开带隙，超胞看色散分支是否穿过带隙，传输仿真看电场云图是否只有一侧有波。三步环环相扣，缺一步都容易让后面现象解释不清。

2. 单胞设计与能带计算：找到谷态位置

2.1 几何结构与材料参数选择

我用的模型是二维三角晶格中的介质柱。晶格常数取 a 作为归一化单位，背景是空气，介质柱折射率按硅取3.4。三角晶格每个原胞包含两个不等价子格位置，我把两个柱子的半径分别记为r1和r2。当 r1 = r2 时，结构具有空间反演对称性，能带在K点附近会形成狄拉克锥；当 r1 不等于 r2 时，空间反演对称性被破坏，狄拉克锥打开，K和K'谷处的模式频率发生分裂。

这个选择比改变折射率更推荐，因为半径是纯几何参数，在COMSOL里做成扫参变量非常方便，后续切换拓扑相也只需要交换半径值。我自己用的典型参数是 a=1，r1=0.45a，r2=0.28a。这个参数组合在介质柱折射率3.4、背景空气的条件下，禁带中心落在归一化频率0.26c/a附近，比较适合做边界态传输。

为什么要用三角晶格而不是正方晶格？因为三角晶格的第一布里渊区是正六边形，存在不等价的两个角顶点，也就是K和K'谷。正方晶格没有这种天然谷结构，硬凑出来的边界态物理图像会差很多。如果你只是想快速看现象，也可以先不管具体频率绝对值，所有长度单位都归一化到a，这样结果可以任意缩放。

2.2 COMSOL单胞建模与Floquet边界条件

在COMSOL中，我使用的是电磁波、频域接口，二维模型，选择面外电场模式，也就是电场极化方向沿柱轴方向。介质柱体系分析这种极化模式最直观，带隙也容易打开。

几何构建要让原胞是一个平行四边形。两个布拉维基矢建议定义为 a1=(a,0)，a2=(a/2, sqrt(3)a/2)。两个不等价子格位置可取 (0,0) 和 (a/3, a/3)，分别放置半径 r1 和 r2 的介质柱。平行四边形四个顶点是原点、a1、a2、a1+a2。这样构造出来的单胞满足三角晶格周期性和空间群条件。

边界条件是最容易出错的地方。必须给左右一对边界、上下另一对边界分别指定Floquet周期性边界条件，并输入布洛赫波矢对应的相位偏移。不要用普通的周期性边界条件，那是等同周期边界，相位差固定为零，算不出来能带。在COMSOL的特征频率研究里，我可以把波矢分量定义为参数 kx、ky，然后通过参数化扫描来跑能带路径。注意COMSOL中Floquet条件需要的相位输入通常以弧度为单位，需要换算成 k·R 的形式。

最简单的做法是在全局参数里定义 kx、ky，扫描路径取三角晶格第一布里渊区的边界折线：Γ–M–K–Γ。如果用参数扫描去扫连续的ky/kx组合，需要先把折线路径离散成点，然后用辅助扫描变量代替。我习惯用“参数化扫描”中的列表，把路径上的25个左右采样点列出来，保证每条能带分支连续。

2.3 能带图判读与谷态验证

扫描完特征频率后，把频率结果用 c/a 归一化，也就是无量纲的 f·a/c。当 r1=0.45a、r2=0.28a 时，原来的狄拉克锥消失，K点附近会出现两条分裂的带，频率差就是谷带隙宽度。典型值约在 0.02c/a 到 0.04c/a 之间，具体取决于半径差大小。

我一般会在K点停下仔细看模式场图。谷模式有一个非常明显的特征：动量空间K点附近的简并模式，其面外电场呈涡旋状，一个顺时针，一个逆时针，对应不同的谷赝自旋。如果你在COMSOL里画箭号图或等值线图，能看到K点的两个模式分别拥有相反的相位绕数。这是确认谷态“确实存在”而不是随便一条平板模式的最好方法。

还有一件事要提前做：确定禁带两侧的本征频率，记下来作为后面超胞分析的频率基准。如果色散曲线在某段出现锯齿状跳变，通常是网格太粗或者特征值搜索范围太窄，需要先调整网格和特征频率基准，不要急着读数据。

3. 超胞构建与边界态色散：界面产生在哪里

3.1 为什么单胞算完还不够

单胞能带只能告诉你体态在哪里开带隙，但拓扑边界态不是单胞几何本身能直接解出来的。谷霍尔边界态必须出现在两种拓扑性质不同的晶格界面上，所以要人为构造一个异质界面。把两种单胞拼在一起，再沿界面切向施加周期边界条件，就构成了超胞问题。超胞里包含所有可能的体态和界面态，通过比较能隙内的模式，就能确认边界态频率和群速度方向。

这种做法的好处是既能拿到色散关系，又能通过场图判断模式是否局域在界面。比直接做有限尺寸传输仿真更省资源，也更容易解释物理。

3.2 两种拓扑相的单胞拼接方式

两种拓扑相的具体定义很简单：A类单胞取 r1=0.45a、r2=0.28a，B类单胞取 r1=0.28a、r2=0.45a。它们的体带结构几乎完全相同，但K和K'谷的模式频率次序发生交换，所以拓扑性相反，拼接后才能在界面上出现谷边界态。

在COMSOL里，我建议先把A类和B类分别建好，然后用数组功能沿界面方向复制若干周期，再在界面法向堆叠若干排。界面沿x方向保留平移对称性，超胞沿x方向长度为一个或几个界面原胞周期，沿y方向包含上下各6到10排原胞。

这里有个关键操作技巧：复制出来的多个域一定要执行合并或连续性设置，否则求解器默认把它们当成互相绝缘的独立域，算出来的模式完全不正确。在几何装配模式下，用“创建并集”把这些域合成一个整体，最省事。如果因为不同圆柱半径导致界面两侧网格节点不匹配，需要开启连续性条件，让场跨越内部边界时自动满足连续条件。

界面构型我采用常见的zigzag型界面，超胞沿x方向的周期长度取 sqrt(3)*a。界面结构算出的边界态色散会落在禁带里，具体的群速度符号取决于界面的具体微观形貌。你可以在扫描后直接看场分布来验证，不用纠结于某种特定界面是否“正宗”。

3.3 超胞色散计算与模式判别

超胞仿真仍然用特征频率研究，但只需要沿界面方向扫描一个波矢分量 kx。上下两个法向边界可以用散射边界条件或PML，但要注意界面态必须局域在内部异质界面附近，不能受上下外边界影响。如果外边界产生表面态，很容易和需要找的边界态混淆。我通常把上下块体厚度增加到8排原胞以上，这样目标模式的场在到达外边界前已经衰减到很小，外边界条件对结果影响可以忽略。

扫描完kx之后，将频率画成色散曲线，保留带隙内的分支。这里有个重点：不是所有带隙内的模式都是谷边界态。普通表面态也可能落在带隙内，但它们的场局域在超胞的外边界，而不是结构拼接界面。所以每次算完一个模式，我都会把模式场图画出来，检查最大值位置是否在异质界面附近，并且沿法向两侧都呈现衰减趋势。

谷边界态还有一个标志性特征：色散曲线在带隙内部并非单调只有一支，而是会出现正斜率和负斜率两个分支。正斜率分支在kx正方向携带一种谷自由度，负斜率分支携带另一种。这两条分支相交于kx=0附近，但如果时间反演对称性保持着，它们的频率会关于kx=0对称。看到这样的结构，基本可以确定界面态是谷型的。

4. 谷单向传输仿真：用场图验证方向锁定

4.1 构建有限尺寸器件模型

超胞给出边界态频率后，下一步是把它放到一个有限尺寸的真实器件里看传输。我一般建一个长方形计算区域，左右两侧设置PML吸收边界，中间是拓扑相反的两种光子晶体，界面放在水平方向。上下两个方向保留一定厚度的完整光子晶体，但不能太薄，否则波会从外边界漏出去。

器件宽度不能太小。如果左右各只有3排原胞，边界态与上下体态的耦合会很强，单向性会被破坏。我试过用左右各8排原胞，计算时间大概是小结构的4倍，但结果干净很多。PML区厚度至少一个波长，并且PML与真实结构之间要留一小段缓冲区域，不能直接把PML贴在结构边界上。

网格设置在这一步特别关键。由于介质柱折射率高达3.4，材料内部波长是空气波长的1/3.4，如果网格统一按空气中波长剖分，柱体内场会严重欠采样。我通常对介质柱内部单独指定最大网格尺寸，取 f=0.26c/a 时硅内波长约 a/0.26/3.4 ≈ 1.13a，网格最大尺寸取0.2a左右就够；如果想更保险，取0.1a，不过计算时间会增加很多。

4.2 激励源与单向性验证流程

激励源我用一个点偶极子，等效于电流源，放置在界面附近，坐标为界面中心稍稍偏向一侧。频率取超胞色散中带隙中心的频率，也就是能谷边界态群速度最大、体态分量最少的频率。频域研究直接求解，观察电场模分布。

判定单向传播不能只看云图。我的标准流程是：

• 在源左右两侧对称位置各放一个点探针，记录电场模平方；
• 计算左右探针能量比，用10log10换算成dB；
• 改变激励频率到带隙外，确认左右两侧信号差异显著缩小；
• 交换源的位置或改变源的激发相位，观察主传输方向是否反转。

我实测下来，合理的结构参数能轻松做到一侧比另一侧高20dB以上。如果只能差到5dB，多半是频率没对准边界态带中心，或者源离界面太远。这时候不要急着改结构，先把源位置朝界面移动，再微调频率。如果源直接放在界面正中间，往往因为两个谷同时被对称激活，导致两端都有波，这是正常的对称选择问题，不代表体系失效。

4.3 谷锁定传输容易被误解的地方

要特别强调：谷霍尔边界态不是传统意义上的“绝对单向模式”。体系没破时间反演对称性，所以K和K'两个谷的边界态方向相反。点源之所以能激发出单方向传输，是因为源的空间位置和频率与某个谷的模式匹配，相当于做了一次谷选择性激发。如果你用均匀面光源或大尺度宽频光源，就可能同时激发两个谷，看到双向传播，这不是仿真出错了，而是物理上应得的。

另一个容易踩的坑是边界粗糙度。实际结构的截断边界不一定刚好处在理想的谷匹配构型上，可能出现微弱的反向散射。要减小这种效应，材料柱边界网格要精细，PML距离要足够，界面沿长度方向不能有随机缺失柱。随机缺失几个柱子就够破坏相位匹配，单向传输比会跳水，这是相当合理的结果。

5. COMSOL参数设置细节与常见问题排查

5.1 特征频率扫描的模式对齐技巧

做单胞和超胞能带时，最大的坑是模式排序混乱。COMSOL特征频率求解器按本征值大小输出，但如果扫描参数变化较大，某个模式的频率可能从第3阶跳到第5阶，画能带图就会得到交叉错位。我常用的办法是：

• 搜索基准频率设为预期带的中心，比如0.26c/a；
• 扫描步长尽量小，特别是在K点附近；
• 关闭“自动排序”，改用“特征值跟踪”功能，或者在参数扫描设置里勾选“上一步解”作为初始猜测。

如果这一步不做好，后面所有色散图都是乱的，超胞边界态也找不到。我在第一次扫描时就吃过亏，当时把kx步长放得很大，从0扫到π/a只取了8个点，能带曲线里K点附近分支完全看不清，还以为是网格问题。后来改成40个点，一下就正常了。

5.2 网格收敛性检验的实操经验

网格密度直接决定能不能算出可靠带隙。我的检验方法是：先粗网格跑一版，记带隙上下沿频率；然后把全域最大单元尺寸减半，再跑一版；如果带隙上下沿变化小于1%，就认为收敛了。不要跳到加密就已经开始处理超胞，那样算得慢不说，还可能把数值误差误认为物理带隙。

介质柱边界是曲率比较大的区域，建议在圆柱表面加局部最大单元尺寸，而不是只依赖全局网格。我用 r=0.45a 的柱体时，表面最大单元尺寸取0.08a，内部用自由三角网格，背景区域可以稍粗一些。超胞尺寸大时，要善用“物理场控制网格”和“自适应网格细化”，把计算资源集中在界面附近，否则一个超胞可能几十万个自由度，跑一趟要等很久。

5.3 常见报错与排查速查表

下面是我在实际中遇到的几个典型问题，整理成了一张速查表。

这个表格不是绝对标准，但按顺序排查，绝大多数问题都能收敛到正确方向。

5.4 半径差参数扫描：确定最佳单向传输窗口

半径差 r1-r2 是控制带隙宽度和边界态带宽的主要参数。我做过一组扫描：固定 a=1，背景折射率1，柱体3.4，取 r1=0.45a，r2从0.3a步进到0.42a，记录K点带隙中心和边界态频率窗口。结果是半径差越大，带隙越宽，但边界态色散与体态距离也在变化。最佳工作点往往不是带隙最宽处，而是在带隙中心频率附近，边界态群速度较大、与上下体态色散最近的一侧。

所以在做传输之前，我建议先花半小时做半径差扫描，画出带隙宽度和边界态频率范围的关系，而不是直接选一个好看的结构算传输。用“参数化扫描+特征频率”组合，可以一次把多组半径差的能带都跑出来。跑完以后，后续有限结构传输的频率选择就有明确依据了。

6. 一点仿真经验与可扩展方向

我印象最深的坑，是在超胞阶段没做模式场判别，结果把外边界表面态当成边界态，后面传输仿真怎么调都出不来单向性。后来老老实实每个候选模式都看场分布，才发现问题出在上下面边界条件上。加厚块体、改用散射边界之后，边界态分支才变得清楚。所以我的经验是，拓扑仿真不能只看色散曲线，模式和场图必须对照着看，这是判断模式身份最可靠的依据。

另外一个小技巧是，所有频率都用 c/a 归一化。C OMSOL默认单位是米，如果直接填频率到10的14次方，数值求解会很别扭，还容易出现数值截断误差。用归一化单位长度把a设为1，频率用无量纲量表示，再在需要时换算成真实尺寸，比如要设计近红外波段，就把a设为500nm，所有几何参数按比例乘上去即可。

做完二维模型后，还可以往三个方向扩展：一是把平板结构改成有限厚度，加上上下包层，做准三维版本的谷边界波导；二是在边界上加拐角、分束器，验证谷霍尔边界态在弯折处的抗散射能力；三是引入可调介质材料，实时改变半径差参数，做成可重构的谷选择性开关。每一步都还有大量的参数扫描和模式分析要做，但基础就是这次讲的单胞、超胞和传输三步法。只要这三步的流程顺了，后面换成空气孔、其他晶格甚至声学晶体，思路都是完全一致的。