先声明一下这篇文章不是从教科书里搬概念而是把我在 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 方向长度为一个晶格常数 ay 方向取一个小尺度比如 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.05a50 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[ffreq{i}].values # 纵轴用归一化频率 a/lambda f * a / c0 y f * a / c0 ax.plot(kx_norm, y, -o, markersize2.5, linewidth1.2) ax.set_xlabel(r$k_x a / \pi$) ax.set_ylabel(r$a/\lambda$) ax.set_xlim(0, 1) ax.grid(True, linestyle--, alpha0.4) plt.tight_layout() plt.savefig(bandstructure.png, dpi300)注意“-o”只是辅助看数据点正式发表时去掉标记。如果某个频率点为负数或 0多半是模式没解出来检查研究设置。也有人用 MATLAB 画原理一样。4.3 能带图怎么看带隙在哪、带边怎么找拿到图后第一眼看 Kx 0Gamma 点附近。均匀介质极限下最低带的频率随 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 一维光子晶体仿真的问题速查表现象常见原因优先排查动作能带全是乱线模式排序混乱减小扫描步长或导出场分布手动排序带隙消失折射率填错、厚度比偏差、网格太粗核对材料折射率和 d_H/d_L加密网格低频出现拓扑怪异模式y 方向取了太大高度减小 y 尺寸或用周期边界抑制横模参数扫描没有改变结果Floquet 波矢填成了固定数字改成全局参数并确认参数扫描启用特征频率出现小虚部材料填了损耗或漏填虚部检查材料折射率虚部是否为 0高频皱褶严重网格不足以分辨短波长在最高折射率材料内加密网格这张表是我反复被折磨之后总结的拿到任何异常结果先按表里顺序查一遍基本能解决九成问题。6. 能带算完后的进一步玩法6.1 从无限周期到有限堆叠布拉格反射镜反射谱能带图算的是无限周期结构的本征性质但实际器件总是有限层数。你可以把晶胞复制 5 到 10 个周期在两侧加端口和 PML 边界扫频计算反射率和透射率。反射率在带隙频率范围内会接近 1范围恰好对应能带图中的带隙。这一步能验证能带计算和实际器件的一致性也是把仿真从“好看”推向“有用”的关键一步。我建议选 8 个周期叠厚约 4.7 μm反射谱会非常漂亮带边缘也比较陡。6.2 在带隙里“开窗”缺陷模的引入一维光子晶体最经典的应用是在带隙中制造缺陷模。方法很简单把周期堆叠里某一层的厚度或折射率改掉比如把中间一层 SiO₂ 厚度加倍带隙中就会出现一个极窄的透射峰。这个峰对缺陷层的光学厚度极其敏感环境折射率一变峰位就漂移因此非常适合做传感器或窄带滤波器。COMSOL 里做缺陷模有两种思路一种是直接用有限周期结构加端口激励扫频看透射峰另一种是用超胞法把带缺陷的周期序列当成新的“大晶胞”再用 Floquet 边界条件算超胞能带带隙里的平带就是缺陷模。第二种方法物理意义更清晰推荐在一维结构练熟后再上手。6.3 从一维到二维三维思路通用工作量和陷阱翻倍当你把一维流程摸清楚后往二维光子晶体比如三角形排列的介质柱或空气孔迁移会轻松很多因为核心操作没变选取元胞、加 Floquet 周期边界、参数扫描高对称点路径。差别在于网格复杂度剧增模式数量也更多后处理排序变得更难。三维全矢量仿真的自由度会再高一个数量级计算成本陡增初学者容易直接陷入“算到死”的困境。我的建议是务必先把一维场景下的边界条件、参数扫描、模式追踪、后处理脚本全部验证好形成一套可复用的 workflow然后再尝试二维。否则一边调边界一边改网格一边追模式出了错都不知道是哪个环节的问题。个人习惯是每次改参数后都把网格加密一倍重新算一次带隙位置稳定到三位有效数字才敢放心画进论文。这一步验证看起来折腾实际上能帮你在审稿人问“网格收敛性如何”的时候拿出硬数据比一句“默认设置”有说服力得多。