做非线性光学的仿真二次谐波产生SHG是绕不开的第一道坎也是我觉得最适合拿来上手 COMSOL 的模型之一。它物理图像清晰、方程形式简单但真在软件里搭起来又会冒出不少反直觉的坑。我自己从二维波导倍频做到块状晶体聚焦倍频再往后加了周期极化结构前后掉过不少头发。这篇就专门讲 COMSOL 仿真二次谐波产生这件事从非线性极化的物理公式到如何在 COMSOL 里把倍频源项写进方程再到怎么把转换效率算得和理论对得上。目标读者是两类人一类是刚开始用 COMSOL 做非线性光学的同学想找一个能跑通又符合物理的参考模型另一类是已经跑通线性仿真、想往非线性方向扩展的工程师。文章不会只给步骤还会解释每个设置背后的原因读完你应该能照着搭出一个自己的倍频模型而不是只会抄参数。1. 二次谐波产生的物理基础与仿真映射1.1 SHG的物理本质二阶极化率如何成为倍频的源二次谐波产生的本质是介质在强光作用下出现了非线性极化。基频电场E(ω)通过二阶极化率χ(2)在介质内激起极化强度P(2ω) ε0 * χ(2) * E(ω) * E(ω)这个极化强度可以理解成一个频率为2ω的振荡偶极子源它不断向外辐射二次谐波。这里有个容易搞混的点基频光并不是“变成”了倍频光而是非线性极化把能量耦合到了2ω的电磁场模式上。所以在 COMSOL 里建模最常见也最不容易出错的做法是把P(2ω)直接当作一个等效电流源加到倍频方程里而不是去修改材料的介电常数或折射率。往后看转换效率主要看三件事泵浦光强、有效非线性系数、以及作用长度内的相位失配。相位失配Δk k(2ω) - 2k(ω)当Δk * L远大于π时倍频光一边产生一边又往回转换最后输出的功率会周期振荡甚至趋近于零。真实器件里用双折射相位匹配或准相位匹配来补偿Δk这也是后面仿真里必须手动控制的关键变量。1.2 为什么用频域接口而不是瞬态求解不少教程一上来就建议用“波动光学电磁波频域”接口但很少说清楚为什么可以用频域。关键在于连续光SHG是一个单频稳态过程泵浦是单频的输出也是固定频率的2ω整个系统里只有这两个频率分量完全可以用频域求解。瞬态仿真当然也能做倍频但它要处理脉冲包络的时间演变、色散以及数值稳定性模型规模和计算量都会大很多。对于入门和器件设计先不用考虑瞬态。频域建模的思路其实很简单把ω和2ω这两个频率看成两个独立的物理场接口分别求解Maxwell方程再用一个耦合项把基频场的平方作为源引到倍频接口上。一旦建立了这个“双频域接口非线性源项”的框架后面所有设置都顺理成章。1.3 单向耦合还是双向耦合仿真策略的取舍严格来讲SHG是基频和倍频互相耦合的过程倍频增强之后会反过来消耗基频功率这就是泵浦耗尽。但在大多数工程场景里转换效率只有百分之几甚至更低泵浦耗尽可以忽略。这时可以用单向耦合先求基频场再把基频场作为已知源去算倍频场。单向耦合的好处是数值稳定、调试简单结果也足够准。COMSOL里做单向耦合最简单的方式是定义两个电磁波频域接口。一个在频率f1负责基频另一个在频率f2 2*f1负责倍频。基频接口的源是端口边界条件指定的输入功率倍频接口的源是手动添加的体电流密度。电流源表达式里直接引用基频接口的电场分量比如写成ewfd1.Ez。如果转换效率超过10%甚至更高单向耦合会明显高估输出功率。这时才需要考虑双向耦合在基频方程里也加入与倍频场相关的反向耦合项两个接口一起迭代。这种做法的收敛性会差很多我一般只在做高Q微环谐振增强时才会硬着头皮去调。2. 动手建模前先把这几个关键决定做对2.1 几何选型块状晶体、波导还是微环很多人第一反应是“我要仿谐振腔”但第一个模型最好别碰腔体。块状晶体倍频适合模拟入射光强低、作用距离短的实验波导结构因为有模式束缚有效模面积很小单位长度转换效率比块状晶体高好几个量级微环或谐振腔则靠循环功率增强效率可以非常高但设计参数多调试麻烦。我建议从二维波导开始。二维模型几何简单、网格量小几分钟就能算完还能清楚看到基频模场与倍频模场的重叠情况。等二维结果和理论公式对上了再往三维扩展或者改微环都不迟。实际器件里波导截面的宽度、脊高、衬底折射率差都会影响模式分布进而影响有效模面积和重叠积分。拿铌酸锂脊波导来说宽度和脊高的变化会让模场从泄漏到强束缚之间切换效率差距可以到几个数量级。所以几何选择本质上是在选“模式约束强度”。2.2 材料参数色散、非线性系数与单位换算COMSOL内置材料库确实有不少光学材料但它的折射率多数是常温下的静态数据色散拟合精度对非线性仿真不够。更好的办法是查对应晶体的Sellmeier方程自己把折射率随波长的变化写进模型。以铌酸锂为例常温下ne在1064nm大约2.233在532nm大约2.320这0.087的折射率差会让相位失配严重到几乎无法产生倍频所以必须靠温度调谐或周期极化来补偿。非线性系数χ(2)需要自己换算。文献里更多给的是d系数单位是pm/V二者关系是χ(2) 2d比如铌酸锂z向极化的d33大约是25pm/V那么χ(2)就是5×10^-11 m/V。这个单位换算我见过太多人栽跟头数值差一个数量级效率就差两个数量级最后怎么都调不回来。另外要注意张量方向。不同晶向允许的非线性系数通道不同不是E^2随便乘一个标量就行。对于z切铌酸锂如果基频和倍频电场主分量都是z方向走的是d33通道如果是x或y偏振用的可能是d31或其他分量效率会明显低一截。2.3 散射边界和网格的取舍模型边界如果不处理电磁波会在边界反射反射波和入射波混在一起最后算出来的效率全乱。对于一进一出的波导输入截面要用端口边界条件给定模式功率输出端和侧壁要么用散射边界条件SBC要么用完美匹配层PML。SBC对垂直入射吸收比较好但对斜入射会有残余反射PML更稳但会增加网格数量。网格尺寸的基准是每个波长方向至少5到6个网格。注意这里说的波长是介质中的波长不是真空波长。对折射率n的材料实际波长是λ/n所以最大网格尺寸要取到λ/n / 5边界处最好再额外加边界层网格去刻画倏逝场。网格太粗会引入数值色散等于把相位失配量给偷偷改变了。你明明材料参数都设对了最后倍频效率还是不对很可能就是网格搞的鬼。3. 实操全过程1064nm到532nm的波导倍频模型3.1 参数设定与几何构建下面给一个可以直接照抄的参数起点。模型是x方向传播的脊波导z切铌酸锂宽度1.2μm脊高0.6μm衬底厚度取2μm包层是空气。芯片长度先从1mm开始验证没问题再拉长。参数数值说明泵浦波长1064 nm基频波长倍频波长532 nm二次谐波波长波导宽度1.2 μm脊宽脊高0.6 μm脊刻蚀深度衬底厚度2 μm下层铌酸锂泵浦功率100 mW端口输入功率波导长度1 mm 起可后续拉长极化周期约6.1 μm准相位匹配周期几何构建没什么特殊技巧画一个长方形覆盖脊波导截面再画一个衬底矩形材料分别赋给铌酸锂和空气即可。注意COMSOL几何里如果把整个晶体都当作铌酸锂来画基频模式可以跨整个基底传播这和真实刻蚀波导不太一致。实际模型要在脊和衬底之间区分开才能得到接近真实的模式约束。材料部分铌酸锂折射率用Sellmeier公式拟合ne和no分开给。非线性系数不需要在材料节点里设置直接写在后续的体电流源表达式中即可。关于周期极化我的做法不是真的在几何里画一排正负畴区那样网格会爆炸。更高效的做法是在体电流源里加一个符号调制项χ_eff(x) χ(2) * sign(cos(2π x / Λ))这样就能在均匀材料里近似模拟周期极化而不用网格去解析每一个畴壁。3.2 物理场设置端口、体电流源和变量映射两个电磁波频域接口的设置分别是基频接口假设叫ewfd1频率设为f1 c_const / lambda1。输入截面用“端口”边界模式设为“TE”或“TM”取决于你要激励的模式。端口输入功率设为100mW。求解后要能输出基频电场分量方便倍频接口引用。倍频接口假设叫ewfd2频率设为f2 2*f1。无需端口输入但要添加体电流密度源。电流源表达式沿z方向写成J_z i * 4 * pi * f2 * eps0_const * chi_eff * ewfd1.Ez^2这里有两个关键点稍微不注意就会全军覆没。第一ewfd1.Ez是频域复振幅直接平方后时域上对应的正好是2ω分量。如果手滑写成abs(ewfd1.Ez)^2那就等于只取了模的平方物理上对应的是直流分量而不是倍频分量输出会完全错误。这个坑我踩过不止一次每次换模型重新写表达式时都要提醒自己。第二电流源里带i是频域方程的习惯写法。因为时域里P(2ω)对时间求两次导对应频域中要乘上i 2ω再乘i 2ω严格说倍频极化源在无泵浦耗尽近似下等效电流密度是J ∂P_nl / ∂t频域中就是J i * 2ω * P_nl。所以J_z i * 2 * ω2 * ε0 * χ_eff * E_z^2而ω2 2π f22 * ω2 4π f2最终表达式就是我上面写的那条。在COMSOL变量里用ewfd1.Ez引用另一个接口的场前提是两个接口求解顺序正确。我一般用两步研究第一步只算ewfd1第二步再算ewfd2这样倍频求解时基频场已经是已知量不会形成隐式耦合。如果强行把两个接口放在同一个频域研究里一起解一般会得到完全不收敛或者物理上不正确的解。另外倍频模式也有可能泄漏到衬底或包层网格必须能解析倍频电场。如果只按基频波长设计网格倍频波长短一半网格很可能不够。最大网格尺寸设置时直接按倍频介质波长来设最保险。3.3 研究顺序与后处理把转换效率算出来研究设置按顺序做两步第一步只求解ewfd1得到基频场。第二步求解ewfd2利用上一步结果作为源得到倍频场。倍频输出功率在COMSOL后处理里可以用输出截面上的坡印廷矢量积分来算。定义一个积分算子intop施加在输出边界上然后在全局变量里写P_out intop(ewfd2.nPoav)这里ewfd2.nPoav是倍频接口的平均功率流法向分量。如果你用的版本变量名不完全一样可以在结果里手动插入“表面积分”表达式选ewfd2.nPoav即可。转换效率定义为η P_out / P_pump在低转换效率极限下理论上倍频功率正比于泵浦功率的平方、作用长度平方、有效非线性系数平方反比于有效模面积。所以换算成归一化效率时常用单位是%/W/cm²。仿真结果如果和这个标度关系对不上比如扫泵浦功率时η不是线性增长基本可以肯定模型里有问题。我习惯的计算流程是这样的先固定泵浦功率100mW扫波导长度看倍频输出是否随L²增长然后再固定长度扫泵浦功率看η是否线性变化。两个扫描都符合预期才认为这个模型是可信的。不要第一步就直接拿一个长波导去算最高转换效率一旦结果不对你很难判断是物理问题还是数值问题。4. 常见问题与排查技巧实录4.1 倍频效率几乎为零先手动算相位失配如果你照上面的步骤搭完模型发现倍频输出低到小数点后一大堆零第一反应别去怀疑COMSOL先手动算一次相位失配。用前面的铌酸锂例子1064nm处ne1 2.233532nm处ne2 2.320波导长度L 1mm。相位失配近似为Δk 4π * (n2 - n1) / λ1代入数值Δk ≈ 4π * 0.087 / 1.064e-6 ≈ 1.03e6 rad/m再看Δk * L ≈ 1030远远大于π。这说明倍频光在第一个振荡周期都没走完输出几乎归零纯粹是相位匹配没做好不是软件问题。解决办法是准相位匹配。周期取Λ 2π / Δk λ1 / (2 * (n2 - n1)) ≈ 6.1μm在体电流源里加上sign(cos(2π x / Λ))调制后等效每走一个周期就给相位失配补一次π倍频输出会指数级上升。做这一步时你会非常直观地体会到“相位匹配才是SHG的灵魂”。4.2 不收敛、场图振荡、报错按顺序排查倍频仿真最常见的数值问题有三个场图出现棋盘状振荡、求解器报奇异矩阵、或者收敛迭代慢到无法接受。遇到这些问题排查顺序应该是网格太粗把最大单元尺寸减半看倍频输出是否明显变化。如果变化很大继续加密直到结果基本稳定。源项太尖锐sign(cos(2πx/Λ))这种周期符号函数在畴壁处不连续局部电场梯度极大。可以用tanh(x/δ)近似平滑符号函数δ取波长的1/20左右既保留物理效果又避免数值振荡。求解器不合适倍频接口里电场突变和强源项会让迭代求解器很难收敛换用直接求解器比如PARDISO通常会好很多。边界反射如果输出端和侧壁用的是SBC而波的入射角很斜残留反射会影响倍频场。换成PML或者把边界拉远。这套排查顺序基本覆盖了我遇到过的绝大多数情况。先把网格和求解器搞定再去抠物理细节效率会高很多。4.3 结果和文献对不上从这三个方向找原因仿真结果和文献实验值对不上时先别急着怪文献也先别怪自己参数填错。我通常按下面几个方向排查第一材料色散模型是否一致。不同Sellmeier公式拟合出来的折射率在可见光和中红外波段可能差出0.001到0.01对相干长度的影响非常大。你用的色散和文献用的不一定是同一套相位匹配位置就会偏移。第二模式重叠是否一致。二维截面近似和三维真实波导的模场分布不同尤其近截止区。倍频效率正比于基频与倍频模式的重叠积分如果重叠很小效率就低。可以算一下归一化重叠因子和文献对比。第三泵浦功率的定义是否一致。COMSOL端口功率是输入总功率但非线性源用的是局部电场强度。如果波导模场面积很大单位面积光强很低即使总功率有100mW倍频效率也可能低得可怜。反过来如果用了归一化模场却没有正确缩放输出的绝对功率就毫无意义。下面这个速查表是我自己调试时一直用的症状优先检查倍频输出极低相位失配、Λ是否设置正确输出随长度振荡网格不够、边界反射输出随功率非线性增长模式场缩放、归一化错误求解不收敛网格、求解器、符号函数太尖锐5. 进阶方向从静态模型走向真实器件5.1 周期极化、温度调谐和微环增强如果你跑通了均匀波导倍频下一步强烈建议做周期极化扫描。把Λ设成扫描参数从5μm扫到7μm你会看到一个明显的倍频效率峰值。这个峰值位置和Δk计算值吻合时说明你的色散和模式都设对了。再进一步可以把温度也加进变量。铌酸锂的折射率随温度变化明显不同温度下最佳极化周期也不同。这样仿真就能指导真实实验的温度控制。微环增强的倍频则更有挑战性。微环里循环功率比输入功率大很多所以同样的非线性系数可以产生高得多的倍频效率。但要注意微环的谐振频率随温度飘移严重仿真里必须把热效应和相位匹配一起考虑否则结果和实验很难对上。5.2 模式重叠积分与波导截面的优化随着模型跑通你会越来越关注模式重叠。倍频产生的模场分布和基频模场不一定完全匹配。如果基频和倍频模场在截面上的空间重叠很小即使相位匹配完美效率也高不起来。可以在后处理里定义一个重叠积分overlap ∫ E_z(ω) * E_z(2ω) dS / sqrt(∫|E_z(ω)|² dS * ∫|E_z(2ω)|² dS)然后用参数扫描分别改变波导宽度、脊高、刻蚀深度观察重叠因子和有效模面积的变化。这个优化过程很考验参数设计能力但也是做真实器件设计最有价值的部分。5.3 脉冲与瞬态什么时候才需要上瞬态连续光倍频做熟了以后你会碰到脉冲泵浦的问题。这时频域接口已经不够用因为脉冲包含多个频率分量而且时间包络会随位置演化。COMSOL里做瞬态倍频需要用含时间项的波动方程并考虑材料色散在时间域的等价形式。瞬态模型的计算量远大于频域网格和时间步长都要严格验证通常只用来研究超短脉冲倍频的波形变换或者腔体内光场的瞬态建立过程。如果你只是要稳态转换效率踏踏实实用频域模型就好别给自己找麻烦。做完了这些工作我越来越觉得COMSOL仿真二次谐波产生真正有价值的部分不在软件操作而在把物理过程拆成“源项、相位匹配、模式重叠”这三个可以分别量化的环节。每一个环节对应一个参数或一个后处理量只要一个一个地验证模型就能从“能跑”变成“可信”。我个人调试模型时养成的最后一个习惯是每次只改一个参数扫完一个变量之后记录结果再改下一个。这个习惯看着慢实际上帮你省掉大量因为多因素纠缠而浪费的时间。如果你照着前面流程搭好了模型先固定其他参数只扫一遍极化周期看到倍频输出从几乎为零跳到峰值的那一刻你对SHG的理解就算真正落地了。