X切LNOI波导倍频仿真：COMSOL建模与相位匹配实战

最近研究X切型绝缘体上铌酸锂薄膜LNOI的倍频SHG转化效率COMSOL仿真前前后后跑了一个多月越跑越觉得这东西比想象中有意思得多。LNOI这两年几乎是集成光子学里的“顶流”平台几百纳米厚的单晶铌酸锂薄膜配合氧化硅埋层形成高折射率差光可以被死死压在亚波长截面的波导里。这个结构先天适合做非线性频率变换尤其是倍频因为模式体积小光场强度高波导色散又给了你调相位匹配的空间。今天就把我这阵子摸索出来的仿真思路、代码片段和踩坑记录整理出来希望能让刚接触这个方向的朋友少走点弯路。这篇内容我尽量用“工程化”的口吻写不追求论文级别的严谨但求每一步能落地。适合正在做LNOI波导倍频设计、准备用COMSOL做模式分析和效率预估、或者搞不清楚X切型材料方向怎么给的人。你不需要是COMSOL高手但需要大概知道波动光学模块长什么样。下面从物理背景开始一路讲到代码和后处理最后是问题排查。1. 为什么大家都在盯X切LNOI做倍频1.1 材料底子好薄膜平台更占便宜铌酸锂在非线性光学材料里属于“全科生”二阶非线性系数大透光范围覆盖可见到中红外还有很好的电光、声光效应。其中最常用的d33系数大约在27pm/V附近这比磷酸钛氧钾、砷化镓这些常见非线性材料高出不少。但块状铌酸锂倍频有个尴尬的地方——相互作用长度够长可光斑尺寸很难压下来转换效率靠高泵浦功率硬堆。到了LNOI薄膜时代情况完全不同。薄膜厚度通常只有300到700nm脊形波导的模场面积可以做到1μm²以下微瓦量级的泵浦功率就能在波导里产生很高的功率密度。而且LNOI用的是单晶薄膜材料损耗低晶轴取向可以精确控制这对需要具体判定非线性张量方向的仿真非常关键。1.2 X切型到底切出了什么铌酸锂是三方晶系常说的切型由晶轴和薄膜法向的关系决定。X切型表示薄膜表面法线沿晶体的X轴而光轴Z轴躺在薄膜平面内。Z切型则是表面法线沿Z轴光轴垂直于薄膜平面。这个区别在倍频仿真里非常重要。LNOI波导通常是在薄膜上刻出脊形结构光沿着波导方向传播。对X切型来说光轴Z在膜面内你可以设计一个电场沿Z方向偏振的横电模式正好激活最大非线性系数d33。而Z切型的TE模电场基本只在膜面内振荡和光轴垂直更多依赖d31这类较小的系数。所以在需要高效率倍频时X切型是常见选择这也是我做仿真时首选它的原因。1.3 仿真到底要回答什么问题做LNOI倍频仿真核心不是“把COMSOL跑通”而是回答几个具体问题波导截面形状、薄膜厚度和刻蚀深度怎么选才能让泵浦光和倍频光在有效折射率上满足相位匹配模式场分布能产生多大重叠在给定泵浦功率和波导长度下倍频输出功率到底是多少。这些问题如果只靠实验盲调成本极高。仿真能提前把可能参数区域扫出来缩小实验范围。这也是我这段时间一直用COMSOL反复算的原因。2. 仿真模型背后的物理2.1 SHG的起点二阶非线性极化倍频过程本质上是由二阶非线性极化产生的。在频率域里二次谐波极化强度可以写成P_i(2ω) ε0 * χ_ijk(2)(-2ω; ω, ω) * E_j(ω) * E_k(ω)这里的χ(2)是二阶非线性磁化率张量而工程上更习惯用非线性系数dijk两者关系是χ(2) 2d。铌酸锂的d张量有确定的空间取向一旦材料坐标系和模型坐标系不对齐计算结果就会完全错误。仿真中我采用的思路是“泵浦不耗尽近似”假设倍频效率不高泵浦光在传播过程中基本不衰减。先求出泵浦模再由泵浦场平方得到倍频频率的非线性极化把它当作等效源去求解倍频场。这个近似在小信号转换效率低于百分之十几时非常可靠而且计算量远小于全耦合的三波混频求解。2.2 转化效率公式和有效模面积在完美相位匹配且无损耗情况下波导倍频转换效率可以用下面这个形式估计P_SHG (2ω² d_eff² L²)/(ε0 c³ n_p² n_s A_eff) * P_pump²其中ω是泵浦角频率d_eff是有效非线性系数L是相互作用长度n_p和n_s分别是泵浦光和倍频光的模式折射率A_eff是有效模面积。这个公式告诉我们三件事效率正比于泵浦功率的平方正比于长度平方反比于有效模面积。因此波导设计的目的就是把模场压缩到很小同时保证两个频率的模场尽量重叠。有效模面积由泵浦模和倍频模的场分布共同决定不是简单拿波导截面几何面积来算。公式里还会出现模式重叠积分所以在仿真后处理时不能用“一个矩形面积”代替。2.3 相位匹配从Δk到准相位匹配倍频要高效必须满足波矢匹配Δk k_2ω - 2k_ω 0换成有效折射率就是Δk (2ω/c) * (n_eff(2ω) - n_eff(ω))LNOI波导的优势在于泵浦光和倍频光在同一个几何结构里的模式色散可能刚好相等也可能非常接近。只要合理调节波导宽度、薄膜厚度、刻蚀深度、上包层材料等就有希望让n_eff(2ω)等于n_eff(ω)实现真正的“模式色散相位匹配”。如果色散怎么调都不为零还可以用周期性极化进行准相位匹配。周期Λ满足Λ 2π / |Δk|仿真时可以通过在非线性系数d_eff上引入周期性的符号反转来实现。手动建模周期极化有点麻烦但用参数化几何或定义空间依赖的材料参数也能做。2.4 COMSOL里怎么把非线性源“装”进去COMSOL的波动光学模块默认是线性系统不会自动算二阶极化。我的做法是在倍频频率的电磁波频域研究中添加外部电流密度作为源。因为时谐场里极化强度P和电流密度的关系是J ∂P/∂t iωP所以对倍频频率外部电流密度为J(2ω) i2ω * P_NL(2ω)在COMSOL域条件里把Jx、Jy、Jz的表达式写成由泵浦模式场分量平方组合成的非线性源项。这样求解出来的就是倍频场。方法不算复杂但要注意源项的表达式和材料坐标系的映射关系写错了效率曲线会非常奇怪。2.5 归一化效率怎么算才可靠做参数扫描时我会统一用“归一化倍频效率”来横向比较不同结构单位通常写成%/W·cm²或者%/W·cm。仿真流程里先固定泵浦功率为1W计算出倍频输出功率后再除以1W和长度平方。但COMSOL模式分析解出来的场是任意归一化的必须先把泵浦模场幅度缩放使其携带的真实功率等于1W再去平方构造非线性源。很多人结果离谱往往就是漏了这一步。3. COMSOL实操从建几何到出效率3.1 几何与材料参数设置以X切LNOI脊形波导为例我通常建二维截面模型传播方向用有效折射率来描述。几何从上到下依次是空气/二氧化硅覆盖层、铌酸锂脊形区、薄膜残留层、二氧化硅埋层、硅衬底。空气和二氧化硅覆盖层在实际器件中常见不能随手省略因为上包层会直接影响模式色散。常用起始参数LNOI薄膜厚度500nm脊宽800nm脊刻蚀深度300nm侧壁倾角稍微留几度模拟实际工艺BOX层厚度2μmPML吸收层放在最外侧。这些参数不是标准答案但适合作为初值扫描的原点。材料设置要非常小心。铌酸锂是单轴晶体折射率需要区分寻常光和非常光对应折射率no和ne。X切型意味着晶轴Z在膜面内材料坐标和模型坐标存在旋转关系。最简单的方式是直接在COMSOL材料节点里定义各向异性介电常数张量根据晶轴方向把主值转动到模型坐标系。SiO2、空气和Si在倍频波段吸收很小折射率设为常数即可。3.2 模式分析要把两个频率的模式都找到我习惯建立两个独立研究。第一个研究做模式分析求解泵浦频率下的本征模第二个研究做倍频频率的模式分析。也可以在同一模型里设置两个“电磁波频域”接口分别绑定不同频率但物理上它们并不会自动耦合。模式分析求解器用特征值求解器搜索一个有效折射率范围。比如泵浦波长1550nm薄膜波导有效折射率大概在1.8到2.1之间就搜索这个范围。倍频波长775nm有效折射率可能因为强色散掉到1.6到1.9需要单独设置搜索起点。每个频率我都要求前几个模式都算出来从中找重叠度最高且强度最集中的基模。网格也是关键。波导截面小场变化快我通常用“极细”网格并加上边界层网格。膜面内至少要有6到8个网格单元跨过脊宽否则有效折射率误差会直接影响相位匹配判断。3.3 从泵浦模到倍频源的完整流程求完泵浦模后真正的SHG仿真才开始。流程是这样的第一步把泵浦模场导出。用mphinterp或COMSOL后处理中的“表面最大值”确认场分布第二步对泵浦模做功率归一化。计算通过波导截面的平均功率然后给模式场乘一个系数让积分功率等于1W第三步根据归一化泵浦模场分量在倍频频率的电磁波频域研究中设置外部电流密度。表达式形如Jx 1i * 2 * omega * 2 * eps0 * (d31*Ey*Ez d32*Ex*Ez ...)这里省略号代表你需要根据有效非线性系数具体展开。COMSOL里还可以用变量定义把这一大串表达式封装起来避免在多个边界和域里重复写。第四步求解倍频频率下的受迫波动方程。由于源项已经确定这是一个线性求解问题不需要迭代第五步积分倍频频率处通过截面或某个监视边界的坡印廷矢量得到输出功率P_SHG。3.4 参数扫描与相位匹配曲线我通常把波导宽度、薄膜厚度、刻蚀深度设成参数化扫描变量。每个参数点都执行“模式分析两次倍频求解一次”最后把所有结果汇总。这里有个经验不要只记录最终SHG功率同时要把两个频率有效折射率随参数的变化导出来。画在一张图里相位匹配点往往一目了然——就是n_eff(2ω)与n_eff(ω)交点附近SHG功率出现尖峰。扫描时还要注意模式追踪。有效折射率随宽度变化会发生模式阶次交替某个宽度下你原本关注的基模可能排到第二或第三阶。所以我每次扫描都会输出几个模式的折射率曲线等位后再判断哪条分支才是需要的模式而不是盲目相信“第一阶就是基模”。4. 核心脚本与代码分析4.1 为什么要脚本化COMSOL界面上手动点当然能跑但LNOI倍频仿真有太多重复环节改一个宽度尺寸、重新剖网格、算模式、提取场、设源、求倍频、出图。手动操作不仅慢还容易漏同步。我选择用LiveLink for MATLAB把建模流程脚本化。这样每组参数都能自动跑半夜挂着扫参数第二天起来直接看结果。下面这段代码是整理后的核心流程不是完整工程文件但展示了从新建模型、模式分析到导出泵浦场的关键命令。COMSOL版本不同API略有差异但整体逻辑一致。import com.comsol.model.* import com.comsol.model.util.* model ModelUtil.create(Model); model.component.create(comp1, true); model.component(comp1).geom.create(geom1, 2); % 参数定义 model.param().set(w, 0.8[um]); model.param().set(h_film, 0.5[um]); model.param().set(h_etch, 0.3[um]); model.param().set(lam0, 1.55[um]); model.param().set(omega, 2*pi*c_const/lam0); % 创建几何脊形波导截面 % 这里只画一个示意实际需要多个矩形合并 model.component(comp1).geom(geom1).create(r_air, Rectangle); model.component(comp1).geom(geom1).create(r_ridge, Rectangle); model.component(comp1).geom(geom1).create(r_slab, Rectangle); model.component(comp1).geom(geom1).create(r_box, Rectangle); model.component(comp1).geom(geom1).run; % 电磁波频域接口 model.component(comp1).physics.create(ewfd, ElectromagneticWaves, geom1); % 材料设置略主要是各向异性折射率张量 LNOI % 模式分析研究 model.study.create(std1); model.study(std1).create(mode, Eigenfrequency); model.study(std1).feature(mode).set(eigenfunctionSearch, 2.0); model.sol.create(sol1); model.study(std1).feature(mode).attach(sol1); model.sol(sol1).run; % 导出泵浦模场数据和有效折射率 coord [0; 0.5e-6]; % 某个采样点或截面坐标 [Ex, Ey, Ez] mphinterp(model, {Ex,Ey,Ez}, coord, coord, dataset, dset1); neff_p mphglobal(model, ewfd.neff, dataset, dset1); % 后续需要把场归一化到1W再作为非线性源这段代码里最关键的是mphinterp的返回值。COMSOL模式分析解出的电场是复数有实部和虚部实际取回后还需要做功率归一化。我不会直接在源项里用原始场因为一旦归一化系数算错效率会差好几个数量级。4.2 功率归一化和非线性源计算归一的思路是先通过后处理得到泵浦模在波导截面上的时间平均功率再算一个缩放系数。这里我通常写一小段MATLAB后处理% 假设之前已经提取了泵浦模的电场和磁场 P_int int_surface_pooynting; % 从COMSOL积分得到 P0 1.0; % 目标泵浦功率单位W scale sqrt(P0 / P_int); Ex_norm Ex * scale; Ey_norm Ey * scale; Ez_norm Ez * scale; % 构造非线性极化源示例分量实际要按张量展开 eps0 8.8541878128e-12; d33 27.0e-12; % 单位m/V Pz_nl eps0 * 2 * d33 * Ez_norm .* Ez_norm; % 简化示意 Jz_src 1i * 2 * omega * Pz_nl;然后把这个Jz_src通过变量或函数写回COMSOL倍频研究的外部电流密度节点。注意这里的表达式是示意性的。真实X切LNOI中d33作用的偏振方向和模式场分量对应关系必须从张量旋转里推导不能想当然复制这个公式。4.3 倍频求解和结果提取新研究不需要再算本征模只需要“电磁波频域”在倍频频率下加上外部电流密度。求解完成后我用mphinterp在输出边界上取坡印廷矢量的法向分量再积分。% 提取倍频场 model.param().set(omega2, 2*pi*c_const/(lam0/2)); [Ex2, Ey2, Ez2, Hx2, Hy2, Hz2] mphinterp(model, ... {Ex,Ey,Ez,Hx,Hy,Hz}, ... coord, coordOnSection, dataset, dset2); % 计算坡印廷矢量并积分得到P_SHG P_SHG real(0.5 * sum(Ex2 .* conj(Hy2) - Ey2 .* conj(Hx2))) * dA; % dA是积分微元的面积这里我用了一个简化近似求坡印廷矢量。实际COMSOL后处理有自动积分表面功率的算子建议直接用它避免手写积分出错。我写这段只是为了说明脚本思路。4.4 脚本防御点跑脚本时我吃过几个亏一是忘了在求解前把网格更新导致几何改了但网格没跟着变二是模式搜索范围太窄某些参数点模式漏掉三是材料坐标系没有跟着参数化角度变化。这些都要在脚本里加入检查逻辑。比如扫描结束后立刻打印每个点的有效折射率如果某个点突然跳变很大赶紧回查模式阶数。5. 常见问题与排查技巧5.1 模式找不到或有效折射率对不上很多人第一步就卡在模式分析。现象是特征值求解器返回很奇怪的折射率或者找不到目标模式。这个问题九成出在网格和搜索范围上。LNOI薄膜波导阶数高模式分布密集网格太粗会把模式“磨”没。建议先用非常细的网格跑一个点确认模式场分布合理再放宽网格做扫描。搜索范围也可以设置大一点比如1.5到2.5等模式出来后再缩小范围追踪特定模式。5.2 倍频效率高得离谱如果算出的SHG功率比泵浦功率还高肯定不是物理。常见原因有三种一是泵浦模没有归一化到1W导致源项强度失真二是倍频频率的吸收或边界反射导致场叠加异常PML没吸收干净三是网格问题导致场奇点处平方项爆炸。我的检查顺序是先看泵浦模和倍频模的功率积分是否与预期一致再加密PML区域网格最后才怀疑公式写错。5.3 材料方向和张量矩阵错乱X切LNOI里最阴间的坑就是坐标旋转。COMSOL全局坐标是笛卡尔坐标但铌酸锂的Z轴可能躺在模型平面任意方向。如果你只在材料节点填一个对角折射率张量那默认主轴完全和全局坐标重合很可能等于用了Z切或Y切。我的做法是先用一个极简单的平板波导做验证让光传播方向沿Y轴电场沿Z轴偏振算出的模式折射率应该接近ne如果接近no说明坐标系旋转没设对。5.4 相位匹配扫描看不到峰值扫完宽度发现效率一片平坦没有尖峰别急着怀疑物理模型。先看看你画的是不是两个频率的有效折射率差。我遇到过扫描范围不够宽根本没跨过零失配点的情况也遇到过模式阶次跳变前后跟踪的不是同一个模。建议先导出n_eff随参数变化找到交点附近再局部加密扫描。5.5 和实验测量对不上仿真里效率很高实验测出来只有几分之一甚至一个数量级差距这很正常。实验里还有波导侧壁粗糙度、刻蚀损伤、端面耦合损耗、实际泵浦功率标定误差等因素。仿真能尽力做到的是把相位匹配位置、波导结构参数的相对趋势算准。如果实验中最佳结构宽度跟仿真差几十纳米先检查是不是工艺侧壁角与仿真不一致这个影响往往比折射率取值还大。写在最后的一个实操习惯跑LNOI倍频仿真给我最大的感受是模型“接线”比求解器更费心思。材料张量方向、泵浦功率归一化、模式追踪任何一个环节出问题结果都能美得离谱或者丑得诡异。我现在每次开新模型前都会先花半小时做一个平板波导验证把X切晶轴方向、d张量作用路径这些基础设定验干净然后再上脊形波导和周期极化。这步看起来慢实际省下的排查时间远不止半小时。再分享一个小技巧参数扫描时别只盯SHG功率把有效折射率和模式电场图一起导出。很多看似反常的现象比如效率峰值偏移、谱线不对称、模式串扰看折射率曲线马上就能解释。仿真做到最后往往就是用一条简单的色散曲线说服自己下一步改哪里。希望这篇语无伦次但全是实操的分享能给你的LNOI倍频仿真省下几天宝贵时间。