我在自己的课题里做了一阵子铌酸锂光子晶体超表面的二次谐波模拟最初也是一头雾水BICs三个字母看着玄乎LN的材料参数一大堆COMSOL里加非线性源又经常不收敛。但把原理一层层拆开就会发现这个课题其实很“顺”——连续域束缚态负责把Q因子拉满铌酸锂提供二阶非线性超表面负责在亚波长尺度把两者咬合在一起COMSOL给你一个把这三个物理故事变成数据和曲线的工作台。这篇文章想从一个实际建模者的角度把“基于连续域束缚态BICs的铌酸锂二次谐波COMSOL光子晶体超表面模拟”这件事完整串一遍为什么会选这样的结构建模前要想清楚哪些参数COMSOL里怎么从找BIC一路算到SHG效率最后再把我在仿真里踩过的坑和排查方法晾出来。适合正在做超表面、非线性纳米光学或者只是单纯想用COMSOL复现一个带点物理深度的周期结构案例的读者。1. 为什么把BICs、铌酸锂和超表面绑在一起1.1 BIC到底是什么为什么能把光“锁”在纳米尺度BICbound states in the continuum连续域束缚态。名字很绕但可以用一个日常画面理解演唱会现场满场都是光可总有一个人站在舞台阴影里不往台下发光身边明明就是铺满光线的“连续谱”他却始终不辐射。在周期性纳米结构里就存在这样的模式它的频率落在辐射连续谱中但因为某种对称性最常见的是面内C2对称与所有出射辐射通道正交导致它完全局域在结构内部既不向外传播也不衰减。这种模式的辐射通道被“焊死”了所以理论Q因子可以做到无限大。把这种几乎不损耗的模式放到非线性过程里意味着什么二次谐波效率不仅依赖材料的二阶非线性系数更依赖局域场强度。共振越强泵浦光在场增强区域积累的能量越高产生的非线性极化越强。BIC恰好提供了最干净的强共振没有辐射损耗只有材料吸收、制造误差这些“外来者”才会限制Q值。在无耗散仿真模型里BIC本征频率的虚部会小到接近机器精度这在工程上就相当于白送增强倍数。1.2 铌酸锂的第二个身份二阶非线性大户铌酸锂LiNbO3在很多人印象里是电光晶体、压电晶体其实它一直在非线性光学里当主力。透明窗口覆盖可见到中红外二阶非线性系数里d33能到27 pm/V左右在已经成熟量产的光子材料中算一流水平。另一个优势是薄膜化现在LNOIlithium niobate on insulator铌酸锂薄膜技术成熟300到600纳米厚的单晶薄膜可以直接被刻蚀出柱体和孔洞做超表面不需要像BBO晶体那样费劲做镀膜结构本身就是标准光子集成工艺。它的压电特性看起来和SHG无关但意味着这套结构以后可以和声学、电学调控结合起来玩。简单说BIC负责把Q因子推向极致铌酸锂负责把场增强转成实实在在的非线性信号两者配合就是11大于2。1.3 超表面加准BIC等于把“增强”做成可调旋钮如果BIC完全不辐射那怎么把光“接”进去这就需要准BICquasi-BIC。通过微调结构比如把圆柱改成椭圆或者让柱体稍旋转一点原本被禁闭的对称性被打破模式获得一个可调的小辐射率。一面漏出来与外部入射光耦合一面又保持极高的Q因子这种“可控泄漏”正是超表面非线性增强需要的状态。和传统相位匹配方案比超表面的非线性过程发生在亚波长厚度上不依赖毫米级晶体、不依赖严格相位匹配角度只要泵浦频率与共振频率对齐、模式重叠足够好就能输出二次谐波。BICs和铌酸锂的组合本质上是在回答一个问题如何在纳米尺度把非线性过程做“响”。2. 建模之前先把物理模型想清楚2.1 结构选型椭圆柱、LNOI、周期怎么定在COMSOL里动手的第一步不是建几何而是把结构和材料的“故事”定下来。我的模型采用了最常见的方案二维方晶格排列的椭圆LN纳米柱柱子长轴沿x方向短轴沿y方向衬底是SiO2上方是空气。LNOI薄膜厚度就是柱高我习惯取300到500纳米周期p取900纳米左右。对应1064纳米泵浦波长时空气侧因为波长大于周期只有零级衍射传播而SiO2衬底侧因为介质中波长变短、周期相对变大会出现多级衍射通道——这些通道就是BIC所寄居的“连续谱”。单位胞内只画一个椭圆柱周期边界把它延拓成无限阵列。一个容易忽略的点是占空比。柱子长轴和短轴不能乱给要保证目标频率附近存在可识别的导模。我的经验是先用几组粗扫描看透射谱哪个频率出现尖锐的Fano线型再回来精细调整周期和占空比。整个过程我都会把几何尺寸设为全局参数lam0、P、h、a、b这样后面做参数扫描时不用改模型结构只要改数据集里的参数值。2.2 LN的折射率与二阶张量别在坐标系上翻车LN是单轴晶体寻常光折射率n_o和非寻常光折射率n_e都要考虑色散。1064纳米附近的典型数值大约是n_o约等于2.232、n_e约等于2.156到了532纳米数值会升到n_o约等于2.324、n_e约等于2.234。COMSOL材料库不一定会内置准确色散我会手工输入Sellmeier色散公式或者把这两个波长的折射率做成插值表。更关键的是晶轴方向Z-cut的LN薄膜光轴垂直膜面而大多数超表面共振的电场集中在面内d33这个大分量根本用不上。我更推荐X-cut或Y-cut LNOI让极轴躺在面内让泵浦偏振方向与极轴对准。接下来是二阶非线性张量。LN属于3m点群在缩并记法下d矩阵的非零项包括d33约27 pm/V、d22约2.1 pm/V、d31约4.6 pm/V、d15同数量级。如果光轴沿着x基频电场E_x被共振放大后主要贡献来自d33路径P_x(2ω)约等于2ε0 d33 E_x平方。要在COMSOL里把这个关系加进去必须把各向异性折射率和非线性张量同时放到正确坐标系下。常见错误是材料默认z轴就是光轴而结构里光轴已经被旋转了算出来的SHG效率能整整差两个量级。这一步要是漏了后面全白做。2.3 我们到底要找出哪种模式BIC不是随便一个驻波而是能带上某个点的特殊模式。在二维方晶格Gamma点kxky0处C2对称保护的BIC最常见的身份是面内偶极矩为零的混合模式。它有一个特点特征频率实部落在辐射连续谱中虚部在理论上为零。COMSOL本征模计算要找的就是这样的一个特征频率实部落在目标波段附近虚部小到接近数值噪声把它的模场画出来能看到能量主要窝在柱体内几乎没有向外辐射的“触角”。动手建几何之前先把这个物理图像装进脑子里后面读结果会少走很多弯路。3. COMSOL实操从本征模到二次谐波3.1 模型框架边界条件、端口、PML一次配齐我用的是“电磁波频域ewfd”接口三维单胞模型。周期边界用Floquet周期边界条件x和y方向各一对分别指定Floquet波矢分量kx和ky。底部SiO2衬底不能无限延伸切出一段厚度以后加完美匹配层PML顶部空气层同样处理。PML的厚度至少要达到自由空间波长的一半并且与结构保持安全距离否则倏逝场会被PML表面反射回来。做频域平面波激励时我推荐用“背景场加散射场”的组合背景平面波设成沿z传播、偏振沿x计算出来的散射场直接反映了结构辐射行为。这样比单纯用端口设置入射波更直观后期提取透射谱也更干净。模型里LN的介电属性一定要指定为各向异性折射率矩阵用对角形式再配合坐标旋转转换到光轴沿x的X-cut方向。3.2 第一步先在本征模里找到BIC这一步不追求绝对数值而是找“候选者”。添加研究步骤选择“特征频率”目标频率设置在目标波长附近比如1064纳米对应约282THz让求解器算出20到30个模态。我会先做一次粗网格的频域透射谱把透射谷的位置找出来然后把特征频率的搜索区间缩小到几个THz以内这样求解器不会到处乱找。在Gamma点扫描时BIC的特征频率虚部会突然变小。如果结构完全对称计算出的特征频率虚部是浮点误差量级一旦把椭圆率参数α设成0.01虚部立刻变大若干个数量级这基本就是准BIC无疑。Q因子用Q等于Re(f)除以两倍Im(f)的绝对值来估算。对称情况下这个值会高到几百万以上对称破缺后降到几百到几万。如果这里找不到虚部特别小的模式先检查Floquet周期边界有没有按正确晶格加约束。哪怕一个方向的周期边界漏掉BIC的对称保护就会失效。3.3 第二步破开对称性用频域扫描把准BIC叫出来结构完全对称时BIC与外部光场的耦合是零透射谱里看不到共振峰。要让入射光真正“接”进去把椭圆率α设成非零值。我一般定义长轴直径为d0乘以(1α)短轴直径为d0除以(1α)这样面积基本不变只是打破了柱子的C2对称。保留参数扫描频率从共振实频附近往两边扩展看透射谱是否出现一个尖锐的非对称线型。注意α等于0时没有峰α越大峰越宽Q越低。很多刚上手的同学看到谱线不对称就以为模型错了其实那是Fano共振的典型特征说明准BIC和辐射连续谱发生了干涉。扫描频率间隔不能太大。准BIC的谱线很窄如果频点间隔大于峰宽的五分之一很容易漏掉共振。我一般先用较粗间隔定位再在共振区局部加密扫点。在共振频率处后处理画出模场模值除以入射场模值的分布增强可以达到几十甚至上百倍。把泵浦偏振、k点、模式对称性都记录下来后面SHG计算会用到。3.4 第三步把基频场变成2ω源算SHG这一步让我绕了不少弯子。标准做法是两步法也叫小信号近似先不考虑SHG对基频的反作用在泵浦频率ω下求解线性问题得到基频场E_ω(x,y,z)再把非线性极化P_NL(2ω)等于ε0乘以χ(2)张量对E_ωE_ω的缩并当成一个已知源项重新在2ω频率下求解麦克斯韦方程。具体到COMSOL我在同一个组件里建两个电磁波频域求解步骤第一个研究求解f0第二个研究求解2f0。在第二个研究里加入“外部电流密度”域节点把电流密度J_2ω写成负的2jω乘以P_NL的形式。P_NL的分量根据2.2节里整理的LN二阶张量写如果光轴沿x且只保留d33主导表达式就是P_x等于2ε0 d33乘E_x平方。跨研究引用基频场这一步不同版本语法略有差异。6.x版本里可以直接用withsol函数把第一个研究的解插值到第二个研究里也可以通过“数据集”手动选择解。如果版本实在不顺手就先把基频电场导出成表格再用插值函数导入笨方法但绝对可靠。算2ω时网格要比基频更细因为倍频之后有效波长变短空间分辨率要求更高。后处理时在输出端面积分2ω频率下的平均能流法向分量除以基频输入功率就得到SHG转换效率。我始终建议用基频和倍频两步独立求解不要一上来就做双向耦合自洽迭代那个对初值极敏感日常科研里多数场景根本不需要那么高的“完整度”。4. 结果分析怎么看Q怎么扫参数4.1 Q因子、谱线形状、场增强怎么读本征模给Q透射谱给“实际耦合状态”两者要对照着看。下面是一组典型的趋势示例具体数值依赖你的周期、厚度和材料参数别拿它当绝对值椭圆率αQ因子本征模估算透射峰线宽局域场增强倍数010的6次方量级几乎无法分辨最大0.011万到10万很窄显著0.05几千明显可测高0.1几百较宽中等有意思的点在于Q最大时不代表SHG效率最高。α增大虽然拉低Q却也提高了准BIC和入射光的耦合效率。实际做参数扫描时SHG转换效率通常会随α出现一个非单调的峰值。找到这个峰比单纯追求最高Q更重要。后处理时也要注意场增强的积分对象不要只盯某个点的最大值要看共振体积内的有效模场因为SHG是一个体积效应模场局域在很小的热点和弥散在更大区域后者对效率贡献往往更可观。4.2 椭圆率α一扭Q从几千调到几十对称破缺准BIC最典型的行为是Q大致与α平方成反比。你可以把α设成0.005、0.01、0.02、0.05、0.1这样一组值分别跑本征模把Q取对数后画出来。如果数据点在log-log坐标下基本落成一条斜率约为负2的直线那恭喜你控制这个BIC的对称性确实被破对了。如果斜率明显偏离说明你破缺的对称性可能不是决定该模式辐射的那个对称性或者有其他泄漏通道在捣乱。还有一种情况当α比较小时网格噪声会被当成辐射损耗导致Q饱和在某个值附近这时需要加密网格重新验证。4.3 批量扫参用MATLAB控制COMSOL省出整个午饭时间单点扫描跑不了几次但要做α、周期、厚度三参数联合扫描手工改参数能累死人。我是用MATLAB来控制COMSOL的。先把模型保存成参数化良好的mph文件把几何尺寸、边界条件、研究频率全部挂到全局参数上。然后在MATLAB里用mphopen打开模型循环里改参数、求解、提取结果再用mphsave存到指定目录。COMSOL 6.4在Linux服务器上支持headless模式配合MATLAB接口几百个参数点一个晚上就能跑完。关键是把参数扫描逻辑写进循环而不是用COMSOL自带的参数扫描研究因为MATLAB循环能同时控制后处理算完一批Q因子顺手把透射谱、SHG效率一起提取了省一次打开GUI的等待。跑批量任务的时候建议先在局部区域用粗网格快速验证趋势再逐步加密。服务器上求解时把日志重定向到文件一旦失败了能看清楚是哪一步、哪个参数组合出的问题。5. 踩坑实录与排查速查表5.1 我在COMSOL里翻过车的五个瞬间第一个坑是特征频率虚部不收敛。粗网格下虚部忽大忽小看起来像噪音其实说明网格没有解析模式在柱体边缘的快速变化。解决方法是使用自适应网格或者在电场能量集中的位置手动加密。第二个坑是透射谱一直平没有尖峰。我排查后发现是背景场偏振方向设成了y而模式响应主要在x方向偏振没有激起共振。第三个坑是SHG效率小到离谱。检查半天发现非线性张量坐标没有从默认z轴旋转到x方向实际激发的根本不是d33分量而是d22那种偏小通道。第四个坑是2ω求解时内存爆了。倍频网格比基频加密之后自由度涨得非常快正确姿势是换直接求解器或者把研究拆开单独求解先算基频并保存再在2ω研究里只加载已知源。第五个坑是PML反射污染远场。PML离结构太近贴得太紧就会产生伪反射。让PML和结构之间至少留出半个自由空间波长的距离基本能解决。5.2 参数化扫描失败或内存炸了的应急思路如果参数扫描中途报错先看错误发生在几何重建、网格剖分还是求解器阶段。几何重建错误通常是参数组合让椭圆相交或变成零厚度需要在参数定义里做范围约束。网格剖分失败常见于把柱体尺寸缩得太小最小单元尺寸跟不上几何变化。求解器不收敛则要回到研究设置将频率扫描改成严格递增并把容差适当放宽一个量级。内存爆炸的应急办法是把“参数扫描”拆成多个独立研究任务每次只算一个点算完用完清除解或者去掉全波模场存储只保存边界上的S参数和积分量。还有一种被很多人忽略的方式把模型降维先用均匀LN平板验证SHG源项流程再切换回纳米柱阵列物理逻辑相同调试速度却快一个量级。我在群里见过有人用COMSOL做移动网格、激光焊接单元活化、沸腾两相流也见过用压电陶瓷做换能器的底层建模习惯和咱们完全一致先弄清楚物理量守恒关系再选模块最后才谈网格和求解器。光学超表面这块COMSOL的波动光学模块把频域、本征模、周期边界整合得很好不用像其他软件那样在网格和方程之间手动来回折腾。5.3 给零基础复现者的行动清单真想把这个课题复现出来按下面的顺序走会顺很多先建一个无图案的均匀LN平板输入Sellmeier色散跑一次基频和一次倍频验证SHG功率计算流程是对的。再在平板上加周期纳米柱线性条件下扫描频率找透射谷确认共振频率和你预想的目标一致。在本征模研究里扫Gamma点确认对称结构的BIC虚部趋近于零破缺对称之后Q按1/α平方下降。在共振频率处加非线性源算SHG先只取一个α值确认量级。最后再做α、周期、厚度联合参数扫描把SHG效率随α的非单调关系画出来。每完成一个阶段就存档一个独立版本别在一个mph文件里把所有东西全堆上。阶段检查比一口气跑完要省力得多。我个人在实际操作中的体会是这种仿真最危险的部分不是COMSOL操作而是你以为自己已经把物理搞懂了结果发现算出来的SHG效率来自一个不知名的杂散模式。所以跑完主流程之后我每次都会回到本征模场分布里再看一遍模式的电场极化确认它确实贡献的是想要的d33分量。这步确认过后面的参数扫描才有意义。如果你也准备复现这个课题建议从2.2节的坐标系问题开始解决先让基频共振场和光轴方向真正对准再做后面的一切。等你看到SHG效率随α呈现非单调变化、而Q又乖乖按1/α平方走的时候这个课题基本就通了。到那时再回头看那些BIC文献里的曲线你会觉得它们并不玄无非是几个参数在COMSOL里来回跑出来的物理故事。