拿光子晶体线缺陷波导做仿真最容易遇到的一个现象是打开COMSOL能带图也画出来了但自己心里并不踏实——不知道算出来的模式是波导模式还是边界引入的杂散模式不知道k点扫得对不对也不知道“线缺陷”到底该怎么在超胞里建模。我这次用二维六角晶格空气孔光子晶体做例子把W1线缺陷波导的能带计算从头到尾捋了一遍这篇就把整个建模思路、COMSOL里的关键设置、能带后处理常用技巧和个人踩过的坑都写出来。光子晶体线缺陷波导这个名字可以拆成两部分光子晶体提供带隙线缺陷在带隙里制造一个可以局域传输的模式通道。它最有价值的点在于光束可以通过缺陷态绕过90度弯、慢光增强、微型光谱分析等场景里都能见到。对仿真来说我们要的不是“给一张好看的场图”而是把带隙范围、导带模式数、群速度这些定量信息算准后面器件的设计和优化才有依据。适合来读这篇的人是那种已经会一点COMSOL基础操作想认真做周期性结构能带仿真但还没完全把一个“超胞布洛赫边界条件本征频率扫描”流程跑通的人。1. 先搞清楚我们到底在算什么线缺陷波导与能带的基本逻辑1.1 光子晶体带隙与缺陷态一笔简单的账在二维光子晶体里介电常数周期性排列会让电磁波在某些频率范围内没有传播模式这就是光子带隙。带隙不是凭空出现的它由填充比、介电常数对比度、晶格对称性共同决定。空气孔型光子晶体通常选高折射率背景材料比如硅、氮化镓、二氧化钛然后把空气孔排列成三角晶格。三角晶格之所以常用是因为它的Brillouin区更接近圆形带隙更容易打开而且沿着Γ-K和Γ-M两个方向都能得到比较完整的禁带。如果我们在一个完整的周期性结构里抽掉一排空气孔原本严格的周期性被局部破坏带隙里就会冒出一支或多支局域态色散曲线。这些曲线落在带隙内意味着它们不能在完整晶体的体态里传播只能沿着缺陷通道走。线缺陷波导的名字就是这么来的。这里我建议你先把“完整光子晶体的带隙”算清楚再动手加缺陷。原因很简单缺陷模的能带曲线是“填补”在带隙里的如果你连完整晶体的带隙范围都不知道后面算出来的孤立的色散曲线是真是假很难判断。完整晶体仿真耗时很短几行胞、几个k点几分钟就能确认带隙位置这笔时间值得花。1.2 W1波导的几何构造超胞与晶格方向在分析线缺陷波导时不能用常规的单胞。因为缺陷破坏了y方向的周期性但沿导波方向x方向仍然具有平移周期性。正确的做法是取一个“超胞”x方向只取一个晶格常数y方向取足够多的晶格周期中间挖掉一排孔。以六角晶格空气孔为例空气孔半径设为r晶格常数设为a。W1波导的标准定义是抽掉沿某个高对称方向的一整排空气孔。在COMSOL的2D模型里我在x方向给的周期是一个晶格常数ay方向取了9排孔中间一排去掉也就是左右各留4排完整孔。超胞y方向两侧使用周期性边界条件时要特别小心相邻超胞之间的缺陷模式是否会互相串扰y方向周期数至少要取7排以上我通常用9排算完后再对比7排与9排的结果如果模式频率变化很小说明超胞尺寸已经收敛。超胞里面那块缺失的孔就是线缺陷区域。从结构上看它是一条沿着x方向连续的高折射率脊。模式的大部分能量会集中在这条脊附近以倏逝波的形式横向衰减。你在后处理看场图时如果发现某个模式能量几乎布满了整个超胞那多半不是缺陷模而是能带边缘附近的体态需要过滤掉。1.3 为什么用COMSOL而不是纯脚本很多人做这类问题会用平面波展开法或者直接写一段Python调MPBMIT Photonic Bands。平面波展开方法在处理纯周期体带结构时非常快但一旦涉及线缺陷、有限高度、材料色散或者损耗就会变得非常拧巴。COMSOL的优势在于RF模块的电磁波频域接口直接支持本征频率求解配合Floquet周期性边界条件可以非常方便地把布洛赫波矢k作为扫描参数一次性得到完整能带。另外COMSOL里面可以比较自然地加入材料吸收、非线性、热效应、结构变形等额外物理场耦合这在做实际器件仿真时是刚需。后面的步骤我会以COMSOL RF模块的“Electromagnetic Waves, Frequency Domain”接口为例。2. COMSOL模型搭建从几何、材料到边界条件的完整设置2.1 参数表与几何构建打开COMSOL选择2D空间维度物理场选择“RF Module Electromagnetic Waves, Frequency Domain”。组件定义里先建全局参数表方便后面批量修改模型。我这里给出一组常用参数你可以直接抄参数名表达式说明a500[nm]晶格常数r0.3*a空气孔半径eps_si12.25硅的相对介电常数n_sisqrt(eps_si)硅的折射率约3.5ky_scan0垂直导波方向的波矢分量kx_scan0沿导波方向的波矢分量几何构建的时候我的做法是先建立一个大矩形作为超胞背景然后用“阵列”功能按六角晶格把空气圆孔铺进去。六角晶格的孔心坐标不是简单的横平竖直而是两排之间有一个横向偏移。更省事的办法是自己写一个“阵列”定义设置两个阵发向量一个沿着x方向长度a另一个沿着60度方向也是长度a这样自动生成的就是三角晶格排列。COMSOL里布尔操作做减法的时候把背景矩形减去所有圆孔得到一个带周期性气孔的背景结构。千万记得几何单元不要用“形成联合体”默认选项和“形成装配体”混在一起如果后面要加上缺陷或者修改其中某个孔用联合体会更灵活。减法完以后在中间的孔上设置一个“隐藏”或直接在参数里加一个“defect_holes 1”控制是否减去这一排这样你还能顺便算一个完好晶体的带隙做对照。2.2 Floquet周期性边界条件的物理意义能带计算的核心公式就是布洛赫定理[ \mathbf{E}(\mathbf{r}) \mathbf{u}_k(\mathbf{r}) e^{i \mathbf{k} \cdot \mathbf{r}} ]其中(\mathbf{u}_k(\mathbf{r}))是一个与晶格周期相同的函数。在COMSOL里你不需要手动把场写成这种形式只需要在超胞的左右以及上下边界上设置周期性边界条件并且指定Floquet周期性的相位关系。左右边界是一对因为它对应真实晶格的周期方向上下边界也要设置周期条件原因不是这个方向真的周期而是它在超胞方法里代表假想的周期重复。这里有一个非常容易搞错的地方。很多人以为设置Floquet周期条件只需要勾选“Floquet periodicity”却不知道还要给“k”向量指定分量。对于沿x方向导波的线缺陷波导k向量写为((k_x, k_y))在计算沿Γ-K方向的能带时一般取(k_y0)只扫描(k_x)范围从0到(\pi/a)。我之前遇到过一种情况漏掉(k_y)默认设置里两个分量都为零最后画出来的曲线只有Γ点附近的一小段其余全是重复模式。在实际软件操作时找到“周期性条件”特征把边界类型设为Floquet然后在“k向量”栏里填入变量比如kx_scan 0把kx_scan定义成参数或者是后面参数扫描的扫描变量。相位关系由COMSOL自动生成不需要自己写表达式。要注意的是不同版本的COMSOL里这个输入框叫法不太一样有时候是“Bloch wave vector”有时候是“Floquet periodicity”你看到带有(e^{ikx})这种提示的地方就对了。2.3 网格划分与求解器选择线缺陷波导的仿真有个特点完整晶体的体态对网格相对宽容但缺陷模的能量集中在脊附近而且倏逝场延伸到周围的几个孔里网格太粗会把模式算飘掉。我的经验是空气孔边界处至少要保证一个波长里有15到20个网格单元缺陷通道区域再用更细的边界层网格。实际操作上我会先对整个超胞用“自由三角形网格”做一个较粗的划分然后选择空气孔的内边界添加“边界层网格”层数设置4到6层拉伸因子1.2。这样既不会让网格数量爆炸也不会在孔边界附近丢失模式细节。如果你用更高阶的单元比如三阶拉格朗日单元模式频率精度会明显提升代价是内存上升。求解器方面研究类型直接选“特征频率”。COMSOL会在后台求解广义特征值问题。关键是要设置一个搜索频率范围也就是“期望特征频率”。我一般先根据先验知识估计带隙归一化频率所在区域比如空气孔半径0.3、背景折射率3.5的三角晶格完整的价带和导带之间的带隙大约在归一化频率(a/\lambda0.25)到0.35附近。这里面的归一化频率换成角频率就是(2\pi c/\lambda)。在COMSOL的特征频率搜索栏里我会填一个包含这个范围的中心值和宽度比如1e15*0.3然后让求解器在这个中心附近找模式。如果搜索范围太小可能漏掉缺陷模太大则可能一次找出几十个体态让后处理变得很乱。一般先跑一次较宽范围找到缺陷模的大致频率后再收窄范围重扫。3. 扫描k点、提取能带与后处理实操3.1 参数化扫描k波矢能带曲线本质上是“特征频率随k的变化关系”。所以我们在COMSOL里要做的不是算单个k而是一个k序列。具体做法是给研究添加一个“参数化扫描”扫描变量设为kx_scan范围从0到pi/a步长可以根据精度需求选择20或30个点。如果步长太稀能带上的拐点、带边位置看不清楚太密计算时间成倍增加所以我个人一般先用21个点摸清形状再对感兴趣的区间加密。这里有个物理细节要注意布里渊区高对称点之间的路径不是简单“从0到π/a”。如果你要画完整的色散关系通常要沿Γ-K-M-Γ走一圈。二维三角晶格的第一布里渊区是六边形高对称路径一般写作Γ → K → M → Γ其中K点的坐标是((2\pi/3a, 2\pi/\sqrt{3}a))之类具体取决于坐标系。每次扫描一条路径段然后把多段扫描结果组合成一条横坐标单调增加的能带图。在COMSOL里一个“参数化扫描”可以扫描多个值序列也可以在“扫描参数”里把多条路径合并但个人觉得后者管理起来比较混乱我更倾向于分开三次扫描然后导出数据到绘图工具里拼接。需要提醒一下很多初学者会觉得“扫到K点就够了”于是直接设kx_scan到π/a。这在沿Γ-K方向时是对的但如果路径包含K转M的过程k方向发生变化单参数扫描就解决不了需要引入一个组合表达式把不同的路径段映射到不同的kx、ky值。如果只是关注线缺陷波导的导模通常扫Γ-K这一段就已经能覆盖主要的工作频率范围。3.2 特征频率求解与模式筛选参数扫描每换一个k点COMSOL都会对当前k设置下的超胞做一次特征频率求解。求解结果列表里会同时出现体态模式和缺陷模式。麻烦的是它们混在一起没有自动标号告诉你哪个是缺陷模。我的筛模式方法有两个。第一观察频率是否落在完整晶体带隙内。带隙内的模式基本可以认定是缺陷态。第二看电场能量分布图。正常的线缺陷导模能量分布在缺陷通道附近沿y方向呈指数衰减体态模式的能量则周期性地布满整个超胞。后处理时按电场模值或者能量密度做表面图一眼就能分辨。再细一点你可以把模式列表按特征频率排序然后和完整晶体的带边频率对比。完整晶体在对应k点的导带底频率差不多就是缺陷模连续分支的上限参考。这里要养成一个习惯不要只看一两个k点要在整个扫描范围内追踪同一条模式分支。有时候COMSOL在相邻k点的模式排序会跳前一k点排在第三个的模式下一k点可能排到第五个如果不做模式追踪画出来的能带曲线会来回跳看着像折线图实际上不是物理问题是排序问题。我的做法是导出每个k点全部分支的频率和标识再用脚本按照“频率接近场分布相似”的原则做连续追踪。你可以用MATLAB把数据读进去按k排序再做一个简单的最近邻匹配。不要直接在COMSOL的绘图里按默认顺序把曲线连起来那很容易画出伪断点。3.3 绘制能带曲线与群速度提取计算完成后在“结果”里新建一维绘图组横轴设为kx_scan纵轴设为特征频率freq。显示成多条曲线时默认会把所有模式同时画出来。这样是对的能带图本来就应该是多条分支重叠在一起。为了让图更能用于论文或工程判断我建议横轴改成归一化波矢(k_x a /2\pi)纵轴改成归一化频率(a/\lambda)或者保留频率本身并标注对应的晶格常数。如果你关心慢光效应就需要从能带图上提取群速度[ v_g \frac{d\omega}{dk} ]在数值上对每一条模式分支做数值微分。COMSOL里可以直接在结果表达式里写导数但平滑性通常不够。我更推荐把数据导出到外部工具先用样条拟合再求导。慢光波导的设计指标里群速度常常用(c/v_g)表达越大越慢。从这个角度看能量色散曲线越平坦群速度越小慢光效果越强。这里要特别小心能带曲线“看起来平坦”的区域并不一定等于慢光。曲线平坦且频率间隔很小说明模式密度高但实际驱动带宽也要看色散和损耗。工程上通常把带宽定义为群速度指数(n_g20)的频率范围而不是简单找斜率最小那个点。3.4 用场图验证模式性质算能带不能只看曲线场图是验证物理图像最重要的手段。我在每个关键k点比如带边、Γ点、K点附近把该k点下缺陷模的电场模绘制出来。如果模式能量紧紧贴着缺陷通道且在通道两侧的几个孔内快速衰减说明超胞尺寸足够、模式描述合理。再补充一个高级一点的操作对二维超胞模型用ewfd.normE做表面图同时在缺陷通道中心线上提取一条沿y方向的归一化场强曲线你就能看到场的横向衰减规律。这条衰减曲线可以拟合成指数衰减形式衰减长度和模式的有效折射率、垂直方向限制能力直接相关。很多人问“这个模式到底是真的导模还是泄漏模”看这个曲线最直观导模的横向场是衰减的如果衰减不明显或者衰减到某个最小值后反弹说明超胞周期不够或者模式接近连续谱。4. 常见问题、排查技巧与应用延伸4.1 模式缺失、色散曲线断点与杂散模做这类仿真最常见的三个报错和异常我按出现频率从高到低排列一下。第一扫描到某些k点时COMSOL提示“找不到特征频率”或者“搜索半径内没有特征值”。这个原因通常是搜索频率范围太小或者是网格问题导致高频模式被过量阻尼。解决办法是把特征频率搜索范围放大先把体态找出来再慢慢缩小范围排查哪一支是缺失的缺陷模。第二能带曲线出现了莫名其妙的不连续。这个问题我想重点讲一下。它不一定是物理上的带隙更多时候是相邻k点之间“模式串线”。求解器在每一个k点都是独立求解的它不会自动帮你按物理逻辑把同一条模式分支延续下去。模式在k点附近排序发生交换后连接曲线就会发生跳跃。处理办法就是我前面说的模式追踪比较相邻k点的归一化场分布重叠度或者比较特征向量之间的相关性。第三算完发现某些“模式”的频率落在带隙之外却依然显示强局域。这多半不是缺陷模而是数值伪模。产生原因通常是网格不对称导致的数值微扰或者周期性边界条件相位设置错误。检查方式就是做收敛性测试把网格加密一倍这类伪模频率会剧烈变化而真正的物理模几乎不动。4.2 超胞尺寸不够对缺陷模的影响超胞尺寸是线缺陷仿真里最容易被低估的参数。原理上缺陷态波函数在横向是指数衰减的但如果超胞y方向不够宽周期性边界条件会让相邻超胞里的缺陷模产生耦合从而抬高或压低模式频率。更麻烦的是这种耦合在不同k点上影响程度不一样会让本来应该平直的能带产生无规律的波动。所以我不建议一上来就拿一个大超胞猛算。一个更聪明的流程是先用7排孔的超胞扫描感兴趣频率范围得到初步能带然后把超胞加到9排、11排重新计算同一批k点对比缺陷模频率的变化量。如果频率变化小于1%这个超胞尺寸就可以接受。如果变化有2%甚至5%就别急着优化结构参数先把超胞尺寸改大。另外超胞改大时计算量上升但有一个省钱的办法。因为在超胞里真正需要高精度网格的地方只有缺陷通道附近远离缺陷的完整晶体区域网格相对粗糙一些也不会显著影响缺陷模频率。你可以利用“自适应网格”功能或者手动把远离缺陷的大块区域设置成较粗的网格。4.3 从能带到器件的延伸透射谱、慢光与调谐能带曲线算准之后下一步就进入实际应用验证。很多人问我能不能只用能带仿真来说明一个波导“能传光”严格说不能。能带仿真给出的是模式存在性和色散关系但真实器件里的耦合、弯折、端面反射都需要另外用频域仿真加端口激励来算透射率。最简单有效的验证方式是把超胞扩展成一段有限长度的直波导一端加端口激励另一端测量透射场。在入射频率处于导模色散曲线上时透射率会出现一个峰被完整晶体带隙挡住的频率透射率会迅速下降。把这个透射谱和能带曲线放在一起对照你会看到非常好的对应关系这是判断能带计算是否正确的一种“外部验证”。再往后延伸线缺陷波导的能带特性直接决定慢光器件的性能。如果你想调节慢光频率常用手段是修改缺陷孔的半径或位置形成“锥形慢光波导”。在COMSOL里你可以把缺陷孔半径设为一个独立的参数重新扫描能带就能看到色散曲线逐渐变平、群速度指数升高的过程。这个过程很有价值但也非常依赖能带计算的准确性因为你追求的那段平坦色带往往只占据很窄的频率窗口一个网格误差就可能把它抹平。我在实际使用中的一个体会是光子晶体器件的仿真表面上看是建个模、扫个参、画张图真正的门槛在于“你明不明白每条曲线对应哪种物理状态”。能带曲线上的一根平线、一个交叉点背后都是模式空间分布和对称性的故事。把COMSOL当成一个能打可验证的计算工具把物理图像放在前面仿真才不会越算越乱。最后分享一个小技巧每次扫完能带把所有缺陷模在缺陷中心线上的场分布快照存成一组图片按k点顺序排列成动画。你会在动画里清晰地看到模式从带边“慢慢挤进”带隙再从带隙另一端消失的完整过程这比任何一张静态场图都能帮你建立直觉。下一次再遇到疑似杂散模调出这个动画和当前模式对比判断速度会快得多。