☰
COMSOL手性超表面模拟:从本征手性BIC到三次谐波的全流程解析
2026/10/24 3:01:14 网站建设 项目流程

做纳米光子学模拟这些年,我手头出现频率最高的一组关键词,就是“comsol三次谐波、本征手性BIC、远场偏振图、手性透射曲线、二维能带图、Q因子图、电场图”。这七个词其实不是七个孤立任务,而是一整套手性超表面表征流程:先用本征频率计算拿到能带和Q因子,再在Floquet周期边界下做频域扫描,回收手性透射谱和远场偏振,在共振点取电场剖面,最后把基频场当作三次谐波的驱动源,输出非线性信号。这篇内容就是把我搭模型、出这些图、排查异常的经验全部摊开,适合正在做超表面、手性光学或非线性光学模拟的人参考,也适合刚上手COMSOL波动光学模块的初学者。

1. 先想清楚物理,再按鼠标:三个核心概念的逻辑关系

1.1 BIC与“本征手性”到底指什么

BIC,即连续区束缚态,听起来很反直觉:它的频率落在辐射连续谱里,却不会向外部辐射能量。一般微纳结构里的模式只要频率允许,就能耦合到自由空间变成辐射损耗;而BIC像是一个“漏不进水池”的漩涡,能量被局域在结构内部,理论上Q因子可以趋近无穷大。对实际超表面来说,这意味着极强的近场增强和极窄的共振线宽。

“本征手性BIC”则是在BIC基础上叠了一只手性属性。结构本身如果缺少所有镜像对称面,就具备本征手性,这种手性不依赖斜入射,也不是晶格旋转造成的假手性。把它和BIC结合后,模式既能像BIC一样强束缚、高Q,又能对手性光场产生强烈选择响应。模拟时最直观的现象就是圆偏振透射不对称、远场偏振图出现奇点、能带里某支模式的线宽突然消失。

做三维超表面模拟时,很多人第一反应是赶紧画几何、加端口、跑频域。但若课题目标确实是本征手性BIC,我建议先确认几何对称性:结构里是否还残留镜面反射对称、旋转对称、滑移对称?这些对称性会直接影响BIC存在的位置和类型。若结构手性不够“本征”,最后透射谱和远场图会变得平庸,甚至查不出任何奇异点。

1.2 三次谐波为什么绕不开高Q共振

三次谐波是非线性频率转换里的经典过程,信号频率等于基频的三倍。效率通常非常低,因为它依赖三阶非线性极化率与基频场强的三次方。提高转换效率的核心思路之一是增强基频电场,高Q共振刚好能做这件事。

BIC把能量锁在腔内,基频场强比普通传播波高出几个数量级;即便三阶非线性系数很小,(I_{3\omega}) 也能被明显放大。再加上本征手性BIC对圆偏振的选择性,可以让某个手性通道非线性信号很强,而另一通道几乎没有,这就是手性非线性光学里最想要的效果。

不过这也会给模拟带来麻烦:高Q意味着模式对网格误差极度敏感,Q越高,数值上需要越精细的计算才能收敛。很多新手一上来就把BIC共振当普通透射峰去跑,结果Q因子忽高忽低,透射曲线里Fano线型肉眼可见但数值不连续。后面我会专门讲如何规避这些坑。

1.3 线性模型与非线性模型的两次求解思路

一个常见误区是希望COMSOL一次性把基频、三次谐波、手性透射全部耦合起来求解。这在理论上可行,实际算起来代价非常大,而且收敛难度高。更稳妥的做法是拆成两代模型:

  • 第一代只算线性问题:频率扫描、能带、Q因子、透射和远场偏振。目的是确定BIC在什么频率、什么倒空间位置出现,基频场分布长什么样。
  • 第二代负责非线性输出:把第一代得到的基频复电场作为已知源,在材料域中定义三阶极化矢量,再在 (3\omega) 频率下求解非线性电磁场。

这相当于在时域谐波框架里做“泵浦不耗尽近似”,对大多数超表面THG研究是足够的。若泵浦功率很高、结构Q因数极大,非线性反作用可能显著,那就需要更严格的自洽耦合,计算量会成倍上涨。但多数论文场景下,两代法已经能很好解释实验趋势。

我把两代法涉及的接口和研究类型整理成一张表,方便规划任务:

工作内容研究类型关键设置主要输出
二维能带图特征频率Floquet边界+倒空间扫描频率-波矢关系
Q因子图特征频率复本征值换算Q-波矢关系
手性透射曲线频域端口+LCP/RCP叠加透射率、CD谱
远场偏振图频域+远场远场计算节点琼斯矢量、偏振椭圆
电场图特征频率或频域共振点附近场存储近场分布
三次谐波频域+弱贡献基频场作为非线性源(I_{3\omega}) 随波长或功率变化

2. 建模准备:模块选择、几何参数与周期边界

2.1 最小可用模型的几何和材料

我通常会选一个最简单的介质超表面单元作为起点:方形周期晶格,周期 (P) 约600 nm,用高折射率介质纳米柱做共振单元,柱高约300 nm,半径约120 nm。模拟时先用一段范围扫描几何参数,找到可能存在BIC的位置,再细调结构细节,比如把圆柱改成椭圆、加一个缺口、调整刻蚀角度,从而打破镜面对称,引入本征手性。

材料参数要区分线性与非线性两部分。在线性模拟里,折射率设为无损耗介质折射率,例如 (n=2.4) 左右,虚部设为0。这个简化能避免吸收损耗干扰Q因子判断,也能让本征频率虚部只反映辐射损耗,便于体系检查。到了非线性阶段,再给介质域补充三阶非线性极化率 (\chi^{(3)})。注意COMSOL内置材料库不一定包含你需要的那组非线性参数,使用时要自己写表达式,并确认单位是国际单位制下的 (\mathrm{m^2/V^2}) 或类似换算形式。很多时候非线性结果数量级不对,就是单位换算出了问题。

想给结构赋予本征手性,几何上最简单的做法是让单胞不再存在任何一个镜面对称面。比如椭圆纳米柱加上一个非对称切角,或者在圆盘上开一个偏心孔,再配合衬底形成垂直方向的破缺。验证手段很简单:把结构镜像后看能否通过旋转和平移重合,如果不能重合,就是手性结构。

2.2 周期性条件与倒空间路径

对于无限周期超表面,单胞四面必须用周期边界条件连接相邻晶胞。COMSOL波动光学模块提供了Floquet周期边界条件,需要指定波矢分量。注意周期晶格基矢不一定是直角坐标,如果用六角晶格,要把晶格基矢对应的波矢关系写对。

能带图扫描是在倒空间沿高对称路径完成的。方形晶格最常用路径是:

倒空间点坐标(以倒格子为单位)物理关注点
Γ((0, 0))高对称BIC、对称保护模式
X((0.5, 0))带边、镜面方向模式
M((0.5, 0.5))对角方向极大倒格矢位置
Γ((0, 0))回到原点验证连续性

COMSOL里通常用一个参数 (t) 在0到1之间扫描,然后把 (k_x, k_y) 写成 (t) 的分段函数。不要直接在派生值里手动输入几十个点,那样后续画能带图很不方便。我会把扫描数据导出成CSV,再用绘图脚本画图,横轴按路径长度重新构造,这样能带曲线才能连成一条条清晰的带。

2.3 求解器与网格的第一个经验阈值

特征频率研究本身不算昂贵,但高Q模式对网格很敏感。我常用的经验是:高折射率介质单元内部至少保证每方向20到30个网格节点,最大网格尺寸控制在工作波长的十分之一以内。如果Q因子准备报道到五位数以上,最好做一次网格收敛测试:把最大单元尺寸从 (\lambda/10) 加密到 (\lambda/20),观察Q因子变化量。若变化超过百分之十几,那这个Q值是不可信的。

求解器方面,特征频率计算时需要指定搜索范围。通常我不直接搜所有低频模式,而是在目标频率附近设定一个窄区间,比如800 THz到900 THz之间,避免低频杂散模式干扰。本征频率扫描时,还要注意把Floquet波矢参数化后,每个参数点都沿用同一个网格。一旦网格随参数点改变,模式连续性会被破坏,能带图就会像“打碎的玻璃”一样到处跳线。

3. 取出能带图、Q因子图与手性透射曲线:线性步骤全流程

3.1 二维能带图:从一个本征频率扫描说起

COMSOL里做能带图,本质是沿着倒空间路径重复求解本征频率。每个波矢点得到若干本征频率,把频率和波矢画出来,就形成能带。实现起来有几种方式,最省事的是新建一个“特征频率”研究,把Floquet波矢分量定义为扫参变量 (t),在参数化扫描中覆盖整条 (\Gamma-X-M-\Gamma) 路径。

要注意,能带图的横坐标不能简单用扫描序号代替,因为倒空间路径上各点之间的实际距离不相等。为了画图正确,我在导出的数据里额外增加一列“路径长度”:假设 (\Gamma X) 段长度是 (0.5\pi/a),(XM) 段长度是 (\sqrt{0.5^2+0.5^2}\pi/a),依次累加。把所有计算结果按这个横坐标画成散点加连线,才能得到连续的带结构。

对于本征手性BIC,能带图上的显著特征是某个模式下存在一条“暗带”或“断开带”,在某个波矢处线宽突然消失。只看实部频率不够,还需要结合Q因子图判断这条带是否真的与连续谱解耦。

3.2 Q因子图:复本征值的正确换算

Q因子最简单的定义是存储在腔体内的能量与每震荡周期消耗能量的比值。从复频率角度,本征频率可以写成 (f = f_r + i f_i),实部对应频率位置,虚部对应衰减。Q因子公式为:

[ Q = \frac{f_r}{2|f_i|} ]

很多仿真初学者会把COMSOL输出的复本征频率直接带入公式,忘记虚部本身已经是频率,前面还要乘2。如果本征频率以角频率 (\omega) 输出,公式依旧适用,只是用 (\omega_r) 和 (\omega_i) 替换。

计算路径扫描的Q因子时,我习惯在全局表达式里定义 (Q),然后在参数化扫描后生成“派生值-全局计算”,把每个波矢点对应的Q值一次性导出。对于理想BIC,理论Q趋向无穷大,实际数值会随着网格细化不断抬升。如果Q值始终上不去,多半不是物理问题,而是网格分辨率不够或边界条件离散化引入损耗。

Q因子图经常采用对数坐标,因为BIC附近Q可以从几百跳到几十万甚至更高。绘制时要把纵轴设为对数标度,才能清楚看到“Q因子急剧上升”的BIC特征。

3.3 手性透射曲线:圆偏振光的端口设置

频域扫描的模型需要添加入射端口和出射端口。在超表面模拟里,Floquet周期边界配合周期性端口,可以直接把平面波分解为衍射级次。零阶透射是最关心的通道。

COMSOL端口默认的偏振基通常是线偏振。要算圆偏振,最直接的办法是在端口激励里叠加两个线偏振分量,相位差为 (\pm \pi/2)。比如入射LCP光可以在 (x) 和 (y) 两个线偏振分量上分别加1和 (-i) 的复系数。注意不同版本COMSOL端口设置方式有差异,有些需要自定义“端口激励”字段,有些需要在“边界设置”里填偏振向量。关键是相位符号要一致,不然LCP和RCP方向反了,手性符号全乱。

射出端的偏振也要投影到LCP和RCP基上。如果端口矩阵给出了 (S_{xx}, S_{xy}, S_{yx}, S_{yy}),那么忽略公共相位因子后,LCP入射到LCP出射的同向透射可近似为:

[ T_{LL} = \frac{1}{4}|S_{xx} + S_{yy} + i(S_{xy} - S_{yx})|^2 ]

类似地,RCP入射后的同向透射 (T_{RR}) 为:

[ T_{RR} = \frac{1}{4}|S_{xx} + S_{yy} - i(S_{xy} - S_{yx})|^2 ]

具体系数还取决于端口场归一化方式,所以我并不会死记这套公式,而是先用一个对称非手性结构验证:如果LCP和RCP透射完全相等,说明端口设置和公式没有搞反。这个验证步骤强烈建议保留,否则后面所有手性结果都可能是符号错误。

3.4 透射曲线的Fano特征与CD计算

BIC本身无法被平面波激励,因为它在远场没有辐射通道。但实验和模拟里看到的“BIC共振”通常不是无穷大Q,而是有限但很大的Q,体现在透射谱上就是一条非对称Fano线型。接近BIC时,直接透射通道与共振通道之间发生干涉,产生尖锐的峰-谷结构。

手性透射曲线的核心结果之一是圆二色性(CD),定义有很多版本。我在文章里统一采用:

[ CD = T_{LL} - T_{RR} ]

也就是同向圆偏振透射率之差。如果结构有本征手性,CD谱会在BIC共振附近出现明显尖锐峰。做频域扫描时,扫描点数要足够密,分点间隔不能大于线宽的一半。高Q共振线宽可能只有零点几纳米,若按5 nm间隔扫描,很可能整个共振峰被夹在两个采样点之间,白白错过。

4. 远场偏振图、电场图和三次谐波:从共振场到非线性场

4.1 远场偏振图:琼斯矢量与Stokes参数

远场偏振图反映的是超表面衍射光在不同出射方向的偏振状态。COMSOL的“远场计算”节点可以把近场结果投影到任意方向的远场方向图,得到复振幅在远场球面上的分布。周期性结构通常只关心零阶衍射,但也可以扩展分析高阶衍射通道。

导出远场复振幅后,我通常用Stokes参数来描述偏振。定义如下:

[ S_0 = |E_x|^2 + |E_y|^2 ]

[ S_1 = |E_x|^2 - |E_y|^2 ]

[ S_2 = 2\mathrm{Re}(E_x^* E_y) ]

[ S_3 = 2\mathrm{Im}(E_x^* E_y) ]

归一化后,(S_1, S_2, S_3) 可以画在一个三维球面上,也就是庞加莱球。本征手性BIC的远场偏振图会出现典型的偏振奇点:某个出射方向上 (S_3) 接近±1,另一方向偏振突然翻转,或者在某个点附近偏振椭圆收缩成圆偏振。这些特征不是随便画张强度分布图能看出来的,必须把振幅和相位都保留下来。

实操中,我会先导出一组随横向波矢变化的远场琼斯矢量,再写一个小脚本计算Stokes参数。COMSOL内置绘图可以看 (S_0, S_1, S_2, S_3) 的二维分布,但想画偏振椭圆方向图,还是导出数据后自己画更方便。

4.2 电场图:共振模式到底长什么样

电场图是所有模拟里最直观的结果,也是判断模式是否物理合理的重要依据。在特征频率研究中,可以直接在模式解里绘制电场分布;在频域扫描里,要在共振波长处先找到透射谷或透射峰,再保存该频率点下的场解。

画电场图时,我喜欢同时看 (|E|) 的等值面图和某切面上的电场分量分布。对于超表面单元,通常要分开画出俯视截面和侧视截面,才能看出模式在柱内还是柱表面集中。若模式场强集中在几何突变的位置,说明该模式对制造误差敏感;若场强均匀分布在共振单元内部,则更像是体共振模式。

注意一个关键细节:特征频率解的绝对场幅值没有物理意义,因为本征值问题不规定激励强度。如果你从特征频率研究里直接导出场图,然后拿这个场强去算三次谐波,必须先确定归一化方式。通常用输入平面波幅度来归一化,或者把特征模式场归一化到每单位晶胞能量为1,再做非线性计算,不然THG结果的绝对数值会失真。

4.3 三次谐波:把基频场转成3ω源

三次谐波模拟是我认为COMSOL工作流里最需要组织逻辑的一步,也是最容易翻车的环节。我采用的步骤大致如下:

第一步,先在基频 (\omega) 处做频域扫描,取共振波长对应的近场复电场分布 (E_\omega(\mathbf{r}))。这个场是复向量,需要完整保存三个分量。

第二步,建立一个三阶极化源。在材料域中定义非线性极化强度:

[ P_i^{(3)}(\mathbf{r}) = \epsilon_0 \chi_{\text{eff}}^{(3)} E_j(\mathbf{r}) E_k(\mathbf{r}) E_l(\mathbf{r}) ]

为简化,可以假设标量有效非线性系数,且忽略张量各向异性。若需要严格对称性,则应完整写出 (\chi^{(3)}) 张量分量。

第三步,把这个极化源作为激励塞进 (3\omega) 频率下的电磁方程。在COMSOL里可以用“弱贡献”加进去,也可以直接在材料极化表达式中增加一个显式依赖基频场大小的额外极化项。前者更通用,适合做自由控制;后者在某些材料模型里更方便,但要注意不能影响 (\omega) 频率的本构关系。

第四步,求解 (3\omega) 场,再计算出射方向的透射功率,得到THG效率随波长或泵浦强度的变化曲线。THG信号通常远小于基频,画对数坐标。

一个让我原以为省事但实际更费时间的做法是直接在一个模型中同时求解 (\omega) 和 (3\omega) 两个频率,用多物理场耦合把基频方程与谐波方程绑在一起。这样做的优点是自洽,但计算矩阵规模成倍增加,且非线性源项容易导致频域迭代不收敛。学术界复制结果时,两代法足够,我也建议初学者从两代法入门。

THG结果的横坐标设计很有讲究。常见做法是固定泵浦功率,扫描基频波长,观察远离和靠近BIC时 (I_{3\omega}) 的差异。由于BIC共振让基频场增强好几个数量级,THG信号会在共振峰附近出现尖锐增强。若扫描波长范围过窄或步长过大,这条增强峰同样会被漏掉。

4.4 计算资源和批量导出建议

高Q手性BIC模拟经常要跑很多点,特别是圆偏振透射和远场偏振需要扫描波长,能带和Q因子又要扫描波矢。批量计算时,参数化扫描的效率直接决定项目周期。

我常用的策略是先把能带和Q因子计算做粗扫,确定BIC位置;再用细扫对准共振频率附近。不要所有任务都用相同的扫描分辨率。比如能带图用50个 (t) 点足够看趋势,Q因子图在BIC附近需要加密到200个点,透射谱则要在共振线宽范围内再细分。三层扫描密度不同,耗时可省下不少。

结果导出时,我会给文件命名带上完整条件,例如“kpt_t120_eps2_4_thickness300_LCP_scan_wl”。这看起来繁琐,但当你需要回溯某个异常结果时,文件命名能救你的命。另一个建议是关闭不必要的结果存储,只保存目标频率点的场快照,避免把整个参数空间的全场都灌进内存。

5. 常见踩坑与排查实录:新手和进阶都能用

5.1 本征模式里混入非物理模式怎么办

特征频率研究经常会在目标模式旁边冒出大量非物理模式,尤其是当搜索范围过宽时。这类模式常见征兆包括:频率非常低但形态复杂、场分布散落在整个背景介质中、Q因子异常高却没有对应的能带结构。

我通常先检查边界条件:Floquet周期边界的波矢表达式是否写对?如果是2D周期,两个方向都必须严格指定,漏掉一个就会出现周期性重复的错误模式。其次是检查搜索频率范围,把它压缩在目标频率附近,可以快速过滤低频伪模式。

还有一个很隐蔽的来源是网格不对称。即使几何完全对称,非结构化网格也可能造成轻微的边界条件破缺,从而在对称性保护的BIC位置打开微小辐射通道。这种情况下本征频率虚部不为零,但收敛行为不稳定。解决方法是采用对称性对齐的网格,或者在几何里设置对称面来强制对称性。

5.2 Q因子突然掉下来,先检查网格而不是物理

我遇到过很多次:参数扫描时Q因子在某一段突然从五位数掉到三位数,看起来很像是出现了新的损耗通道。但仔细排查后,多数原因是网格单元在那个结构位置发生了突变,或者几何参数扫描改变了某个细薄层的厚度,导致网格质量下降。

排查方法很简单:在相同几何下做一次网格收敛测试。如果细化网格后Q因子大幅回升,那就不是物理现象,而是数值色散或离散化误差。不要在物理讨论里把这种数值假象当结论。

高Q模式还有一个数值陷阱:本征频率虚部本身很小,可能和求解器收敛容差在同一量级。如果虚部结果随求解器选项轻微变化,那就说明该点Q因子已经超出该网格能分辨的范围。此时报告“Q大于某个值”比报告具体数字更稳妥。

5.3 透射率与CD的符号约定不清

手性相关模拟里,LCP和RCP的定义、入射方向、坐标手性、相位因子,每一个环节都可能让CD符号翻转。我见过同一个结构在不同文献里一个报 (CD>0),一个报 (CD<0),双方都“算对了”,只是约定不同。

所以我在自己结果里会明确写清楚:光沿哪个方向传播、电场旋转方向怎么定义、左手圆偏振是按 (E_y/E_x = -i) 还是 (+i)。更稳妥的做法是先用一个已知手性方向的结构做校准,例如从几何上观察旋向,再看模拟CD方向是否匹配。把这一条写进笔记本,能省下一整天的自我怀疑。

5.4 非线性结果中出现“无限”场强奇怪

做THG时,如果直接从特征频率场取基频场源,会出现一个容易忽略的问题:特征模式的绝对幅度是任意的,可能很大也可能很小。若直接用未归一化特征场计算 (P^{(3)}),THG结果会随求解器归一化方式而变,甚至出现数值溢出的“无限大”信号。

解决方式是在THG之前先把基频场规整到与入射功率对应的幅度。假设入射平面波幅度为 (E_0),近场增强因子为 (|E_\omega|/\max|E_\omega|),则把特征模式场乘以 (E_0/(\max|E_\omega|)),再用于非线性源计算。这样得到的THG效率才有意义。

另外,三阶非线性极化率的单位错误也会导致离谱结果。很多文献用静电单位制,COMSOL默认国际单位制,两者换算很容易出错。建议在全局表达式中显式写上 (\frac{1}{4\pi\epsilon_0}) 等换算因子,并通过一个已知小信号解析模型自检。

5.5 节约内存的几个小习惯

高Q模拟对内存的消耗主要来自全场存储。参数化扫描如果每个参数点都保存所有域的场,数据量会爆炸。我习惯把存储频率点设为“仅选中的实体”,比如只保存超表面单元内和出射端口区域的数据。远场计算则可以在结果求值阶段再算,不必把远场数据全部存下来。

网格优化方面,高折射率介质内部需要加密,但衬底内部规则网格即可。如果模型包含空气层,用映射或扫掠网格划分,比自由四面体网格更能减少单元数量。对周期结构,先验证完整单胞结果,再根据对称性选择合适简化的建模域。注意手性结构可能破坏部分对称性,简化时要确认简化后的模型仍然能复现完整结果。

我每次带新手做这个课题,都会让他们先建一个最简单的对称破缺单胞,把透射谱、能带、Q因子三条曲线全部跑通,再往里面加手性和非线性细节。这个过程看起来慢,实际上是最快的。因为七类输出之间的衔接关系,只有亲手把线性到非线性、频域到本征频率的链路走一遍,才能真正理解哪些参数影响什么结果。模拟软件只是工具,真正撑住结论的,是每一步数据背后是否清楚自洽的物理判断。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询