拿光子晶体线缺陷波导做仿真,最容易遇到的一个现象是:打开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方向给的周期是一个晶格常数a,y方向取了9排孔,中间一排去掉,也就是左右各留4排完整孔。超胞y方向两侧使用周期性边界条件时,要特别小心相邻超胞之间的缺陷模式是否会互相串扰,y方向周期数至少要取7排以上,我通常用9排,算完后再对比7排与9排的结果,如果模式频率变化很小,说明超胞尺寸已经收敛。
超胞里面那块缺失的孔就是线缺陷区域。从结构上看,它是一条沿着x方向连续的高折射率脊。模式的大部分能量会集中在这条脊附近,以倏逝波的形式横向衰减。你在后处理看场图时,如果发现某个模式能量几乎布满了整个超胞,那多半不是缺陷模,而是能带边缘附近的体态,需要过滤掉。
1.3 为什么用COMSOL而不是纯脚本
很多人做这类问题会用平面波展开法,或者直接写一段Python调MPB(MIT 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”。组件定义里先建全局参数表,方便后面批量修改模型。我这里给出一组常用参数,你可以直接抄:
| 参数名 | 表达式 | 说明 |
|---|---|---|
| a | 500[nm] | 晶格常数 |
| r | 0.3*a | 空气孔半径 |
| eps_si | 12.25 | 硅的相对介电常数 |
| n_si | sqrt(eps_si) | 硅的折射率,约3.5 |
| ky_scan | 0 | 垂直导波方向的波矢分量 |
| kx_scan | 0 | 沿导波方向的波矢分量 |
几何构建的时候,我的做法是先建立一个大矩形作为超胞背景,然后用“阵列”功能按六角晶格把空气圆孔铺进去。六角晶格的孔心坐标不是简单的横平竖直,而是两排之间有一个横向偏移。更省事的办法是自己写一个“阵列”定义,设置两个阵发向量:一个沿着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_y=0),只扫描(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/\lambda=0.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_g>20)的频率范围,而不是简单找斜率最小那个点。
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点顺序排列成动画。你会在动画里清晰地看到模式从带边“慢慢挤进”带隙,再从带隙另一端消失的完整过程,这比任何一张静态场图都能帮你建立直觉。下一次再遇到疑似杂散模,调出这个动画和当前模式对比,判断速度会快得多。