☰
COMSOL仿真一维光子晶体能带:原理、参数与禁带调控
2026/10/2 22:29:50 网站建设 项目流程

一维光子晶体这东西,名字挺唬人,其实就是把两种折射率不同的介质按固定“节拍”交替排列,光走进去就不再是闷头直冲,而是会被这排周期结构反复筛选——有些频率顺利放行,有些频率直接被弹回来。我习惯拿高速公路收费站来打比方:一排收费杆按固定间距立着,不同频率的光子就像不同车速的汽车,杆间距就是晶格常数,杆子之间的空隙就是介质占空比,能不能抬杆放行,全看这套周期性介电结构是怎么设计的。这篇文章要干的事很具体:用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. 建模之前的参数设计,为什么先算一步再动COMSOL

2.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) ≈ 111nm,d_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 = 400nm,f = 0.35,n_H = 3.48,n_L = 1.45。然后厚度用表达式定义:d_H = f*a,d_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向量扫描、特征频率筛选,这些概念听起来很理工男,但本质上和你在高速上判断“这个收费杆抬不抬”没太大区别。模型建对了,算法自己会给你漂亮的带隙图;模型建错了,再强的求解器也救不了结果。建议你跑通这个案例后,把参数记录表留好,然后把晶格常数、占空比、材料折射率都换成自己的目标值重新扫一遍,看看带隙怎么跟着参数走。这个“自己动手换一次”的过程,远比我写上十页教程都管用。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询