一维光子晶体的能带计算说难不难说简单也有一堆细节等着你踩。先把“高速公路收费站”这个比喻记在脑子里一维光子晶体就是给光子设卡收税的周期性介质结构折射率每隔一段距离变化一次频率合适的光子一路畅通频率不对的直接被反射回去——这不就是收费站嘛。这篇文章我打算带你用 COMSOL 在硅基底上搭一个周期性介电结构也就是 SiO2/Si 交替叠层的一维光子晶体用特征值法把光子能带图算出来再从这个能带图里读出光子带隙。适合正在学光子晶体、微纳光学、COMSOL RF 模块的朋友参考尤其是那种“看了教材但不知道模型从哪下手”的初学者。我会把建模思路、参数设定、边界条件、网格划分、求解器配置、后处理出图到问题排查全部过一遍尽量让你照着做就能复现结果。1. 一维光子晶体能带到底在算什么1.1 收费站的比喻怎么理解先把这个比喻拆开讲清楚。一维光子晶体最常见的形式就是高低折射率介质交替堆叠比如 SiO2折射率 1.46和 Si折射率 3.5一层一层交替排列。光在里面传播时每一层界面都会发生菲涅耳反射无数个界面的反射光互相干涉最后在某些频率范围内形成强烈的反射带这就是光子带隙。能带计算做的事情本质上是求解光在这种周期性折射率分布中的本征模。和电子的能带理论一样布洛赫定理告诉我们周期介质中的电磁场可以写成平面波与周期函数的乘积于是我们只需要在一个晶胞里求解再通过周期性边界条件推广到整个结构。COMSOL 里的特征值求解器就是干这个活的给定一个布洛赫波矢 k它就会把对应的特征频率算出来k 从 0 扫到布里渊区边界一条条能带就出来了。1.2 布拉格条件带隙位置从哪来一维光子晶体带隙的中心频率完全可以先用手估出来这就是布拉格条件。对折射率为 n1、n2、厚度为 d1、d2 的两层材料组成的周期结构一阶带隙中心满足[ 2(n_1 d_1 n_2 d_2) m\lambda ]m 取 1 时带隙中心波长就是 (2(n_1 d_1 n_2 d_2))。我这次用的结构是 Si 层 200 nm、SiO2 层 400 nm那么[ \lambda_{gap} 2(3.5\times 200 1.46\times 400) 2(700 584) 2568 \text{ nm} ]换算成频率大概是 116.8 THz。这个估算很有用它能帮你设定特征值搜索的频率范围让你不必在 COMSOL 里盲扫。但要强调一下这只是中心频率的位置带隙到底多宽、带边在哪布拉格条件给不出来必须靠数值仿真。1.3 折射率对比度决定带隙宽度带隙宽度直接取决于两种材料的折射率对比度。对比度越大带隙越宽。Si 和 SiO2 的对比度大概是 2.4 倍作为入门案例很合适带隙明显但又不至于像 Si/空气这种大对比度结构那样出现复杂的表面态和漏模。而且硅基底在近红外波段损耗低、工艺成熟这也是它被大量用于滤波器、反射镜、传感器设计的核心原因。如果你后面想做宽带反射镜可以考虑用更大折射率差的结构或者做渐变/准周期结构。但入门阶段我建议先用 Si/SiO2原因很现实材料参数稳定、模式少、结果容易解释。2. 仿真之前模型选型与参数设计2.1 用 2D 模型模拟一维周期结构很多人会纠结一维光子晶体是不是用一维模型算实际上 COMSOL RF 模块里最顺手的是 2D 模型。为什么因为一维周期结构虽然在 x 方向周期变化但 y 方向可以看成无限均匀延伸我们取一个 2D 截面就能完整描述整个物理图像而且还能后处理看场分布。模型边界取一个晶胞就够了x 方向宽度就是晶格常数 a d1 d2 600 nm。y 方向高度我一般取 300 nm 左右但要注意这个方向不能随随便便设成 PEC 边界完事。后面的物理场设置里我会专门讲这个问题这里先记住上下边界也要用周期性边界处理否则会出现一堆假模式。2.2 材料参数怎么设才不翻车COMSOL 的材料库里确实有 Si 和 SiO2但我建议你在入门阶段直接手填常数介电常数不要用材料库里的色散模型。原因有两点一是库里的硅模型往往覆盖很宽的波长范围表达式复杂低频或高频边界容易出数值警告二是能带计算本质上关心折射率比和厚度用常数已经足够反映物理。具体设置很简单在“材料”节点里新建空材料把相对介电常数设为常数Si 取 12.25即 3.5 的平方SiO2 取 2.13即 1.46 的平方相对磁导率都是 1电导率 0。这组参数在近红外到中红外区间是很合理的近似。2.3 几何尺寸和占空比怎么给定我用的一组参数如下参数取值说明晶格常数 a600 nmx 方向一个周期的长度Si 层厚度 d1200 nm高折射率层SiO2 层厚度 d2400 nm低折射率层模型高度 h300 nmy 方向截断高度占空比 f1/3d1/a占空比直接决定带隙中心频率和高阶带结构。一般来说高低折射率层的光学厚度接近四分之一波长时带隙最强。我们这里 n1d1 700 nmn2d2 584 nm差别不大所以带隙效果很明显适合做教程案例。想优化的话可以用 COMSOL 的参数扫描直接扫占空比看哪个位置的带隙最宽。3. COMSOL 实操从建模型到出能带图3.1 模型向导与全局参数打开 COMSOL模型向导里选择“二维”物理场选择 RF 模块下的“电磁波频域ewfd”研究选择“特征值”。这一步很容易漏掉如果你安装的时候没选 RF 模块授权模型向导里根本看不到这个选项所以先检查许可证。进入模型后先在“全局定义”里把参数填好。我会写一个参数列表名称表达式描述a600e-9晶格常数单位 md1200e-9Si 层厚度d2400e-9SiO2 层厚度h300e-9模型高度kx0布洛赫波矢 x 分量eps_Si12.25硅介电常数eps_SiO22.13二氧化硅介电常数单位一定要统一COMSOL 默认国际制几何尺寸用米。别在这里图省事写 nm后面算波矢和频率的时候单位换算会让你疯掉。3.2 几何和材料分配几何很简单两个矩形并排放。第一个矩形从原点开始宽 400 nm、高 300 nm作为 SiO2 层第二个矩形从 x 400 nm 开始宽 200 nm、高 300 nm作为 Si 层。总共就是一个 600 nm 宽的晶胞。材料分配时把第一个矩形指给定义好的 SiO2 材料第二个矩形指给 Si 材料。如果你不想手动建材料也可以直接在“域”上右键“指定”但用材料节点更规范后面改参数也方便。3.3 物理场与 Floquet 周期边界这一步是能带计算的核心也是最容易翻车的地方。首先要选偏振方向在一维光子晶体里我们通常关心 TE 偏振也就是电场方向沿面外 z 轴。所以物理场设置里把“面外电场”选上这样求解变量就是 Ez模型可以被大幅简化。然后是边界条件。左右两个边界必须是一对 Floquet 周期边界用来施加布洛赫条件。在 COMSOL 里添加“周期边界条件”节点选择左边界为源、右边界为目标并选择“Floquet 周期”类型。在波矢分量里x 分量填全局参数 kxy 分量填 0。这里有个关键细节上下两个边界也要处理。很多人图省事直接给上下边界加完美电导体结果能带图上出现一堆奇怪的平直伪模式。正确做法是给上边界和下边界也加一对 Floquet 周期边界波矢 y 分量设 0x 分量设 0。你可能觉得两个方向的周期边界有点奇怪但这等价于假设 y 方向无限均匀并且完全不会引入人工边界反射。我实测下来用这个方案得到的能带最干净。3.4 网格划分的匹配细节网格是另一个容易埋雷的地方。Floquet 周期边界要求源边界和目标边界的网格节点一一对应否则边界条件匹配会出现数值误差严重时直接报错。最好的办法是用映射网格对两个矩形统一划分。网格尺寸方面x 方向要能分辨最短波长。如果你扫描到 300 THz光在硅中的波长大约是 285 nm按 1/10 网格规则最大尺寸应该控制在 30 nm 以内所以 x 方向我设最大单元尺寸 30 nm。y 方向因为场是均匀的可以松一点设 50 nm 也没问题。这样整个模型网格数量只有几百个特征值求解几乎瞬间完成。如果你发现 Floquet 边界报“源和目标网格不一致”多半是因为两个矩形被分开划分导致边界节点错位。解决办法就是把两个域同时选进映射网格或者干脆用自由三角形网格后检查边界节点数。3.5 特征值研究与波矢扫描研究设置里特征值求解器有几个参数需要注意。首先是“期望特征值数”我一般填 10 到 15 个太少会漏带太多会增加计算量。然后是“搜索特征值范围”下限设 1e12 Hz上限设 3e14 Hz也就是 1 THz 到 300 THz。这个范围覆盖了我们关心的 116.8 THz 带隙附近的几条带。最关键的是参数化扫描。要获得完整能带图需要让 kx 从 0 扫到 π/a。在“研究”里添加“参数化扫描”扫描参数选 kx表达式写 range(0, pi/(30*a), pi/a)也就是分成 30 步共 31 个点。先这样跑通流程后面想细化曲线可以把 30 改成 60 或 100。求解器用默认配置就行网格这么小根本不会遇到内存问题。等计算结束COMSOL 会输出每个 kx 点上的多个特征频率。3.6 后处理导出数据画能带图计算完成后结果里创建一个“一维绘图组”添加“全局”绘图。数据选择“研究 1/参数化解”表达式填 freq/1[THz]x 轴数据选表达式 kx。这样 COMSOL 就能把所有特征频率随 kx 变化的曲线画出来每条曲线就是一条能带。不过 COMSOL 原生的绘图在论文里通常没法直接用。我一般是导出数据到表格然后拿到 Origin 或者 Python 里重画。具体操作是结果里右键“派生值”选“全局计算”表达式写 freq/1[THz]然后点计算得到一张表再把表导出成 txt 或 csv。这个过程虽然有点繁琐但数据到手后就能随便排版了。4. 结果怎么看能带图、带隙和模式场4.1 从能带图中圈出光子带隙能带图横轴是波矢 kx纵轴是频率或归一化频率 f*a/c。你会在图里看到几条从左下到右上弯曲的曲线最低那条从 0 附近逐渐上升这就是第一条能带。往上还会有第二条、第三条。如果没有伪造模式这些曲线之间会出现一个明显的空白区域——这个区域里没有任何一条曲线穿过它就是光子带隙。对我们这个参数带隙大致在 110 到 130 THz 附近具体边界以你的仿真结果为准。想量化带隙宽度就看带顶和带底的频率差。你会发现一个有趣的现象在 kx π/a 这个布里渊区边界上能带往往不连续或有明显的弯曲转折这正是布拉格反射最强的波矢点。4.2 带边缘模式场分布验证能带图只能告诉你“哪里有带隙”但要说清楚“为什么”必须看模式场分布。比如在 kx π/a 处的第一条带顶模式电场能量倾向于集中在低折射率的 SiO2 层第二条带底模式则相反能量集中在高折射率的 Si 层。这是因为带边缘的波是驻波驻波的波腹会按照折射率不同的层重新分配能量从而拉大了两个模式的频率差。在 COMSOL 里你可以在结果中选择某个特征频率然后画 Ez 的面分布图。记得把“数据集”切换到参数化解里对应的 kx 点否则你看到的场不是想要的模式。这个验证过程强烈建议做一遍它比单纯画能带图更能帮你建立物理直觉。4.3 从能带到反射谱带基底的实器件怎么处理能带计算有个前提假设整个空间都被周期介质填满。但实际场景是“硅基底上长周期结构”真实器件上方是空气、下方是硅基底周期结构只有有限层数。这时候反射率和透射率才是你真正关心的量。想算这种实际结构不建议再用特征值研究而是用频域研究。把几何建成长条从上到下依次是空气层、周期叠层比如 8 到 10 对 SiO2/Si、硅基底再往下是空气层或 PML 吸收边界。x 方向依然加 Floquet 周期边界用端口激发平面波计算 S 参数反射谱。你会看到反射谱在能带隙对应的频率范围内出现高反射平台平台高度和周期对数直接相关层数越多反射越接近 1。这个流程很简单但很实用工程上做布拉格反射镜基本就是这么算的。我建议你先算完能带图再顺手算一个反射谱对比着看两边的带隙位置会对得很齐那种“理论和仿真对上了”的感觉非常爽。5. 常见问题与排查技巧实录5.1 伪模式与上下边界条件问题我最开始算一维光子晶体能带时上下边界顺手设了 PEC结果能带图里冒出一堆几乎平直的假曲线。这些平直模式并不是真正的布洛赫模式而是上下边界把电磁场约束成波导模式造成的。判断方法很简单真能带曲线随 kx 有明显变化伪模式则是一条条水平的或近似水平的线而且改变模型高度 h 时伪模式会移动真模式不动。解决办法就是前面说的上下边界改用 Floquet 周期边界ky 设 0。改完后你会发现能带图立刻干净了。这个经验算是我踩过的最深的坑之一分享出来希望你不走弯路。5.2 Floquet 周期边界常见报错实际操作中Floquet 边界最常见的问题是“源边界和目标边界网格不一致”。原因往往是两个域分开划分网格导致左右边界节点不对齐。解决办法是强制用映射网格并把两个域一起选中或者先划分边界单元再同步到域。还有个容易忽略的点Floquet 波矢的单位。kx 的单位是 rad/m很多人会惯性写成 pi/a但如果你把 a 设成了 600以为单位是 nm结果就是 pi/600波矢小了 1e9 倍能带图完全不对。必须写 pi/600e-9。这种问题不报错但结果一看就离谱排查起来还挺费时间。5.3 特征值太少或搜索不到带隙模式有时候你只设了 6 个期望特征值扫到高频段时某些带缺失能带图上断了一条。这不是结构没带隙而是你要求算的特征值数量不够。解决方法是增加期望特征值数比如 15 个同时把搜索范围上限往下调只关心低频的前几条带避免求解器在无用模式下浪费资源。另外如果你在带隙内部看到了孤立的特征频率先别急着高兴它很可能是网格太粗导致的伪模式。加密网格后再看如果这条频率消失了就是数值假象。网格越密特征值越收敛伪模式越少。下面是常见问题的速查表现象可能原因解决方法能带图出现大量平直线上下边界条件不匹配改用 Floquet 周期边界ky0Floquet 边界报错网格不一致左右边界网格节点不对齐用映射网格统一划分能带缺条或断带期望特征值数太少增加到 10~15 个带隙内出现孤立频率网格过粗加密 x 方向网格至 30 nm结果完全不对但无报错单位写错如 nm 当 m 用检查 kx 表达式和几何单位6. 最后再分享两个小技巧第一个技巧如果你拿到了 COMSOL 导出的表格用 Python 画能带图其实只要几行代码。比如导出的 csv 里有两列kx 和 freq直接读取然后 scatter 就行这样图表排版更自由投稿有截稿压力时效率极高。import pandas as pd import matplotlib.pyplot as plt df pd.read_csv(banddata.csv) plt.plot(df[kx], df[freq], ., markersize2) plt.xlabel(kx (rad/m)) plt.ylabel(Frequency (THz)) plt.ylim(0, 300) plt.show()第二个技巧想快速预估带隙位置先用手算布拉格条件定出中心频率再设 COMSOL 的搜索范围通常一次就能拿到完整能带。我带学生做这个案例时发现 90% 的时间其实花在调试边界条件上真正物理建模很快。后面如果你想进阶可以试试把结构从一维叠层改成二维光子晶体板也就是 x、y 两个方向都有周期z 方向厚度有限。思路是一样的区别在于布里渊区变成了二维需要沿着 Γ-X-M 这样的高对称路径扫描波矢网格量和计算量会明显上涨。但底层逻辑和今天讲的一模一样懂了这套进阶只是时间问题。