简介这是一份关于激光器被动调Q速率方程的MATLAB模拟资源面向激光物理、光电子技术方向的初学者或研究人员用于理解饱和吸收体实现被动Q开关的动态过程。压缩包内共2个文件均为.m脚本分别对应速率方程求解与调Q过程仿真整体仅2KB代码精简、便于快速运行和二次修改。目前已有639人学习下载。通过该资源读者可结合速率方程模型观察粒子数反转、腔内损耗及脉冲形成过程掌握可饱和吸收体在高光强下吸收率降低的机制进而为激光器脉冲输出设计与参数优化提供参考。1. Qswitch 系列做到第二篇该把被动调Q单独拿出来讲了被动调Q和主动调Q最大的区别是腔内没有电光、声光调制器也没有高压驱动只在谐振腔里放一块可饱和吸收体——Cr:YAG 晶体、SESAM、染料片都算。低光强时吸收体不透明腔损耗高能量在增益介质里憋着腔内光强一旦超过饱和光强吸收体在几个纳秒内被漂白Q 值瞬间抬高存储的能量以短脉冲放出来。结果就是脉宽、峰值功率、重复频率这三样东西几乎全部由速率方程里的截面和损耗项决定。这篇把被动调Q的速率方程从建立、数值求解讲到参数优化和校核适合要写仿真、要压低脉宽或者给现有谐振腔调饱和吸收体参数的激光工程师。2. 被动调Q速率方程建模饱和吸收体的损耗项与三方程组2.1 可饱和吸收与漂白饱和光强 Is 怎么算可饱和吸收saturable absorption是吸收系数随入射光强增大而下降的现象。对二能级模型吸收系数写成 α(I) α0/(1 I/Is)其中 α0 是小信号吸收系数Is 是饱和光强。饱和的物理机制是强光把基态粒子大量抽到激发态基态布居数下降吸收随之减弱也就是常说的“漂白”。饱和光强的表达式是Is hν / (σ_gs · τ_a)其中 σ_gs 是基态吸收截面τ_a 是基态恢复时间。对 1064 nm 的 Cr:YAGσ_gs 取 4.3e-18 cm²、τ_a 取 3.5 μs算出来 Is 在 MW/cm² 量级。仿真里更顺手的是换成饱和光子数密度phi_sat 1.0 / (SIGMA_GS * C * TAU_A) # 光子数密度饱和值cm^-3 I_sat phi_sat * H * NU * C # 对应的光强W/cm^2phi_sat的量级大约在 2.2e12 cm⁻³代码里用SIGMA_GS * TAU_A的乘积形式而不是直接写数字是为了换吸收体材料时只改两个参数就行。这段计算的意义在于被动调Q的“被动”二字就体现在吸收体透过率随腔内光子数密度自动变化速率方程里其他项可以照搬主动调Q的写法唯独损耗项要换成可饱和吸收体的动态损耗。把可饱和吸收现象和可饱和吸收体器件分开理解后面的建模才不会混。2.2 被动调Q速率方程三方程组光子数、反转粒子数、基态布居完整的被动调Q速率方程一般写成三个耦合的一阶常微分方程dφ/dt (φ/tr)·[ 2σnl − 2σ_gs·n_a·l_a − 2σ_es·(N_a−n_a)·l_a − (L_par − ln R) ]dn/dt −γσcφn − n/τ_cdn_a/dt −σ_gs·c·φ·n_a (N_a−n_a)/τ_a第一个方程是光子数密度 φ 的演化方括号里是单次往返的净增益减净损耗。2σnl 是双程增益——光子每往返一次穿过增益介质两次所以增益和吸收项的“2”都来自双程。σ 是增益介质受激发射截面n 是反转粒子数密度l 是晶体长度。后面四项是损耗2σ_gs·n_a·l_a 是吸收体基态吸收在小信号时最大这正是憋能量的来源2σ_es·(N_a−n_a)·l_a 是激发态吸收造成的残余损耗吸收体完全漂白后也消不掉L_par 是散射、衍射等寄生损耗−ln R 是输出镜透过损耗。tr 2L_cav/c 是光子往返时间L_cav 必须用光学腔长不能拿几何腔长直接算晶体折射率要折算进去。第二个方程是反转粒子数演化−γσcφn 是受激辐射消耗γ 是能级简并因子四能级取 1三能级取 2−n/τ_c 是自发辐射弛豫。脉冲只有几个纳秒而 τ_c 是几百微秒所以脉冲期间弛豫项可以忽略——这是单脉冲仿真能成立的前提也是后面第 4 章讨论重复频率时要把泵浦项重新加回来的原因。第三个方程是吸收体基态布居 n_a 的演化−σ_gs·c·φ·n_a 是漂白项(N_a−n_a)/τ_a 是基态恢复项N_a 是吸收体总掺杂密度。脉冲期间恢复项相对漂白项小到可以忽略但做脉冲序列仿真时不能省否则吸收体永远回不到初始状态。2.3 快吸收体与慢吸收体Cr4:YAG 该用哪个模型选模型看吸收体恢复时间 τ_a 和输出脉宽 τ_p 的相对大小这是一个经常被忽略的前提。类型代表器件τ_a 与 τ_p 关系建模方式快吸收体染料、SESAMτ_a 远小于 τ_pn_a 取准平衡方程组降为两维慢吸收体Cr4:YAGτ_a 远大于 τ_p保留三方程脉冲期间恢复项可忽略Cr4:YAG 的 τ_a 约 3.5 μs输出脉宽通常在 ns 量级差三个数量级按快吸收体近似会有明显偏差直接解三方程最稳。慢吸收体模型还有一个实际推论脉冲结束后基态恢复需要几个 τ_a如果泵浦过强导致脉冲间隔小于恢复时间吸收体来不及复位会出现多脉冲甚至完全不调Q。这个坑在实验上比仿真更容易遇到参数上表现为重复频率不是随泵浦线性上升而是突然掉进不稳定区。3. Python 数值求解被动调Q速率方程最小可复现代码3.1 单位制与参数表cm-s 体系下的截面和寿命激光仿真最容易出错的不是方程是单位。全文统一用 cm-s 体系长度用 cm面积用 cm²密度用 cm⁻³截面用 cm²时间用 s。波长 1064 nm 写作 1064e-7 cm光速 3.0e10 cm/s。下面是一组能直接出单脉冲的参数参数符号数值单位Nd:YAG 受激发射截面σ2.8e-19cm²增益介质长度l0.5cm上能级寿命τ_c230e-6s光学腔长L_cav10.0cm输出镜反射率R0.71寄生往返损耗L_par0.051Cr:YAG 基态吸收截面σ_gs4.3e-18cm²Cr:YAG 激发态吸收截面σ_es8.2e-19cm²吸收体恢复时间τ_a3.5e-6s吸收体初始透过率T00.801模场面积A1.0e-2cm²T0 0.8 意味着单程小信号吸收 −ln(0.8) ≈ 0.223往返可饱和损耗约 0.446与输出耦合损耗 0.357 同量级是工程上“调制深度与输出耦合匹配”的常见配置。这个搭配下脉冲能稳定建立吸收体也能被充分漂白。3.2 用 solve_ivp 求解被动调Q速率方程的完整代码import numpy as np from scipy.integrate import solve_ivp C 3.0e10 # 光速cm/s H 6.626e-34 # 普朗克常数J·s LAMBDA 1064e-7 # 波长cm1064 nm NU C / LAMBDA # 频率Hz # 增益介质 Nd:YAG SIGMA_G 2.8e-19 # 受激发射截面cm^2 L_G 0.5 # 晶体长度cm TAU_C 230e-6 # 上能级寿命s GAMMA 1.0 # 四能级系统取 1 # 谐振腔 L_CAV 10.0 # 光学腔长cm TR 2.0 * L_CAV / C # 光子往返时间s R_MIR 0.7 # 输出镜反射率 L_PAR 0.05 # 往返寄生损耗 # 可饱和吸收体 Cr4:YAG SIGMA_GS 4.3e-18 # 基态吸收截面cm^2 SIGMA_ES 8.2e-19 # 激发态吸收截面cm^2 TAU_A 3.5e-6 # 基态恢复时间s T0 0.80 # 小信号单程透过率 N_A_L_A -np.log(T0) / SIGMA_GS # 等效面掺杂密度cm^-2 L_A 0.2 # 吸收体厚度cm N_A_TOT N_A_L_A / L_A # 掺杂密度cm^-3 BEAM_AREA 1.0e-2 # 腔内模场面积cm^2 def qswitch_rhs(t, y): phi, n_inv, n_a y gain 2.0 * SIGMA_G * n_inv * L_G # 双程增益 loss_gs 2.0 * SIGMA_GS * n_a * L_A # 吸收体基态损耗 loss_es 2.0 * SIGMA_ES * (N_A_TOT - n_a) * L_A # 激发态残余损耗 loss_total loss_gs loss_es L_PAR - np.log(R_MIR) dphi phi / TR * (gain - loss_total) dn_inv -GAMMA * SIGMA_G * C * phi * n_inv - n_inv / TAU_C dn_a -SIGMA_GS * C * phi * n_a (N_A_TOT - n_a) / TAU_A return [dphi, dn_inv, dn_a] # 初始状态吸收体全在基态反转粒子数取阈值之上的裕度 thresh_inv (2.0 * SIGMA_GS * N_A_L_A L_PAR - np.log(R_MIR)) \ / (2.0 * SIGMA_G * L_G) ALPHA 1.3 n_inv0 ALPHA * thresh_inv sol solve_ivp( qswitch_rhs, (0.0, 1.0e-6), [1.0, n_inv0, N_A_TOT], methodLSODA, rtol1e-8, atol1e-12, max_step1e-10 )thresh_inv是初始高损耗态的阈值反转粒子数分母 2σl 来自双程增益。ALPHA 1.3 表示初始反转比阈值高 30%这是被动调Q单脉冲比较稳的区间。积分方法用 LSODA它在非刚阶段自动用显式格式、进入脉冲上升沿后切到刚求解器比固定步长 RK4 省心得多。max_step1e-10即 0.1 ns保证脉冲前沿至少有几个采样点但这同时意味着 1 μs 窗口最多一万步计算量完全可以接受。3.3 从解中提取脉冲能量、脉宽与峰值功率t sol.t phi sol.y[0] P_out -np.log(R_MIR) * H * NU * BEAM_AREA * L_CAV * phi / TR E_pulse np.trapz(P_out, t) # 脉冲能量J P_peak P_out.max() # 峰值功率W tau_p E_pulse / P_peak # 等效脉宽s print(f脉冲能量 : {E_pulse * 1e3:.3f} mJ) print(f峰值功率 : {P_peak / 1e3:.1f} kW) print(f等效脉宽 : {tau_p * 1e9:.2f} ns)P_out的表达式是腔内光子能量密度乘以输出镜逃逸速率。hν·φ·A·L_cav 是腔内总光子能量除以 tr 得到能量流逝速率再乘以 −ln R 就是真正从输出镜出去的功率。用 −ln R 而不是 (1−R) 是因为反射率越偏离 1二者的差异越大工程上习惯用前者当等效损耗。等效脉宽用能量除以峰值功率定义对单脉冲比较稳健比半高宽好算也不怕波形不对称。这个参数组下的典型量级如下具体数值会随容差设置略有浮动指标参考量级脉冲能量12 mJ峰值功率300500 kW等效脉宽48 ns腔内峰值光子密度约 1e16 cm⁻³3.4 初值与积分窗口的选择依据φ0 取 1.0 cm⁻³对应自发辐射噪声的量级只影响脉冲建立时间不影响脉冲能量和脉宽。窗口取 1 μs以 ALPHA 1.3 算光子从种子涨到饱和光子数密度大约需要 7080 ns脉冲本身几纳秒1 μs 足够让脉冲尾部完全收干净。一个常见错误是窗口只开到脉冲峰值附近导致 n_f 没降到最终值能量少算判断标准是把窗口加倍若脉冲能量变化小于 0.1%说明尾部截断没有影响结果。4. 被动调Q参数扫描饱和吸收体透过率、模场比与重复频率4.1 初始透过率 T0 与调制深度的折算初始透过率 T0 exp(−σ_gs·N_a·l_a) 是最容易直接测量的参数也是买 Cr:YAG 时厂家唯一会写在规格书上的东西。它对应的小信号损耗和调制深度需要自己换算T0σ_gs·N_a·l_a往返小信号损耗 2σ_gs·N_a·l_a往返调制深度 2·N_a·l_a·(σ_gs−σ_es)0.900.1050.210.170.800.2230.450.360.700.3570.710.580.600.5111.020.83往返调制深度等于往返小信号损耗减去全漂白后的残余损耗它才是真正参与调Q的“损耗变化量”。T0 越低调制深度越大越容易压短脉宽但代价是阈值反转粒子数升高、吸收体上的能量密度变大。工程经验是 T0 不要低于 0.6否则镀膜损伤和热效应会在高重频下先于调Q性能出问题。4.2 反转粒子数裕度 α 决定脉宽与稳定性α n_i/n_th 是泵浦给出来的自由参数也是被动调Q里最敏感的一个。α 越接近 1脉冲建立越慢脉宽越宽α 提高脉宽变窄、峰值功率升高但超过一定值后容易出现多脉冲和不稳定。被动调Q的稳定单脉冲区间大致在 1.25 ≤ α ≤ 1.6低于下限会退化成弛豫振荡高于上限会出现双脉冲。调整 α 的工程手段是调泵浦能量或泵浦脉宽而不是换吸收体。对给定 α最终反转粒子数 n_f 由超越方程 2σl·(n_i−n_f) L_f·ln(n_i/n_f) 决定其中 L_f 是吸收体完全漂白后的往返损耗。这个式子直接说明能量提取效率取决于 n_i/n_f 的对数而 n_f 与 L_f 成正比。所以降低残余损耗比如选 σ_es 更小的吸收体比单纯增加 n_i 更能提高效率。4.3 增益模与吸收体模的面积比 C 怎么修正上一章代码假设饱和吸收体处的光子数密度和增益介质处相同。实际布腔时吸收体往往放在焦点附近模场面积 A_a 和增益介质处 A_g 不相等。定义 C A_g/A_a代入时只改吸收体的漂白项C_AREA 1.5 # A_g / A_a吸收体处模更小则大于 1 phi_a phi * C_AREA # 吸收体实际感受到的光子数密度 dn_a -SIGMA_GS * C * phi_a * n_a (N_A_TOT - n_a) / TAU_A光子方程里的损耗项不变因为损耗是按往返路径积分出来的只和 n_a 与 l_a 有关。C 大于 1 意味着吸收体处光强更集中同样的腔内能量就能把吸收体漂白得更彻底脉宽更短、建立更快。但 C 不能无限大吸收体表面峰值光强与损伤阈值挂钩一般取 12 之间。腔内光强 I hν·c·φ可以直接换算成功率密度和损伤风险。4.4 重复频率的解析估算与泵浦功率单脉冲仿真不包含泵浦项重复频率要用恢复过程单独估算。脉冲结束后反转粒子数从 n_f 向泵浦平衡值 n_p 指数恢复当恢复到 n_i 时触发下一个脉冲所以周期近似为T_rep τ_c · ln[(n_p − n_f) / (n_p − n_i)]用前面的参数τ_c 230 μsn_p 取 2 倍阈值时算出来 T_rep 约 160 μs对应约 6 kHz。这个公式的条件是 n_p 大于 n_i且吸收体在两次脉冲之间已经基本恢复后者要求 T_rep 远大于几个 τ_a代入 Cr:YAG 的 3.5 μs 是成立的。泵浦进一步提高时n_p 增大使周期对数下降但过高会进入多脉冲区此时吸收体恢复不充分实验上表现为脉冲序列抖动变大——这个区域速率方程也能复现只要把 n_a 的初值改为上次脉冲结束时的残余值就能判断吸收体复位是否到位。5. 被动调Q速率方程解的校核与数值排错5.1 用解析能量式交叉验证仿真跑通后先别急着信用一个与数值解独立的结果做交叉验证。退化分析给出脉冲能量的解析式E_out (hν·A / (2γσ)) · ln(1/R) · ln(n_i/n_f)其中 n_f 用超越方程 2σl·(n_i−n_f) L_f·ln(n_i/n_f) 求根注意 L_f 必须用吸收体漂白后的最终损耗from scipy.optimize import brentq L_final L_PAR - np.log(R_MIR) 2.0 * SIGMA_ES * N_A_TOT * L_A f_root lambda nf: 2.0 * SIGMA_G * L_G * (n_inv0 - nf) \ - L_final * np.log(n_inv0 / nf) n_f_theory brentq(f_root, n_inv0 * 1e-3, n_inv0 * 0.99) E_theory H * NU * BEAM_AREA / (2.0 * GAMMA * SIGMA_G) \ * (-np.log(R_MIR)) * np.log(n_inv0 / n_f_theory)n_f_theory 应该与仿真末尾的 n_inv 接近偏差超过 5% 时优先怀疑积分窗口截断或吸收体没有完全漂白而不是解析式的问题。E_theory 与数值积分得到的 E_pulse 的偏差一般在几个百分点以内这就把数值解、解析解和物理模型三者互相锁住了。5.2 数值上不稳时的四个信号与对策第一φ 出现负值或小幅振荡通常是 tol 太松或输出点跨越了脉冲前沿把 rtol 降到 1e-10、max_step 降到 tr/5 即可。第二脉宽恰好等于 max_step说明峰值被步长削平脉冲峰值功率偏低需要进一步缩小步长。第三脉冲尾部出现非物理的二次鼓包多半是窗口末端反射粒子数方程里恢复项开始起作用检查是否错把单脉冲窗口当成了连续序列仿真。第四LSODA 报错或步长收敛极慢先查单位——最常见的错误是把波长写成 1064e-9 m 混进 cm 体系频率差了 100 倍所有饱和项全部失真。5.3 从实测 T0 反推掺杂浓度注意表面反射偏差拿到一片 Cr:YAG实测单程透过率后掺杂密度反推公式是 N_a −ln(T0) / (σ_gs · d)d 为厚度。常用厚度下的换算实测 T0d 1 mm 时 N_ad 2 mm 时 N_a0.902.5e17 cm⁻³1.2e17 cm⁻³0.805.2e17 cm⁻³2.6e17 cm⁻³0.708.3e17 cm⁻³4.1e17 cm⁻³这里有一个容易踩的偏差未镀增透膜的 Cr:YAG 两个表面各反射约 7%直接用分光光度计测的透过率会把表面反射算进吸收导致反推的掺杂密度系统性偏高。正确的做法是用镀增透膜样品或者先测两个表面的反射基线再扣除否则 N_a 会差 10% 以上代入速率方程后阈值反转粒子数跟着偏仿真结果和实验对不上时先查这一项。本文还有配套的精品资源点击获取
