PPLN波导倍频仿真：从周期计算到COMSOL建模全流程解析

做PPLN波导倍频仿真这件事很多学光学的人一上来就想直接把整个3D波导扔进COMSOL然后吃满内存等一晚上出来一片乱码。其实真正能落地、能复现、能帮你写论文和改实验方案的跑法是先把原理掰开再一层层把模型降维拆开。我这篇就按我实际调通的思路从PPLN周期怎么算、COMSOL模型怎么建、非线性极化怎么加、周期怎么扫描一直写到常见问题排查全串一遍。适合正在做集成光子学、激光频率变换或者刚被导师安排做薄膜铌酸锂波导仿真的人。1. 倍频原理与PPLN周期估算先算对周期再谈建模1.1 倍频的物理过程与相位匹配为什么重要倍频也叫二次谐波产生SHG本质上是两个基频光子在一个非线性介质里湮灭产生一个频率为两倍、波长减半的光子。这个过程的强度完全取决于三件事非线性系数的有效大小、光场的功率密度、以及相位匹配条件。前两个好理解最后一个才是仿真里最容易翻车的点。相位匹配的物理含义很直白。光沿着波导传播时基频光和倍频光的速度不一样。因为材料有色散两个频率对应的折射率不同所以它们在传播过程中会产生一个不断变化的相对相位。非线性极化产生的倍频光理论上应该随着传播距离不断叠加增强但实际上不同位置发出的倍频光之间如果相位对不上就会互相抵消结果倍频功率在传播方向上振荡而不是单调增长。所以相位匹配决定了你能不能在一个合理的长度内得到可用的倍频输出。完全相位匹配时倍频功率随波导长度平方增长完全失配时功率在相干长度 (L_c\pi/\Delta k) 内涨落长度再长也没有意义。这也是为什么很多新手把模型建得漂漂亮亮结果倍频功率长得跟心电图一样说白了就是相位没对上。1.2 为什么选PPLN做准相位匹配相位匹配的实现方式主要有三类双折射角度匹配、温度匹配和准相位匹配QPM。双折射匹配是最传统的方式利用晶体两个偏振方向折射率不同的特性让基频光和倍频光在特定角度入射时实现折射率相等。问题在于角度匹配对方向特别敏感而且只能用正交偏振往往用不上晶体最大的非线性系数。温度匹配则是利用材料折射率随温度漂移的特点来精细调谐但在波导结构里作用有限因为波导中的色散曲线还受几何尺寸影响。PPLN走的是另一条路它的非线性系数符号沿传播方向周期性翻转每隔半个相干长度就翻转一次让原本因为色散而累积的相位失配被周期性补偿回来。这相当于主动给相位“续命”好处非常明显可以用同偏振态相互作用直接用铌酸锂最大的非线性系数d33约27 pm/V不需要苛刻的入射角度降低对准难度只要把周期做对理论上可以在任意目标波长上实现相位匹配。所以做波导倍频PPLN基本是绕不开的选择。1.3 周期计算公式与手算示例PPLN周期的核心公式不复杂[ \Lambda \frac{\lambda_1}{2(n_2 - n_1)} ]其中 (\lambda_1) 是基频波长(n_1) 是基频光对应的有效折射率(n_2) 是倍频光对应的有效折射率。这里我特意写“有效折射率”而不是材料折射率因为实际波导里算的是模式折射率而不是体材料折射率。拿最常见的1550 nm基频光倍频到775 nm为例。铌酸锂e光在1550 nm附近的材料折射率约2.137在775 nm附近约2.211代入公式[ \Lambda \frac{1550}{2 \times (2.211 - 2.137)} \approx 10.47 \text{ μm} ]这就是一个非常典型的PPLN波导周期量级。实际薄膜铌酸锂波导因为模式色散的影响最优周期通常会在9.8到11 μm之间浮动所以仿真里一定要把周期设成可扫描参数。建模的第一步不是画图而是先在笔记本上把这个数算出来记下来。后面所有模型参数、扫描范围都围绕这个数值展开。2. 模型简化与总体设计别一上来就建3D2.1 全三维模型为什么难做如果你老老实实建一个真实的3D波导模型几何尺寸大概是这样波导截面宽度0.5到1.5 μm厚度0.3到0.8 μm波导长度却要几百微米甚至几毫米。而PPLN周期才10 μm左右典型长度500 μm就有50个周期每个周期里还要解析倍频波长约350 nm的电磁场振荡。有限元网格的思维是每个波长至少5到10个单元倍频在材料里的波长约0.35 μm算下来沿传播方向要几千上万个单元。再加上横向剖面和PML边界3D模型自由度动辄千万级普通工作站跑起来基本是煎熬。更坑的是非线性耦合本身让方程之间强耦合网格一旦不均匀或者周期边界处分辨率不够求解器会频繁报错修一次模型要等很久。2.2 分层建模策略先2D模式分析再做2D纵向耦合我的建议永远是把问题拆成两层第一层是二维模式分析。只看波导横截面求解基频和倍频的模式有效折射率确定模式轮廓和传播常数。这一步计算量小几分钟就能出结果主要用来修正周期并确认模式匹配度。第二层才是核心的倍频耦合仿真。用二维纵向剖面模型模拟光的传播方向波导的横向限制用一个理想化的模式分布代替重点是研究周期极化结构下基频到倍频的能量转移。这种思路牺牲了一部分三维波导的精确度但保留了对倍频物理过程最关键的相位关系和非线性耦合计算量降低一两个数量级非常适合前期的参数扫描和工艺容差分析。等这个模型跑通了再用波束包络法或者3D模型做最终验证效率会高很多。2.3 二维纵向剖面模型的几何搭建具体来说我常用的几何结构是在x-z平面上建立纵向剖面z方向是光的传播方向波导长度L设为500 μm左右约50个PPLN周期x方向从下到上依次是衬底层、铌酸锂薄膜波导层、空气包层衬底材料用SiO₂折射率约1.45铌酸锂芯层厚度约0.6 μm用来支持近红外的单模传输上包层是空气。有一个容易被忽略的关键设置在COMSOL的2D电磁模型中默认面外方向是无限延伸的但实际波导在y方向有有限宽度。所以在计算功率和转换效率时需要在物理场设置里指定“面外厚度”比如设为0.5 μm。如果不设置系统默认1 m你算出来的倍频功率会凭空大出一百万倍。2.4 COMSOL模块与物理场如何选做倍频仿真最理想的情况是装了非线性光学模块Nonlinear Optics Module其中有一个专门的“电磁波二次谐波产生”Electromagnetic Waves, Second Harmonic Generation物理场接口它在一个模型里同时求解基频和倍频两个频率的电场分布并在非线性极化项中自动耦合。如果没有这个模块也不是完全没法做。可以退而求其次用波动光学模块的“电磁波频域”接口做两步法先求基频场再把基频场结果作为倍频频率下的非线性极化电流源。但这种方法公式推起来比较繁琐而且在COMSOL里手动加入非线性电流项需要改成PDE形式门槛反而更高。所以我建议有条件的直接用SHG接口。文章后面的操作流程也是围绕SHG接口展开但周期计算、几何建模、网格策略这些思路换到其他工具同样适用。3. 材料参数与周期极化结构设定仿真成败的关键3.1 铌酸锂各向异性与色散参数怎么填铌酸锂是单轴晶体折射率分为寻常光折射率no和异常光折射率ne。倍频过程用哪个折射率取决于你的偏振方案。最常见的PPLN方案是e光e光相互作用也就是基频和倍频都用异常光偏振这样能用上d33系数。那么COMSOL里材料怎么填在SHG接口中基频场和倍频场各自需要一套折射率。以1550 nm基频、775 nm倍频为例基频ne为2.137倍频ne为2.211。注意这两个折射率是体材料折射率。如果是薄膜铌酸锂波导在横向模式分析中得到的是模式有效折射率材料折射率并不会变变的只是模式传播常数。这些值要分别填到两个频率场的材料属性中千万别只填一个频率的属性。更严谨的做法是把Sellmeier色散公式直接写成耦合变量这样扫波长时就不用一个个手动改参数。不过Sellmeier公式版本很多不同文献的拟合系数略有差别实际项目中建议根据你晶体的供应商测试数据来校准。3.2 非线性系数d33与等效deff的关系嵌锂酸锂最大的二阶非线性系数是d33典型值约27 pm/V。在COMSOL的SHG接口里非线性极化强度与电场的关系通过二阶非线性张量建立可以直接按d矩阵形式输入也可以转成二阶电极化率张量χ(2)。这里有个新手容易犯的错误有人看到文献里的PPLN有效非线性系数是2/π乘以d33就直接把这个等效值填进去。但仿真软件里的非线性材料参数是局部物理参数如果你已经把周期性畴结构建出来或者用变量做了符号翻转那每个位置的非线性系数就应该填真实的d33等效系数是仿真结果自动浮现出来的不需要手动除π。如果你填了等效值又把符号翻转做了一遍相当于多除了一个π效率会偏低。至于其他非线性系数分量d31约4.3 pm/Vd22约2.1 pm/V同偏振e-e方案里其他分量贡献很小可以简化忽略。但如果做的是全矢量仿真建议把整个d张量都填进去避免人为遗漏一些非主项产生的干扰。3.3 周期极化结构的两种实现方式PPLN的周期极化结构在几何上有两种做法。第一种是把波导沿z方向切成几十个小矩形域每个域代表一个铁电畴奇数域和偶数域的非线性系数设为相反符号。这种做法的优点是每个域的材料属性一目了然缺点是几何对象数量太多。50个周期就要100个域后期扫参数时几何重建会变慢还容易因为相邻域的网格不匹配导致求解器报错。第二种是我更推荐的做法不建周期几何直接用一个坐标相关的变量来控制非线性系数的符号。在全局定义里写一个表达式d_pol d33 * sign(sin(2*pi*z/Lambda))然后在非线性极化特征里把d矩阵的有效值指向这个变量。这样几何里只有一个连续波导域网格干干净净周期调整也只需要改Lambda这个参数而且扫描周期时不用重建几何速度快得多。这个做法背后没有物理损失因为铁电畴反转本身就是一个空间上的符号切换用符号函数来建模是非常自然的数学描述。3.4 波导损耗与温度效应要不要加仿真初期建议先不加损耗把问题简化。等你已经在无损耗模型里拿到了稳定的倍频曲线再尝试加入材料吸收损耗和波导散射损耗观察效率下降幅度。温度效应主要通过折射率变化影响相位匹配条件实际中热效应很关键因为激光功率稍微大点PPLN晶体的温度变化就会让周期失配。如果做的是实验室级仿真建议把温度设为全局参数用温度相关的Sellmeier公式替换固定折射率然后做温度扫描看最佳周期向哪边漂移。4. 实操步骤从模式分析到倍频仿真完整流程4.1 步骤一二维模式分析找模式有效折射率先不要碰纵向模型单独建一个波导横截面的2D几何。这一步目的是算出基频和倍频在这个波导结构里的模式轮廓和有效折射率。物理场用“电磁波频域”频域设置为基频1550 THz对应的角频率研究步骤选“模式分析”。模式分析的目标是找到支持哪些传播模式以及对应的模式有效折射率。对一片0.6 μm厚的铌酸锂薄膜加上上下包层折射率对比基频条件下通常能支持一个或两个TE模式。跑完模式分析把基频TE0模式的neff记录下来然后把频率改成倍频775 nm再跑一次模式分析得到倍频模式的neff。这个倍频模式可能与基频模式形状差异很大如果两个模式的重叠积分太小就得调整波导厚度或宽度。最后用这两个neff重新计算周期[ \Lambda \frac{\lambda_1}{2(n_{\text{eff},2} - n_{\text{eff},1})} ]这一步得到的周期才是真正匹配这个波导结构的周期比用材料折射率算出来的要准。4.2 步骤二建立纵向SHG模型与端口设置横向模式分析完成后新建一个2D纵向剖面模型。几何与前面2.3节描述一致坐标z从0到Lx从下到上覆盖衬底、波导层、空气。物理场选择“电磁波二次谐波产生”接口。默认情况下它会生成两个耦合的场变量基频场和倍频场。我这边COMSOL里变量的命名通常是e1和e2不同版本可能略有差别如果你用的是5.6以前的版本建议在帮助文档里确认一下变量前缀。在物理场设置中左边界z0设置端口端口模式类型选择“端口”基频场侧输入入射功率P1比如1 mW倍频场同样在左边界或右边界设置端口但入射功率设为0作为倍频输出的测量端口其他边界设置散射边界条件模拟开放边界如果模型长度较短建议在左右两端加PML层吸收残余反射。关键点是倍频端口必须耦合正确的模式。如果端口模式没设置对你会得到“端口模式对象不匹配”的警告或者倍频输出全是数值噪声。4.3 步骤三非线性极化项与周期变量的接入在LN芯层域中添加“非线性极化”特征。这个特征的作用是告诉COMSOL这个区域存在二阶非线性响应。在非线性极化设置中指定二阶非线性张量的分量。如果使用d形式需要把d33对应的张量元素写成d33 * sign(sin(2*pi*z/Lambda))同时把材料参数中的d33设为27 pm/V。如果COMSOL的版本支持直接在表达式中引用变量我可以定义一个全局变量d33 27e-12 [m/V] Lambda 10.47[um]然后在非线性张量里直接引用。如果你使用的是3D或更复杂的波导结构d张量的方向可能与全局坐标系不同需要做坐标旋转。这时候建议先查清楚晶体切割方向和光轴相对全局坐标系的方位。4.4 步骤四网格划分的核心思路纵向模型的网格划分要分区域处理波导芯层是最关键的网格最大尺寸取倍频波长的六分之一到八分之一即50到60 nm左右空气包层和衬底可以稍大但靠近波导界面的部分也要过渡加密沿z方向建议用映射网格保证每个PPLN周期内网格数量一致避免周期边界因为网格疏密不同引入数值反射。如果模型长度是500 μm这个网格规模在2D模型下大概几十万个自由度普通工作站几分钟就能算完不会太夸张。但如果把长度加到2 mm网格量会成倍增长建议前期用短波导验证物理过程后期再拉长。求解器方面SHG接口的非线性耦合比较强建议先用直接求解器PARDISO。如果出现不收敛优先检查是不是泵浦功率设置过高导致非线性迭代发散。可以用参数化扫描把泵浦功率从0.01 mW逐步加到目标值获得更稳定的收敛路径。4.5 步骤五参数化扫描周期与功率模型跑通后最重要的两个扫描就是周期扫描和功率扫描。周期扫描范围设计在理论周期的正负5%左右。以10.47 μm为例扫描范围设为9.8到11.1 μm步长0.02 μm。每次计算都重新求解完整非线性问题最后得到倍频功率对周期的响应曲线。曲线上应该出现一个明显的峰值峰值位置就是当前波导的最优周期。泵浦功率扫描相对简单。从0.1 mW到100 mW按倍数步进扫。低功率区间转换效率随功率线性增长高功率区间会出现基频耗尽效率增长放缓。这一步除了看效率还能帮你判断泵浦耗尽是否需要在设计实验时考虑。5. 结果处理与细节分析从功率数值到物理规律5.1 如何正确提取倍频输出功率COMSOL里提取端口功率的方式比较灵活。最简单的办法是在“派生值”里选择“全局计算”然后利用端口变量计算功率。对于SHG接口倍频场在输出端的功率可以通过对zL边界上的坡印廷矢量积分来获得。坡印廷矢量的z分量表达式是[ P_z \mathrm{Re}(0.5 \cdot E_{2x} \cdot H_{2y}^* - 0.5 \cdot E_{2y} \cdot H_{2x}^*) ]注意乘以面外厚度。如果这部分觉得复杂也可以直接在端口特征的结果里读取S参数倍频输出功率近似等于入射基频功率乘以倍频端口的透射系数但非线性情况下S参数的意义需要小心理解我更推荐直接用坡印廷积分。转换效率就是倍频功率除以基频入射功率通常记作(\eta P_2/P_1)。为了方便对比不同波导长度和功率还可以归一化成功率归一化效率(\eta_{\text{norm}} \eta/(P_1 \cdot L^2))单位%/(W·cm²)。5.2 周期失配曲线与sinc²分布验证扫周期结束后画一条以周期为横轴、倍频功率为纵轴的曲线。如果是认真的仿真这条曲线应该长得很像sinc²函数中心峰值对应最优周期两侧有递减的旁瓣。如果扫出来的曲线不是sinc²样而是多个杂峰或者干脆平得没有峰第一反应应该是看周期变量有没有真正影响到非线性系数。检查一下符号函数表达式确认z坐标在模型中确实沿波导传播方向。另一个可能的坑是波导长度太短导致相干长度差异不明显。比如只有5个周期时周期失配造成的效率变化可能还没体现出来曲线自然不漂亮。至少要跑到20个周期以上再考虑曲线形态。5.3 模式重叠对转换效率的影响模式重叠度是波导倍频的一个隐藏变量。基频模式和倍频模式在横截面上的重叠积分越大非线性转换效率越高。如果两者模式形状差异很大即使周期完美匹配效率也会低得可怜。在2D纵向模型中模式匹配体现在端口激励的时候你是否选对了模式。如果基频端口激励的是TE0模式倍频端口也应该对应TE0模式的输出。如果波导设计得不好倍频TE0模式截止了那就什么也输出不了。需要在模式分析阶段就检查倍频波长下你的波导结构是否存在至少一个导模。如果倍频模式损耗太大或截止就调整波导芯层厚度或者换成更大折射率差的包层材料。5.4 泵浦耗尽带来的基频功率下降小信号假设下基频功率沿波导基本不变倍频功率按抛物线增长。但当泵浦功率升高基频场被明显消耗倍频增长就不再是理想抛物线。仿真中观察基频输出功率随泵浦的增加而下降就是泵浦耗尽出现的证据。这个效应在连续波倍频实验里要到几百毫瓦甚至瓦级才明显但波导里光场被紧束缚功率密度极高几毫瓦就可能看到。这也是为什么波导倍频仿真尤其需要关注泵浦耗尽它决定了实验中最大可用转换效率。6. 常见问题与排查技巧实录6.1 倍频输出功率为零或数值异常小这是我见过最多的问题。产生这种情况的原因排在前面的有非线性极化特征没有添加上去或者添加到了错误的域d33数值单位写错pm/V和m/V之间差了12个数量级周期符号函数表达式里的变量名和全局参数名不一致导致计算出的d_pol恒定不变端口模式没有设置对倍频场在端口处匹配不了导模。排查方法比较实用的是把周期表达式临时改成常数d33然后跑一个光强特别低、长度特别短的模型看倍频输出是否还有微弱但非零的响应。如果是说明非线性项正常问题出在周期结构或相位匹配如果还是零说明非线性项本身没生效。6.2 扫描周期时看不到明显峰值扫周期不出峰首先要确认扫描的变量真正进入了非线性系数表达式。一个常见错误是符号函数表达式直接调用了物理场坐标z但模型坐标系里传播方向不是z而是x或y这样符号根本不随传播方向变化。另一个问题是步长选得太大把峰值跳过去了。周期响应峰的半高宽通常很窄和波导长度直接相关长度越长峰越窄。比如500 μm的波导最优周期的半高宽可能只有0.1 μm量级扫描步长如果设为0.5 μm很容易错过峰值。一开始扫粗扫描找到大致范围以后再加密。6.3 求解器不收敛或内存不足不收敛的原因通常是网格质量或者初始值问题。先在低泵浦功率下求解比如0.01 mW让非线性项对解的扰动足够弱。然后再打开辅助扫描从低功率逐步增加到目标功率每一步以上一步的解作为初始值收敛稳定性会大幅提升。内存不足时优先考虑减短波导长度不要从500 μm开始。先用100 μm验证物理过程确认所有设置正确后再拉长。6.4 仿真最优周期与实验数据对不上仿真算出来的最优周期和实验最优周期存在偏差非常正常主要来源有三个。一是材料色散参数有误差不同批次铌酸锂晶体的Sellmeier系数会有微小差异温度、掺镁浓度都会改变折射率最终导致最优周期移动零点几微米。二是模式色散与材料色散的差别。仿真如果只用了材料折射率算周期没有把波导模式有效折射率代进去误差会很大。三是制造误差。PPLN工艺的畴周期精度通常在几十纳米到几百纳米对倍频相位匹配的影响不容忽视。所以实验上做倍频几乎都要扫周期就是为了覆盖这些误差。6.5 二维纵向模型和三维仿真、实验的结果能差多少务必对二维纵向模型的定位有清醒认识。它能准确反映周期、相位失配、泵浦耗尽这些核心物理趋势但给出的转换效率绝对值只能作为数量级参考。原因在于它忽略了侧向模式的完整分布面外厚度只是一个等效参数。实际三维波导的模式重叠积分可能和二维模型的假设差不少所以我在对比实验数据时会用二维模型得到物理趋势再用三维模型或实验来校准绝对效率。个人建议的流程是二维模式分析定周期二维纵向SHG模型扫趋势、选参数最后三维模型出头版数据实验数据再回来修正材料色散参数。这套链路我在多个芯片设计里都验证过比一股脑建3D模型靠谱得多。最后再分享一个小技巧。PPLN周期扫描是最便宜的实验不光在COMSOL里便宜在流片之前它几乎能帮你把所有参数风险摸一遍。你只需要让周期以0.02 μm的步长扫过理论值附近就会看到倍频功率从零爬到峰再跌回零这个峰就是你的工艺容差和设计中心的直接证据。所以不管模型多复杂周期扫描永远值得先跑。