用COMSOL算一维光子晶体能带：建模、边界条件与带隙分析

先声明一下，这篇文章不是从教科书里搬概念，而是把我在 COMSOL 里跑一维光子晶体能带的全过程摊开来说，包括中间踩过的坑和反复试错之后留下的经验。一维光子晶体听着唬人，其实说白了就是两句话：把两种折射率不同的介质按周期交替叠起来，光在里面走的时候，某些频率会被周期性界面反射到“进不来”，另一些频率却能正常穿过去。这种“能不能过”的规则，由晶格排列决定，最后落在光子能带图上。本文用 COMSOL 在硅基底上搭一个由硅和二氧化硅交替构成的周期性介电结构，从几何建模、材料参数到 Floquet 周期边界条件，一步步算出一维光子晶体能带，并教你怎么把带隙位置和宽度读出来。适合刚接触光子晶体仿真、正纠结如何用 COMSOL 做周期结构计算的人，也适合已经跑过案例但被能带乱线和边界条件设置折磨过的人。

1. 先想清楚：一维光子晶体能带计算到底在算什么

1.1 用“高速公路收费站”理解带隙的本质

我给朋友讲光子晶体的时候，最喜欢用“收费站”这个比喻。高速公路单向只有一条道，不同频率的光子就像不同型号的车。收费站就是个周期设置的关卡：如果所有关卡都敞开，所有车都能走过去；如果关卡隔一段就设一个，某些速度的车会因为每个关卡都在同一位置配合默契地被拦下，另一些速度的车反而能踩着间隙冲过去。周期性介电结构对光子做的事情完全类似：折射率周期变化，相当于在光路上周期性地设置“折射率关卡”，光每经过一个周期都会在两个界面上产生反射。当某一频率的反射光在各周期界面处相位相同、相互加强，入射光就会被整体弹回去，这个频率就“禁行”；当各周期反射相互抵消，光就畅通无阻地通过。

用术语说，这个“禁行区”就是光子带隙，英文 photonic band gap。一维结构中它也叫布拉格带隙，因为核心物理就是布拉格反射。带隙的位置主要由晶格常数、两种介质的折射率差和厚度比决定。很多人问：是不是只要折射率差大就能出现带隙？不一定。周期、厚度比必须满足一定条件，带隙才会在某个频率范围稳定打开，而且带隙宽度随折射率对比增大而增大。硅对二氧化硅的折射率比接近 3.45 比 1.44，对比度相当可观，用来做演示模型再合适不过。

1.2 能带图在说什么：频率和波矢的“契约”

一维光子晶体能带图横轴是波矢 k，纵轴是频率 ω（或者归一化频率）。均匀介质中，光子频率与波矢是正比关系，ω = ck/n，所以画出来是一条直线。加入周期结构后，色散关系在布里渊区边界处断开，原来那条直线被“劈”成若干段，每一段叫一条能带，能带之间没有实频率解的区域就是带隙。

为什么要限制波矢范围？因为周期结构具有平移对称性，Bloch 定理告诉我们，不同布里渊区的波矢本质上可以折叠到第一布里渊区。一维晶格的第一布里渊区是 -π/a 到 π/a，又因为系统中心对称，通常只画 0 到 π/a 这一段，右边端点在不同文献里叫 X 点或边界点。能带图的价值在于：只要给定光频率，你就能判断这个频率的光在这个周期性介质里能不能传播、传播速度多快（能带斜率对应群速度），带边频率还能用来设计反射镜、滤波器和传感器。说到底，能带计算就是求解一个含周期参数的色散本征问题。

1.3 为什么选 COMSOL，而不是自己写传输矩阵

我见过很多人第一反应是写传输矩阵法，一维堆叠解析公式几分钟就能出结果。那确实快，但它只能处理无损耗、线性、规则分层的问题。如果后面想加缺陷层、想模拟有限尺寸结构、想把热效应或应力耦合进来、想从一维扩展成二维三维，脚本的复杂度会瞬间膨胀。COMSOL 的优势在于模块化：几何画好、材料填好、周期边界条件设好，剩下的求解器替你干活，而且建模思路能从一维平滑迁移到二维和三维。尤其 6.4 版本把 Floquet 周期边界条件和模式跟踪做得越来越顺手，配合参数化扫描和 LiveLink，批量跑能带数据很方便。

当然 COMSOL 也不是没有代价：每个 k 点都要解一次特征频率，整体计算量比传输矩阵大不少；模式排序和后处理还得靠手动经验。所以我的建议是：一维定性判断用解析公式，成体系的能带扫描和后续扩展用 COMSOL。两条腿走路，效率最高。

2. 几何建模与材料参数：在硅基底上把周期性结构搭出来

2.1 先把晶格参数定死：晶格常数、厚度比、折射率

开始建模之前，先把手算的预期值算好，这样后面拿到仿真结果才有对照。我这次用的结构参数如下：

• 晶格常数 a = 1 μm；
• 硅层厚度 d_H = 0.45 μm；
• 二氧化硅层厚度 d_L = 0.55 μm；
• 硅折射率 n_H = 3.45（近红外取 3.45，实际随波长略有变化）；
• 二氧化硅折射率 n_L = 1.44。

为什么选这个组合？首先硅在近红外损耗很小，且折射率高；二氧化硅折射率低、材料库自带，工艺也成熟。其次，厚度比不取 1:1 而是 0.45:0.55，主要是想让带隙中心位置落在一个比较好记的频段。布拉格带隙中心频率近似满足

f₀ = c / [2(n_H d_H + n_L d_L)]

把数字代进去：n_H d_H + n_L d_L = 3.45×0.45 + 1.44×0.55 = 1.5525 + 0.792 = 2.3445 μm，于是 2 倍为 4.689 μm，对应 f₀约 64 THz，换算成自由空间波长约 4.7 μm。这个波长在红外中红外区，实际仿真时你完全可以按比例缩放：把 a 缩到 0.5 μm，带隙中心就跑到 9.4 μm 附近；把 a 缩到 0.2 μm，就到近红外 2.35 μm 左右。归一化频率 a/λ 才是真正重要的量，它约为 0.213。

2.2 材料参数设置：折射率、相对介电常数、虚部

COMSOL 里做光学仿真最省事的方式是直接用波动光学模块的“电磁波，频域”接口，这个接口认的是相对介电常数 ε_r 或折射率。非磁性材料默认 μ_r = 1，所以只填折射率就行。COMSOL 中新建材料时，可以在“折射率”栏填实部和虚部，也可以直接在“相对介电常数”里填 n²。

这里有一个非常容易踩的坑：很多人看到材料库里有“Silicon (Si)”就直接选，结果算出来能带全在低频漂移，因为材料库里的折射率是随波长变化的一组离散数据。如果是初学，建议先用常数折射率跑通流程，后续再换成色散材料。硅层填 n = 3.45，虚部暂时填 0，代表无损介质；二氧化硅填 n = 1.44，虚部同样为 0。如果你以后要模拟吸收、金属或损耗介质，特征频率会变成复数，虚部对应损耗，那种情况下要注意区分“损耗模式”和“传播模式”，前期不建议碰。

2.3 模型维度与几何搭建：一个晶胞就够了

严格的一维光子晶体可以在一条线段上建模，但 COMSOL 的电磁波频域接口在二维下用得更顺手，可视化也更清楚。我的做法是建一个二维矩形，x 方向长度为一个晶格常数 a，y 方向取一个小尺度，比如 0.1a，然后在模型里明确：x 方向是周期方向，y 方向是“无限延伸”的垂直方向。

这个 y 方向的小矩形会不会引入额外的模式？会。所以上下两条边要设置周期边界条件，且相位差为 0，等效于把 y 方向无穷堆叠，模仿无限大平面膜；这样一来，主要关注的模式就是沿 x 方向传播的平面波模。实际建模时也可以把 y 方向缩得非常小，比如 a/20，这样可以减少横向高阶模式，低频段就干净很多。总之，几何越简单越好，一个晶胞、两到三个域，别画多余结构。

画几何时，将矩形分成左右两块：左边宽 0.45 μm 设为硅，右边宽 0.55 μm 设为二氧化硅。画完记得用布尔“并集”或者“形成联合体”把它们合成一个装配体，否则后续网格和边界处理容易出现内部边界不连续的问题。很多人卡在这一步：两个域没有合并，求解时界面条件不对，能带算出来一团糟。

2.4 网格划分：到底要多细才算够

光子晶体仿真里网格尺寸是精度和速度的平衡。均匀介质中波长是 λ₀/n，硅的折射率最高，所以硅里的波长最短。如果你关心的是归一化频率 a/λ 在 0.2 附近的模式，对应自由空间波长约 5 μm，在硅里的波长只有 5/3.45 ≈ 1.45 μm，也就是大约 1.45 a。理论上每个波长放 8 到 10 个网格单元就能得到不错的精度，所以最大单元尺寸设在 0.05a（50 nm）左右比较稳。

我建议用映射网格：x 方向把硅层和二氧化硅层各剖分 30 到 40 个单元，y 方向剖 3 到 5 个单元。这样得到的网格规整、自由度少、求解快。如果直接用三角形自由网格，也能算，但同一个算例的离散误差会大一些，而且高频模式容易出“皱褶”。等模型跑通之后，再像模像样地做一次网格收敛性验证：把网格加密一倍，看前几条能带频率是否稳定到第三位小数，这比什么都管用。

3. 核心实操：特征频率 + 周期边界条件求解能带

3.1 最关键的一步：Floquet 周期边界条件怎么设

能带计算和普通本征模问题的最大区别在于边界条件。计算一个晶胞的特征频率时，左右边界不能是简单的完美电壁或完美磁壁，因为那样会把波矢强制成 0 或某些离散值，得不到完整的色散关系。正确做法是在左右边界上加 Floquet 周期边界条件，也叫 Bloch 周期边界条件，它要求电场满足

E(r + a) = E(r) exp(i k_x a)

这句式表示：周期方向相隔一个晶格常数处的场，除了一个相位因子外完全相同。相位因子里的 k_x 就是我们要扫描的波矢。COMSOL 操作时，在“电磁波，频域”节点下添加“周期边界条件”，先选左侧边界为源边界，右侧边界为目标边界，然后选择 Floquet 模式，并把波矢分量填成全局参数，比如 kx = Kx*a，其中 Kx 是归一化波矢。注意设置里周期矢量的方向要和几何一致，别把 x 方向填成 y 方向。

这一步是最容易出问题的。常见错误是只选了单侧边界，或者源目标边界选反，导致相位方向反了，能带图不对称。还有人在 Floquet 参数里直接填数字而不是全局参数，导致扫描不生效。建议把波矢分量定义成“Kx * pi / a”，然后扫描归一化参数 Kx 从 0 到 1，这样横轴就是 k_x a/π，非常直观。

3.2 参数化扫描波矢：从 Gamma 点到布里渊区边界

波矢扫描是能带计算的例行过程。在特征频率研究的设置里加一个参数化扫描，把参数 Kx 从 0 扫到 1，步长取 1/40 或 1/60。40 步已经能看出趋势，60 步会让能带更顺滑，但耗时接近线性上升。如果是第一次跑，建议先用 20 步，快速验证流程，确认没报错再加密。

这里有个需要注意的细节：传统能带图通常还包括负 k 方向，即 k_x 从 -π/a 扫到 π/a。但一维体系是中心对称的，正负 k 方向的能带完全一样，所以从 0 扫到 π/a 就够用了。如果你以后做二维三角晶格，扫描路径就要绕高对称点走折线，那另当别论。

扫描过程中，每个 Kx 都会对应一组特征频率。如果扫描点太少，后续模式追踪会很痛苦；扫描点太多则计算时间猛涨。我的习惯是先扫 30 点看结构，再针对带隙附近加密到 60 或 100 点，兼顾效率和精度。

3.3 特征频率的数量：不是越多越好

特征频率研究默认只求 6 个特征值，但光子晶体能带我们至少要拿到前 10 到 20 条带。否则图太稀疏，看不出带隙。COMSOL 里可以设置“所需特征频率数”为 20 甚至 40，也可以限定频率范围，比如只计算 0 到 200 THz 之间的特征模式。

特征值数量设得越多，每次求解的时间越长。计算量有点像爬楼梯：低阶模式快，高阶模式慢，而且高频模式对网格更敏感。我的经验是：第一次跑 15 个特征频率，如果带隙区域落在其中，就不要再盲加；如果目标是高频带，则必须同时加密网格。别贪心想一次算 100 条带，那会让单次求解慢到你怀疑人生。合理的目标是前面几条能带，一维光子晶体研究里最常用的通常是最低带和第一带隙。

3.4 模式追踪：能带图“串线”的土办法

这是整个流程里最烦也最耗时的部分。特征频率求解器在参数扫描时，输出的特征频率顺序是乱的：某个特征频率在 Kx = 0.1 时是第 3 条带，到 Kx = 0.2 时可能变成第 5 条带，因为不同模式的场形态在变化，求解器只按频率大小排序。如果直接按“每个 Kx 的第 n 个解”连接绘制能带图，你会得到一坨乱麻。

解决思路有三种。第一种最简单：把扫描步长收紧，让相邻 Kx 点之间同一模式的频率变化很小，然后取最低的几条模式，顺序自然就对了；但高频段依然容易乱。第二种是用场分布识别：在几个关键 Kx 点（0、半程、边界）导出电场或磁场分布，数清楚节点数，然后按模式物理连续性手动排序。第三种是正规做法：用 COMSOL 的模式跟踪功能，或者导出结果到 MATLAB/Python，用特征向量重叠积分自动重排，代码量也不大。

我个人在实际项目中通常会写一个半小时级的后处理脚本：每个 Kx 点导出所有特征向量的复振幅，计算相邻 k 点特征向量的内积绝对值，绝对值最大的一对视为同一模式的延续。对一维模型来说，这个办法又快又稳，而且顺便能发现是否有模式交叉。手工检查至少要做到：能带图里没有明显的“断线跳线”，带边连续平滑。

4. 后处理：从一堆频率点画出能带图

4.1 数据导出：把参数化扫描结果整理成表格

COMSOL 里解完参数化扫描后，结果节点下会生成一组“特征频率”解，每个 Kx 点对应一个解编号。你可以直接导出“全局计算”里的特征频率值，但默认导出格式是每个解一张表，不好直接画图。更推荐的做法是使用“结果 > 导出 > 数据”，在导出设置里勾选“所有参数值”，并把特征频率写入一行，列头可以设为频率 1 到频率 15。导出的 CSV 结构大致是：

Kx, freq1, freq2, freq3, ...

这就是后面画图的数据源。如果在 COMSOL 里看到某些频率是复数（比如填了损耗），画图时取实部。

4.2 用 Python 画能带图：小代码解决大麻烦

我习惯导出 CSV 后用 Python 画图，因为能带图经常需要反复修改样式、叠加实验数据或解析理论线，脚本比 COMSOL 后处理灵活得多。下面这段代码足够完成基本绘制：

import pandas as pd import numpy as np import matplotlib.pyplot as plt # 导出文件假设包含：Kx, freq1, freq2, ..., freqN df = pd.read_csv("band_data.csv") kx_norm = df["Kx"].values # 归一化波矢 k*a/pi c0 = 299792458.0 # m/s a = 1e-6 # 晶格常数 1 um fig, ax = plt.subplots(figsize=(6, 5)) for i in range(1, 16): # 画前 15 条能带 f = df[f"freq{i}"].values # 纵轴用归一化频率 a/lambda = f * a / c0 y = f * a / c0 ax.plot(kx_norm, y, "-o", markersize=2.5, linewidth=1.2) ax.set_xlabel(r"$k_x a / \pi$") ax.set_ylabel(r"$a/\lambda$") ax.set_xlim(0, 1) ax.grid(True, linestyle="--", alpha=0.4) plt.tight_layout() plt.savefig("bandstructure.png", dpi=300)

注意，“-o”只是辅助看数据点，正式发表时去掉标记。如果某个频率点为负数或 0，多半是模式没解出来，检查研究设置。也有人用 MATLAB 画，原理一样。

4.3 能带图怎么看：带隙在哪、带边怎么找

拿到图后，第一眼看 Kx = 0（Gamma 点）附近。均匀介质极限下，最低带的频率随 Kx 线性上升，斜率对应有效折射率；如果斜率很平，说明模式群速度慢，这就是慢光效应的基础。第二眼看 Kx = 1（布里渊区边界）处：第一条带在边界处频率达到某个值后不再上升，第二条带从更高频率开始，两条带之间没有模式覆盖的区域就是光子带隙。带边频率分别记为 f₁ 和 f₂，带隙宽度是 Δf = f₂ - f₁，相对带宽是 Δf / f₀，其中 f₀ 是带隙中心频率。

按我们前面算的布拉格中心频率 a/λ ≈ 0.213，仿真出来的带隙会明显落在 0.18 到 0.25 这个区间附近。因为硅和二氧化硅折射率差大，这个一维带隙通常很宽，相对带宽可能达到 30% 甚至更宽。如果看到 Kx = 0 处出现重频、频率简并，那是正常的；如果没有带隙，先回去检查是不是折射率没填对、周期厚度比偏了、或者网格太粗。

4.4 与解析解对比验证：仿真结果靠不靠谱

算完能带之后，至少用两个解析判断验证一下。一是带隙中心频率。前面用布拉格条件算出的 f₀ = c/[2(n_H d_H + n_L d_L)] 应当落在仿真带隙的中心附近。代入具体数值得 a/λ₀ ≈ 0.213，仿真结果如果偏离超过 5%，优先怀疑材料折射率或几何尺寸填错。二是带隙宽度趋势。对于两介质交替堆叠，折射率对比度越大带隙越宽；如果你把 SiO₂ 换成空气，n_L = 1，带隙会显著变宽。这个趋势不需要精确解析，看一眼就能心里有数。

如果连这个都对不上，就做网格收敛性测试：把最大网格尺寸从 0.05a 缩到 0.025a，重新算 10 个 Kx 点，比较带边频率变化。变化小于 1% 就说明网格已经够用，不用再烧时间。

5. 常见问题与排查技巧实录

5.1 能带曲线乱跳、交叉全是毛刺

这是出现频率最高的问题，本质就是模式排序混乱。如果你已经在后处理时把每个 Kx 点的特征频率手工排过序，还是乱，那么可能的原因有二：一是扫描步长太大，相邻 Kx 之间同一模式的频率变化超过了模式间距，导致追踪失效；二是特征模式数设得太多，高频模式在网格不够密时本身就不准，排序自然更乱。

排查办法说起来很简单：先从 Kx = 0 到 1 之间取 5 个点，每个点导出场分布，观察电场节点数，把模式按物理解编号，再看曲线是否连续。如果中间有模式交叉，物理上也是允许的，但在无耦合二维模型中交叉通常是模式“绕过”而不是真正交叉，需要看场分布确认。

我自己的土办法是：在脚本里对每个 Kx 解，计算相邻 Kx 的特征向量投影矩阵，把最大投影行视为同一模式。这套方法对一维结构成功率接近百分之百，但要求你导出的是完整复特征向量，而不只是频率。COMSOL 支持导出特征向量分量，值得花时间设置一次。

5.2 Floquet 边界条件设了没反应或直接报错

最常见的报错是“边界条件参数不一致”或“未识别边界”。检查三项：第一，左右边界是否都选上了，并且源边界和目标边界是成对的；第二，周期矢量的方向是否与几何方向一致，一维模型里周期矢量就是 (a, 0)，如果填成 (0, a) 那就什么都对不上；第三，波矢分量是否引用全局参数，如果直接填数字，参数扫描就不会更新边界条件，结果和 Kx = 0 完全一样。

还有一个隐蔽问题：当几何由多个域组成但没有“形成联合体”时，左右边界的点可能有重复边或内部边，COMSOL 会把内部边也当成实体边界，导致周期边界选到错误的边上。处理方式是在建立几何时统一用“形成联合体”，并且网格划分时启用“自动移除内部边界”，这样干净很多。

5.3 高频模式不收敛，能带背部出现褶皱

这个问题几乎都是网格太粗。频率越高，波长越短，相同网格尺寸对应的每波长单元数越少，数值色散就会以“非物理褶皱”的形式出现。还有一个容易被忽视的因素：y 方向尺寸如果取得太大，会出现高阶横向模式，这些模式混在低频区域，很容易让人误判成光子晶体的新能带。

解决方法有两个方向：一个是加密网格，尤其是硅层内部，用映射网格把每层分成 40 段以上；另一个是减小 y 方向高度，比如从 0.1a 减到 0.02a，把横向模式推到更高频段，低频区自然只剩平面波模。如果减了高度之后最低几条带不变，说明之前的横模干扰确实存在。

5.4 计算时间太长：怎么压性价比

参数扫描 + 特征频率本身就是时间大户。提升效率的办法按优先级排：先限制特征频率数量，从 40 减到 15 或 20，通常前几条带足够用了；再减扫描点数，先用 20 个 Kx 点跑通，后期再在带隙附近加密；最后才是加内存、开多核。

COMSOL 6.4 在多核并行上表现不错，参数化扫描的各个 Kx 点也可以多进程并行。在集群或工作站上用 Floquet 边界条件做能带计算时，建议把参数扫描拆成多个独立任务，每个任务负责一段 Kx 范围，最后合并结果，这样能显著缩短墙钟时间。个人电脑上如果实在慢，检查是不是网格剖得太密：0.025a 以下实际上很少必要。

5.5 一维光子晶体仿真的问题速查表

这张表是我反复被折磨之后总结的，拿到任何异常结果，先按表里顺序查一遍，基本能解决九成问题。

6. 能带算完后的进一步玩法

6.1 从无限周期到有限堆叠：布拉格反射镜反射谱

能带图算的是无限周期结构的本征性质，但实际器件总是有限层数。你可以把晶胞复制 5 到 10 个周期，在两侧加端口和 PML 边界，扫频计算反射率和透射率。反射率在带隙频率范围内会接近 1，范围恰好对应能带图中的带隙。这一步能验证能带计算和实际器件的一致性，也是把仿真从“好看”推向“有用”的关键一步。我建议选 8 个周期，叠厚约 4.7 μm，反射谱会非常漂亮，带边缘也比较陡。

6.2 在带隙里“开窗”：缺陷模的引入

一维光子晶体最经典的应用是在带隙中制造缺陷模。方法很简单：把周期堆叠里某一层的厚度或折射率改掉，比如把中间一层 SiO₂ 厚度加倍，带隙中就会出现一个极窄的透射峰。这个峰对缺陷层的光学厚度极其敏感，环境折射率一变，峰位就漂移，因此非常适合做传感器或窄带滤波器。

COMSOL 里做缺陷模有两种思路：一种是直接用有限周期结构加端口激励，扫频看透射峰；另一种是用超胞法，把带缺陷的周期序列当成新的“大晶胞”，再用 Floquet 边界条件算超胞能带，带隙里的平带就是缺陷模。第二种方法物理意义更清晰，推荐在一维结构练熟后再上手。

6.3 从一维到二维三维：思路通用，工作量和陷阱翻倍

当你把一维流程摸清楚后，往二维光子晶体（比如三角形排列的介质柱或空气孔）迁移会轻松很多，因为核心操作没变：选取元胞、加 Floquet 周期边界、参数扫描高对称点路径。差别在于网格复杂度剧增，模式数量也更多，后处理排序变得更难。三维全矢量仿真的自由度会再高一个数量级，计算成本陡增，初学者容易直接陷入“算到死”的困境。

我的建议是务必先把一维场景下的边界条件、参数扫描、模式追踪、后处理脚本全部验证好，形成一套可复用的 workflow，然后再尝试二维。否则一边调边界一边改网格一边追模式，出了错都不知道是哪个环节的问题。个人习惯是每次改参数后都把网格加密一倍重新算一次，带隙位置稳定到三位有效数字才敢放心画进论文。这一步验证看起来折腾，实际上能帮你在审稿人问“网格收敛性如何”的时候拿出硬数据，比一句“默认设置”有说服力得多。