一维光子晶体这东西名字挺唬人其实就是把两种折射率不同的介质按固定“节拍”交替排列光走进去就不再是闷头直冲而是会被这排周期结构反复筛选——有些频率顺利放行有些频率直接被弹回来。我习惯拿高速公路收费站来打比方一排收费杆按固定间距立着不同频率的光子就像不同车速的汽车杆间距就是晶格常数杆子之间的空隙就是介质占空比能不能抬杆放行全看这套周期性介电结构是怎么设计的。这篇文章要干的事很具体用COMSOL在硅基底上搭一维光子晶体模型把光子能带算出来找到禁带位置再讲讲怎么通过调整结构参数移动带隙。适合正在做硅光栅、多层膜反射镜、或者建了模型却不知道能带图怎么稳定算出来的朋友。后面不堆大公式但关键原理、参数设置和实操步骤都会一步步拆开讲包括Floquet周期边界条件怎么加、特征频率扫描怎么避免扫出大量假模式、硅基底怎么处理才能不污染能带结果。可以提前说一句光能带不玄只要周期没错、边界没错、网格没烂COMSOL会给出一张非常干净的带隙图。1. 能带计算之前先清楚一维光子晶体到底在“拦”什么1.1 收费站比喻背后的光学原理我们要把“高速公路收费站”这个故事翻译成光学语言。两种介质折射率分别是 n_H 和 n_L厚度为 d_H 和 d_L它们组成一个周期单元晶格常数 a d_H d_L。光打进去每经过一层界面就会发生部分反射和透射这些反射光之间会相互干涉。当波长满足布拉格条件也就是 2(n_H d_H n_L d_L) 等于波长的整数倍时所有界面的反射光会在入射方向同相叠加形成强反射。这个强反射区间在能带图上就表现为“禁带”。简单说收费站的杆子不是随机摆放的而是每隔一段固定距离放一根。车辆光子按固定节奏通过这些杆界面符合这个节奏的频率被规则接纳不符合的就被干涉消掉。这里的节奏就是布拉格条件。仿真里我们算的不是某一条光线的透射或反射而是这个周期结构里所有可能存在的电磁波本征模式。本征模式存在的频率区域叫通带不存在的区域就是禁带。所以“能带计算”本质上是在求解周期介质中的麦克斯韦方程组得到频率和波矢之间的关系也就是色散关系 ω(k)。这段关系为什么重要因为一个光子晶体的全部光学行为几乎都能从能带图中读出来。禁带位置决定它能反射哪个波段带边缘的色散曲率决定光的群速度甚至慢光效应、负折射这些现象也都能从能带图里追踪。所以我一直建议哪怕你只是想做一层反射镜也最好先把能带图算出来因为能带图是“全视角”而单波长透射率只是“单视角”。1.2 为什么选硅基底和周期性介电结构硅在近红外波段折射率约为3.48和二氧化硅约1.45形成大约0.85的折射率差。这个对比度决定了禁带的宽度。一维光子晶体的带隙宽度大致正比于折射率对比度对比度越大禁带越宽器件对工艺误差的容忍度也越高。如果换用聚合物或者玻璃材料折射率差只有0.3甚至更低带隙就会窄很多仿真和实测都容易被加工误差打穿。硅基底的另一个好处是工艺成熟。SOI晶圆、深紫外光刻、干法刻蚀这些都是半导体行业玩了几十年的东西。你仿真里画出来的硅条流片厂能直接按设计加工出来这对后续做实验验证至关重要。而且硅本身在通讯波段是透明的1300nm和1550nm这两个光通信主力波段硅的吸收损耗都非常低。所以“Si/SiO₂周期性结构”是入坑光子晶体最合理的第一站。周期性介电结构在一维情形下就是交替的层堆叠或者在硅膜上刻出周期性的光栅条。前者是一维层状结构光线垂直入射则问题退化成一维后者实际上是二维截面但只要沿着光传播方向保持周期性同样满足一维光子晶体的物理模型。我在这篇里用层状堆叠来做演示因为几何最简单、物理最直观、参数最少带隙规律也最容易验证。1.3 能带图里的禁带到底意味着什么能带图的横轴是波矢 k纵轴是频率 ω或者归一化频率 a/λ。对自由空间来说ω 和 k 是线性关系一条光锥线走到底。加入周期结构后由于布拉格散射这条线在布里渊区边界处会被“掰开”形成一个禁止区间。这个区间里没有任何实数波矢对应的传播模式意味着光子走不进结构或者穿不过结构。这样的频段安在反射镜上就是镜子的工作波段安在滤波器上就是阻带。在实际算能带时有一件特别容易忽略的事情能带图的横轴波矢是有范围的。对一维晶格第一布里渊区就是 k 从 -π/a 到 π/a。因为周期性能带结构和 k 是周期重复的所以我们只需要计算 k 在 0 到 π/a 这个不可约布里渊区内的色散关系。很多新手第一次画能带图时发现曲线一堆纠缠原因往往是 k 扫出了范围把重复的高阶模式也画进去了。我还要提一个归一化问题。COMSOL里特征频率求解得到的频率单位一般是Hz而能带图横轴是波矢两者直接画出来数值大小会受结构绝对尺寸影响。做参数研究时我喜欢先归一化成 a/λ。这样横轴从0变到0.5对应k从0到π/a纵轴就是归一化频率无论晶格常数是500nm还是1μm物理规律完全一致。这个习惯能帮你快速对照文献、迁移设计。2. 建模之前的参数设计为什么先算一步再动COMSOL2.1 从目标波长反推晶格常数和占空比打开COMSOL之前最好先在纸上把几何参数定下来否则直接在软件里乱调参数最后会陷入“反复重画—算不出合适带隙—再重画”的循环。设计一维光子晶体的关键参数只有四个晶格常数 a、占空比 f d_H / a、高折射率层厚度 d_H、低折射率层厚度 d_L。一阶布拉格条件的中心波长近似是 λ_c 2(n_H d_H n_L d_L)。如果我的目标是把禁带中心设计在1550nm附近且采用四分之一波长堆叠也就是让 n_H d_H n_L d_L λ_c / 4那就可以反推厚度。设 d_H λ_c / (4 n_H) ≈ 1550 / (4×3.48) ≈ 111nmd_L λ_c / (4 n_L) ≈ 1550 / (4×1.45) ≈ 267nm晶格常数 a d_H d_L ≈ 378nm。这个尺寸在半导体工艺里完全做得出。但要注意这个反推公式只对垂直入射、一维无限周期堆叠成立。你仿真里如果用周期性边界条件模拟无限周期这个估算值就比较准。如果只堆叠五六层带隙中心会轻微移动我在实操部分会展示怎么校准。也可以反过来先定一个整数周期的工艺尺寸比如 a 400nm、f 0.35再扫频看带隙落在哪里。两种做法都行关键是有个起始依据而不是随机填数。占空比的选择也有讲究。f 太小高折射率层太薄折射率调制作用不强带隙窄f 太大低折射率层被挤得太薄同样会削弱调制。经验上 f 在0.3到0.5之间比较合适。如果希望带隙更宽又不做特殊优化就取0.5高低折射率层等厚几何上最好画工艺上也最好刻。2.2 用离散模型还是连续模型有人会问硅基底上刻光栅那个基底要不要在整个模型里画成一个大方块这里有个建模的重要决策。算光子能带默认的前提是结构无限周期性延展。一维光子晶体既然是无限周期的那么在周期方向上我们只需要取一个晶胞配合Floquet周期边界条件就能代表无限结构。这个思想来自Bloch定理它能让我们用极小的计算量获得无限周期结构的完整物理信息。硅基底如果也参与周期性排列那问题简单晶胞里带上基底切片整个模型就是一个周期单元。但如果基底是一整块厚硅片不是周期排列的那就要区分算能带时基底最好当作“衬底”处理在半导体器件里它就是一个物理支撑不影响光在周期堆叠中的本征模式。实际操作上绝大多数人能带计算的目的是分析周期介质本身的色散特性所以基底可以不画或者画得很薄只作为几何说明不参与周期调制。我在COMSOL中默认画法是一个二维矩形代表低折射率层上面叠一个高折射率层再叠低折射率层组成一个周期单元然后在单元外面加一个硅基底薄片。这个薄片只在垂直方向加厚周期方向上仍然和单元边界对齐。这样做有两个好处一是后续如果要算反射率或者透射率可以直接在基底边缘加端口二是如果只想看本征模式可以把基底厚度设置成可调参数扫一下看看基底厚度对带隙数值的影响小不小以此判断模型是否需要包含它。2.3 Floquet周期边界条件和特征频率研究的配置思路COMSOL里做能带计算最核心的两个设置是“周期性条件”和“特征频率研究”。周期性条件选择Floquet周期性边界后软件会要求你指定波矢 k 的方向和大小。对一维光子晶体波矢方向就是晶格周期方向k 的大小从0标到 π/a。这里有坑COMSOL的Floquet边界条件输入的波矢分量需要是实际数值不是归一化值不同软件约定不一致切换软件时特别容易搞错。特征频率研究出来的是本征频率。理论上对一个k点做特征频率求解能得到这个波矢下所有允许的本征模式频率。这些模式包括了电磁波在周期结构中的所有传播模式其中一部分是我们要的物理模式还有一些可能是数值模式或者集中在边界的模式。为了区分它们我通常先算波长范围内的前几个模式再在能带图上观察连续性。真正的物理模式在相邻k点之间应该平滑变化而数值模式往往出现得很突兀。这里补充一个求解的小细节如果想在一条能带曲线上采样20个点不需要手动切20次用参数化扫描就行。在研究中把 k 设成参数从0到π/a扫过去每步都做一次特征频率求解最后把所有结果汇总就能一次画出完整的能带图。我一般扫描步长取21个点也就是步进π/20a精度和计算量比较平衡。多取一些点能看清带边缘的细节少取一些则响应更快显然后续工作都是以快速迭代为先。3. 手把手实操COMSOL 6.x里从零搭建一维光子晶体3.1 几何搭建和材料指派下面所述全是基于我自己的操作经验版本是COMSOL 6.2/6.4不同版本菜单名称会有些差异但逻辑基本一致。打开COMSOL新建模型选择二维空间维度。物理场选择“电磁波频域”Electromagnetic Waves, Frequency Domain这属于RF模块或波动光学模块。如果你只有AC/DC模块那做不了光学频段仿真这一点要先确认。全局参数里填上a 400nmf 0.35n_H 3.48n_L 1.45。然后厚度用表达式定义d_H f*ad_L (1-f)*a。材料部分不需要从材料库选复杂的色散模型直接在“空材料”里手动输入相对介电常数和相对磁导率。对无损近似相对介电常数就是折射率的平方ε_r_H n_H²ε_r_L n_L²。硅在1550nm附近的损耗很小初始计算先用无损模型等后续需要严格模拟吸收损耗时再加入虚部。几何就画三个矩形按顺序上下排列第一个矩形宽a、高d_L放在底部第二个矩形宽a、高d_H叠在第一个上面第三个矩形再叠一层d_L。这就是一个完整的周期单元。为什么用三层而不是两层因为一维光子晶体的一个周期严格说包含高折射率层和低折射率层各一层两层矩形就够了但COMSOL里画三层在视觉上更好看出周期性而且如果你后面想做缺陷腔中间插入缺陷层也方便。不过在建模里两层矩形就足以定义一个周期单元三层也是一样的物理。我倾向于用“低—高—低”三层结构这样顶部和底部边界条件都对称提取模态更直观。3.2 边界条件设置和网格策略边界条件有两个层次。一是周期性条件给模型的左右两个边界面设置Floquet周期边界k向量设为(kx, 0)其中kx就是扫描参数。二是上下边界在严格无限周期假设里上下边界也可以用周期条件但这样模型等价于无限厚度周期堆叠物理上没问题如果你希望模拟“有限厚度堆叠两侧自由空间”那就不能把上下边界设成周期而要设置散射边界条件或完美匹配层。这里其实藏着初学者的一个大坑。很多教程里画一层晶胞四个边全设成周期条件算出来的能带确实干净但这个模型只代表“无限大周期材料”不是“硅基底上有限厚度的光子晶体”。两种模型算出来的禁带中心位置几乎一样但禁带深度、带边缘形状会有细微差别。我建议你第一次做能带计算时采用“四周全周期”的最简单模型把能带结构本身搞清楚后面做传输仿真时再换包含基底的有限模型。分清这两个模型的差别能帮你省掉不少排查问题的时间。网格设置上高折射率硅层内的光波波长约为 λ0/3.48如果目标是1550nm硅内波长约445nm。想解析电磁模式的振荡网格尺寸至少要在介质内波长的十分之一到八分之一也就是大约40~50nm。低折射率层因为波长更长网格可以稍微松一点。二维模型用映射网格把周期单元切成规整的矩形网格计算效率和精度都最好。我一般先在每条边上手工设网格数高折射率层竖直方向设4个单元低折射率层竖直方向设6个单元水平方向整个宽度设16个单元。这个密度对于一维能带计算来说已经很保守结果出来会非常平顺。3.3 特征值扫描和k向量参数化研究类型选“特征频率”。在研究中开启“参数扫描”扫描对象就是我们定义的k_x。这里再强调一次k_x的单位是1/m在COMSOL默认米制单位下450nm的周期对应k_x最大值是π/(450e-9) ≈ 6.98e6 rad/m。扫描范围写0到这个值扫描步数设置成20。有些版本允许直接把扫描范围写成0到pi/a的表达式a是全局参数这样可以随时改周期不用重设扫描范围。特征频率求解有“指定特征值搜索”的选项。我不建议直接空着让求解器瞎找因为容易找到一堆高频数值模式。更稳的做法是先估算目标频率范围。以1550nm目标波长频率约为193THz。在搜索范围里把目标值设成193e12 Hz左右再设置一个合适的搜索半径。这样求解器会优先收敛在物理模式附近。之后对每个k点基本都能扫到前6~8条能带足够看清第一带隙了。求解器内部设置里我习惯把“特征值数”稍微调大一点比如一次算10个特征值而不是默认的6个。原因很简单有些模式是简并的或者某些模式在带边缘变化剧烈数量设太少可能漏掉关键模式画出来的能带图断断续续。多算几个后面筛选时再删数值模式比漏算后重新扫要快得多。3.4 能带图后处理和带隙判断扫描完成后把所有特征频率结果汇集到一个图上。画图时要注意坐标转换横轴是波矢但COMSOL直接画的话不会自动变成很好看的一条条曲线你要用“一维绘图组”横轴选择k_x纵轴选择特征频率的绝对值。如果模型里有两个相似模式它们的频率几乎一样画出来就是两条几乎重合的曲线这样的“双带结构”在一维光子晶体中是正常现象不用紧张。归一化频率的处理方法是建立衍生值。我通常定义 f_norm freq * a / c0这样纵轴无量纲化。归一化之后的第一带隙中心大约在0.5附近取决于布拉格条件。你在能带图上看到某一段频率区间没有任何曲线穿过那就是禁带。斜着看这条带隙的上下边缘下边缘的频率是多少上边缘是多少它们的差除以中心频率就是相对带隙宽度。我强烈建议每次算出能带图顺手把带隙上下边缘记录下来放进一个表格再改晶格常数或占空比多跑几组对比。你会发现带隙中心几乎跟着λ_c的公式走而带隙宽度由折射率对比度主导。这个记录过程会让你对参数规律产生直觉以后设计新结构时根本不用盲扫参数。4. 常见问题与排查记录4.1 为什么特征频率扫出来一堆“幽灵模式”第一次做能带计算最常见的现象是能带图里冒出大量杂乱的曲线而且这些曲线在相邻k点之间完全不连续相位跳来跳去。这大概率不是物理问题而是求解器收敛到了非物理数值模式。原因主要有三个。一是搜索频段太宽。如果你把特征值搜索范围设置成0到1000THz那解出来会有大量高频模式其中很多是边界模式或数值寄生模式。解决方法是把搜索中心设置到目标频率附近搜索半径控制在100THz左右。二是模型里有多余的自由度。如果你画了几何域但没指认材料软件默认空气填充结构就变成“空气包围硅条”的三维散射问题模式形态和纯周期介质完全不同。三是网格粗细差别太大导致某些高梯度场没有解析。按我前面说的映射网格一般就不会有这个问题。如果确实出现了疑似幽灵模式一个快速的判断技巧是看电场分布。物理模式会呈现布洛赫波的空间调制场分布沿周期方向呈现规律性强弱变化幽灵模式往往出现场集中在小区域、高幅值振荡的情况。在COMSOL后处理里把几个模式挨个画出来看一遍一目了然。4.2 Floquet边界条件报错或结果异常Floquet边界条件的报错最常见的是波矢方向和边界组不匹配。一维光子晶体周期方向是x你要确保施加Floquet条件的左右两条边界正好代表的是一对周期对应的边。如果模型是从CAD导入的边界可能存在碎线导致周期条件没有完整覆盖这时候要手动合并边界。另一个常见问题是k向量写错。COMSOL中k向量要写成实数波矢分量(kx, ky)单位是1/m不是归一化单位。假设有人从文献抄了个带隙图横轴写的是0到0.5直接把这个数值填进COMSOL那算出来的能带是完全错的因为实际kx应该比这个小好几个数量级。解决办法很简单用表达式 π/a 作为扫描上限而不是写死数字。周期性条件下还有一个细节如果结构中包含硅基底而基底超出周期单元边界Floquet边界条件会把不需要周期延展的基底也周期化导致结果变差。处理方式是确保“周期域只包含周期堆叠部分”基底要么不画要么放在模型一侧并在相对面上也设置合理的边界条件。我说过能带计算本身就是无限周期体系的数学抽象加基底只会多出不必要的模式第一次建模不加就完事了。4.3 算出来的带隙和理论透射率对不上这种情况我遇到太多次了明明能带图显示有大带隙但用同一结构做透射率仿真却看不到明显的阻带。原因一般是两种模型的条件不同。能带模型是无限周期、光子从任意波矢方向传播透射率模型是有限层数、光垂直入射。如果周期数太少带隙处的反射率可能只有50%看起来当然不够“禁”。要验证这两者的一致性建议先保证周期数足够。一维光子晶体的反射率随着周期数增加而提高一般8到10个周期才能让阻带透射率降到-20dB以下。能带图的带隙宽度是从无限周期理论推导的透射率谱的阻带宽度在周期数不足时会“缩水”。不要用3周期结构去验证理论带隙那必然对不上。还有一点透射率仿真里要用端口激励且端口尺寸、模式阶数都要一致。如果端口设置的是TEM模或平面波垂直入射而能带图里禁带对应的是布洛赫模两边的本征模式不完全匹配透射率阻带的边缘也就会有轻微偏移。这个偏移通常在几个纳米到几十纳米取决于折射率对比度。4.4 我建议你避开的三个坏习惯做一维光子晶体能带计算折腾了多年之后我总结出三个最容易让人返工的坏习惯分享出来给你避坑。第一个是“参数全写死”。几何尺寸、材料折射率全部硬编码在几何里改一次参数要重新画图。第一次跑通可能觉得省事但后面的对比扫描会非常痛苦。一定要把晶格常数、占空比、厚度都定义成全局参数几何直接用表达式。这样扫占空比时只需要改一个数COMSOL自动重建几何和网格。第二个是“不记录探索过程”。跑了几十组参数最后连哪组参数对应哪张能带图都分不清这是最常见的科研事故。每次参数扫描把结果导出成CSV或者截图标上参数值。我自己习惯在文件命名里写全关键参数比如“a400_f035_dH109_dL26”三个月后回看依然清清楚楚。第三个是“只看能带不验证”。能带图算完就结束了这是很多教程的常态但这不是工程的做法。至少应该再做一次透射率仿真或者文献对比确认带隙边缘和理论值一致。如果不一致就回头检查模型。能带图的数值和其他方法算出来的结果完全一致才能说明你的模型没有隐含错误。一步验证能省下的时间远超多跑一次仿真的花费。5. 从能带计算到器件设计一维光子晶体还能怎么玩5.1 引入缺陷层在禁带里“凿”出一个谐振峰一维光子晶体最大的可玩性在于缺陷。你在完整的周期堆叠中间插入一层厚度异常或折射率异常的缺陷层就会在禁带里形成局域模式。这个模式的频率落在原来禁带的范围内光子不能穿过完美周期区却可以局域在缺陷层附近振荡。一维缺陷态就是一台微腔谐振器在硅基光子学里可以用来做窄带滤波器、激光器腔体或传感器。要在能带仿真中看到缺陷态最简单的方法是把几何扩展成“缺陷层夹在两个周期堆叠之间”的超晶胞。超晶胞的横向周期变大了第一布里渊区变小禁带中会出现一条或者几条平坦的缺陷模能带。这条能带的群速度趋近于零对应强烈的慢光效应是传感器件的核心工作点。我这里只提醒一点超晶胞的能带图和单晶胞能带图积分后是一样的只是折叠方式不同画图时不要被曲线数量增多吓到。5.2 用参数扫描和优化自动找最佳带隙手动扫参数只能试探几个点。当你需要系统性地寻找最大带隙时建议把COMSOL和MATLAB或Python联起来跑。COMSOL支持通过LiveLink for MATLAB或者Java脚本批量修改参数、执行研究、提取结果。把这套自动化跑起来后你可以一次性扫上百组占空比和晶格常数组合把相对带宽和带隙中心做成热图肉眼就能找出最优区域。我做过多组扫描后发现一维Si/SiO₂结构的相对带隙通常在20%到30%之间想做得更宽就得引入更高对比度的材料比如硅和空气、锗和空气或者改变结构维度。一维能带计算的价值就在这里它能用极低成本快速筛选材料组合和结构参数筛完再用三维全波电磁仿真验证最终器件整个设计流程的效率比盲目建完整模型高太多了。5.3 能带知识迁移到二维/三维光子晶体一维光子晶体的能带计算学会后迁移到二维光子晶体比如三角晶格空气孔平板基本只差两步一是把模型从二维换成三维或对二维晶格使用面外传播模型二是布里渊区扫描路径从一条线变成多边形——典型的是Γ-M-K-Γ闭合路径。Floquet周期边界条件和特征频率扫描的方法完全一样。这正是一维结构作为入门训练的价值。你在这个项目里搞懂了Bloch定理在软件里是怎么落地的搞懂了k向量扫描与布里渊区之间的关系也搞懂了如何区分物理模式和数值模式那后续所有周期结构仿真都只是几何和物理场替换的问题核心流程你已经完全掌握了。我见过不少一开始直接挑战二维光子晶体结果卡在Floquet边界条件和k路径设置上好几周的人。如果先花一天把一维跑通二维项目至少能砍掉一半的调试时间。我个人在实际操作中还有个体会一维光子晶体能带计算这类仿真真正的难点从来不是软件操作而是模型抽象。周期性边界条件、k向量扫描、特征频率筛选这些概念听起来很理工男但本质上和你在高速上判断“这个收费杆抬不抬”没太大区别。模型建对了算法自己会给你漂亮的带隙图模型建错了再强的求解器也救不了结果。建议你跑通这个案例后把参数记录表留好然后把晶格常数、占空比、材料折射率都换成自己的目标值重新扫一遍看看带隙怎么跟着参数走。这个“自己动手换一次”的过程远比我写上十页教程都管用。