最近和一个做非线性超表面的朋友聊天他说自己把孔径、周期、厚度轮换着扫了几轮二次谐波效率始终停留在 10⁻⁶ 量级上不去。我第一反应是问他你确定自己在追踪的模是连续域束缚态BIC准模吗如果只是拿透射谱里的某个谷当高 Q 谐振那很可能从一开始就追错对象了。这类光子晶体超表面模拟的真正门槛不在于会点几个菜单命令而在于搞清楚 BIC 的能量到底被锁在哪里以及 COMSOL 里哪些设置能把这个高 Q 模老老实实提取出来再用于二次谐波计算。这篇文章就把这条链路拆开讲为什么 BIC 值得用来做倍频增强、铌酸锂材料参数怎么进模型、COMSOL 中周期结构和边界条件的搭建要点、Q 因子的准确提取方式以及二次谐波模拟里最容易让人翻车的非线性源处理。内容主要面向正在搭这类模型的研究生和工程师默认你用过波动光学模块但没专门做过 BIC 课题。1. BIC 为什么能成为二次谐波效率的破局点1.1 高 Q 与非线性效率的定量关系先看一个最基本的事实二次谐波效率不是简单正比于泵浦功率而是受局域电场强度强烈调控。在小信号近似下转换效率大致可以写成η_SHG ∝ |χ⁽²⁾|² · |E_ω|⁴ · L_eff² · sinc²(ΔkL/2)这里的 E_ω 是基频处的局域电场振幅。对超表面这种亚波长结构有效作用长度 L_eff 很短相位匹配项又很难做到完美剩下的最有效杠杆就是 E_ω 的四次方。换句话说基频电场振幅只要增强三倍效率就能提升两个数量级如果增强一个数量级效率就是四个数量级的差距。谐振腔的作用正在于此。当结构支持一个高 Q 光学模式时能量会在腔内有很长的寿命局域场被反复积累最终电场振幅显著高于入射场。这就像荡秋千每次只推一小下但如果空气阻力和摩擦很小振幅也能越荡越大。Q 因子描述的就是这个能量耗散速度Q 越高能量存得越久场增强越可观。但注意效率对基频场是四次方关系对倍频模式的场是平方关系。所以如果只增强基频不增强倍频收益会打折扣如果两个模式都被谐振增强效果才最明显。这也是 BIC 类结构这两年特别受关注的原因——它能把两个波段的模式同时拉到很高的 Q 值而不需要复杂的腔体设计。1.2 连续域束缚态到底束缚了什么普通的光学谐振腔比如微环、法布里-珀罗腔都是靠反射镜或全内反射把光关在一个区域。这种束缚态容易理解光出去了就没法再回来但在两个界面之间还能反复弹射。而 BIC 从名字上就显得反直觉——它处在连续辐射谱之中频率和许多可辐射模式重叠理论和直觉都说它应该往外漏但它偏偏不辐射。对称保护型 BIC 是超表面中最常见的一类。想象一个正方晶格中的圆孔结构具有面内镜面对称性。在 Γ 点布里渊区中心附近某个布洛赫模的两个辐射通道因为对称性而产生了相位上的抵消远场辐射的净分量消失能量被彻底约束在结构内部。如果把圆孔改成椭圆孔或者圆孔稍微错位对称性被破坏辐射通道重新打开原本理想的 BIC 就退化成一个准 BICQ 值不再无穷大但仍然可以做到几千甚至上万。这里有一个很重要的直觉BIC 不是无损耗的同义词。材料吸收、加工粗糙度、衬底折射率失配都会引入额外损耗实际 Q 值会被这些非辐射因素限制。模拟中如果完全不设置材料损耗得到的 Q 值就是纯辐射 Q可以非常高但这种理想值只能用于理解物理不能直接和实验测到的数据比。2. 铌酸锂平台与超表面结构参数设计2.1 为什么是 LNOI 而不是别的材料做二次谐波超表面材料选择的本质是比三样东西二阶非线性系数、透明窗口、微纳加工成熟度。铌酸锂的二阶非线性系数 d33 一般在 27 pm/V 上下部分掺杂和特定波长下会更高一些这个数值在薄膜非线性光子学平台里属于第一梯队。它从可见到中红外都有不错的透过率1550 nm 泵浦到 775 nm 倍频这个常见组合完全落在透明窗口内两个波长都不会有明显吸收。另一个容易被忽略的点是折射率对比度。LNOI 薄膜折射率在 2.1 到 2.2 左右与空气和二氧化硅衬底都有足够差距。光子晶体超表面要做成高 Q 谐振必须靠足够强的折射率调制把模式束缚在结构层里折射率对比太小模式会很快泄漏到衬底或空气中很难形成 BIC。我在实际建模时的一个体会是不要在材料参数上省事。铌酸锂是各向异性晶体折射率不是一个数而是两个寻常光折射率 no 和非常光折射率 ne。z 切薄膜的 z 轴是光轴方向模型里需要把介电常数写成对角张量x 和 y 分量用 noz 分量用 ne两个波长分别对应不同的数值。想当然地用一个各向同性折射率替代基频模的偏振分布会算错后续非线性耦合强度的判断也就失真了。2.2 初始几何参数怎么定典型的 BIC 超表面结构很简单一块铌酸锂薄膜上周期性排列圆孔或方孔晶格是正方形。结构几何只需要三个参数晶格周期 a、圆孔直径 d、薄膜厚度 t。我的建议起点如下表后续再根据模式分析结果微调。参数建议初始范围说明晶格周期 a950–1100 nm保证基频布洛赫模落在 1550 nm 附近圆孔直径 d400–650 nmd/a 在 0.4–0.6 之间比较稳妥薄膜厚度 t300–500 nm太薄模式约束弱太厚会引入高阶垂直模式衬底折射率1.45SiO₂如果不做衬底模拟结果会和实际样品明显偏离这些数值不是拍脑袋。周期 a 决定了布洛赫模在 Γ 点的等效折射率位置d 和 t 一起调节模式的场分布和辐射泄漏最终要落到目标波长上。更严谨的做法是先不做圆孔只做一块平板薄膜的模式分析看无孔结构的导模在哪个波长相位匹配然后开孔微调把目标模式拉回 1550 nm。2.3 材料色散与损耗设定的取舍COMSOL 里定义铌酸锂材料时我建议直接把两个波长下的折射率作为参数分别填进去而不是让软件在频域自动插值。原因是超表面模拟往往只关心基频和倍频这两个离散频率点中间频率的值对结果没有实际影响自动插值反而可能引入来源不明的色散误差。损耗方面要分情况讨论。如果目标是研究 BIC 的辐射 Q 极限材料损耗可以设为零这样特征值分析得到的 Q 值反映的纯粹是辐射泄漏如果目标是和实验对比那需要在折射率虚部里加一点损耗。但注意铌酸锂在 1550 和 775 nm 的吸收很低体吸收对高 Q 模的影响往往小于散射损耗模拟中过度强调吸收意义不大。我个人的做法是先做无损耗模型摸清物理再在最后阶段加入一个小虚部观察 Q 值变化趋势。3. COMSOL 建模全过程几何、边界与激发设置3.1 物理场选择与研究类型的匹配COMSOL 中做 BIC 超表面模拟主战场是电磁波、频域接口。但同是这个接口不同研究步骤干的事完全不同很多人一开始就混淆了。要找到 BIC 模式并提取 Q 因子应该用特征值研究。特征值求解直接给出模式复频率从虚部可以算出辐射率进而得到 Q 值这比扫描入射波长再拟合透射谱效率高得多也直观得多。要得到二次谐波效率需要先做一个频域研究来获得基频稳态场然后再用一个频域研究求解倍频场。这两个研究可以在同一个组件里用两个研究步骤串在一个文件里也可以拆成两个文件。习惯上我更倾向拆成两个文件一个专门做基频模式和 Q 因子分析另一个做非线性源加载和倍频求解。原因后面会细说主要是非线性源的表达式比较复杂混在一个模型里容易把研究步骤的耦合搞乱。3.2 几何建模一个晶胞就够了超表面是周期结构COMSOL 里不需要建一大片只建一个包含圆孔的晶胞周围用周期性边界条件模拟无限阵列即可。几何上我习惯是一个长方体代表铌酸锂薄膜在顶面画一个圆向下拉伸并做布尔减法扣除圆孔然后在上方和下方分别加一段空气或二氧化硅背景域。这里有几个容易出问题的细节。第一布尔运算后要确保所有域做一次并集避免相邻域之间存在虚假的内部界面否则网格划分时会在这些界面上产生多余的单元影响高 Q 模的收敛。第二圆孔壁上的网格质量直接决定 Q 值的可靠性圆孔边缘最好在几何建模时就把倒角或切向边界留出来让网格生成器能在壁面附近自然加密。第三背景域的高度要足够通常上方至少留一个工作波长的长度给端口和可能出现的衍射场展开。如果需要完整模拟衬底的影响衬底域也要加进去但衬底底部的边界条件要处理好。实际样品中衬底是半无限的模拟中要么用足够厚的衬底加吸收边界要么用周期性端口把衬底侧的辐射通道也考虑进去。只放一层薄衬底然后什么都不接会让能量错误地堆积在模型边界上。3.3 周期性边界与端口激发和泄漏通道要分清楚晶胞的四个侧面设置 Floquet 周期性边界条件x 方向两对面和 y 方向两对面分别成对绑定并给出 Bloch 波矢 kx、ky。在 Γ 点kx ky 0如果要看带结构可以扫描 kx。上下两个面根据用途不同有两种处理方式。如果做特征值分析只需要在外围加一层完美匹配层吸收辐射模式本身不需要外部激发算出来的是结构本征模。如果做透射谱或倍频效率就需要端口边界。上端口的周期性端口设置为入射平面波功率归一化为 1 W下端口的周期性端口作为透射探测如果需要反射谱上端口也要同时探测反射。这里提醒一个关键点端口必须设置为周期性端口类型而不是普通的均匀端口。普通端口假设场在横截面内是常数这对超表面的高阶衍射模是错误的描述。周期性端口会依据 Floquet 波矢自动判断各衍射级次这样散射到各阶衍射通道的能量才能被正确记录。3.4 网格划分的底线网格是超表面模拟里最没有技术含量但又最容易翻车的环节。我的经验规则是铌酸锂薄膜内部最大网格尺寸不超过 λ/(8n)空气区域不超过 λ/8圆孔壁面至少加三层边界层网格。这里 λ 取基频波长n 是铌酸锂折射率。如果你发现一个模式算出来的 Q 值超过预期一个量级以上不要高兴太早先做一次网格减半对比。特征值随网格加密会逐渐收敛如果加密后 Q 值显著变化说明之前的结果是伪收敛。反过来如果透射谱中的谐振峰位置对网格不敏感但 Q 值敏感问题也出在网格。BIC 模的场主要集中在孔洞边缘和薄膜内部能量泄漏通道又很细网格不够密时会把辐射泄漏算得偏大或偏小Q 值自然不可信。4. 从特征值到 Q 因子BIC 的定位与对称性扰动4.1 特征值模式分析和透射谱各有各的用处我习惯先用特征值研究找到目标模。具体操作上在特征值求解器中设置一个搜索区域让求解器只在 180–220 THz 范围内找模式对应 1360–1660 nm然后观察每个模式的复频率和电场分布。特征值的结果是一组复数值。假设复角频率为 ω ω_R jω_I那么 Q 因子定义为 Q |ω_R| / (2|ω_I|)。但 COMSOL 输出的特征值单位可能是角频率也可能是频率符号约定也因版本而异千万不要只看虚部绝对值就套公式一定要先用一个低 Q 的简单模式做验算确认提取方式正确后再批量处理。实际定位 BIC 时单纯看 Q 值还不够。对称保护 BIC 在 Γ 点的特征值虚部理论上趋于零但在数值实现中会得到一个有限的但异常大的 Q 值比如 10⁶ 甚至 10⁸。除了看 Q还要看电场分布BIC 模的场集中在圆孔边缘或薄膜内部而不是均匀分布在空气域。如果你追踪的模式主要落在空气背景里那大概率是某个数值伪模而不是物理 BIC。透射谱的用途则是在你确认模式的物理机制之后。对结构加一个斜坡入射端口激发扫描波长找到透射谷或反射峰拟合洛伦兹线形可以得到 Q 值。这个 Q 值和特征值提取的 Q 应该相互印证。我个人的流程是特征值分析粗筛候选模 → 电场分布做物理确认 → 透射谱验证 → 微调参数。4.2 对称性破缺从一个不可激发的束缚态到一个可用的准 BIC理想 BIC 的最大问题是它不辐射也不被外部入射波激发属于好得没法用。要让 BIC 变成实际可以被平面波激发、又能把能量耦合出去的高 Q 谐振必须引入微弱扰动打破对称性。最常用的做法是把圆孔改成椭圆孔。假设椭圆短轴和长轴原本相等即圆形扰动参数 δ 0现在设定长轴和短轴的差值为 δ那么当 δ 非零时面内对称性被破坏模式获得一个辐射泄漏通道透射谱中出现一个谐振峰。这个峰的 Q 值随 δ 增大而迅速下降在小扰动区间近似满足 Q ∝ 1/δ²。所以设计上你要做一个权衡δ 太小Q 值极高但透射耦合太弱谐波输出也弱δ 太大Q 值降下来场增强不够倍频效率上不去。实际设计中存在一个最优 δ这只能通过参数扫描找到。建议 δ 从 0.02 扫到 0.15步长 0.01 左右把每个 δ 对应的 Q 值和透射峰深度都记下来在 log-log 坐标里画 Q-δ 曲线斜率接近 -2 的地方就是准 BIC 行为区。4.3 参数扫描与模式追踪的实用脚本思路做参数扫描时有一个麻烦随着几何参数变化模式频率也在移动追踪同一个物理模式需要额外技巧。我通常的做法是固定特征值搜索范围提取每个参数点下的所有模式然后依据电场分布的重叠积分来判断哪个模式是之前追踪的那个。这个思路在参数变化后尤其重要。孔径变大、厚度改变可能某个模式的顺序发生交换如果你只按频率排序去对号入座很容易追到另一个模式上。代码层面可以用 COMSOL 的 Java API 或 LiveLink for MATLAB 批量提取特征值然后按场分布相似度自动匹配模式。手动做虽然可行但扫上几十组参数后肯定累到怀疑人生。把模式追踪逻辑自动化省下的时间足够多做几轮网格收敛测试。5. 二次谐波模拟链路从基频场到倍频源5.1 非线性极化的物理表达与偏振匹配问题非线性超表面的倍频模拟本质上是求解一个包含非线性极化源的电磁问题。二阶非线性极化可以写成P_i(2ω) ε₀ · χ_ijk⁽²⁾(2ω;ω,ω) · E_j(ω) · E_k(ω)对铌酸锂来说最大的非线性系数是 d33即当参与作用的所有场分量都沿光轴方向z 方向时耦合最强。这意味着如果你的基频模和倍频模的主要偏振分量不在 z 方向那么即使 Q 值算得再高非线性效率也会被 d31、d22 这些小系数的数值限制住。所以基频模式设计要让电场偏振尽可能与晶体光轴产生有效重叠。对 z 切铌酸锂薄膜这意味着我们需要模的电场在 z 方向有显著分量。这个要求在波导模式里往往对应 TM 模但在超表面这种强约束结构中模式的偏振成分是混合的不能简单说是 TE 或 TM必须从模式场图里实际判断。我见过不少模拟报告把基频模算出来之后完全不检查偏振分布就直接算倍频最后效率上不去还不知道卡在哪。5.2 两步法而不是一次非线性求解直接在 COMSOL 里做全非线性频域求解把 χ⁽²⁾ 放进材料本构关系然后一次解出基频和倍频听起来很优雅实际做起来非常痛苦。高 Q 模式本身的数值条件数就不好再引入非线性耦合求解器的收敛性和内存开销都会急剧恶化。对绝大多数 BIC 超表面倍频问题两步微扰法足够准确步骤也更可控。第一步在基频处做无源频域求解得到稳态电场 E_ω。把几何参数、边界条件和基频波长固定好端口注入归一化泵浦功率。这一步和普通透射谱模拟没有区别。第二步建立第二个频域研究频率设为 2ω。在这个研究中把非线性极化 P(2ω) 作为外加源项加入波动方程。COMSOL 的具体实现方式是把 P 写成变量表达式通过外部电流密度或直接修改方程的源项施加到求解域中。第二步的难点在于把基频场的 E_x、E_y、E_z 表达式正确映射到倍频域的变量里。我的做法是先用全局定义里的变量公式写出 P 的三个分量调用步骤一的结果场变量然后把这个源项赋给倍频域。注意这里有个小坑COMSOL 中不同研究的场变量名不同直接在同一模型里跨研究引用时要确认变量引用指向的是研究 1而不是当前研究否则倍频求解时 COMSOL 会重新计算基频场造成循环依赖。5.3 倍频效率的计算口径得到倍频场之后效率计算方式决定你的结果能不能和实验对比。效率公式很简单倍频输出功率除以泵浦输入功率。泵浦输入功率在端口定义里归一化为 1 W所以关键是准确提取倍频输出功率。倍频输出功率的积分面要放在合理的位置。我的建议是放在倍频端口面之前、PML 或吸收边界之前同时统计向衬底侧和向空气侧两个方向的辐射功率。如果在整个模拟域内做体积分数值噪声会被放大尤其是高 Q 模产生的强局域场会把积分结果污染得面目全非。计算得到的效率通常是 10⁻⁷ 到 10⁻⁵ 量级别嫌低这是亚波长结构加非理想相位匹配的正常水平。要提升效率除了提高 Q 值还要让基频模和倍频模在空间上有最大重叠。可以计算一个模式重叠因子η_overlap |∫ E_ω² · E_2ω* dV|² / (∫|E_ω|⁴ dV · ∫|E_2ω|² dV)这个式子虽然不能直接预测绝对效率但能告诉你结构参数改变时两个模式的耦合趋势是好是坏。5.4 相位匹配在超表面语境下的现实含义传统非线性晶体中相位匹配条件 Δk 0 非常重要因为作用长度长。超表面属于亚波长器件有效作用长度只有几百纳米到一两个微米相位失配对效率的惩罚没有体块晶体那么致命但也不是可以完全忽略。在 BIC 超表面里基频模和倍频模通常都是离散的谐振模式它们的频率由结构参数决定。要让倍频高效产生需要倍频模式的频率恰好落在二次谐波频率附近这本质上是一个模式匹配而不是波矢匹配问题。实际操作中我会在几何参数扫描时同时追踪基频模波长和倍频模波长当二者的比值接近 2:1 时倍频效率往往会出现峰值。这个设计思路意味着你的参数扫描不能只看基频的 Q 值要同时记录倍频模的谐振位置。我一般会建一个二维参数表横轴是扰动参数 δ纵轴是晶格周期 a颜色标注是倍频效率效率峰值区域往往对应基频和倍频模式同时满足谐振条件的一小片参数区间。找到这个区间后再细化扫描。6. 高 Q 模拟中的数值稳定性与常见坑位6.1 特征值求解的搜索策略高 Q 模式的特征值虚部极小如果直接用默认求法很容易漏模或者在大量低 Q 泄漏模中把目标模淹没。我习惯把求解器设置为在目标频率附近搜索指定数量的模式比如在 190–210 THz 范围内搜索 20 个模式然后按虚部绝对值升序排列虚部最小的那个通常就是最接近 BIC 的模式。COMSOL 特征值求解器有一个容易让人踩坑的地方默认可能会输出大量空间振荡的伪模。这些伪模往往分布在网格的高频极限区域电场分布呈棋盘状虚部也不小。判断方法很简单看电场分布是否光滑如果相邻网格点电场剧烈跳变那几乎可以肯定不是物理模式。6.2 PML 和端口边界的互相干扰很多人在高 Q 超表面模拟里把 PML 直接加在上下两端同时又在端口面上设置入射波。这两个东西的功能有重叠配置不对时会出现额外反射污染谐振峰。我的建议是透射谱计算中上下端面用周期性端口PML 放在端口后方的延伸域末端特征值分析中才用 PML 直接包裹结构上下方并延伸到周期边界的上下范围。PML 的厚度至少取一个基频波长折射率渐变曲线用默认的多项式即可。我发现很多人把 PML 厚度设成 λ/4这在低 Q 问题里可能够用但高 Q BIC 模式的远场拖尾很长PML 太薄会产生数值反射反过头来抬高 Q 值给你一个虚假的高 Q。6.3 色散数据单位与坐标轴方向导入铌酸锂的折射率色散数据时单位是一个经典坑。COMSOL 中材料属性对话框把波长单位默认设为微米而多数 Sellmeier 公式用的是微米但也有文献用纳米或森特米。我习惯在导入前统一把波长换算成微米并把折射率数据列单独存成 CSV。如果发现算出来的透射谱谐振峰整体偏移了几十纳米先检查色散数据单位而不是急着改几何参数。光轴方向同样要反复确认。z 切铌酸锂的介电常量张量是 diag(no², no², ne²)但如果你建模时把薄膜法向选成了 x 方向张量就要对应旋转。这个错误在模式电场分布里容易暴露如果算出来的 TE 模和 TM 模频率差得离谱多半是材料方向设置不对。6.4 与实验对比时的归一化口径最后提醒一个从模拟走向实验时最容易产生误解的问题。模拟中的效率是单胞端口效率对应一个无限周期超表面在平面波照射下的响应实验中的效率则要考虑聚焦光斑尺寸、入射角展宽、加工误差造成的 Q 值退化。所以模拟效率普遍比实验高一到两个量级这不一定是模拟错了而可能是理想周期结构和实际样品之间本来就存在差距。如果要把模拟和实验放在一起比较至少要把 1550 nm 波长的光斑面积与超表面有效面积的比值做一个修正再把实验测到的 Q 值代入模拟的损耗模型里而不是拿无损耗模拟的 Q 值直接和实验谱线对比。我自己做这类 BIC 超表面模拟最大的体会是整个项目里COMSOL 操作只占很小一部分大部分时间都花在确认自己算的确实是物理上那个模式这一件事上。特征值分析、透射谱、场分布、参数追踪每一步都在交叉验证同一个目标。先把 Q 因子做实再谈倍频效率这条路径比一开始就直接往全非线性求解器里灌参数要稳妥得多。最后再分享一个小技巧所有几何参数的扫描结果都记录在同一套命名规则下比如 d050_t400_delta003 这种后期回溯数据时你会感谢当初的自己。