1. 内容整体设计与思路拆解1.1 光子晶体到底在算什么从“高速公路收费站”说起把一维光子晶体比喻成“光子的高速公路收费站”这个类比其实挺传神的。你可以想象一束光在两种不同折射率的介质交替堆叠的结构里穿行每经过一个界面就相当于过了一次收费站——一部分光被反射回去一部分光继续前行。收费站设得密不密集、每个站收多少“过路费”取决于两种介质的折射率差和每层介质层的厚度。而整个收费站的布局规则——也就是晶格排列的方式决定了哪些频率的光子能顺利通过、哪些频率的光子会被彻底拦下。咱们要做的这件事就是在COMSOL Multiphysics里把这个“收费站”给搭出来然后用特征频率求解器去扫描波矢空间得到结构在不同波矢下的本征频率画出一张“频率—波矢”关系图也就是平时说的光子能带图。这张图最直观的作用就是用来找带隙band gap某一频率范围内的光无法在该结构中传播这对应着光子晶体最核心的应用——光的反射镜、微腔、波导、滤波器等等。挺好的一点是一维光子晶体是三种维度一维、二维、三维里最容易通过手算验证的。为什么这么说因为一维结构沿着一个方向周期性变化麦克斯韦方程组可以退化成简单的传输矩阵问题甚至用MATLAB几十行代码就能解出来。所以用COMSOL算一维结构不只是为了拿一张好看的能带图更重要的目的是理解整个建模仿真流程——从几何构建、材料参数、边界条件到求解器配置然后把这一套流程平移到二维、三维光子晶体上去时就不会发懵。1.2 为什么选COMSOL而不是其他工具说到光子晶体能带计算行业里常见的工具还有Lumerical FDTD、CST、RSoft等。FDTD靠的是时域有限差分通过宽带脉冲入射然后做傅里叶变换来提取透射谱CST侧重高频电磁仿真RSoft里集成了专门的能带计算模块使用起来也相当顺手。那我为什么在这里选COMSOL原因有几个。第一COMSOL把几何建模、物理场定义、网格划分、求解和后处理都整合在同一个图形化界面里对新手来说上手门槛比写脚本低得多。第二COMSOL的“特征频率”求解器加“周期性边界条件”的组合天然适配能带计算这种本征值问题——你不需要自己去实现布洛赫定理只需要把波矢分量作为参数扫过去就行。第三COMSOL的“App开发器”和“参数化扫描”功能很强扫一条能带曲线基本上就是设置一个参数就能自动跑完非常适合做参数化研究。它也有它自己的问题后面对比的时候我再细说。总体来说如果你需要明确看到模式场分布电场、磁场、能量密度、需要反复调整结构尺寸并观察能带变化COMSOL是一个非常顺手的选择。1.3 能带计算背后的物理布洛赫定理和简约布里渊区这里必须稍微提一下布洛赫定理因为整个建模方式就是围绕它展开的不理解这个你后面设置边界条件时会很摸不着头脑。在周期性结构里电磁场的本征模式可以写成布洛赫波的形式场 周期函数 × 平面波因子。这个平面波因子里有一个波矢 k它的物理含义是光子在晶格中传播时的“拟动量”。k 不是一个连续的任意值而是落在第一布里渊区这里就是 [-π/a, π/a]a 是晶格常数内。由于能带结构具有周期性和反演对称性实际计算时通常只扫第一布里渊区一半也就是 [0, π/a]这个区间对应一维晶格里的简约布里渊区边界通常标记为 X 点。扫描过程说起来很直接固定一个 k 值求解该 k 值下结构的所有本征频率然后让 k 从 0 逐渐增加到 π/a把每个 k 下求得的频率点连成一条条曲线。每条曲线就是一个能带相邻能带之间如果存在一个没有任何模式覆盖的频率区间那就是带隙。你可能会问为什么这个计算要用“特征频率求解器”而不是普通的高频电磁波仿真道理在于能带图关心的是“结构里存在哪些本征模式以及它们对应的频率”这本质上是一个本征值问题不是输入激励然后看响应的散射问题。COMSOL里的特征频率求解器就是专门干这个的它求解的是无源情况下的谐振模式。1.4 场景选择频域还是特征频率在COMSOL的RF模块里建这个仿真核心物理场接口选“电磁波频域ewfd”。特征频率求解模式下它实际解的是无源波动方程的本征模式也就是亥姆霍兹方程的特征值问题。这里面有一个细微的区别需要说清楚特征频率求解器求解出的“特征频率”对应的是结构“纵向谐振”的频率。但对于一维光子晶体能带计算我们关心的其实是布洛赫波矢 k 和频率 ω 的色散关系——周期性边界条件已经引入了 k 的相位关系所以求解得到的模式频率就是这个 k 下允许传播的模式频率。这里的“谐振”是周期性单元Virtual内部场的自洽解并不代表物理上存在一个封闭腔体的谐振。说个容易混淆的点有些初学者会把能带计算和超材料等效参数提取搞混。超材料等效参数提取需要用S参数反演本质上是散射问题而能带计算是纯本征问题核心设置就是“周期性边界条件 特征频率扫描 波矢参数扫描”这点请先把逻辑盘清楚再动手建模。2. 核心细节解析与实操要点2.1 从物理模型到几何构建硅基底上搭一维光子晶体咱们的标题里已经给了一个很明确的应用场景硅基底上搭周期性介电结构。为什么选硅基底因为硅在近红外波段1.1μm~1.5μm附近是透明的折射率约为3.45和空气n1或二氧化硅n1.45的折射率对比很强正是做光子晶体最经典的材料体系之一。实际建模中可以有两种路径一种是严格建模硅基底一个厚硅层 上方周期性介电结构另一种是只建周期性结构本身硅和空气交替层叠不考虑基底。第一种更贴近实际的实验样品——比如你在SOI硅片上干法刻蚀出一排排周期性凹槽那自然就有基底存在第二种更干净便于理论验证和能带分析基底的存在并不会改变带隙的中心频率和宽度只要基底足够厚。两种路径我都走过要说仿真效率去掉基底速度能快不少要说实用价值保留基底更贴近你最终要去加工验证的结构。在教程这个层面我建议的方式是先把无基底模型跑通等你会读能带图、会调参数了再在模型里加一个硅基底层看看有什么变化。这样拆解下来每一步的误差来源都更可控。2.2 材料参数折射率的实部与虚部在光学波段COMSOL里定义介质材料最常用的方式就是直接指定折射率 n 和消光系数 k注意这里的 k 不是波矢别混淆。在频域仿真中这种处理完全够用因为咱们关心的是定性规律——能带的形状、带隙位置和宽度。不过有一点我说下我自己的习惯在初步设计阶段把材料当成无损耗介质也就是 k0来算。为什么因为损耗会使得本征频率变成复数特征频率求解器解出的“虚部”会变大你就得区分“这个模式是真的衰减了还是数值误差”同时求解器的收敛性也会变差。先无损耗跑通拿到干净的能带图再在后续优化阶段引入真实材料色散看看损耗对带隙边缘的影响。硅在红光波段以下1.1μm开始明显吸收做可见光波段光子晶体时最好不要用硅基底、改用二氧化钛TiO2n≈2.4~2.8或氮化硅Si3N4n≈2.0更合适。紫色标题里的“硅基底”应该理解为近红外应用场景。2.3 几何参数周期 a、填充比 f、厚度 h——它们怎么定一维光子晶体的核心几何参数就三个晶格常数 a一个周期内的总长度、填充比 f高折射率材料在一个周期内占的体积比例一维情况下就是厚度占比、以及层数 N如果建模对象是有限周期结构。选择参数有个经典的参考起点对于四分之一波长堆栈即每层光学厚度都是 λ/4带隙中心频率可以粗略估算成 f_center c / (2 × (n1×d1 n2×d2))。在近红外中心波长1.55μm处硅n3.45厚度取112nm左右、空气n1厚度取388nm左右总周期a约500nm这就是一个很常见的起点。有一个可以快速手算验证的办法用传输矩阵法TMM算一个10周期结构的光谱找到反射带位置再和COMSOL算的能带图对比。两者应该能对得上。这也是我强烈推荐各位实践的第一步——别急着直接上COMSOL先用TMM把理论值摸清楚再去仿真你会对自己搭的模型有多大的信心。下面给一个参考参数表后续实操演示都用这组参数数值说明a500 nm晶格常数d1112 nm硅层厚度d2388 nm空气层厚度n13.45硅折射率无损耗n21.0空气折射率基底不建模/后续可选硅厚层2.4 周期性边界条件和布洛赫周期条件在COMSOL的“电磁波频域”接口里设置周期性边界条件关键是理解“周期性”的类型。能带计算必须选Floquet周期性边界条件在COMSOL里叫“周期性条件”周期类型选择“Floquet”而不是简单的周期性条件。Floquet边界条件的数学表达式是目标边界的场 源边界的场 × exp(-i·k_F·r)。其中 k_F 就是布洛赫波矢。在一维情况下只有 x 方向周期方向需要设置 Floquet 条件y 面、z 面则需要设置其他边界条件。2.5 边界条件的坑一维模型的YZ平面怎么设这是个非常关键的细节。一维光子晶体只在x方向周期变化在y和z方向结构是均匀无限的。仿真时单元胞在x方向取一个周期a在y和z方向取多长如果y、z方向取的尺寸很小比如就取纳米宽那y和z方向的边界条件就要搞对。这里建议用“周期性条件”或“完美磁导体/完美电导体”组合视偏振情况而定。如果你的关注点是TE偏振电场平行于界面那就用PEC边界TM偏振磁场平行于界面用PMC边界。如果y、z方向取足够宽比如几个微米边界影响会变小但计算量会上升。我实测下来对于能带计算这种本征值问题建议y、z方向各取100nm就够配合PEC/PMC边界处理计算量极小、结果也准。如果你想要结果更干净直接在“特征频率”设置里把y、z方向的“面外波矢分量”设为0就完事了。2.6 网格划分光子晶体仿真的一碗水端平网格是这类仿真的分水岭。光子晶体里的电场分布对介质界面非常敏感光学波段折射率高硅n3.45波长在介质里会被压缩成真空波长的1/3.45因此网格必须细。我的经验法则最大网格尺寸控制在介质内最小波长的1/10到1/12。对近红外1.55μm硅里波长约449nm网格最大尺寸45nm就够了。空气里波长1.55μm网格可以稍微放松到100nm。但COMSOL的自动网格未必能满足建议在介质层内部加“大小”节点手动指定最大单元尺寸。还有一个容易被忽略的点一维光子晶体结构非常薄112nm如果网格剖分得不好会有高度方向分辨率不足的问题导致模式频率漂移。我建议在硅层内至少剖分2~3层六面体或超薄四面体单元这直接关系到高频带数据的准确性。3. 实操过程与核心环节实现3.1 建立模型的八步流程总览拿到这个仿真任务操作流程可以拆成八步设置模型空间维度为二维或三维后面我会说两者的差异。添加“电磁波频域”物理场接口。构建几何一个周期内的高折射率层低折射率层。指定材料折射率。设置Floquet周期性边界条件。划分网格。配置特征频率求解器并设置波矢参数扫描。后处理绘制能带图检查带隙。3.2 几何建模一个周期就够别傻乎乎的建一堆层很多新手一开始容易习惯性建很多重复单元比如一次性画10个周期。除非你是要算透射谱或反射谱否则能带计算根本不需要。周期结构的本征问题只要一个单元胞就够了这与布洛赫定理直接对应——整个无限周期结构的模式信息全部压缩在一个周期内通过不同的布洛赫波矢 k 来区分。在COMSOL里模型维度和几何尺寸怎么定得看你的偏振需求。如果只做二维近似直接建一个 500nm × 200nm 的矩形里面用两条线把x方向切成三段左边硅层112nm、中间空气层388nm、右边硅层112nm但这右边那段其实属于下一个周期的起始你可以这样理解左侧和右侧分别与相邻单元胞共享界面。还有一个更直观的建法把单元胞建在 [-a/2, a/2] 区间内中心是空气层两侧各是半个硅层厚度56nm。这种布局会让Floquet边界条件的两侧完全对称周期性的相位差刚好落在边界条件里物理上更清晰。3.3 参数设置和物理场设定在“全局定义”里定义参数a 500[nm] d1 112[nm] d2 388[nm] n_Si 3.45 n_air 1.0 k_x 0 # 布洛赫波矢后续扫描用在“电磁波频域”接口里选择“特征频率”研究在“电磁波频域”节点设置中把“电场”分量设为面外方向即Ez这对应TM偏振习惯定义略有差异关键在于模型里只有一个分量算的是标量亥姆霍兹方程设置相对介电常数硅层为 n_Si^2空气层为 n_air^23.4 Floquet周期性边界条件的精确配置在“周期性条件”节点选择两个相对的面x-a/2 和 xa/2周期类型选“Floquet”。注意k矢量的分量写法。一维模型里要指定k_Floquet_x k_x k_Floquet_y 0 k_Floquet_z 0其中 k_x 在后续参数扫描里从 0 变化到 π/a。有一个比较隐蔽但重要的地方COMSOL 在 Floquet 边界条件里使用的 k_Floquet 有正负号约定exp(i(k·r))还是exp(-i(k·r))如果设置方向反了你会发现能带曲线出现负群速度模式或者模式个数不对。需要检查模型的边界坐标系方向来确认正负号。我这边的经验是在单元胞左侧边界选“源边界”右侧边界选“目标边界”目标边界上的相位设置为 k_Floquet × (x_goal - x_source)正负号在目标边界上加“-”。3.5 特征频率扫描的设置方式这是整个操作里最能体现COMSOL技巧的一环。我们想要的是k_x从0到π/a每个k值解出一组特征频率。在COMSOL里实现这个有两种做法方式A辅助扫描推荐在研究中新建“辅助扫描”选择全局参数 k_x扫描值用范围range(0 π/a/20 π/a)也就是21个点。勾选“扫描时更新几何或物理场”不用开因为k_x只影响边界条件不影响几何。方式B手动批量扫描不建辅助扫描手动把k值设置成多个研究这在参数特别多时会比较繁琐。我推荐方式A因为后期如果要在Energies里提取数据辅助扫描的结果会多出一个维度导出表格时顺理成章。在特征频率求解器中有两点必须注意所需特征频率数要设置的比期望的带数多。例如一个一维二组分结构低k处至少前5条能带里会出现带隙那么你至少设置“所需特征频率数8”防止重点频段被漏掉。搜索频率基准建议围绕预期中心频率设置比如1.55μm对应的频率约193THz则把搜索基准设在193THz这样才能快速收敛到光学模式而不是零频附近的数值噪声。3.6 求解和后处理把数据变成能带图求解完成后结果里会有一堆特征频率值它们是 k_x 扫描参数的函数后处理里需要把数据重新整理一下。在COMSOL中常用做法是在“结果”中添加一维绘图组。绘图数据选“特征频率”数据集。x轴设为表达式参数即k_xy轴设为特征频率值单位为THz或归一化频率 a/λ。更丝滑的做法是直接导出表格文件→导出→数据然后放到Origin、Python或者Excel里画图。我个人的偏好是用Python的matplotlib因为能同时对比不同结构的能带批量处理起来更灵活。归一化频率 a/λ 也是个常用横纵坐标即用 a/λ 做y轴这样不同周期长度的结构可以直接对比。换算关系freq_THz c/(a×归一化频率) 反过来归一化频率 a/λ freq_THz × a / c。3.7 配一个完整的操作验收流程拿前面那组参数我给出一个完整的验收实验扫描k_x 0到π/a共21个点每个点求解8个特征频率得到能带图后找出带隙位置用传输矩阵法手算反射带范围对比以 d1112nmd2388nmn13.45n21.0 为例TMM手算预期反射带带隙中心约在1.55μm反射带边缘大致覆盖1.3μm~1.85μm。COMSOL能带图里应该能看到带隙中心频率在193THz附近。这个验证一过你的模型基本就靠谱了。4. 常见问题与排查技巧实录4.1 能带图里出现“断带”——跨波段频率跳变跑完辅助扫描之后把能带图画出来最常见的现象是曲线不是连续的有的点突然跳到很高的频率连不成一条平滑的能带。这个问题的根源在于特征频率求解器按频率升序返回结果但不同 k 点处第N个解可能对应的是不同的模式分支。解决思路是后处理时做“模式追踪”mode tracing。COMSOL从5.x开始有模式追踪功能它可以根据前一个k点的场分布来匹配当前k点最相似的模式。开启方式是在特征频率节点下勾选“模式追踪”并指定一个跟踪量通常是电场积分或最大场位置。如果版本没有模式追踪还有个朴素的办法把k扫描步长缩小特别是靠近能带边缘的地方模式之间靠得很近步长不够就会出现模式跳变。把扫描从21个点加到41个点能解决大部分问题。4.2 带隙位置和手算值差很多这是新手最容易懵的问题。我排查过的案例里绝大多数原因是边界条件没设对。情况一周期性边界条件没有设成Floquet而是默认的周期条件这会在布里渊区边界X点处给出错误的本征频率。情况二Floquet波矢符号和方向反了导致带结构在1/2布里渊区处出现不对称的假象。情况三少设了面外方向的波矢分量比如二维模型里面外方向也需要指定 k0缺失会引入额外自由度。还有一个容易被忽视的用相对介电常数还是折射率来定义材料。COMSOL里RF模块的相对介电常数是复数折射率的平方等于相对介电常数无损耗时这个换算不能错。用“折射率”材料模型的时候要注意频率单位一致。4.3 高频模式不收敛或收敛慢一维光子晶体在COMSOL中网格剖分相对简单但高频带若出现不收敛或特征频率虚部大最常见的原因是网格不够细尤其在高折射率层硅层。硅中的波长是真空的1/3.45如果自动网格最大单元是 λ/5相当于硅内只有不到一个单元当然不准。解决方法是给硅层单独加一个“大小”节点最大单元尺寸设为40nm同时在高度方向至少剖分2层。实测这样设置后前10条能带的精度都能控制在1%以内。4.4 特征频率出现“零频模式”和“伪模式”特征频率求解器除了物理上有意义的模式还会解出一些数值噪声模式。通常它们的特点是频率接近0或者虚部异常大场分布呈棋盘状高频振荡对应网格尺度上的伪振荡。处理办法设置“搜索频率基准”为预期中心频率比如193THz这样可以压制零频附近的无用解在求解器配置中增大“所需特征频率数”丢掉最低几个接近0的模式用“电场模”的绘图检查可疑模式如果场分布密密麻麻、没有清晰的空间结构多半是伪模式4.5 常见问题速查表问题现象可能原因解决办法能带曲线断带k扫描步长过大缩小扫描间隔开启模式追踪带隙位置偏离理论值超过10%Floquet边界条件符号错误检查源/目标边界设置及相位符号高频能带抖动严重硅层网格太粗硅层最大网格单元设为40nm出现一堆0附近模式未设置搜索频率基准设置基准为预估中心频率模式和文献差距大材料折射率单位或频率单位不一致统一用THz确认相对介电常数 n^2求解器不收敛几何上存在尖锐角或重复面检查几何布尔运算清理多余面4.6 关于二维和三维扩展的几点建议二维光子晶体的能带计算比如六角排列空气孔和三维光子晶体比如Yablonovite结构在COMSOL里的建模思路完全一样的区别在于二维结构需要扫描布里渊区的高对称点路径Γ-M-K-Γ涉及面外波矢分量的参数变化三维结构的Floquet条件要设置两个方向的相位差或三个方向布里渊区高对称点更多计算量呈指数增长如果你已经能跑通一维模型注意你是把“波矢分量”当成参数扫描的——一维是一个参数二维就是两个参数的组合扫描三维就是三个参数。COMSOL的参数化扫描器本身可以处理多维参数组合但计算量要提前评估。我在实际使用中的体会是一维光子晶体虽然结构简单但它是检验你对周期性边界条件、布洛赫定理、特征频率求解器理解程度的黄金试金石。把一维模型跑得滚瓜烂熟之后再去碰二维、三维你会有一种“不过如此”的底气。最后再分享一个小技巧COMSOL的模型测量里可以看到“周期性条件”的面方向方向设反是新手最容易踩的坑。每次建模后先画一个本征模的电场分布看看x方向是否存在预期的布洛赫相位调制——如果边界处的场不连续回查边界相位符号就没错了。 COMSOL光子晶体能带计算：从布洛赫定理到带隙仿真
1. 内容整体设计与思路拆解1.1 光子晶体到底在算什么从“高速公路收费站”说起把一维光子晶体比喻成“光子的高速公路收费站”这个类比其实挺传神的。你可以想象一束光在两种不同折射率的介质交替堆叠的结构里穿行每经过一个界面就相当于过了一次收费站——一部分光被反射回去一部分光继续前行。收费站设得密不密集、每个站收多少“过路费”取决于两种介质的折射率差和每层介质层的厚度。而整个收费站的布局规则——也就是晶格排列的方式决定了哪些频率的光子能顺利通过、哪些频率的光子会被彻底拦下。咱们要做的这件事就是在COMSOL Multiphysics里把这个“收费站”给搭出来然后用特征频率求解器去扫描波矢空间得到结构在不同波矢下的本征频率画出一张“频率—波矢”关系图也就是平时说的光子能带图。这张图最直观的作用就是用来找带隙band gap某一频率范围内的光无法在该结构中传播这对应着光子晶体最核心的应用——光的反射镜、微腔、波导、滤波器等等。挺好的一点是一维光子晶体是三种维度一维、二维、三维里最容易通过手算验证的。为什么这么说因为一维结构沿着一个方向周期性变化麦克斯韦方程组可以退化成简单的传输矩阵问题甚至用MATLAB几十行代码就能解出来。所以用COMSOL算一维结构不只是为了拿一张好看的能带图更重要的目的是理解整个建模仿真流程——从几何构建、材料参数、边界条件到求解器配置然后把这一套流程平移到二维、三维光子晶体上去时就不会发懵。1.2 为什么选COMSOL而不是其他工具说到光子晶体能带计算行业里常见的工具还有Lumerical FDTD、CST、RSoft等。FDTD靠的是时域有限差分通过宽带脉冲入射然后做傅里叶变换来提取透射谱CST侧重高频电磁仿真RSoft里集成了专门的能带计算模块使用起来也相当顺手。那我为什么在这里选COMSOL原因有几个。第一COMSOL把几何建模、物理场定义、网格划分、求解和后处理都整合在同一个图形化界面里对新手来说上手门槛比写脚本低得多。第二COMSOL的“特征频率”求解器加“周期性边界条件”的组合天然适配能带计算这种本征值问题——你不需要自己去实现布洛赫定理只需要把波矢分量作为参数扫过去就行。第三COMSOL的“App开发器”和“参数化扫描”功能很强扫一条能带曲线基本上就是设置一个参数就能自动跑完非常适合做参数化研究。它也有它自己的问题后面对比的时候我再细说。总体来说如果你需要明确看到模式场分布电场、磁场、能量密度、需要反复调整结构尺寸并观察能带变化COMSOL是一个非常顺手的选择。1.3 能带计算背后的物理布洛赫定理和简约布里渊区这里必须稍微提一下布洛赫定理因为整个建模方式就是围绕它展开的不理解这个你后面设置边界条件时会很摸不着头脑。在周期性结构里电磁场的本征模式可以写成布洛赫波的形式场 周期函数 × 平面波因子。这个平面波因子里有一个波矢 k它的物理含义是光子在晶格中传播时的“拟动量”。k 不是一个连续的任意值而是落在第一布里渊区这里就是 [-π/a, π/a]a 是晶格常数内。由于能带结构具有周期性和反演对称性实际计算时通常只扫第一布里渊区一半也就是 [0, π/a]这个区间对应一维晶格里的简约布里渊区边界通常标记为 X 点。扫描过程说起来很直接固定一个 k 值求解该 k 值下结构的所有本征频率然后让 k 从 0 逐渐增加到 π/a把每个 k 下求得的频率点连成一条条曲线。每条曲线就是一个能带相邻能带之间如果存在一个没有任何模式覆盖的频率区间那就是带隙。你可能会问为什么这个计算要用“特征频率求解器”而不是普通的高频电磁波仿真道理在于能带图关心的是“结构里存在哪些本征模式以及它们对应的频率”这本质上是一个本征值问题不是输入激励然后看响应的散射问题。COMSOL里的特征频率求解器就是专门干这个的它求解的是无源情况下的谐振模式。1.4 场景选择频域还是特征频率在COMSOL的RF模块里建这个仿真核心物理场接口选“电磁波频域ewfd”。特征频率求解模式下它实际解的是无源波动方程的本征模式也就是亥姆霍兹方程的特征值问题。这里面有一个细微的区别需要说清楚特征频率求解器求解出的“特征频率”对应的是结构“纵向谐振”的频率。但对于一维光子晶体能带计算我们关心的其实是布洛赫波矢 k 和频率 ω 的色散关系——周期性边界条件已经引入了 k 的相位关系所以求解得到的模式频率就是这个 k 下允许传播的模式频率。这里的“谐振”是周期性单元Virtual内部场的自洽解并不代表物理上存在一个封闭腔体的谐振。说个容易混淆的点有些初学者会把能带计算和超材料等效参数提取搞混。超材料等效参数提取需要用S参数反演本质上是散射问题而能带计算是纯本征问题核心设置就是“周期性边界条件 特征频率扫描 波矢参数扫描”这点请先把逻辑盘清楚再动手建模。2. 核心细节解析与实操要点2.1 从物理模型到几何构建硅基底上搭一维光子晶体咱们的标题里已经给了一个很明确的应用场景硅基底上搭周期性介电结构。为什么选硅基底因为硅在近红外波段1.1μm~1.5μm附近是透明的折射率约为3.45和空气n1或二氧化硅n1.45的折射率对比很强正是做光子晶体最经典的材料体系之一。实际建模中可以有两种路径一种是严格建模硅基底一个厚硅层 上方周期性介电结构另一种是只建周期性结构本身硅和空气交替层叠不考虑基底。第一种更贴近实际的实验样品——比如你在SOI硅片上干法刻蚀出一排排周期性凹槽那自然就有基底存在第二种更干净便于理论验证和能带分析基底的存在并不会改变带隙的中心频率和宽度只要基底足够厚。两种路径我都走过要说仿真效率去掉基底速度能快不少要说实用价值保留基底更贴近你最终要去加工验证的结构。在教程这个层面我建议的方式是先把无基底模型跑通等你会读能带图、会调参数了再在模型里加一个硅基底层看看有什么变化。这样拆解下来每一步的误差来源都更可控。2.2 材料参数折射率的实部与虚部在光学波段COMSOL里定义介质材料最常用的方式就是直接指定折射率 n 和消光系数 k注意这里的 k 不是波矢别混淆。在频域仿真中这种处理完全够用因为咱们关心的是定性规律——能带的形状、带隙位置和宽度。不过有一点我说下我自己的习惯在初步设计阶段把材料当成无损耗介质也就是 k0来算。为什么因为损耗会使得本征频率变成复数特征频率求解器解出的“虚部”会变大你就得区分“这个模式是真的衰减了还是数值误差”同时求解器的收敛性也会变差。先无损耗跑通拿到干净的能带图再在后续优化阶段引入真实材料色散看看损耗对带隙边缘的影响。硅在红光波段以下1.1μm开始明显吸收做可见光波段光子晶体时最好不要用硅基底、改用二氧化钛TiO2n≈2.4~2.8或氮化硅Si3N4n≈2.0更合适。紫色标题里的“硅基底”应该理解为近红外应用场景。2.3 几何参数周期 a、填充比 f、厚度 h——它们怎么定一维光子晶体的核心几何参数就三个晶格常数 a一个周期内的总长度、填充比 f高折射率材料在一个周期内占的体积比例一维情况下就是厚度占比、以及层数 N如果建模对象是有限周期结构。选择参数有个经典的参考起点对于四分之一波长堆栈即每层光学厚度都是 λ/4带隙中心频率可以粗略估算成 f_center c / (2 × (n1×d1 n2×d2))。在近红外中心波长1.55μm处硅n3.45厚度取112nm左右、空气n1厚度取388nm左右总周期a约500nm这就是一个很常见的起点。有一个可以快速手算验证的办法用传输矩阵法TMM算一个10周期结构的光谱找到反射带位置再和COMSOL算的能带图对比。两者应该能对得上。这也是我强烈推荐各位实践的第一步——别急着直接上COMSOL先用TMM把理论值摸清楚再去仿真你会对自己搭的模型有多大的信心。下面给一个参考参数表后续实操演示都用这组参数数值说明a500 nm晶格常数d1112 nm硅层厚度d2388 nm空气层厚度n13.45硅折射率无损耗n21.0空气折射率基底不建模/后续可选硅厚层2.4 周期性边界条件和布洛赫周期条件在COMSOL的“电磁波频域”接口里设置周期性边界条件关键是理解“周期性”的类型。能带计算必须选Floquet周期性边界条件在COMSOL里叫“周期性条件”周期类型选择“Floquet”而不是简单的周期性条件。Floquet边界条件的数学表达式是目标边界的场 源边界的场 × exp(-i·k_F·r)。其中 k_F 就是布洛赫波矢。在一维情况下只有 x 方向周期方向需要设置 Floquet 条件y 面、z 面则需要设置其他边界条件。2.5 边界条件的坑一维模型的YZ平面怎么设这是个非常关键的细节。一维光子晶体只在x方向周期变化在y和z方向结构是均匀无限的。仿真时单元胞在x方向取一个周期a在y和z方向取多长如果y、z方向取的尺寸很小比如就取纳米宽那y和z方向的边界条件就要搞对。这里建议用“周期性条件”或“完美磁导体/完美电导体”组合视偏振情况而定。如果你的关注点是TE偏振电场平行于界面那就用PEC边界TM偏振磁场平行于界面用PMC边界。如果y、z方向取足够宽比如几个微米边界影响会变小但计算量会上升。我实测下来对于能带计算这种本征值问题建议y、z方向各取100nm就够配合PEC/PMC边界处理计算量极小、结果也准。如果你想要结果更干净直接在“特征频率”设置里把y、z方向的“面外波矢分量”设为0就完事了。2.6 网格划分光子晶体仿真的一碗水端平网格是这类仿真的分水岭。光子晶体里的电场分布对介质界面非常敏感光学波段折射率高硅n3.45波长在介质里会被压缩成真空波长的1/3.45因此网格必须细。我的经验法则最大网格尺寸控制在介质内最小波长的1/10到1/12。对近红外1.55μm硅里波长约449nm网格最大尺寸45nm就够了。空气里波长1.55μm网格可以稍微放松到100nm。但COMSOL的自动网格未必能满足建议在介质层内部加“大小”节点手动指定最大单元尺寸。还有一个容易被忽略的点一维光子晶体结构非常薄112nm如果网格剖分得不好会有高度方向分辨率不足的问题导致模式频率漂移。我建议在硅层内至少剖分2~3层六面体或超薄四面体单元这直接关系到高频带数据的准确性。3. 实操过程与核心环节实现3.1 建立模型的八步流程总览拿到这个仿真任务操作流程可以拆成八步设置模型空间维度为二维或三维后面我会说两者的差异。添加“电磁波频域”物理场接口。构建几何一个周期内的高折射率层低折射率层。指定材料折射率。设置Floquet周期性边界条件。划分网格。配置特征频率求解器并设置波矢参数扫描。后处理绘制能带图检查带隙。3.2 几何建模一个周期就够别傻乎乎的建一堆层很多新手一开始容易习惯性建很多重复单元比如一次性画10个周期。除非你是要算透射谱或反射谱否则能带计算根本不需要。周期结构的本征问题只要一个单元胞就够了这与布洛赫定理直接对应——整个无限周期结构的模式信息全部压缩在一个周期内通过不同的布洛赫波矢 k 来区分。在COMSOL里模型维度和几何尺寸怎么定得看你的偏振需求。如果只做二维近似直接建一个 500nm × 200nm 的矩形里面用两条线把x方向切成三段左边硅层112nm、中间空气层388nm、右边硅层112nm但这右边那段其实属于下一个周期的起始你可以这样理解左侧和右侧分别与相邻单元胞共享界面。还有一个更直观的建法把单元胞建在 [-a/2, a/2] 区间内中心是空气层两侧各是半个硅层厚度56nm。这种布局会让Floquet边界条件的两侧完全对称周期性的相位差刚好落在边界条件里物理上更清晰。3.3 参数设置和物理场设定在“全局定义”里定义参数a 500[nm] d1 112[nm] d2 388[nm] n_Si 3.45 n_air 1.0 k_x 0 # 布洛赫波矢后续扫描用在“电磁波频域”接口里选择“特征频率”研究在“电磁波频域”节点设置中把“电场”分量设为面外方向即Ez这对应TM偏振习惯定义略有差异关键在于模型里只有一个分量算的是标量亥姆霍兹方程设置相对介电常数硅层为 n_Si^2空气层为 n_air^23.4 Floquet周期性边界条件的精确配置在“周期性条件”节点选择两个相对的面x-a/2 和 xa/2周期类型选“Floquet”。注意k矢量的分量写法。一维模型里要指定k_Floquet_x k_x k_Floquet_y 0 k_Floquet_z 0其中 k_x 在后续参数扫描里从 0 变化到 π/a。有一个比较隐蔽但重要的地方COMSOL 在 Floquet 边界条件里使用的 k_Floquet 有正负号约定exp(i(k·r))还是exp(-i(k·r))如果设置方向反了你会发现能带曲线出现负群速度模式或者模式个数不对。需要检查模型的边界坐标系方向来确认正负号。我这边的经验是在单元胞左侧边界选“源边界”右侧边界选“目标边界”目标边界上的相位设置为 k_Floquet × (x_goal - x_source)正负号在目标边界上加“-”。3.5 特征频率扫描的设置方式这是整个操作里最能体现COMSOL技巧的一环。我们想要的是k_x从0到π/a每个k值解出一组特征频率。在COMSOL里实现这个有两种做法方式A辅助扫描推荐在研究中新建“辅助扫描”选择全局参数 k_x扫描值用范围range(0 π/a/20 π/a)也就是21个点。勾选“扫描时更新几何或物理场”不用开因为k_x只影响边界条件不影响几何。方式B手动批量扫描不建辅助扫描手动把k值设置成多个研究这在参数特别多时会比较繁琐。我推荐方式A因为后期如果要在Energies里提取数据辅助扫描的结果会多出一个维度导出表格时顺理成章。在特征频率求解器中有两点必须注意所需特征频率数要设置的比期望的带数多。例如一个一维二组分结构低k处至少前5条能带里会出现带隙那么你至少设置“所需特征频率数8”防止重点频段被漏掉。搜索频率基准建议围绕预期中心频率设置比如1.55μm对应的频率约193THz则把搜索基准设在193THz这样才能快速收敛到光学模式而不是零频附近的数值噪声。3.6 求解和后处理把数据变成能带图求解完成后结果里会有一堆特征频率值它们是 k_x 扫描参数的函数后处理里需要把数据重新整理一下。在COMSOL中常用做法是在“结果”中添加一维绘图组。绘图数据选“特征频率”数据集。x轴设为表达式参数即k_xy轴设为特征频率值单位为THz或归一化频率 a/λ。更丝滑的做法是直接导出表格文件→导出→数据然后放到Origin、Python或者Excel里画图。我个人的偏好是用Python的matplotlib因为能同时对比不同结构的能带批量处理起来更灵活。归一化频率 a/λ 也是个常用横纵坐标即用 a/λ 做y轴这样不同周期长度的结构可以直接对比。换算关系freq_THz c/(a×归一化频率) 反过来归一化频率 a/λ freq_THz × a / c。3.7 配一个完整的操作验收流程拿前面那组参数我给出一个完整的验收实验扫描k_x 0到π/a共21个点每个点求解8个特征频率得到能带图后找出带隙位置用传输矩阵法手算反射带范围对比以 d1112nmd2388nmn13.45n21.0 为例TMM手算预期反射带带隙中心约在1.55μm反射带边缘大致覆盖1.3μm~1.85μm。COMSOL能带图里应该能看到带隙中心频率在193THz附近。这个验证一过你的模型基本就靠谱了。4. 常见问题与排查技巧实录4.1 能带图里出现“断带”——跨波段频率跳变跑完辅助扫描之后把能带图画出来最常见的现象是曲线不是连续的有的点突然跳到很高的频率连不成一条平滑的能带。这个问题的根源在于特征频率求解器按频率升序返回结果但不同 k 点处第N个解可能对应的是不同的模式分支。解决思路是后处理时做“模式追踪”mode tracing。COMSOL从5.x开始有模式追踪功能它可以根据前一个k点的场分布来匹配当前k点最相似的模式。开启方式是在特征频率节点下勾选“模式追踪”并指定一个跟踪量通常是电场积分或最大场位置。如果版本没有模式追踪还有个朴素的办法把k扫描步长缩小特别是靠近能带边缘的地方模式之间靠得很近步长不够就会出现模式跳变。把扫描从21个点加到41个点能解决大部分问题。4.2 带隙位置和手算值差很多这是新手最容易懵的问题。我排查过的案例里绝大多数原因是边界条件没设对。情况一周期性边界条件没有设成Floquet而是默认的周期条件这会在布里渊区边界X点处给出错误的本征频率。情况二Floquet波矢符号和方向反了导致带结构在1/2布里渊区处出现不对称的假象。情况三少设了面外方向的波矢分量比如二维模型里面外方向也需要指定 k0缺失会引入额外自由度。还有一个容易被忽视的用相对介电常数还是折射率来定义材料。COMSOL里RF模块的相对介电常数是复数折射率的平方等于相对介电常数无损耗时这个换算不能错。用“折射率”材料模型的时候要注意频率单位一致。4.3 高频模式不收敛或收敛慢一维光子晶体在COMSOL中网格剖分相对简单但高频带若出现不收敛或特征频率虚部大最常见的原因是网格不够细尤其在高折射率层硅层。硅中的波长是真空的1/3.45如果自动网格最大单元是 λ/5相当于硅内只有不到一个单元当然不准。解决方法是给硅层单独加一个“大小”节点最大单元尺寸设为40nm同时在高度方向至少剖分2层。实测这样设置后前10条能带的精度都能控制在1%以内。4.4 特征频率出现“零频模式”和“伪模式”特征频率求解器除了物理上有意义的模式还会解出一些数值噪声模式。通常它们的特点是频率接近0或者虚部异常大场分布呈棋盘状高频振荡对应网格尺度上的伪振荡。处理办法设置“搜索频率基准”为预期中心频率比如193THz这样可以压制零频附近的无用解在求解器配置中增大“所需特征频率数”丢掉最低几个接近0的模式用“电场模”的绘图检查可疑模式如果场分布密密麻麻、没有清晰的空间结构多半是伪模式4.5 常见问题速查表问题现象可能原因解决办法能带曲线断带k扫描步长过大缩小扫描间隔开启模式追踪带隙位置偏离理论值超过10%Floquet边界条件符号错误检查源/目标边界设置及相位符号高频能带抖动严重硅层网格太粗硅层最大网格单元设为40nm出现一堆0附近模式未设置搜索频率基准设置基准为预估中心频率模式和文献差距大材料折射率单位或频率单位不一致统一用THz确认相对介电常数 n^2求解器不收敛几何上存在尖锐角或重复面检查几何布尔运算清理多余面4.6 关于二维和三维扩展的几点建议二维光子晶体的能带计算比如六角排列空气孔和三维光子晶体比如Yablonovite结构在COMSOL里的建模思路完全一样的区别在于二维结构需要扫描布里渊区的高对称点路径Γ-M-K-Γ涉及面外波矢分量的参数变化三维结构的Floquet条件要设置两个方向的相位差或三个方向布里渊区高对称点更多计算量呈指数增长如果你已经能跑通一维模型注意你是把“波矢分量”当成参数扫描的——一维是一个参数二维就是两个参数的组合扫描三维就是三个参数。COMSOL的参数化扫描器本身可以处理多维参数组合但计算量要提前评估。我在实际使用中的体会是一维光子晶体虽然结构简单但它是检验你对周期性边界条件、布洛赫定理、特征频率求解器理解程度的黄金试金石。把一维模型跑得滚瓜烂熟之后再去碰二维、三维你会有一种“不过如此”的底气。最后再分享一个小技巧COMSOL的模型测量里可以看到“周期性条件”的面方向方向设反是新手最容易踩的坑。每次建模后先画一个本征模的电场分布看看x方向是否存在预期的布洛赫相位调制——如果边界处的场不连续回查边界相位符号就没错了。