★ FEATURED ARTICLE 简介面向非线性光学仿真与数值计算的完整工程资料包适合光学工程、物理及相关专业的学生与研究人员侧重二次谐波产生、参量下转换、四波混频等典型非线性效应覆盖光场传播、相位匹配、材料极化率与转换效率计算等核心问题。压缩包内共包含五百七十七个文件核心文件类型包括MATLAB脚本、C源文件与头文件、预编译的MEX内核以及Visual Studio工程配置、CHM帮助文档、数据文件与图示文件全部内容压缩后约五点二三MB目录结构清晰便于按功能模块检索和学习。当前已有八百五十三人学习或下载。压缩包内提供了信号计算、最小误差内核修正、二次谐波求解等仿真代码与MEX工具并配有工程配置和说明文档可帮助学习者将非线性光学原理转化为可运行的程序实践也适合在此基础上扩展新的数值实验深入理解高功率光场与材料相互作用时的各种非线性效应。1. 拿到 Nonlinear-Optics 代码包先搞清楚它在算什么相位匹配比非线性方程更要命做光学仿真的人第一次打开 Nonlinear-Optics 这类代码包第一反应多半是去找非线性薛定谔方程的求解器。我的经验是反过来的这类项目里最容易翻车的从来不是方程本身解不动而是你根本没把相位匹配算对。相位匹配没对准再高阶的积分器、再细的网格都白搭。这篇文章把非线性光学仿真最常用的一条技术路线讲透从 I 类倍频切入用耦合波方程加 RK4 把转换效率算出来再逐一扫角度、温度这些调谐参数最后落在代码和事务怎么组织上。适合正在做倍频、OPO、OPA 或者激光频率变换的工程师和研究生照着改参数就能复现。2. 先把物理模型立住从 I 类倍频到三波互作用的耦合波方程2.1 为什么先啃 I 类倍频三波互作用里最省事的试金石倍频SHG是非线性光学里最经典、也最容易做数值验证的过程一束 1064 nm 的基频光打进晶体出来一束 532 nm 的倍频光。它背后是三波互作用最简形式——两个基频光子湮灭产生一个倍频光子信道只有一个极化率只看 χ⁽²⁾。我一般建议新手或者新项目先从 I 类倍频入手理由有三个。第一实验数据遍地都是BBO、LBO、KTP 的倍频效率曲线在文献里能查到一大堆仿真对得上对不上一眼就能看出来。第二它的耦合波方程只有两个数值上是最低门槛却已经包含了后面 OPO、OPA 的几乎所有核心要素三个波、一个相位失配量、一组耦合系数。第三从 I 类倍频出发把方程从两个扩到三个就是完整的光参量放大OPA再加一个谐振腔边界条件就是 OPO。所以这一章把倍频的物理模型啃下来后面全是机械扩展。很多人会问为什么不用 FDTD 或者非线性薛定谔方程来做。FDTD 确实能算非线性过程但代价是三维网格加时间步进一个毫米级晶体的计算量非常可观而且它给的是场分布不是工程师最关心的转换效率和匹配带宽。连续波或者纳秒脉冲场景下慢变包络近似下的耦合波方程是行业里最常用的做法几行 RK4 就能解决没必要上重型求解器。至于非线性薛定谔方程那是处理自聚焦、自相位调制这类效应时才需要的倍频场景里如果泵浦功率还没到损伤阈值这些效应通常可以放一边。2.2 耦合波方程怎么落到可数值求解的形式从麦克斯韦方程组出发做三个假设平面波、无损耗晶体、慢变包络近似。慢变包络的意思是场的振幅包络在波长尺度上变化很慢所以可以把快振荡的相位项单独拎出来处理。做完这些假设二阶非线性极化项会以驱动项的形式出现在波动方程里最后整理成一组一阶常微分方程也就是耦合波方程。对 I 类倍频两个基频光子合成一个倍频光子场振幅的演化可以写成dE1/dz i (ω1·deff)/(n1·c) · E2·E1*·exp(−iΔkz)dE2/dz i (ω2·deff)/(2·n2·c) · E1²·exp(iΔkz)E1 是基频场E2 是倍频场deff 是有效非线性系数Δk k2 − 2k1 是相位失配量。第一个方程描述泵浦光在转换过程中的损耗第二个方程描述倍频光的增长。注意第二个方程前面那个 1/2 因子它来源于 E² 展开后 2ω 分量前面的系数丢了它能量守恒直接破功效率算出超过 100% 就是从这里来的。这套方程本身没有解析解除了小信号近似下所以它的落地形式必然是数值积分。把实部和虚部分开或者直接用复数场推进沿着晶体长度 z 一步一步算每一小步里场和场之间的耦合按上面的右端项更新。参数上的关键点在于deff 的量纲是 m/Vn1、n2 是折射率ω1、ω2 是圆频率这几个量一旦单位混了结果会差出几个数量级后面第 5 章专门讲这个坑。2.3 相位匹配条件与失配量的物理直觉为什么一定要相位匹配因为介质有正常色散基频光和倍频光在晶体里的折射率不一样走的速度不一样。两列波在前进过程中不断累积相位差转换出来的倍频光和后面新产生的倍频光相位对不上就会发生能量倒流——先转过去的又转回来。这就是所谓 back-conversion效率曲线长到一定长度后不升反降的原因。所以 Δk 是整件事的灵魂。Δk 为零时转换效率随长度近似二次方增长Δk 不为零时效率变成振荡的峰值被压低。真实晶体里 Δk 几乎不可能天然为零因为任何材料的色散都让 k2 大于 2k1。工程师做的事就是通过角度、温度、周期极化这三种手段把 Δk 压回零附近。相位匹配不是一个点而是一个窗口。窗口的宽容度由晶体长度决定长度越长允许的 Δk 偏差越小角度或温度的容忍度就越窄。这个直觉在调参时非常重要——你扫角度扫出来的峰宽度往往只有零点几度不是程序写错了是物理本来就这个脾气。3. 用 Python 把倍频效率算出来一套可直接改的 RK4 脚本3.1 最小脚本RK4 积分两个耦合方程先定义参数区把所有物理量统一到 SI避免中途换算import numpy as np # 物理常数 c 2.99792458e8 # 光速, m/s # 晶体与光束参数BBO I类基频 1064 nm - 倍频 532 nm wl1 1.064e-6 # 基频波长, m wl2 wl1 / 2.0 # 倍频波长, m n1 1.655 # 基频 o 光折射率近似精确值查 Sellmeier 方程 n2 1.556 # 倍频 e 光在匹配角附近的折射率近似 deff 2.0e-12 # 有效非线性系数, m/V数量级近似 I0 500e4 # 峰值功率密度, W/m^2即 500 MW/cm^2 Lc 3.0e-3 # 晶体长度, m w1 2 * np.pi * c / wl1 w2 2 * np.pi * c / wl2 k1 2 * np.pi * n1 / wl1 k2 2 * np.pi * n2 / wl2 dk k2 - 2 * k1 # 相位失配量, 1/m # 振幅归一化让 |A|^2 直接对应功率密度W/m^2 A10 np.sqrt(I0) A20 0.0 # 两个方程的耦合系数 k1_cpl w1 * deff / (n1 * c) k2_cpl w2 * deff / (2 * n2 * c)这套初始化的思路是不把场振幅定义成电场强度 V/m而是把 |A|² 直接映射到功率密度。好处是末端效率可以直接用 |A2|²/(|A1|²|A2|²) 算不用再转换单位。坏处是耦合系数的表达式跟着变但这个变是透明的所有量都在 SI 下不容易错。然后是 RK4 积分器def shg_rk4(z_grid, A1_0, A2_0, dk, k1_cpl, k2_cpl): # 恒定步长 RK4每一步推进两个复振幅 A1, A2 A1_0, A2_0 for i in range(len(z_grid) - 1): h z_grid[i1] - z_grid[i] def f(z, y): a1, a2 y[0], y[1] return np.array([ -1j * k1_cpl * a2 * np.conj(a1) * np.exp(-1j * dk * z), -1j * k2_cpl * a1 * a1 * np.exp(1j * dk * z), ]) y np.array([A1, A2], dtypecomplex) k1v f(z_grid[i], y) k2v f(z_grid[i] h/2, y h*k1v/2) k3v f(z_grid[i] h/2, y h*k2v/2) k4v f(z_grid[i] h, y h*k3v) y y h/6 * (k1v 2*k2v 2*k3v k4v) A1, A2 y[0], y[1] return A1, A2主流程z_grid np.linspace(0, Lc, 4000) A1_end, A2_end shg_rk4(z_grid, A10, A20, dk, k1_cpl, k2_cpl) eta abs(A2_end)**2 / (abs(A1_end)**2 abs(A2_end)**2) print(f晶体末端倍频效率: {eta*100:.2f}%)逻辑说明右端项里 exp(−iΔkz) 这个快振荡项不能忽略它代表两波之间的相位滑移RK4 的优势就在处理这种带振荡相位的耦合方程。为什么不用欧拉法因为欧拉法的局部截断误差是 O(h²)对 3 mm 晶体、4000 步这种尺度误差累积会让效率曲线出现锯齿。RK4 的局部误差是 O(h⁵)同样步数下精度高得多。参数方面步数 4000 不是拍脑袋。基频在晶体里的波长约 0.64 μm3 mm 晶体约 4700 个波长步长取 0.75 μm 左右能保证每个波长周期内有接近一个采样点包络近似才站得住。如果你要扫大 Δk步长上限要用 π/|Δk| 来约束后面第 5 章细说。3.2 从小信号到泵浦耗尽效率曲线怎么看把这段代码跑起来你会看到效率随晶体长度并不是一路涨到 100% 就停下来。小信号区转换效率低于 10%效率近似正比于 L²这是解析解能覆盖的范围。继续涨到 50% 以上泵浦开始被明显消耗曲线斜率变缓。再往后可能出现效率回落这就是 back-conversion——倍频光又差频回去了能量倒流。很多人第一次看到效率回落以为自己 bug 了。不是是 Δk 不为零或者长度超过相干长度时三波之间的相位关系反转了。这个现象在真实实验里也看得到不是数值伪影。判断标准很简单把 Δk 设成零看效率是否单调饱和并趋近 100%如果这个行为对了说明程序逻辑没问题回落现象是物理不用修。另一个验证手段是能量守恒检验。每走一步|A1|² |A2|² 应该保持常数误差在 1e-6 量级以内。如果这条不成立优先怀疑耦合系数写错了。这是我每次跑完新参数必做的一件事比盯着效率曲线本身更有用。3.3 参数表哪些量要调调了会怎样参数示例取值作用调大之后的效果wl11064 nm基频波长改变光子能量和 k 矢量一般固定I0500 MW/cm²峰值功率密度越高转换越快但过高会触发其他非线性效应Lc3 mm晶体长度越长小信号效率越高但超过相干长度后回落dk由折射率算出相位失配量往零调效率显著上升是核心调参对象deff2e-12 m/V 量级有效非线性系数整体缩放效率不改变峰的位置和宽度步数4000积分精度太少出锯齿太多只增加计算量注意 dk 很少是自由参数。它由角度、温度、周期极化决定所以实际项目里你调的是晶体角度或者烘烤温度而不是直接改 dk。下一章就讲这两个最常用的调谐手段以及它们各自带来的副作用。4. 相位匹配调谐的落地参数角度、温度与走离效应4.1 角度调谐把 θ 扫一遍找到 I 类匹配峰与带宽BBO 是负单轴晶体I 类匹配配置是 oo→e也就是两个基频光子都是寻常光倍频光子是非寻常光。非寻常光折射率随传播方向与光轴夹角 θ 变化关系是n_e(θ) n_o·n_eo / sqrt(n_o²·sin²θ n_eo²·cos²θ)n_eo 是主轴方向上的非寻常光折射率。θ 一变k2 跟着变Δk 就能扫过零。所以角度调谐的本质是用几何关系补偿材料色散。下面这段代码把角度扫描接进上一章的求解器# 角度扫描在 20°~30° 范围内扫 BBO 的相位匹配角 n_o 1.655 n_eo 1.556 theta_deg np.linspace(20, 30, 2001) theta np.deg2rad(theta_deg) n2_theta (n_o * n_eo) / np.sqrt( n_o**2 * np.sin(theta)**2 n_eo**2 * np.cos(theta)**2 ) # 每个角度对应一个新的 dk重跑 RK4 eta_scan [] for n2_i in n2_theta: dk_i 2*np.pi*n2_i/wl2 - 2*(2*np.pi*n1/wl1) A1_end, A2_end shg_rk4( z_grid, A10, A20, dk_i, k1_cpl, k2_cpl ) eta_scan.append(abs(A2_end)**2 / (abs(A1_end)**2 abs(A2_end)**2)) eta_scan np.array(eta_scan) best_idx np.argmax(eta_scan) print(f最佳匹配角: {theta_deg[best_idx]:.2f}°) print(f峰值效率: {eta_scan[best_idx]*100:.2f}%)跑出来的峰值角度文献上 BBO 对 1064 nm 的 I 类匹配角一般在 22.8° 上下你的 Sellmeier 常数版本不同会差个零点几度这很正常。真正要注意的是峰的宽度半高全宽通常只有零点几度稍微偏一点效率就掉一半。这是因为 3 mm 晶体的相干长度有限Δk 稍大整个晶体长度内累积的相位滑移就超过 π。角度容忍度这个参数直接决定实验里要用什么精度的旋转台。如果仿真显示容忍度只有 0.1°那手动调节基本没戏得上带刻度的步进电机。这一条信息比效率绝对值更有工程价值。4.2 温度调谐与 90° 非临界相位匹配角度调谐有个天然副作用就是倍频光相对基频光会走离下一节细说。如果想躲开走离就得让传播方向垂直于光轴也就是 θ 90°这时非寻常光和寻常光的传播方向一致但 Δk 通常不为零。怎么办换温度。以 LiNbO₃ 为例它的折射率随温度变化很明显升温可以把 Δk 压回零这就是所谓 90° 非临界相位匹配。在脚本里做温度扫描比角度扫描更简单找一个温度依赖的 Sellmeier 方程算出每个温度下的 n1(T)、n2(T)然后套同一套 RK4 循环。折射率对温度的变化率大约在 10⁻⁵/K 到 10⁻⁴/K 量级换算成温度带宽通常也是几度的量级比角度容忍度宽松得多。另一个完全不同的思路是周期极化比如 PPLN。它不靠温度也不靠角度而是让晶体的自发极化方向周期翻转每过一个相干长度就把相位误差掰回来。极化周期 Λ 2π/Δk设计好周期之后可以在很宽的温度范围里做到非临界匹配同时用上晶体最大的有效非线性系数。仿真上周期极化只是在右端项里把 dk 换成 dk − 2π/Λ代码改动极小但物理效果完全不同。4.3 走离效应为什么临界匹配的晶体不能开大口径走离是角度调谐绕不开的代价。非寻常光在晶体里传播时能量传播方向坡印廷矢量和波矢方向不一致夹角叫走离角。BBO 的走离角大约 3°看起来不大但基频光和倍频光在横向上会慢慢分开重叠区域变小有效作用长度变短。数值上怎么处理一种常见做法是把晶体长度 L 换成有效长度 L_eff 光束直径 / tan(ρ)然后再跑同一套方程。比如光束直径 1 mmtan(3°) 约 0.052有效长度约 19 mm比晶体长走离可以忽略但光束直径 100 μm 时有效长度只有约 1.9 mm这时候再不建模走离效率会被高估不少。所以选型的时候有一个基本权衡匹配方式走离角角度容忍度温度容忍度典型晶体临界角度匹配有BBO 约 3°窄零点几度较宽BBO、LBO非临界温度匹配几乎为零宽窄可控LiNbO₃周期极化准相位匹配几乎为零宽由周期设计决定PPLN如果你的应用要求大光束、高效率、长晶体优先考虑非临界或者周期极化方案如果只是小口径、短晶体、追求简单BBO 角度调谐完全够用。这个取舍在做系统设计时比盲目上高非线性系数晶体重要得多。5. 非线性光学仿真常踩的 5 个坑与排查路径5.1 坑一效率曲线在中段出现锯齿或振荡现象eta 随晶体长度不是光滑上升而是出现周期性的锯齿甚至效率在某个位置跳变。原因步长不够。当 Δk 很大时右端项里的 exp(−iΔkz) 振荡非常快步长必须小于 π/|Δk|否则采样点落在振荡周期的错误位置上RK4 会把振荡混叠成低频误差。解决先算出 π/|Δk|然后把步长设成它的 1/5 到 1/10。比如 |Δk| 10⁶ /m 时π/|Δk| 约 3 μm步长取 0.3 μm对 3 mm 晶体就是 10000 步。不要嫌多这个场景下步数不是性能瓶颈。5.2 坑二泵浦耗尽后效率超过 100%现象|A2|²/(|A1|²|A2|²) 算出来 1.2能量无中生有。原因第二个方程的 1/2 因子丢了。倍频方程里 E1² 本身带有 2ω 分量系数很多人在推导时只记了形式把 ω2·deff/(n2·c) 直接写上漏掉分母里的 2。这个错误在小信号区效率 10%几乎看不出来一旦进入泵浦耗尽区就爆发。解决在 RK4 循环里加一行断言abs(A1)**2 abs(A2)**2 的极差小于 1e-6不满足就直接打印错误退出。这个校验应该成为所有耦合波方程的标配换晶体、换波长之后都先跑一遍。5.3 坑三角度扫描整个区间都找不到峰值现象把角度从 10° 扫到 40°效率全程接近 0或者峰值出现在明显不合理的位置。原因折射率数据用错了。最常见的是把文献里某一温度下的 Sellmeier 常数直接拿来算 1064 nm没注意那个常数只覆盖某个波长范围另一个常见问题是偏振配置搞反oo→e 写成了 oe→o计算出来的 Δk 永远不为零。解决先在脚本里打印 Δk 的扫描范围。如果 Δk 全程没有跨过零那不管扫多少个角度都不可能找到峰值。然后回去核对 Sellmeier 常数的波长适用范围和温度条件。我在换晶体参数时总是先用小信号解析解算一个粗略匹配角再上 RK4匹配角差太远就先修数据不急着看图。5.4 坑四单位制混乱结果差几个数量级现象同样的物理参数不同人跑出来的效率差一万倍而且谁都觉得自己没错。原因cm 和 m 混用是最常见的。功率密度有人用 W/cm²有人用 W/m²折射率倒是无量纲但波长一会儿 μm 一会儿 m最后 k 的量级差出 4 个次方。pm/V 和 m/V 的有效非线性系数也经常被直接套用。解决脚本头部单独做一个换算区所有物理量一次性转到 SI后续变量只引用不换算。比如功率密度统一写成 W/m²在参数区用 500e4 表示 500 MW/cm²但注释里写明原始工程单位。这个习惯能省掉大量排错时间属于血泪经验。5.5 坑五仿真结果和 SNLO 对不上现象同样的晶体、波长、功率密度自己的脚本算出来效率 60%SNLO 给 45%带宽也有偏差。原因SNLO 这类工程软件内部做了平面波近似之外的修正包括走离平均、孔径效应、聚焦参数、高斯光束的重叠积分。你的耦合波方程是理想平面波两者算的本来就不是同一个物理对象。这不叫 bug叫模型边界差异。解决先把场景设成小信号、薄晶体、低效率这个区间里解析解、你的数值解、SNLO 三者应该高度一致。对不上就先查自己的数据对得上再逐步提高功率、加长晶体观察两者分叉。分叉本身就是信息它告诉你什么时候平面波模型失效什么时候必须上更完整的模型。盲目追求和 SNLO 完全一致很多时候是在浪费时间调一个本来就不该一致的参数。6. 从倍频走向 OPO、OPA脚本扩展与 master 分支代码组织倍频的方程组只有两个扩展到 OPA 就是再加一个闲频光方程结构完全一样# OPA 三波耦合方程骨架信号光、闲频光、泵浦光 def rhs_opa(z, y, kappa_s, kappa_i, kappa_p, dk): As, Ai, Ap y[0], y[1], y[2] return np.array([ -1j * kappa_s * Ap * np.conj(Ai) * np.exp(-1j * dk * z), -1j * kappa_i * Ap * np.conj(As) * np.exp(-1j * dk * z), -1j * kappa_p * As * Ai * np.exp(1j * dk * z), ])这里的 dk kp − ks − ki符号定义和倍频场景一致。OPA 的物理行为和倍频完全不同泵浦功率超过阈值后信号光和闲频光会指数增长增益曲线对波长非常敏感所以它的仿真重点从效率计算变成了增益谱和阈值扫描。脚本骨架不变RK4 循环不用动只换右端项和初值这就是我推荐先做倍频的原因——后面全是复用。代码组织上我现在的习惯是把这类 Nonlinear-Optics 项目按 src、sims、results 三块分master 分支只保留可复现的版本。有人问 master 上的代码怎么切到 devgit checkout -b dev 一条命令的事但在科学计算项目里真正该做的不是频繁开分支而是每次跑批前在 config 文件里固定参数、跑完打 tag。tag 里记的是“这个结果对应哪一组晶体参数、哪个温度、哪个功率密度”比分支名可靠得多。我自己曾经吃过没打 tag 的亏同一份代码改了 Sellmeier 常数忘了记录两周后对着旧结果怎么都对不上。现在每个仿真任务新建一个 config 文件task 名和 tag 同步后悔药算是吃上了。现在拿到一个非线性频率变换题目我第一件事不是找求解器而是先查这批晶体的折射率数据对不对、匹配角落在哪个区间、走离影响大不大这些确认了再写 RK4。这个习惯救了我好几次希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
