光子晶体这名字听起来唬人说白了就是把手头能拿到的几种不同折射率材料按固定规律一层层码起来让光在里面走的时候像过收费站——能放行的频段一路畅通不能放行的频段直接劝返。一维的情况最简单也最贴近工程DBR反射镜、光纤光栅、各种高低折射率交替膜层底层都是同一套物理。今天我用COMSOL把硅基底上一维周期性介电结构的能带计算完整跑一遍从边界条件到后处理全部复盘。整个过程不写一行代码但每个设置为什么这么做我会尽量讲透。想入门光子晶体仿真或者正在做多层膜、光栅设计却不知道怎么出能带图的朋友这篇可以直接照着抄。1. 光子能带到底在算什么先搞懂“收费规则”1.1 一维光子晶体的物理本质一维光子晶体最典型的形态就是多层膜两种折射率不同的介质交替堆叠形成一个周期性介电结构。光在其中传播时会在每个界面上发生反射这些反射光之间会产生干涉。当每个周期的光学厚度接近入射波长的四分之一时反射光同相叠加形成很强的反射——这就是布拉格反射。反射强到一定程度就会形成一个频率范围让光完全无法在结构里传播这个频率范围就是光子带隙。这个现象和电子在晶体中运动的能带理论可以一比电子受周期性原子势场调制形成能带和禁带光子则受周期性介电常数调制形成光子能带和光子禁带。所以“光子晶体”这个名字并不是碰瓷它背后的数学结构和电子能带理论是高度类似的。在COMSOL里做能带计算本质上不是把一堆频点挨个扫一遍而是在一个周期单元上求解本征值问题。把求解频率作为本征值把周期方向上的波矢k作为输入参数每给定一个k就求出一组可传播的模式频率。把k扫完把所有模式频率画在一起就是能带图。1.2 布拉格条件与带隙中心频率的快速估算动手建模之前最好先拿笔算一个大概的频率位置后面好用来验证COMSOL的结果靠不靠谱。对于一维交替多层膜如果两种材料的光学厚度满足n1 × d1 n2 × d2 λ0 / 4那么带隙中心波长λ0和中心频率f0分别为λ0 2 × (n1 × d1 n2 × d2) f0 c / λ0我这次选的模型参数是晶格常数a 500 nm其中硅层厚度d1 200 nm二氧化硅层厚度d2 300 nm。硅在近红外的折射率取3.45二氧化硅取1.44。先算一下光学厚度n1 × d1 3.45 × 200 690 nm n2 × d2 1.44 × 300 432 nm所以布拉格中心波长λ0 2 × (690 432) 2244 nm对应的中心频率f0 c / λ0 ≈ 3 × 10^8 / 2.244 × 10^-6 ≈ 1.34 × 10^14 Hz这个值一会儿用来对照带隙位置非常有用。如果你的模型频率单位是THz或rad/s记得先换算统一这是新手最容易栽的坑。1.3 为什么能带图要沿布里渊区边界扫波矢能带图的横轴是波矢k纵轴是频率f。一维光子晶体的周期方向只有一个第一布里渊区就是从-π/a到π/a。由于能带图是偶函数实际计算时只扫0到π/a就够。布里渊区边界k π/a处布拉格条件满足带隙一定在这里打开。理论上无限周期的多层膜在任何位置、任何频段都可能存在特征解我们只需要关注一个周期单元内、满足布洛赫周期条件的特征模式。COMSOL里用Floquet周期边界来施加这个布洛赫相位条件这是整个仿真中最关键的一步。2. 动手前的方案设计2D模型、物理场和单胞选择2.1 空间维度与物理场接口怎么选一维光子晶体虽然是周期方向只有一维但我不建议真的用1D模型去搭。COMSOL的“电磁波频域”接口在2D空间维度里用起来最顺手原因有两个一是2D模型能直观看到电场模式分布方便后续判断模式性质二是如果要扩展成带基底的超胞模型或多周期透射计算2D几何直接就能改。物理场接口选“电磁波频域ewfd”研究类型选“特征频率”。这里要注意一个很多人问过的问题能带计算不是用频域研究去扫描频率而是用特征频率研究去解本征值。你把波矢k作为输入参数特征频率作为输出这才是能带图的正确打开方式。在2D模型里COMSOL默认的求解模式是面外电场EzTE极化也支持面内电场TM。一维多层膜对TE和TM的反应不同带隙位置会有差异。新手入门建议先算TE极化也就是默认的电场垂直纸面方向物理图像更简单。2.2 单胞几何怎么切分才不会出错单胞的选择是能带仿真最容易翻车的地方。原则只有一条沿周期方向一个单胞平移复制之后必须能无缝拼回完整结构。我的模型这样设计矩形单元宽度W 300 nm沿x方向高度a 500 nm沿y方向即周期堆叠方向。单胞内包含完整的硅层和完整的二氧化硅层各一层。注意一个问题上下边界处的材料必须是同一层介质这样Floquet周期边界连接上下边界时几何和材料才能自然吻合。我上面这个参数上边界是硅层的顶面下边界也是硅层的底面周期延拓后正好是层叠结构。如果你要做一个周期内包含三种材料的结构比如高-低-高单胞高度就要取一个完整周期的总厚度边界两侧也必须是同种材料。切边界的时候多留一个心眼后面建模会省很多事。2.3 材料参数的取法直接用手输还是用内置材料库硅和二氧化硅在近红外波段的折射率COMSOL内置材料库里可以找到但不同版本、不同波长下数值略有差异。我对能带计算的建议是用用户自定义材料直接手动输入折射率最稳妥后面调参也方便。内置材料“Si (Crystalline)”给出的折射率在1550 nm附近约3.48我实际用手输入的常数值n 3.45这样和手算的布拉格条件能对上。做能带计算阶段材料损耗折射率虚部可以不填填了反而会把本征频率变成复数后处理时会多出很多麻烦。等研究透射率、损耗特性时再加虚部也不迟。在建立全局参数时我会把折射率、厚度、晶格常数都设为参数后续优化周期结构尺寸时只改参数表不需要动几何和材料。3. COMSOL完整实操从空模型到能带图3.1 全局参数定义打开COMSOL新建一个2D模型物理场选择“电磁波频域”研究选择“特征频率”。进入模型开发器后第一步是把全局参数表建好。我习惯把参数分成三类几何参数、材料参数、扫描参数。具体如下参数名表达式说明a500[nm]周期厚度d1200[nm]硅层厚度d2300[nm]二氧化硅层厚度W300[nm]单胞宽度n_Si3.45硅折射率n_SiO21.44二氧化硅折射率eps0_Sin_Si^2硅相对介电常数eps0_SiO2n_SiO2^2二氧化硅相对介电常数ky0周期方向布洛赫波矢扫描用ky这个参数先设为0后面在“研究”里做参数扫描时它会作为扫描变量被替换。3.2 几何搭建两个矩形拼一个晶胞在几何节点下创建两个矩形矩形1是硅层宽度W高度d1左下角坐标(0, d2)这样它正好坐在二氧化硅层上面。 矩形2是二氧化硅层宽度W高度d2左下角坐标(0, 0)。两个矩形拼起来之后整个单胞高度正好是a 500 nm。宽度W取300 nm纯属为了演示因为结构在x方向是无限均匀的W不影响垂直入射模式的频率但会影响横向高阶模的位置。W取得太大低频段会出现一堆面内高阶横模能带图会变得乱七八糟取得太小比如几十纳米网格数量又上去了。300 nm是一个比较平衡的选择。创建完几何后在“形成联合体”节点里做一次布尔并集把两个矩形合成一个域。这个步骤看起来多余但如果不合并后续材料分配和边界条件选择时容易误选到内部界面。3.3 材料分配直接给域指定折射率在“材料”节点下新建一个空材料命名为“Si”然后在“折射率”子节点中把相对介电常数设置为eps0_Si也就是n_Si的平方。同样再建一个“SiO2”材料相对介电常数设为eps0_SiO2。之后把Si材料赋给上半部分的域也就是硅层把SiO2材料赋给下半部分的域。如果你不确定哪个域对应哪层可以在几何节点里给每个矩形起一个有意义的名字比如“SiLayer”和“SiO2Layer”。COMSOL的电磁波频域接口默认求解的是相对介电常数所以直接填折射率的平方就行。如果你填的是折射率实部求解器也能工作但电场和功率的物理意义会有点别扭。我一般还是填介电常数少一层换算的麻烦。3.4 周期边界条件两个方向分开处理边界条件是能带计算的灵魂。我先说结论左右边界x方向用“周期条件”类型选普通周期条件模拟面内无限均匀。上下边界y方向用“Floquet周期条件”布洛赫波矢设为(0, ky)ky为扫描参数。为什么要分开因为物理上一维光子晶体的周期只存在于堆叠方向即y方向。我们算能带时关心的是沿堆叠方向能不能传所以y方向必须加布洛赫相位。x方向是面内的连续平移对称方向使用普通周期条件就够表示电场在x方向严格周期重复这也正好对应平面波垂直入射的情况。在COMSOL里的具体操作是选中“电磁波频域”节点右键添加“周期条件”。然后分别选择左右边界和上下边界。对于上下边界在周期条件的设置窗口中把“类型”选为“Floquet周期”然后在“Floquet波矢”一栏输入kx 0ky ky。注意这里输入的ky会引用全局参数ky扫描时它会被参数扫描替换。有一个细节很多人会忽略Floquet周期条件的源边界和目标边界必须对应正确的几何位置。COMSOL会自动判断但如果你发现解出来模式在边界处突变先检查一下是不是源和目标选反了。3.5 特征频率研究设置与参数扫描在研究节点下找到“特征频率”研究步骤。默认情况下COMSOL只搜索前几个特征频率。我建议把“所需特征频率数”设为8搜索范围的下限设为0上限设为3×10^14 Hz。这个范围覆盖了前面手算的中心频率1.34×10^14 Hz应该能看到至少两个带。接着在研究中添加“参数扫描”扫描变量选择全局参数ky扫描范围从0到pi/a。步长建议取pi/a的20等分也就是0.05 × (pi/a)大概能画出平滑的能带曲线。扫描24个点、每个点解8个特征频率这个计算量在普通笔记本电脑上也就一两分钟。如果发现速度慢先检查是不是把网格建得太细了。一维结构用映射网格就能很好处理。计算完成后COMSOL会把所有ky值对应的特征频率都列出来。这里要提醒一句不同ky下特征频率会自动按频率大小排序但在带交叉区域这种排序可能会造成能带分支看起来不连续。后面画图时如果发现曲线跳变不要太惊讶这是正常现象。3.6 后处理画能带图的正确姿势在“结果”节点下新建一个“一维绘图组”数据源选择“特征频率”研究。横轴表达式填ky / (pi/a)纵轴填f。如果你希望纵轴用归一化频率a/λ那就是f × a / c_const。这里c_const是COMSOL内置的光速常量。画出来之后你会看到几条从左往右上升的曲线并在某个频率区间出现空白区域这个空白区域就是光子带隙。为了让图更专业建议把横轴标签改成“ky (π/a)”纵轴改成“Frequency (Hz)”或“a/λ”。然后在绘图组属性里设置“x轴标签”和“y轴标签”坐标轴范围也可以根据结果稍微裁剪一下。如果适合可以在同一张图里加一个“全局计算”节点把带隙上下边缘的频率用水平线标出来方便直接读数。我个人习惯把带隙上边缘和下边缘的频率值用表格形式复制出来跟手算值对照。4. 结果解读与交叉验证能带图不骗人但别急着信4.1 能带图长什么样才算正确按我给的参数跑完你应该能看到类似这样的趋势第一带和第二带在k π/a处出现频率分离形成带隙带隙中心频率大约在1.34×10^14 Hz附近。这就是布拉格反射条件得到满足的位置。能带图一个常见的误读是只要你看到有空白区域就以为是带隙。不对。有些空白区域是模式密度低造成的伪带隙有些只对特定极化方向存在。判断真正的带隙要看它是否在所有波矢方向上都存在频率空白。对于一维光子晶体只要k从0扫到π/a带隙内没有实频模式就可以认为存在完全带隙。如果你把TE和TM两种极化都算了两种极化下带隙位置有重叠那在实验上才能实现真正全极化的截止频率。4.2 用布拉格公式验证数值结果我现在最喜欢做的事就是用解析公式去“打假”仿真结果。把COMSOL算出的带隙中心频率和手算的布拉格频率对照一下。刚才手算f0 ≈ 1.34×10^14 Hz。如果COMSOL给出的带隙大致在这个频率附近上下张开说明模型没有问题。偏差在百分之几以内都属于正常因为布拉格公式是中心频率的估算带隙并不完全对称。如果偏差很大先检查折射率有没有填对再做一次几何检查确认硅层和二氧化硅层的厚度不是反的。这个低级错误我踩过一次跑出来的带隙位置差了快一倍排查了半天才发现是材料赋反了。4.3 交叉验证用频域透射谱看带隙能带图只是告诉你理论上有带隙但实际“光能不能透过去”最好再做一个频域透射仿真来验证。我习惯另外建一个2D模型硅基底放在最下面厚度取1到2微米基底上交替堆叠8个周期的硅/二氧化硅层最上面是空气区域。这个模型更贴近“硅基底上的周期性介电结构”这个实际场景。物理场同样用“电磁波频域”但研究类型选“频域”。上下边界用散射边界条件在底部施加一个平面波激励从基底一侧入射扫频范围比如1×10^14到1.7×10^14 Hz然后计算顶部边界的透射功率。透射谱上你会非常清晰地看到在能带图预测的带隙频率范围内透射率断崖式下跌可能下降好几个数量级。带隙外的频率透射率则保持高位。这一步验证做完你才算真正把能带计算落地到一个可以跟实验对比的结果。4.4 如果非要带上硅基底算能带超胞玩法有时候你确实想直接看“硅基底上周期结构”的向导模式或表面态而不是理想周期堆叠的体能带。这时可以用超胞法把几何改成一大块硅基底加上几个周期层作为一个超大周期单元。超胞的y方向高度包含基底高度加周期结构高度然后y方向仍然用Floquet周期边界。但由于基底不是周期性的这种方法得到的模式会包含基底导模能带图里会出现很多近似平坦的曲线这些对应基底的波导模式。能带图里落在带隙区域内、但电场分布明显局域在周期结构表面的模式才是你要找的表面态或局域界面态。超胞法对模式数目的要求更高一般要把特征频率数扩大到20到30个计算时间也会线性增加。新手建议先把无基底的体能带跑熟再碰超胞。5. 常见问题与排查记录这五个坑我最常踩5.1 特征频率全是0或者重复这个一般出在周期边界条件没设置对或者几何域没有正确划分。检查一下是否所有边界都被正确选中尤其注意内部边界不能多选否则求解时会出现零能模式。还有一种情况是全局参数ky没有进入Floquet波矢导致所有扫描点的边界条件一模一样特征频率自然全部重复。去Floquet周期条件设置里确认波矢表达式引用的参数名是“ky”而不是一个写死的数字。5.2 能带图曲线严重跳变、分支交叉混乱特征频率研究在每个ky点都是单独求解的COMSOL会自动按频率升序排序但在模式交叉区域排序会让分支看起来突然断裂。解决办法有几种一是降低特征频率搜索范围让曲线上只有你关心的那些模式避免高频杂散模式干扰排序。 二是利用COMSOL的“模态扩展”或“特征值追踪”功能不过这个功能在不同版本里位置不太一样需要自己找一下。 三是把数据导出去在MATLAB或者Python里按场分布相似度做模式追踪和排序。我后来处理高折射率对比度的结构时都是直接导出数据用脚本排序。5.3 网格怎么画才能又快又准一维周期结构强烈建议用映射网格。在网格节点下选择“映射”然后把网格尺寸设为手动最大单元尺寸设为10 nm左右。为什么是10 nm硅和二氧化硅的折射率大约3.45和1.44当频率上限为3×10^14 Hz时对应的材料内最短波长约为λ_min c / (f_max × n_max) ≈ 3×10^8 / (3×10^14 × 3.45) ≈ 290 nm经验法则是一个波长内至少要有10个以上网格点所以30 nm以下都能接受。但Floquet周期边界要求源边界和目标边界的网格完全匹配映射网格天然满足这个要求所以我一般直接压到10 nm算得也很快整个模型不过几千个自由度。5.4 模式数不够导致带隙不明显如果你只解了2到3个特征频率很可能画出来的能带图只有一两条曲线带隙根本看不出来。我建议起步8个如果结构折射率对比度高或者你加了基底再多加几个。搜索频率上限也要跟上。如果你只搜到2×10^14 Hz而第二带在2.5×10^14 Hz处才打开那你看到的只是第一带带隙自然不存在。别偷懒先按我给的3×10^14 Hz试一次再把范围缩小到你真正关注的频段。5.5 归一化频率单位换算搞混画图时如果选择纵轴为a/λ那是一个无量纲数。我算的例子中a 500 nm中心波长2244 nma/λ约等于0.223。如果看到能带图的纵轴数值在0.2附近说明结果合理。如果你用的是频率f直接作图注意COMSOL里特征频率的默认单位是Hz但在某些版本界面里会显示为THz。这个单位显示问题曾经让我一度以为结果差了1000倍。整理成一张速查表问题可能原因排查方向特征频率为0周期边界未生效检查Floquet波矢参数能带图跳变模式排序不稳定减少搜索范围/外部重排带隙位置偏差大材料参数或几何厚度填错对照布拉格公式手算验证曲线看起来乱横向高阶模介入缩小单胞宽度W计算极慢网格过细先用映射网格粗网格试跑这几条基本覆盖了我自己从零开始算光子晶体能带时遇到的绝大部分问题。每次换新材料、换参数我都会先在脑子里过一遍这张表能省下大量调试时间。6. 一个容易被忽略但很实用的收尾技巧最后分享一个我自己的小习惯每跑完一个能带模型我会顺手在全局参数里改一下占空比也就是硅层厚度占周期厚度的比例比如把d1从200 nm改成150 nm再重新算一遍能带观察带隙宽度怎么变化。这个操作几分钟就出结果但对理解和优化光子晶体结构非常值。带隙宽度通常不是占空比的线性函数存在一个最优值。手算布拉格条件只能告诉你带隙中心在哪带隙到底有多宽必须靠数值扫描。我算过一组硅/二氧化硅一维结构占空比从0.2变到0.5时带隙宽度先增后减在0.4附近达到最大。这种规律不实际跑一遍参数扫描很难凭直觉猜准。如果你后续要做更复杂的结构比如二维光子晶体、带缺陷的波导或者想引入非线性材料这个能带计算流程完全可以平滑迁移过去。核心逻辑不变单胞、Floquet周期边界、特征频率扫描只是几何和网格需要跟着升级。希望这份实操记录能帮你少绕几次路一次把能带图画明白。
COMSOL一维光子晶体能带计算实操：从多层膜结构到带隙分析
光子晶体这名字听起来唬人说白了就是把手头能拿到的几种不同折射率材料按固定规律一层层码起来让光在里面走的时候像过收费站——能放行的频段一路畅通不能放行的频段直接劝返。一维的情况最简单也最贴近工程DBR反射镜、光纤光栅、各种高低折射率交替膜层底层都是同一套物理。今天我用COMSOL把硅基底上一维周期性介电结构的能带计算完整跑一遍从边界条件到后处理全部复盘。整个过程不写一行代码但每个设置为什么这么做我会尽量讲透。想入门光子晶体仿真或者正在做多层膜、光栅设计却不知道怎么出能带图的朋友这篇可以直接照着抄。1. 光子能带到底在算什么先搞懂“收费规则”1.1 一维光子晶体的物理本质一维光子晶体最典型的形态就是多层膜两种折射率不同的介质交替堆叠形成一个周期性介电结构。光在其中传播时会在每个界面上发生反射这些反射光之间会产生干涉。当每个周期的光学厚度接近入射波长的四分之一时反射光同相叠加形成很强的反射——这就是布拉格反射。反射强到一定程度就会形成一个频率范围让光完全无法在结构里传播这个频率范围就是光子带隙。这个现象和电子在晶体中运动的能带理论可以一比电子受周期性原子势场调制形成能带和禁带光子则受周期性介电常数调制形成光子能带和光子禁带。所以“光子晶体”这个名字并不是碰瓷它背后的数学结构和电子能带理论是高度类似的。在COMSOL里做能带计算本质上不是把一堆频点挨个扫一遍而是在一个周期单元上求解本征值问题。把求解频率作为本征值把周期方向上的波矢k作为输入参数每给定一个k就求出一组可传播的模式频率。把k扫完把所有模式频率画在一起就是能带图。1.2 布拉格条件与带隙中心频率的快速估算动手建模之前最好先拿笔算一个大概的频率位置后面好用来验证COMSOL的结果靠不靠谱。对于一维交替多层膜如果两种材料的光学厚度满足n1 × d1 n2 × d2 λ0 / 4那么带隙中心波长λ0和中心频率f0分别为λ0 2 × (n1 × d1 n2 × d2) f0 c / λ0我这次选的模型参数是晶格常数a 500 nm其中硅层厚度d1 200 nm二氧化硅层厚度d2 300 nm。硅在近红外的折射率取3.45二氧化硅取1.44。先算一下光学厚度n1 × d1 3.45 × 200 690 nm n2 × d2 1.44 × 300 432 nm所以布拉格中心波长λ0 2 × (690 432) 2244 nm对应的中心频率f0 c / λ0 ≈ 3 × 10^8 / 2.244 × 10^-6 ≈ 1.34 × 10^14 Hz这个值一会儿用来对照带隙位置非常有用。如果你的模型频率单位是THz或rad/s记得先换算统一这是新手最容易栽的坑。1.3 为什么能带图要沿布里渊区边界扫波矢能带图的横轴是波矢k纵轴是频率f。一维光子晶体的周期方向只有一个第一布里渊区就是从-π/a到π/a。由于能带图是偶函数实际计算时只扫0到π/a就够。布里渊区边界k π/a处布拉格条件满足带隙一定在这里打开。理论上无限周期的多层膜在任何位置、任何频段都可能存在特征解我们只需要关注一个周期单元内、满足布洛赫周期条件的特征模式。COMSOL里用Floquet周期边界来施加这个布洛赫相位条件这是整个仿真中最关键的一步。2. 动手前的方案设计2D模型、物理场和单胞选择2.1 空间维度与物理场接口怎么选一维光子晶体虽然是周期方向只有一维但我不建议真的用1D模型去搭。COMSOL的“电磁波频域”接口在2D空间维度里用起来最顺手原因有两个一是2D模型能直观看到电场模式分布方便后续判断模式性质二是如果要扩展成带基底的超胞模型或多周期透射计算2D几何直接就能改。物理场接口选“电磁波频域ewfd”研究类型选“特征频率”。这里要注意一个很多人问过的问题能带计算不是用频域研究去扫描频率而是用特征频率研究去解本征值。你把波矢k作为输入参数特征频率作为输出这才是能带图的正确打开方式。在2D模型里COMSOL默认的求解模式是面外电场EzTE极化也支持面内电场TM。一维多层膜对TE和TM的反应不同带隙位置会有差异。新手入门建议先算TE极化也就是默认的电场垂直纸面方向物理图像更简单。2.2 单胞几何怎么切分才不会出错单胞的选择是能带仿真最容易翻车的地方。原则只有一条沿周期方向一个单胞平移复制之后必须能无缝拼回完整结构。我的模型这样设计矩形单元宽度W 300 nm沿x方向高度a 500 nm沿y方向即周期堆叠方向。单胞内包含完整的硅层和完整的二氧化硅层各一层。注意一个问题上下边界处的材料必须是同一层介质这样Floquet周期边界连接上下边界时几何和材料才能自然吻合。我上面这个参数上边界是硅层的顶面下边界也是硅层的底面周期延拓后正好是层叠结构。如果你要做一个周期内包含三种材料的结构比如高-低-高单胞高度就要取一个完整周期的总厚度边界两侧也必须是同种材料。切边界的时候多留一个心眼后面建模会省很多事。2.3 材料参数的取法直接用手输还是用内置材料库硅和二氧化硅在近红外波段的折射率COMSOL内置材料库里可以找到但不同版本、不同波长下数值略有差异。我对能带计算的建议是用用户自定义材料直接手动输入折射率最稳妥后面调参也方便。内置材料“Si (Crystalline)”给出的折射率在1550 nm附近约3.48我实际用手输入的常数值n 3.45这样和手算的布拉格条件能对上。做能带计算阶段材料损耗折射率虚部可以不填填了反而会把本征频率变成复数后处理时会多出很多麻烦。等研究透射率、损耗特性时再加虚部也不迟。在建立全局参数时我会把折射率、厚度、晶格常数都设为参数后续优化周期结构尺寸时只改参数表不需要动几何和材料。3. COMSOL完整实操从空模型到能带图3.1 全局参数定义打开COMSOL新建一个2D模型物理场选择“电磁波频域”研究选择“特征频率”。进入模型开发器后第一步是把全局参数表建好。我习惯把参数分成三类几何参数、材料参数、扫描参数。具体如下参数名表达式说明a500[nm]周期厚度d1200[nm]硅层厚度d2300[nm]二氧化硅层厚度W300[nm]单胞宽度n_Si3.45硅折射率n_SiO21.44二氧化硅折射率eps0_Sin_Si^2硅相对介电常数eps0_SiO2n_SiO2^2二氧化硅相对介电常数ky0周期方向布洛赫波矢扫描用ky这个参数先设为0后面在“研究”里做参数扫描时它会作为扫描变量被替换。3.2 几何搭建两个矩形拼一个晶胞在几何节点下创建两个矩形矩形1是硅层宽度W高度d1左下角坐标(0, d2)这样它正好坐在二氧化硅层上面。 矩形2是二氧化硅层宽度W高度d2左下角坐标(0, 0)。两个矩形拼起来之后整个单胞高度正好是a 500 nm。宽度W取300 nm纯属为了演示因为结构在x方向是无限均匀的W不影响垂直入射模式的频率但会影响横向高阶模的位置。W取得太大低频段会出现一堆面内高阶横模能带图会变得乱七八糟取得太小比如几十纳米网格数量又上去了。300 nm是一个比较平衡的选择。创建完几何后在“形成联合体”节点里做一次布尔并集把两个矩形合成一个域。这个步骤看起来多余但如果不合并后续材料分配和边界条件选择时容易误选到内部界面。3.3 材料分配直接给域指定折射率在“材料”节点下新建一个空材料命名为“Si”然后在“折射率”子节点中把相对介电常数设置为eps0_Si也就是n_Si的平方。同样再建一个“SiO2”材料相对介电常数设为eps0_SiO2。之后把Si材料赋给上半部分的域也就是硅层把SiO2材料赋给下半部分的域。如果你不确定哪个域对应哪层可以在几何节点里给每个矩形起一个有意义的名字比如“SiLayer”和“SiO2Layer”。COMSOL的电磁波频域接口默认求解的是相对介电常数所以直接填折射率的平方就行。如果你填的是折射率实部求解器也能工作但电场和功率的物理意义会有点别扭。我一般还是填介电常数少一层换算的麻烦。3.4 周期边界条件两个方向分开处理边界条件是能带计算的灵魂。我先说结论左右边界x方向用“周期条件”类型选普通周期条件模拟面内无限均匀。上下边界y方向用“Floquet周期条件”布洛赫波矢设为(0, ky)ky为扫描参数。为什么要分开因为物理上一维光子晶体的周期只存在于堆叠方向即y方向。我们算能带时关心的是沿堆叠方向能不能传所以y方向必须加布洛赫相位。x方向是面内的连续平移对称方向使用普通周期条件就够表示电场在x方向严格周期重复这也正好对应平面波垂直入射的情况。在COMSOL里的具体操作是选中“电磁波频域”节点右键添加“周期条件”。然后分别选择左右边界和上下边界。对于上下边界在周期条件的设置窗口中把“类型”选为“Floquet周期”然后在“Floquet波矢”一栏输入kx 0ky ky。注意这里输入的ky会引用全局参数ky扫描时它会被参数扫描替换。有一个细节很多人会忽略Floquet周期条件的源边界和目标边界必须对应正确的几何位置。COMSOL会自动判断但如果你发现解出来模式在边界处突变先检查一下是不是源和目标选反了。3.5 特征频率研究设置与参数扫描在研究节点下找到“特征频率”研究步骤。默认情况下COMSOL只搜索前几个特征频率。我建议把“所需特征频率数”设为8搜索范围的下限设为0上限设为3×10^14 Hz。这个范围覆盖了前面手算的中心频率1.34×10^14 Hz应该能看到至少两个带。接着在研究中添加“参数扫描”扫描变量选择全局参数ky扫描范围从0到pi/a。步长建议取pi/a的20等分也就是0.05 × (pi/a)大概能画出平滑的能带曲线。扫描24个点、每个点解8个特征频率这个计算量在普通笔记本电脑上也就一两分钟。如果发现速度慢先检查是不是把网格建得太细了。一维结构用映射网格就能很好处理。计算完成后COMSOL会把所有ky值对应的特征频率都列出来。这里要提醒一句不同ky下特征频率会自动按频率大小排序但在带交叉区域这种排序可能会造成能带分支看起来不连续。后面画图时如果发现曲线跳变不要太惊讶这是正常现象。3.6 后处理画能带图的正确姿势在“结果”节点下新建一个“一维绘图组”数据源选择“特征频率”研究。横轴表达式填ky / (pi/a)纵轴填f。如果你希望纵轴用归一化频率a/λ那就是f × a / c_const。这里c_const是COMSOL内置的光速常量。画出来之后你会看到几条从左往右上升的曲线并在某个频率区间出现空白区域这个空白区域就是光子带隙。为了让图更专业建议把横轴标签改成“ky (π/a)”纵轴改成“Frequency (Hz)”或“a/λ”。然后在绘图组属性里设置“x轴标签”和“y轴标签”坐标轴范围也可以根据结果稍微裁剪一下。如果适合可以在同一张图里加一个“全局计算”节点把带隙上下边缘的频率用水平线标出来方便直接读数。我个人习惯把带隙上边缘和下边缘的频率值用表格形式复制出来跟手算值对照。4. 结果解读与交叉验证能带图不骗人但别急着信4.1 能带图长什么样才算正确按我给的参数跑完你应该能看到类似这样的趋势第一带和第二带在k π/a处出现频率分离形成带隙带隙中心频率大约在1.34×10^14 Hz附近。这就是布拉格反射条件得到满足的位置。能带图一个常见的误读是只要你看到有空白区域就以为是带隙。不对。有些空白区域是模式密度低造成的伪带隙有些只对特定极化方向存在。判断真正的带隙要看它是否在所有波矢方向上都存在频率空白。对于一维光子晶体只要k从0扫到π/a带隙内没有实频模式就可以认为存在完全带隙。如果你把TE和TM两种极化都算了两种极化下带隙位置有重叠那在实验上才能实现真正全极化的截止频率。4.2 用布拉格公式验证数值结果我现在最喜欢做的事就是用解析公式去“打假”仿真结果。把COMSOL算出的带隙中心频率和手算的布拉格频率对照一下。刚才手算f0 ≈ 1.34×10^14 Hz。如果COMSOL给出的带隙大致在这个频率附近上下张开说明模型没有问题。偏差在百分之几以内都属于正常因为布拉格公式是中心频率的估算带隙并不完全对称。如果偏差很大先检查折射率有没有填对再做一次几何检查确认硅层和二氧化硅层的厚度不是反的。这个低级错误我踩过一次跑出来的带隙位置差了快一倍排查了半天才发现是材料赋反了。4.3 交叉验证用频域透射谱看带隙能带图只是告诉你理论上有带隙但实际“光能不能透过去”最好再做一个频域透射仿真来验证。我习惯另外建一个2D模型硅基底放在最下面厚度取1到2微米基底上交替堆叠8个周期的硅/二氧化硅层最上面是空气区域。这个模型更贴近“硅基底上的周期性介电结构”这个实际场景。物理场同样用“电磁波频域”但研究类型选“频域”。上下边界用散射边界条件在底部施加一个平面波激励从基底一侧入射扫频范围比如1×10^14到1.7×10^14 Hz然后计算顶部边界的透射功率。透射谱上你会非常清晰地看到在能带图预测的带隙频率范围内透射率断崖式下跌可能下降好几个数量级。带隙外的频率透射率则保持高位。这一步验证做完你才算真正把能带计算落地到一个可以跟实验对比的结果。4.4 如果非要带上硅基底算能带超胞玩法有时候你确实想直接看“硅基底上周期结构”的向导模式或表面态而不是理想周期堆叠的体能带。这时可以用超胞法把几何改成一大块硅基底加上几个周期层作为一个超大周期单元。超胞的y方向高度包含基底高度加周期结构高度然后y方向仍然用Floquet周期边界。但由于基底不是周期性的这种方法得到的模式会包含基底导模能带图里会出现很多近似平坦的曲线这些对应基底的波导模式。能带图里落在带隙区域内、但电场分布明显局域在周期结构表面的模式才是你要找的表面态或局域界面态。超胞法对模式数目的要求更高一般要把特征频率数扩大到20到30个计算时间也会线性增加。新手建议先把无基底的体能带跑熟再碰超胞。5. 常见问题与排查记录这五个坑我最常踩5.1 特征频率全是0或者重复这个一般出在周期边界条件没设置对或者几何域没有正确划分。检查一下是否所有边界都被正确选中尤其注意内部边界不能多选否则求解时会出现零能模式。还有一种情况是全局参数ky没有进入Floquet波矢导致所有扫描点的边界条件一模一样特征频率自然全部重复。去Floquet周期条件设置里确认波矢表达式引用的参数名是“ky”而不是一个写死的数字。5.2 能带图曲线严重跳变、分支交叉混乱特征频率研究在每个ky点都是单独求解的COMSOL会自动按频率升序排序但在模式交叉区域排序会让分支看起来突然断裂。解决办法有几种一是降低特征频率搜索范围让曲线上只有你关心的那些模式避免高频杂散模式干扰排序。 二是利用COMSOL的“模态扩展”或“特征值追踪”功能不过这个功能在不同版本里位置不太一样需要自己找一下。 三是把数据导出去在MATLAB或者Python里按场分布相似度做模式追踪和排序。我后来处理高折射率对比度的结构时都是直接导出数据用脚本排序。5.3 网格怎么画才能又快又准一维周期结构强烈建议用映射网格。在网格节点下选择“映射”然后把网格尺寸设为手动最大单元尺寸设为10 nm左右。为什么是10 nm硅和二氧化硅的折射率大约3.45和1.44当频率上限为3×10^14 Hz时对应的材料内最短波长约为λ_min c / (f_max × n_max) ≈ 3×10^8 / (3×10^14 × 3.45) ≈ 290 nm经验法则是一个波长内至少要有10个以上网格点所以30 nm以下都能接受。但Floquet周期边界要求源边界和目标边界的网格完全匹配映射网格天然满足这个要求所以我一般直接压到10 nm算得也很快整个模型不过几千个自由度。5.4 模式数不够导致带隙不明显如果你只解了2到3个特征频率很可能画出来的能带图只有一两条曲线带隙根本看不出来。我建议起步8个如果结构折射率对比度高或者你加了基底再多加几个。搜索频率上限也要跟上。如果你只搜到2×10^14 Hz而第二带在2.5×10^14 Hz处才打开那你看到的只是第一带带隙自然不存在。别偷懒先按我给的3×10^14 Hz试一次再把范围缩小到你真正关注的频段。5.5 归一化频率单位换算搞混画图时如果选择纵轴为a/λ那是一个无量纲数。我算的例子中a 500 nm中心波长2244 nma/λ约等于0.223。如果看到能带图的纵轴数值在0.2附近说明结果合理。如果你用的是频率f直接作图注意COMSOL里特征频率的默认单位是Hz但在某些版本界面里会显示为THz。这个单位显示问题曾经让我一度以为结果差了1000倍。整理成一张速查表问题可能原因排查方向特征频率为0周期边界未生效检查Floquet波矢参数能带图跳变模式排序不稳定减少搜索范围/外部重排带隙位置偏差大材料参数或几何厚度填错对照布拉格公式手算验证曲线看起来乱横向高阶模介入缩小单胞宽度W计算极慢网格过细先用映射网格粗网格试跑这几条基本覆盖了我自己从零开始算光子晶体能带时遇到的绝大部分问题。每次换新材料、换参数我都会先在脑子里过一遍这张表能省下大量调试时间。6. 一个容易被忽略但很实用的收尾技巧最后分享一个我自己的小习惯每跑完一个能带模型我会顺手在全局参数里改一下占空比也就是硅层厚度占周期厚度的比例比如把d1从200 nm改成150 nm再重新算一遍能带观察带隙宽度怎么变化。这个操作几分钟就出结果但对理解和优化光子晶体结构非常值。带隙宽度通常不是占空比的线性函数存在一个最优值。手算布拉格条件只能告诉你带隙中心在哪带隙到底有多宽必须靠数值扫描。我算过一组硅/二氧化硅一维结构占空比从0.2变到0.5时带隙宽度先增后减在0.4附近达到最大。这种规律不实际跑一遍参数扫描很难凭直觉猜准。如果你后续要做更复杂的结构比如二维光子晶体、带缺陷的波导或者想引入非线性材料这个能带计算流程完全可以平滑迁移过去。核心逻辑不变单胞、Floquet周期边界、特征频率扫描只是几何和网格需要跟着升级。希望这份实操记录能帮你少绕几次路一次把能带图画明白。
