做纳米光子学模拟这些年我手头出现频率最高的一组关键词就是“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的位置再细调结构细节比如把圆柱改成椭圆、加一个缺口、调整刻蚀角度从而打破镜面对称引入本征手性。材料参数要区分线性与非线性两部分。在线性模拟里折射率设为无损耗介质折射率例如 (n2.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^20.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符号翻转。我见过同一个结构在不同文献里一个报 (CD0)一个报 (CD0)双方都“算对了”只是约定不同。所以我在自己结果里会明确写清楚光沿哪个方向传播、电场旋转方向怎么定义、左手圆偏振是按 (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因子三条曲线全部跑通再往里面加手性和非线性细节。这个过程看起来慢实际上是最快的。因为七类输出之间的衔接关系只有亲手把线性到非线性、频域到本征频率的链路走一遍才能真正理解哪些参数影响什么结果。模拟软件只是工具真正撑住结论的是每一步数据背后是否清楚自洽的物理判断。