做光学超表面仿真的朋友大概率绕不过非线性谐波这个坑。最近我把COMSOL里介质超表面三次谐波非线性模型完整跑通了一遍顺手把倍频模型和转换效率计算也一起实现了还特意处理了大家常问的功率依赖问题。这篇分享适合正在做非线性超表面、谐振增强光学或者谐波转换仿真的工程师和学生也适合刚入手COMSOL波动光学模块、想搞明白非线性项怎么加进去的初学者。1. 项目概览与建模思路1.1 这个模型到底要算到什么这个项目的最终诉求很明确在介质超表面结构上用COMSOL模拟三次谐波产生THG和二次谐波产生SHG也就是倍频并且要计算出转换效率。转换效率不能只是一个定值它得随入射泵浦功率变化也就是所谓功率依赖关系。很多文献里只给一个固定功率下的效率但实际实验里效率对泵浦强度非常敏感模拟时必须能复现这种趋势。展开来说模型需要回答三个问题第一基频光在超表面里的共振增强了多少局部电场有多强第二增强之后的电场通过介质材料的非线性极化率生成倍频和三次谐波源这个源在结构里以什么模式传播第三在给定入射功率下输出端能收集到多少谐波能量换算成效率是多少。这个小项目非常适合作为非线性超表面仿真的入门骨架。因为它的物理模型相对干净低泵浦、弱转换、不涉及太复杂的热效应和自由载流子效应主要就是围绕电场极化强度展开。一旦跑通后续再扩展和频、差频、太赫兹产生、甚至是瞬态非线性响应都只是在这个骨架上加内容。1.2 为什么是介质超表面而非金属早期非线性超表面大量使用等离激元结构因为金属纳米天线的局部场增强非常直观非线性响应也强。但金属有两个致命问题一是光热损耗太大结构很容易把能量变成热量损伤阈值低二是在近红外和可见光波段为增强而优化的金属结构往往带宽受限难以稳定输出高功率谐波。介质超表面用硅、氮化镓、氧化钛、砷化镓这类高折射率材料本身损耗低靠几何共振米氏共振、导模共振等把光束缚在亚波长体积内。虽然单点场增强比等离激元低但整体发热小可以承受高功率泵浦而且可以通过调节纳米柱的半径、高度和周期比较方便地控制共振波长的Q值。这个项目选择介质超表面还有一个重要的建模优势材料参数稳定非线性系数通常可以从晶体材料手册或者单晶薄膜文献里查到不像金属表面电荷和等离激元近场那么依赖界面状态。所以在COMSOL里给一个均匀介质柱体设定介电常数和非线性系数得到的仿真结果和实验通常会更接近调试也更容易。1.3 建模路线图整个建模过程可以拆成五步后面每一节我都会详细展开第一步设计超表面单元几何设置衬底、纳米柱周期和尺寸用周期性边界模拟无限阵列。第二步线性仿真基频共振先搞清楚场增强峰在哪个波长这是所有非线性响应的基础。第三步定义非线性极化强度的表达式把它添加到倍频和三次谐波接口的源项中。第四步分别求解倍频和三次谐波频率下的电磁场。第五步统计输入和输出功率计算转换效率扫描入射功率得到效率随泵浦强度变化的曲线。这五步是环环相扣的。我见过不少初学者一上来就想着把三个频率全耦合在一起自洽求解结果非线性源表达式一旦有误整个解都发散最后也不知道是哪一步出了问题。按上面的顺序逐级做每一级都可以单独验证出了问题也好定位。2. 先把非线性物理和参数准备讲清楚2.1 三次谐波、倍频与功率依赖的来龙去脉在强光电场作用下介质内部的电极化强度不再与电场线性递增而是展开成电场幂级数的形式P ε0 χ(1) E ε0 χ(2) E² ε0 χ(3) E³ ...其中 χ(1) 对应线性介电常数χ(2) 对应二次非线性效应倍频、和频、差频等χ(3) 对应三次非线性效应三次谐波、自相位调制、双光子吸收等。在频率为 ω 的基频光入射到材料时χ(2) 项会产生一个频率为 2ω 的振荡偶极子源向外辐射二次谐波χ(3) 项会产生频率为 3ω 的振荡源向外辐射三次谐波。这里有两个容易混淆的地方。第一中心对称材料比如硅、锗的 χ(2) 严格为零因为反演对称性要求二阶极化强度偶极子在空间反演下必须改变符号。所以纯介质硅超表面通常只做三次谐波很少做倍频。要做倍频需要选择氮化镓、砷化镓这类非中心对称材料或者利用硅表面/应变处的对称性破缺但那效率通常低很多。所以在动手建模前先想清楚你是要算倍频还是三次谐波需要哪一阶非线性系数不然材料设置会出大问题。第二功率依赖并不是一个额外的“物理机制”而是谐波辐射强度对基波场强的高次依赖关系。简单看二次谐波源正比于 E(ω)²所以等效输出功率 P(2ω) 正比于 P(ω)²三次谐波源正比于 E(ω)³所以 P(3ω) 正比于 P(ω)³。因此在效率曲线上双对数坐标下倍频效率的斜率约为 1三次谐波效率的斜率约为 2。功率依赖本质上就是这种幂次关系在超表面共振放大下的数值体现。但是这和“拿一个固定线性的基波场去算谐波”还是有区别的。真正的功率依赖要求基波场本身随入射功率线性增大而超表面内部的局域电场受几何共振影响并不是均匀变化。所以仿真时必须显式地把入射功率作为参数重新计算每一功率对应的基波场和谐波场再统计效率才能得到正确的依赖曲线。2.2 非线性系数怎么定、单位怎么换非线性系数的定义是建模时最容易被坑的地方。COMSOL默认采用SI单位制极化强度用C/m²电场用V/m。那么 χ(2) 的单位是 m/Vχ(3) 的单位是 m²/V²。材料手册里的 χ(2) 或者 χ(3) 数值往往会用别的单位系统比如esu单位或者给出的是“有效二阶非线性系数 deff”而不是 χ(2)。以常见的三次谐波材料硅Si为例近红外波段的 χ(3) 通常在 10⁻¹⁸ 到 10⁻¹⁹ m²/V² 量级但这个数值对晶向、泵浦偏振、波长和测量方法都很敏感。所以你最好从同一篇文章中找到与你共振波长接近的测量值不要随手抄一个。如果手头只有 esu 单位下的 χ(3)换算方法记住一句话1 esu 下的 χ(3) 乘以 4.19×10⁻⁸大概可以转成 m²/V² 单位的 χ(3)。具体换算不是简单乘以系数因为 esu 体系中电场单位是 statV/cm体积等定义也不同最稳妥的办法是找到文献里同样目标波段、用SI单位给出的值。实在找不到就用一个相对保守的量级去试算先看谐波能否收敛再调整系数乘以一个量级。倍频材料如果选砷化镓GaAs它的二阶非线性系数 deff 大约在 90 pm/V 左右。注意这里 deff 对应的是一个二阶张量的分量与 χ(2) 的关系是 χ(2) 2 deff。所以在COMSOL里写极化强度时需要乘系数 2否则倍频效率会算错一倍。2.3 在COMSOL中怎么把非线性源项写进方程COMSOL早期版本里非线性光学建模没有一个现成的“非线性极化”按钮需要自己在电磁波接口里加入源项。现在我用的做法是在“电磁波频域”接口下增加一个“外部电流密度”或者“安培定律”节点在那里直接定义谐波频率处的极化电流。以三次谐波为例谐波频率 3ω 处的非线性极化强度表达式为P_NL(3ω) ε0 χ(3) |E(ω)|² E(ω)需要注意的是这个表达式里的 E(ω) 是基频复数电场矢量|E(ω)|² 是它的模平方。如果先用“基波频率的计算结果”作为已知场那么应将基波场变量映射到谐波接口中。COMSOL里常用的做法是给基波和谐波各建一个电磁波接口谐波接口的源项写入基波接口的电场变量例如ew1.Epx * ew1.Epy 这种分量组合。如果直接用右手规则写出张量表达式还需要处理 χ(3) 的张量属性。很多简化模型里默认 χ(3) 是各向同性的标量于是三次极化强度沿某一方向的分量可以直接写成 ε0 * chi3 * (Epx * abs(E1)^2 ... )。但实际上晶体硅的 χ(3) 是四阶张量有多个独立分量。为了工程计算可以先假设一个等效标量并强调这是“极化方向与电场偏振方向对齐”的简化。等结果跑通之后再根据晶向投影填入更复杂的张量表达式。对于倍频非线性极化强度为P_NL(2ω) ε0 χ(2) E(ω) E(ω)同样需要在 2ω 接口中作为外部电流密度源项。如果只考虑单一偏振也可以简化为分量乘积。这里有个关键点直接加入外部电流密度源后谐波接口里的方程会把该源项当作驱动场因此谐波场幅度正比于源项大小。比如基波场如果在共振峰处局部放大 10 倍三次谐波源就放大 10³ 1000 倍输出谐波功率可能放大 10⁶ 量级。这就是为什么谐波对共振增强极其敏感。3. 几何、边界和求解配置实际操作3.1 超表面晶胞与周期性边界仿真对象是二维周期阵列一般只需建一个晶胞。介质超表面最常见的几何是一个矩形衬底层上面立一个圆柱或方形纳米柱。材料给定折射率周围背景要么是空气要么是某种低折射率介质。入射方向一般是垂直照射所以晶胞是规则的方格子边长通常是亚波长尺度比如 600 nm 到 1000 nm 之间。几何构建时我习惯把衬底厚度建得足够厚通常超过 500 nm避免在频域求解时把衬底本身的模式边界效应引入。衬底下表面要设置成端口或者散射边界纳米柱上方需要一层空气间隔再放一个端口接收透射光。如果想进一步消除背向反射可以在空气层上方和衬底下方加 PML完美匹配层。周期性边界在COMSOL的电磁波频域接口中叫 Floquet 周期条件也可以直接使用工程里的“周期性端口”。垂直入射即波矢平行于表面法线所以 Floquet 周期矢量的相位差为 0。如果后续要做斜入射扫描只需要在 Floquet 周期条件里输入各方向的相位模型不需要大改。尺寸设计上可以用一个简单原则纳米柱的高度大致对应目标共振波长的 1/3 到 1/2半径则决定米氏共振的位置。比如想增强基频光在 1300 nm 附近硅柱折射率约 3.4那么柱高大约 500 nm半径在 250-400 nm 之间周期约为 700-900 nm。这只是起点正式做之前必须跑一遍线性透射谱找到共振位置。3.2 多频率接口怎么搭我会在同一个COMSOL模型中添加两个或三个电磁波频域接口。第一个接口的频率设置为基频 ω第二个设置为倍频 2ω第三个设置为三次谐波 3ω。每个接口对应同一套几何但是物理方程域略有不同基波接口的场可以输入端口谐波接口没有入射端口只有非线性极化源项。具体操作是基波接口的边界条件用端口去定义入射平面波端口功率设为 P_in端口类型是“矩形”或“衍射光栅”模式是横电或横磁。谐波接口则不设置任何端口入射而是在“域”上添加“外部电流密度”节点并把电流密度表达式定义为非线性极化源的时谐电流J_NL(2ω) i * 2ω * P_NL(2ω)J_NL(3ω) i * 3ω * P_NL(3ω)这里 i 是虚数单位系数来自 J ∂P/∂t在频域下就是乘积 iω。很多人漏掉这个 i会导致谐波相位不对效率虽不受影响但后续如果要做相位匹配或相干叠加就会出问题。两个接口之间的变量传递很简单只要在谐波接口的表达式里直接引用基波接口的电场分量变量。COMSOL会自动识别不需要额外设置耦合节点。需要注意的是谐波接口里没有基波场所以引用的是 ew1假设基波接口名称是 ew1中的解。如果两个接口采用相同的网格变量映射是完全一致的。这样分开搭的好处是求解顺序可以控制先单独算基波再算谐波。基波结果和实验可以对比能确认共振是否模拟准。如果一上来就把基波和谐波全耦合在一起任何一个方程设置的错误都会淹没在大量未知数里。如果你非要做一个自洽模型那也不是不行但需要把非线性源同时加回基波接口表示高泵浦下的基波消耗。这时源项中既有基波场三次方又有谐波场耦合回基波的项方程会变成高度非线性收敛和内存都会很痛苦。弱非线性场景下没有太大必要先做分步法足矣。3.3 网格、求解器和内存管理网格是非线性超表面仿真的重中之重。基波频率下材料波长是 λ/n比如 1300 nm 入射到硅中波长约 380 nm。如果按照每个波长划分 5 个网格点介质内部网格单元可以放到 80 nm。但是谐波频率更高三次谐波波长是 1300/3 ≈ 433 nm在硅里约 127 nm网格就得加密到 25 nm 左右才能描述谐波场的空间振荡。所以三次谐波模型比倍频模型对内存敏感得多。我一般把纳米柱区域内设置“较细”网格并用“边界层”加密柱体表面因为场增强集中在这里。衬底和空气层则用稍粗的网格尤其是 PML 区域设置扫描网格即可。如果模型是三维总自由度可能达到数百万量级首次跑之前建议先用二维近似试算或者在柱体附近只做局部网格加密然后开“自适应频率扫描”或“直接线性求解器”。COMSOL的求解器选“直接求解器MUMPS”通常比迭代求解器稳对于电介质小模型速度也尚可。但是如果自由度超过几百万内存不够可以换用迭代求解器 GMRES 配合适当预条件或者降低网格密度。实践中计算转换效率对绝对场幅度精度要求不高最关心比值所以网格不用过于极致只要确保谐波场相位折射率分辨率足够即可。内存优化还有一个取巧办法只模拟半个周期结构即利用对称边界条件把模型减半。如果基波是平面波正入射且极化方向沿 x 方向那么在 y 方向结构又是对称的理论上可以沿对称平面切割加完美磁导体或对称条件。不过要注意非线性过程的对称性三次谐波的场分布在某些对称面可能是反对称的用对称边界条件时容易误删模式做之前先把基波线性模式在对称边界下能否正确解出来验证一遍。4. 转换效率计算和功率扫描4.1 谐波功率怎么统计转换效率是仿真最后要交出的数值。在COMSOL中要统计某一频率的透射功率可以在截面边界上做坡印廷向量的面积分。具体操作在“结果”菜单下添加“体/边界派生值”表达式填nx * W/m²? 其实是 0.5 * real(Epx * conj(Hpy) - Epy * conj(Hpx)) · n用界面上的外法向量 n所以实际计算是P_trans 0.5 * Re(∫(E × H*) · n dS)但COMSOL里不好直接写叉乘向量分量我习惯用内置变量ew2.Poavx、ew2.Poavy、ew2.Poavz时间平均坡印廷向量分量去积分这样最直接。截面可以取超表面下方衬底内部的一个横切面也可以取空气层上方出口。如果下方是衬底注意面积分时要把折射率不匹配造成的模式阻抗考虑进去但功率流本身是守恒的所以取哪个截面结果应该一致。对于基波接口同样统计输入功率或者透射功率。输入功率其实可以直接用端口功率COMSOL在端口节点里会显示入射功率、反射功率和传输功率。为了让效率定义更接近实验我一般用“透射侧总谐波功率 / 入射泵浦功率”作为效率物理含义就是有多少泵浦光被转换成了谐波并透射出去。需要留意的是谐波接口里没有入口源因此从端口边界出去的功率是谐波向上/下两个方向传播的总和。如果只关心顺向透射谐波就只对超表面下侧的横截面积分如果关心背向发射对上方截面积分即可。实验里一般用透射方向但超表面对称辐射两方向都有份额。4.2 转换效率定义与归一化倍频转换效率 η_SHG 的定义是η_SHG P_out(2ω) / P_in(ω)三次谐波转换效率 η_THG P_out(3ω) / P_in(ω)。有些文献里会用“光子转换效率”来算考虑的是光子通量而不是功率这样数值会比功率效率高因为谐波光子能量更大。在超表面非线性研究中通常默认功率效率。为了不被误导写报告时要说明清楚。很多仿真新手会觉得效率算出来总是一个很小的数比如10⁻⁵或10⁻⁶就怀疑是不是模型错了。实际上在没有特别复杂相位匹配的介质超表面单晶胞模型里这个量级非常正常。真正的实验能测得百分级别效率往往依赖的是准束缚连续态BIC这类超高Q共振或者特殊的相位匹配方案。普通米氏共振增强下三次谐波效率在10⁻⁵以下不用太焦虑。效率归一化还有一个细节超表面是周期阵列你仿真的是一个晶胞端口功率对应的是单个晶胞面积上的入射功率。如果你用实验里那个总光斑功率去和单个晶胞功率曲线比较差了一个“周期内晶胞数量”的倍数。其实效率是一个无量纲比值单个晶胞仿真和完整阵列在理想周期性阵列下的效率是相同的。所以曲线直接可用只是 P_in 值要理解成每个晶胞的平均输入功率。4.3 功率依赖仿真实操要得到功率依赖曲线方法很简单把基波接口的端口功率设置成一个参数 P_in然后做参数化扫描比如从 0.1 mW 扫到 100 mW。在每一个功率点重新求解基波和谐波并记录谐波功率和效率。为什么不能只算一个功率然后按幂次规律推导因为超表面的非线性源项虽然正比于场强的幂次但入射功率改变时电场幅度在整个频谱中分布几乎不变谐波功率确实会严格按幂次增长前提是材料非线性饱和和热效应不出现。在弱功率下双对数曲线基本就是直线斜率为对应阶数减一。既然线性规律简单为什么还要扫功率因为实际中随着功率增加可能出现共振热漂移、非线性折射率改变共振相位、高泵浦导致的基波损耗这会让曲线偏离理想幂次。功率依赖仿真就是要把这些偏离算出来。另外如果你做的是“包含功率依赖”的项目实验方会给一组不同功率的效率数据你的模型必须能直接复现这个趋势所以参数化扫描是必经步骤。具体在COMSOL里参数的 Scope 分两类全局参数和扫描参数。把 P_in 设成全局参数端口功率赋值给它。扫描时选择“参数扫描”步骤是先求解基波线性问题再求解2ω和3ω接口。如果你还是沿用“每个接口独立求解”的方式最好在每个功率点都重新计算全部接口也就是说让两个或三个研究步骤都落在扫描的范围内。否则如果只在最后一个步骤里扫描功率前面基波接口不会重新计算谐波接口的源项引用到的还是之前固定功率的结果曲线会完全错误。还有一个更省时间的做法先做一次基波共振波长扫描把每个波长的场增强曲线跑出来然后固定在最强的波长上做功率扫描。这样既能看清功率依赖又能确认共振位置。我强烈建议先把第一步“找共振”做扎实再上功率扫描否则一次全扫描会浪费大量时间。5. 问题排查与经验技巧实录5.1 常见报错与排查速查表下面是实际跑这个项目时最常碰到的几个问题和我的排查思路整理成一张表方便大家遇到报错时快速对照。现象可能原因排查与解决方向谐波接口处处为零非线性源项没写对或引用了错误的接口变量先检查源项表达式里是否用到基波接口的电场变量确认谐波接口是否勾选了“外部电流密度”临时把源项设成一个常数看谐波场是否建立谐波效率异常大或异常小非线性系数单位错误或使用了不对的 χ(3)/χ(2) 换算重新核对公式和文献数值先在简单平板中测试常见数量级再回到超表面求解不收敛基波场发散端口功率过大局部场超过材料击穿或网格分辨率不够降低功率一两个量级看看能否收敛加密纳米柱周围网格并尝试改为直接求解器谐波功率有负值坡印廷积分方向不对或者截面取在 PML 内部检查法向方向将截面移到物理介质中不要放在 PML 内部基波共振波长和实验差异很大结构尺寸、衬底折射率或周期边界相位设置不准确重新检查几何尺寸和 Floquet 周期条件尝试扫一个尺寸范围与实验对比内存不足三维模型网格太密减小模型到1/4对称降低外部区域网格等级或采用二维等效截面对比必要时用“有效折射率”法效率曲线不成预期幂次谐波接口源项使用固定场而不是重算基波检查扫描是否覆盖全部研究步骤把基波接口也放到同一个参数扫描研究中5.2 我踩过的坑和调参心得第一个坑是非线性系数的单位。我当时在做一个硅超表面三次谐波模型直接参考材料文章里写着 χ(3) ≈ 2.8×10⁻¹⁸我当时没注意那是 esu 单位直接填进COMSOL导致三次谐波能量大得离谱反射比基波还强。后来核对换算才发现差了大约四个量级。从那以后我每次建模前会先把资料里的非线性系数抄在一个专门的小本子上标注单位来源和适用波长而不是只复制数值。第二个坑是谐波接口的“外部电流密度”一定要放在求解域内不是边界上。有人会把源项放到端口边界上结果谐波在边界被当成入射激励重新反射效率虚高。正确的做法是放在纳米柱所在体积域中让它真实地在超表面内部激发。第三个坑是 FLOQUET 周期边界的相位设置。正入射时看似简单但如果模型旋转了或者坐标方向没对齐Floquet 周期矢量的相位可能不为零需要先用线性解检查透射谱是否和实验一致。我遇到过横空出世的共振峰最后发现是 Floquet 周期相位设错了把零阶模式都搞乱了。第四个坑是“功率依赖”与“自洽耗尽”的混淆。用户说“包含功率依赖”但很多场景下只是需要扫描功率不需要让基波场被谐波过程消耗。分步法里基波场在每个功率下重新求解但谐波场不会反向影响基波。如果实验效率很高比如超过10%基波耗尽不可忽略那你必须用耦合方程组做自洽。这时建议先跑通分步法再逐步打开反向耦合项否则一上来就自洽损坏的不只是心情还有硬盘。最后一个小技巧在参数扫描时把“谐波效率”的表达式直接定义为派生值表达式并勾选“在求解中计算”这样扫描结束后直接得到效率随功率的表格不用每个功率点手动点一次积分。导出曲线时建议用双对数坐标观察斜率如果三次谐波斜率不是接近2就说明源项或单位可能还有问题。这个检查方法比盯着绝对数更可靠。整个模型跑通之后我最大的体会是非线性超表面仿真的关键不在于界面上的参数多炫而在于把每一步的物理量纲和边界条件吃透。线性基波的共振、非线性源的定义、谐波功率的统计这三件事分开做每一步都验证基本不会有大偏差。如果你也是刚开始摸索COMSOL的超表面非线性建议先照着这个骨架跑一遍拿到一条功率依赖效率曲线后再往材料张量、斜入射、高功率自洽这些方向拓展。这样后面遇到新问题你至少知道自己是从哪个地基上加上去的。