COMSOL仿真魔角光子晶体激光器：能带计算与参数化建模实践

直接进入主题。最近一段时间我密集地用COMSOL做了魔角光子晶体激光器的光学模型，从能带扫描到模式分析再到参数化几何建模，来回折腾了将近三个月，终于把一套相对稳定的仿真流程跑通。这篇文章把我在这个项目里的思路、参数设置、关键操作和踩过的坑都交代清楚，后面接手相关仿真工作的朋友可以少走不少弯路。

1. 项目背景与整体思路

1.1 为什么盯着魔角光子晶体激光器

光子晶体激光器不是新概念，它靠周期性介电结构形成光子带隙，把光场限制在有源区里，实现低阈值、高品质因子的激光输出。传统的设计思路是在晶格常数、占空比、折射率对比度这些参数上做文章，结构一旦定下来，光学性质基本就锁死了。魔角光子晶体激光器多了一个自由度——层间转角，这是从石墨烯魔角物理里迁移过来的想法。

石墨烯魔角的精髓在于：两层材料之间转一个特定角度后，电子能带在费米面附近变成极平的平带，电子有效质量变大，关联效应被显著放大。光子晶体里也有类似的机制，通过扭转两层光子晶体，可以让光子色散曲线在某个频率区间变得非常平缓。色散越平，光的群速度越低，光与增益介质的相互作用时间就越长，增益阈值自然下降。这就是魔角光子晶体激光器低阈值、高相干性的物理根源。

我建模型的目标很明确：第一，把转角变化对能带结构的影响定量算出来，找到平带出现的魔角位置；第二，分析转角结构下的光学模式，看模式场分布和Q因子如何随几何参数变化；第三，用参数化几何建模把这个过程自动化，方便批量扫描。

1.2 COMSOL在光学仿真中解决什么问题

COMSOL这套东西的好处是灵活，尤其适合做这种带多物理场耦合和复杂几何参数扫描的问题。它不像FDTD那样纯粹做时域推进，也不像RCWA那样绑定周期结构，而是基于有限元方法解偏微分方程，你可以在同一个平台里同时做特征频率分析、频域响应和参数化扫描。

对魔角光子晶体这个课题，COMSOL的核心优势体现在三块：

• 能带计算可以直接用特征频率求解器配合Floquet周期性边界条件完成，不需要额外写复杂的本征求解代码；
• 几何建模支持参数化，转角、晶格常数、填充因子都可以定义成全局参数，改一个数重算一遍，批量扫描非常顺手；
• 后处理工具完善，能直接提取本征频率、Q因子、电场分布、能带曲线，还能导出数据画图。

当然，COMSOL的有限元方法也有它的脾气，网格不合适会出伪模，周期性边界条件设置不对能谱会断，这些后面会详细展开。

2. 能带计算的物理基础与COMSOL实现

2.1 光子晶体能带到底算的是什么

在动手之前，先要明确能带的物理图像。电磁波在周期性介质中传播时，满足Bloch定理，电场可以写成布洛赫形式：

E(r) = u(r) exp(i k r)

其中u(r)是周期函数，k是波矢。给定一个k，求解麦克斯韦方程组会得到一组离散的本征频率，把k扫过第一布里渊区，得到的频率-波矢关系就是能带图。

平带意味着什么？群速度正比于dω/dk，能带越平，斜率越小，群速度越低。对激光器来说，低群速度意味着光在腔内待的时间更长，增益介质能够更充分地与光场相互作用，阈值增益会降低，Q因子会升高。魔角光子晶体激光器的核心设计目标，就是在目标波长附近造出一个平带。

2.2 COMSOL特征频率求解能带的基本流程

在COMSOL中计算能带，我用的是波动光学模块的特征频率求解。几何上只需要建一个晶胞，不需要建整块晶体。晶胞四周加Floquet周期性边界条件，形如：

E_dest = E_src exp(-i k (r_dest - r_src))

这里k是布洛赫波矢，在COMSOL里通过周期性边界条件中的k向量参数定义。扫描能带时，把k向量从Γ点扫到M点再到K点，再回到Γ点，记录每个k对应的特征频率，然后画成曲线。

具体的操作步骤大概是这样的：

1. 新建模型，选择二维模型，波动光学模块，特征频率研究；
2. 定义全局参数：晶格常数a、转角theta、材料折射率n1、n2、占空比等；
3. 建几何，用两个多边形或者两个圆孔阵列表示上下两层光子晶体，中间留出很薄的间隙层；
4. 添加Floquet周期性边界条件，设置布洛赫波矢分量kx、ky；
5. 设置材料折射率，注意有源区可以选择复折射率来近似增益；
6. 网格划分，默认物理场控制网格先跑通，后续细化；
7. 研究设置里使用特征频率求解器，设置所需最低频率范围，求解；
8. 后处理里提取特征频率，用辅助扫描扫k向量。

2.3 能带扫描时布洛赫波矢的写法

这一步是很多人容易踩坑的地方。COMSOL的Floquet周期性边界条件里有两种设置布洛赫波矢的方式：一种是直接用周期性子域，另一种是定义“周期边界条件”并在参数里写入波矢分量。

我习惯把k向量归一化到第一布里渊区边界。对于二维六角或者四方晶格，首先要算出倒格子基矢，然后把k点写成倒格子坐标。比如四方晶格的倒格子基矢是2π/a沿x和y方向，那么k点写作：

kx = kx_norm * pi / a ky = ky_norm * pi / a

其中kx_norm、ky_norm在0到1之间变化。这样扫描比较直观，布里渊区边界也容易对应上。

2.4 能带仿真参数设置参考

给一组我实际用过的基础参数，结构是双层光子晶体平板，单层为三角晶格空气孔阵：

这些参数启动后先用特征频率求解器算一下，确保网格尺寸在λ/5左右，太粗的网格会把能带算偏。

3. 模式分析与参数化几何建模的关键细节

3.1 模式识别：别被伪模带偏

算完特征频率后，COMSOL会给出一大堆本征频率，有些是真实的物理模式，有些是数值伪模。伪模的特征很明显：场分布没有明显的周期性，常常集中在几何的角上或者网格密集区域，频率值对网格密度非常敏感。

判断真实模式的方法我通常有三个：

• 看电场模分布是否符合晶胞对称性，真实模式通常具有明确的对称性或反称性；
• 对同一结构用两套不同密度的网格各算一遍，频率差小于0.1%才认为是收敛的；
• 检查模式是否满足布洛赫定理的预期——改变k值，频率应该平滑连续变化，如果有跳变点需要警惕。

还有一个更物理的方法：把计算区域放大到超胞，比如5x5个晶胞，算出来的模式频率如果和单胞加周期边界的结果吻合，那模式就是可信的。这个方法虽然慢，但可以作为交叉验证。

3.2 Q因子怎么从仿真里提出来

Q因子表征腔体的储能能力，对激光器来说几乎是核心指标。在COMSOL里提取Q值有三种常见做法：

第一种，直接用复特征频率的实部和虚部算。特征频率求解器可以给出复数值，虚部反映了损耗或者增益。Q = Re(f) / (2 |Im(f)|)，这里的Im(f)在无源结构里来自辐射损耗和材料吸收。这个方式最简单，但前提是你设置的材料折射率里带了虚部，或者通过完美匹配层吸收了辐射。

第二种，在频域做一次扫频，对模式中心的电场能量谱线做洛伦兹拟合，用中心频率除以线宽得到Q值。这种方法更准确，但需要在特定频率附近做精细扫描，计算量稍大。

第三种，用本征模分析之后做瞬态仿真，看场衰减的包络时间常数，Q = 2π f τ。这个计算量最大，通常只在最后验证阶段使用。

我平时先用第一种方法快速筛选，候选结构再用第二种方法复核，效率和精度都比较满意。

3.3 参数化几何建模如何设计

魔角结构的几何建模是整个项目里最费心思的部分。所谓魔角，就是上下两层光子晶体之间有一个相对旋转角度。建模时如果手动移动几何，每改一次角度都要重新装配接触对，非常繁琐，而且旋转后的周期结构并不严格满足单胞的周期边界条件——这是魔角模型最棘手的地方。

我采用的方案是把上下两层分别建成独立的组件，旋转角度设为一个全局参数theta，通过旋转特征完成相对定位。这里有一个关键点：上下层旋转后，整个结构不再具有平移周期性，严格来说不能用单胞加Floquet边界来处理，需要构建一个足够大的超胞来近似。

超胞的大小选择是一个权衡。如果超胞太小，旋转后上下层的晶格失配带来的边缘效应会污染能带；如果太大，计算量爆炸。我试过从3x3、5x5到7x7的超胞，发现对600 nm晶格常数、220 nm平板厚度这个体系，5x5超胞在能带低频区的结果已经基本收敛，高频区还差一些，最终选择了7x7作为平衡点。

3.4 参数化扫描的策略

参数化扫描的目的不仅是找平带，还要看模式演化。我设置的扫描变量是转角theta和气孔半径r，偶尔也会扫层间距d。

这里推荐一个技巧：先用粗扫描确定大的趋势。以0.5度为步长从0度扫到5度，看能带图里哪几个角度出现了明显的平带特征，然后再用0.1度甚至0.05度的细步长在那个角度附近精确找极值。

参数化扫描有个容易忽略的问题：COMSOL默认每个参数点都从头开始求解，前面的解不会作为下一个点的初值。这会导致很多没有物理意义的模式被计算出来，白白耗计算时间。解决办法是在研究的“扫描”设置里启用“从上一个解继续”选项，或者手动把上一个参数点的解设为下一个点的初始猜测，这种做法对连续性比较好的模式能大幅提升求解效率。

4. 实操过程与核心环节实现

4.1 从零搭建魔角光子晶体模型

这里给出一套可以直接照着做的流程，我以7x7超胞、三角晶格空气孔为例。

第一步，参数定义。在全局参数里写好：

• a = 600e-9
• r = 0.3a
• h = 220e-9
• theta = 1.5
• n_slab = 3.45
• n_gap = 1.0
• d_gap = 20e-9

第二步，几何建模。先在xy平面画一个7x7的超胞区域，尺寸是7a x 7a。然后通过阵列画出空气圆柱孔，这些孔代表下层光子晶体。把整个几何复制一份并整体旋转theta角度，放在上层。上下层之间用d_gap分隔。

注意：这里不能采用二维多孔板直接拉伸，需要处理平板上下表面和空气孔的边界层。实际做法是建三维模型，在z方向设置平板厚度h，上层和下层之间加薄空气层。

第三步，材料设置。整个域设置折射率n_slab，空气孔和间隙层设置为n_gap。如果想模拟增益介质，把n_slab的虚部设为负值，例如3.45-0.01i，负虚部代表受激放大。

第四步，边界条件。超胞的侧面需要处理电磁辐射问题。我试过两种方案：一是全部设为完美电导体边界，简单但会让泄漏模式变成寄生模；二是加完美匹配层吸收边界，更符合开放腔的物理情况。推荐用PML包裹超胞四周，PML厚度设为0.5a，与超胞间隔至少0.5a的空隙。

第五步，设置特征频率求解器。目标频率范围先设宽一点，例如100 THz到400 THz，让求解器把所有可能的模式都算出来。扫描完后对模式进行筛选。

第六步，网格划分。超胞尺寸为7a，扣除空气孔，几何细节复杂。建议用四面体网格，最细网格尺寸设为a/10到a/12。空气孔内部网格适当细化，因为场在孔界面附近变化剧烈。

4.2 扫描k向量的自动化实现

能带计算需要在布里渊区边界扫点。COMSOL没有内置的k路径扫描功能，得自己想办法。

我的做法是使用“辅助扫描”功能。先定义一个扫描变量s，范围从0到1，然后计算：

如果扫描路径是Γ-M-K-Γ，就把整个路径分成三段，每一段用参数s映射到对应的k坐标。比如Γ到M段：

kx = s * (pi / (a * 2)) ky = s * (pi / (a * 2))

这里需要根据具体晶格方向调整。每条路径段的端点坐标要提前算好，写成分段函数的形式。

在实际操作里，我通常把这三个段写进同一个参数表达式，用if函数判断s落在哪个区间，然后输出对应的kx、ky。这样一次参数化扫描就能覆盖整条布里渊区边界路径。

这个方案生成的能带图有个特点：横坐标不是线性的，因为Γ-M段和M-K段的长度不同，做图时要把每个段的实际k距离算出来，否则横轴会被压缩拉伸得很难看。

4.3 后处理与数据提取

后处理主要做三件事：提取特征频率、绘制能带曲线、可视化模式场分布。

提取特征频率的方法是先跑完参数化扫描，再在“结果”里用“数据集”选中不同k点的解，导出每个解的频率实部、虚部和Q值。我通常把数据以表格形式导出，再用脚本处理成能带图。

模式场分布的后处理，优先看Ez分量。对横电模平板系统，面外电场分量能够很直观地反映模式的横向分布。得到的模式场加上箭头图可以画出能量流动方向，分析回音壁式的旋转模还是驻波型模式。

参数化扫描得到的数据量大，建议在COMSOL里先做一次筛选，只导出模式下概率大于一定阈值的特征频率，减小数据量。具体阈值要靠对模式场分布的观察确定。

5. 常见问题与排查技巧实录

5.1 特征频率总是出现零模或负频率

刚起步时最容易遇到的一个问题就是一算特征值，出来一堆接近零或者负的很奇怪的频率，物理上完全不成立。排查方向有两个：

第一，检查几何里是否存在独立悬浮的域。比如做旋转操作时，上层和下层之间如果留有微小空隙，但没有接触或者耦合边界条件，这部分可能形成弱连接域，产生数值泄漏模。第二，检查求解器设置，尤其是特征频率数目。COMSOL默认计算6个特征值，如果设置的目标频率范围太小，求解器会返回边缘的低质量解。把目标频率范围扩大，或者把特征频率搜索数量增加到20、30个，问题一般就解决了。

我在建模时吃过一个亏：上下层之间的间隙层用了四面体网格，厚度只有20 nm，但宽度有4.2微米，网格纵横比非常大，导致求解缓慢而且结果不稳定。后来我改用扫掠网格处理间隙层，问题才得到缓解。

5.2 能带曲线在布里渊区边界不连续

Floquet周期性边界条件下，能带应该在布里渊区边界连续，但实际扫描时经常出现曲线断裂或跳变。这个问题大部分时候出在k向量的映射方式上。相位因子在边界两侧应该有连续性，如果kx、ky设置不连贯，边界处自然不连续。

解决办法是确保扫描路径上每个点的k向量都精确落在布里渊区边界上，并且相位因子写法一致。另外，网格的周期必须严格匹配几何的周期。使用超胞后，晶格常数的有效值变了，周期性边界条件的距离也要按超胞尺寸设置，不能再用原来的单胞常数。

5.3 参数化扫描计算量爆炸怎么破

7x7超胞加PML加精细网格，一次特征频率求解大约需要5到10分钟，听起来还好，但如果转角步长0.1度从0扫到3度，那就是31个点，再加上每个点要算20个特征频率，整体耗时几个小时到一天。想提速，我有几个心得：

• 粗网格起步。先用最大的网格尺寸算出能带趋势，锁定目标角度区间后再用细网格精确计算。粗网格结果虽然频率偏一点，但趋势是对的。
• 利用模态渐变。扫描前把上一个转角点的解作为下一个点的初始值，这样求解器收敛更快，而且能确保同一个物理模式在不同转角点的连续性。COMSOL的“从先前步骤继续”功能在特定条件下可以做到这一点。
• 只计算感兴趣频率范围附近的特征值。如果目标波长在某个频段，可以在求解器设置里把搜索范围缩小，计算量能减少将近一半。
• 超胞尺寸可以先从3x3开始。3x3的转角结构算得很快，用来做参数扫描初步筛选用，确认某个角度附近有平带趋势后，再用7x7甚至更大的超胞做精确验证。

5.4 一个容易被忽视的网格连续性问题

魔角结构旋转之后，上下两层的网格很难保证界面处的严丝合缝。如果不做处理，界面上会出现网格节点不匹配，导致无法正常计算。

我采用的方案是设定装配模式为“非连续装配”，然后在上下层的接触界面上添加一致性对或者周期性对，让COMSOL在求解时自动在界面处做映射插值。但需要注意，如果层间间隙很小，界面映射会带来额外的数值误差。后来我把间隙层厚度提高到20 nm以上，误差才降到可接受范围。间隙再小，耦合是更强了，但数值精度很难保证，这也是这一结构的仿真瓶颈。

5.5 优化参数的实用经验

仿真最终要落到激光器设计上。我最关注的是两个量：Q因子和平带宽度。转角变化会同时影响这两个量，但影响方向不完全一致。转角接近魔角时，Q因子通常升高，但带宽会变窄，增益带宽也随之变窄。对激光器性能而言，最好是兼顾高的Q值和一定的工作带宽。

一个可供参考的做法是：先画出Q因子随转角变化的曲线，找到峰值，再回到能带图里看这个转角对应的群速度。如果群速度低于光速的千分之一，说明慢光效应明显，对激光器有益。但群速度太低也会导致线宽过窄，对泵浦光谱的匹配要求更苛刻，这个度要在具体器件设计里取平衡。

6. 工具选型与性能权衡

6.1 为什么不用FDTD或RCWA做这个课题

很多人会问，能带计算用平面波展开或RCWA不是更成熟吗？确实，常规光子晶体能带用平面波展开法最快。但魔角结构有个特点：上下层的旋转破坏了严格周期，整块结构实际上是准周期的。平面波展开法在这种几何下展开项数会非常多，收敛的速度很慢。RCWA擅长多层周期薄膜结构，但处理有源介质和复杂几何时需要调整的细节很多。

COMSOL的性能优势在于做磁角结构的参数化扫描。它能通过统一的接口把几何、材料、边界条件和求解器串起来，改一个转角参数就可以自动重新建几何、重画网格、重新求解。这在工作流效率上比脚本化拼接口要省事。当然，也有缺点——计算量偏高，资源占用大，不太适合大规模倒格矢扫描。

6.2 内存与算力配置建议

7x7超胞建PML的三维模型，自由度大概在300万到500万之间，加上特征频率求解器，建议内存至少32 GB起步，64 GB会更加从容。CPU方面，COMSOL的默认求解器对多核支持不错，8核以上有明显的提速，不过特征是频率分解算法对内存带宽更加敏感，CPU慢一些顶多多等一会儿，内存不够直接报错。

如果实验室机器比较紧张，可以先把超胞降为5x5跑完整流程，结果做预研完全够用，最终报数据时再用7x7。这样既节省了调试时间，又保证了最终结果的质量。

7. 数据解读与设计指导

7.1 从能带图上怎么判断魔角出现

把不同转角下的能带图放在一起对比，重点关注低频区的某一条带。随着转角从0度增大，这条带的斜率会逐渐变小。当误差角的某个值得附近，这条带几乎变成水平线，这就是“魔角”。继续增大转角，能带又开始变斜。

还有一个辅助判据：观察模式场分布。平带对应的模式通常表现出局域化的特征，电场能量集中在小区域内，而不是均匀铺开。用COMSOL后处理里的能量密度图查看，非常直观。

7.2 参数化几何建模对器件的意义

参数化建模不只是为了出几张图。激光器设计中，需要对谐振频率做精确调谐。晶格常数决定了工作波长，转角是调节模式耦合的旋钮，气孔半径控制带隙宽度。在COMSOL里把这几个变量都定义成全局参数后，可以对目标波长做“反向设计”——先定好想要的频率和Q值，然后扫参数空间找最佳组合。这就是参数化几何建模在实际工程中最有价值的应用。

8. 实操过程中的一些个人心得

最后分享一个我反复调整才意识到的问题：不要一上来就用完整的大小配置把所有细节做满。魔角光子晶体激光器的仿真链条长，从几何建模到能带计算再到模式分析，每一步都容易出错。如果一开始就建7x7超胞加精细网格加PML，出问题之后难以定位到底是几何的问题、网格的问题还是边界条件的问题。

我现在的习惯是先在一个大大的简化版模型上把流程跑通——用2x2或者3x3超胞，粗网格，甚至是二维模型，先用二维模型验证能带趋势，确认物理效应符合预期之后，再切换到三维精细模型做定量计算。这样试错成本低，也更容易建立对模型物理行为的直觉。

另外，在转角扫描前务必做一次网格收敛性验证。选取一个中间转角，比如1.5度，分别用粗、中、细三套网格计算同一条带的最低频率，确保结果差异在1%以内。这一步做好了，后面所有扫描才有意义，不然可能把网格误差当成物理趋势，白折腾。

COMSOL的官方案例库里有光子晶体能带的模板，但是魔角结构这种旋转双层的案例很少见。建议把官方模板的结构理解透，尤其要理解Floquet边界条件的相位设置原理，然后自己动手改几何、加转角。理解原理再动手，比照抄案例有效得多。