简介ASE.zip 面向激光物理与光学工程方向的研究者和高年级学生聚焦激光放大与激光放大自发辐射ASE两大核心主题并以氮气离子型空气激光的自发放大辐射过程为具体应用场景开展数值计算与仿真。压缩包体积约3KB仅含1个m文件即 MATLAB 格式脚本借助数值方法模拟辐射强度随时间与空间的变化可用于分析激活介质浓度、能量转移速率及辐射衰减系数等关键参数直观呈现粒子在高能级与低能级之间跃迁产生的增益行为。该小巧仿真脚本目前已有 323 人浏览学习。资源围绕空气激光 ASE 的典型建模思路展开适合作为理解能级跃迁、自发辐射干涉增强及激光增益机制的入门样例同时便于在 MATLAB 环境中快速运行并允许调整泵浦强度、介质浓度等输入条件通过输出结果反推受激放大条件为改善激光器单色性、方向性或抑制不同强度ASE干扰提供依据。1. 激光放大与放大自发辐射ASE.m 到底在算什么做气体激光或者超快光丝实验的人几乎都绕不过一个词ASEAmplified Spontaneous Emission放大自发辐射。你在文献里看到别人把氮气在强场电离后发出的 337 nm 或者 391 nm 信号叫做“氮离子空气激光”其实那里面很大一部分贡献就是 ASE——没有谐振腔靠单程增益把自发辐射放大到可探测的强度。这份ASE.zip里的ASE.m是一个 MATLAB 脚本用速率方程模型模拟氮气离子型空气激光的 ASE 强度随泵浦条件的变化。它能算的是给定氮气密度、泵浦速率、上能级寿命和增益截面输出辐射强度随时间的演化曲线。适合正在做空气激光实验、需要预先估算信号量级的研究生或者想理解 ASE 和真正激光放大本质区别的入门者。这篇文章我按自己拆这个脚本的思路写先讲物理模型再逐段拆代码最后把调参时容易翻车的地方列出来。2. 从速率方程到 ASE 强度先把物理模型立住2.1 为什么用速率方程而不是 Maxwell-Bloch 方程ASE.m的核心不是求解 Maxwell 方程组而是用速率方程近似。原因是氮气离子空气激光的增益介质长度通常在厘米到米量级泵浦是飞秒或者皮秒脉冲电离产生的粒子数反转建立的时间和上能级寿命纳秒量级远小于光在增益介质中的渡越时间。这种情况下光场与介质的相干效应不是主导因素用粒子数密度的时间演化来描述增益过程就足够。速率方程的基本形式是% 上能级粒子数密度 N2 的时间演化 dN2dt Pump * N1 - N2 / tau_upper - sigma_ase * N2 * I_ase;这里的Pump是泵浦速率单位 1/sN1是基态粒子数密度tau_upper是上能级寿命sigma_ase是受激辐射截面I_ase是 ASE 光强。方程的含义很直接泵浦把粒子从基态抽运到上能级自发辐射和受激辐射共同把上能级粒子数消耗掉。需要注意的是这个方程里没有包含谐振腔反馈项。这是 ASE 和激光振荡的本质区别——光子在增益介质里只走一趟靠单程增益放大没有模式选择所以输出光谱宽、方向性差。如果你之后要对比“真正的激光”和“ASE”这个方程就是分水岭。2.2 增益系数的计算小信号增益和大信号增益脚本里有一段计算增益系数的代码这是整个模拟的关键。小信号增益系数 g0 由粒子数反转密度和受激辐射截面决定% 小信号增益系数 g0 sigma_ase * (N2 - N1 * exp(-Delta_E / (kB * T))); % 其中 Delta_E 是上下能级能量差kB 是玻尔兹曼常数这个表达式的物理含义是只有上下能级粒子数满足玻尔兹曼分布偏离时才有净增益。温度越高下能级热布居越多增益越小。我实际跑脚本时习惯先固定sigma_ase和tau_upper然后扫N2氮气密度看 g0 的变化趋势——这样能快速判断实验上该加压还是减压。大信号增益则是光强增长到不可忽略时的表现。脚本里用了一个迭代公式% 大信号增益考虑增益饱和 G g0 * L / (1 I_ase / I_sat); % I_sat 是饱和光强L 是增益介质长度这里的I_sat是关键参数它由hbar * omega / (sigma_ase * tau_upper)给出。饱和光强决定了 ASE 能放大到什么程度——超过它之后增益被粒子数消耗拖住输出光强不再线性增长。我在实验里把这个参数当作“天花板”如果模拟出来的 ASE 强度已经接近饱和光强那实验上再增大泵浦能量也不会带来明显收益了。3. 拆解 ASE.m从参数初始化到数值求解3.1 参数初始化哪些量必须由实验定哪些可以查文献ASE.m开头是一段参数定义我建议你不要直接跑先逐行确认每个数值的物理依据。这里贴一段典型初始化代码%% 参数初始化 % 氮气分子数密度1/cm^3常压 300K 下约 2.4e19 N_total 2.4e19; % 上能级寿命sN2 B2Sigma_u 态约 60 ns tau_upper 60e-9; % 受激辐射截面cm^2取值在 1e-16 到 1e-14 之间 sigma_ase 1e-15; % 泵浦速率1/s由飞秒激光电离速率估算 Pump 1e9; % 增益介质长度cm L 10; % 时间窗口s t_end 100e-9;这里最容易出错的是sigma_ase。不同文献里氮气离子 391 nm 跃迁的受激辐射截面数据差别很大从 1e-17 到 1e-14 cm² 都有报道。原因是截面和气压、温度、碰撞展宽有关你在模拟时应该把它当可调参数而不是固定常数。我的习惯是先用 1e-15 跑一遍看输出强度量级是否合理再根据实验数据反推修正。Pump这个参数也容易误设。它不是你加在介质上的光强而是电离产生的自由电子通过碰撞激发或直接激发把氮分子离子打到上能级的速率。飞秒光丝实验中泵浦速率沿光丝是不均匀的中间高两边低但零维模型只能给一个平均值。如果你需要更精确的空间分布就得升级成多维模型那不是这个脚本能干的活。3.2 数值求解用 ode45 还是自己写递推ASE.m里求解速率方程用的方法是 MATLAB 自带的 ode45。代码大概是这样%% 定义速率方程 ODE function dNdt ase_rate_eq(t, N, params) % N(1) 是上能级粒子数密度 % N(2) 是 ASE 光子数密度 Pump params.Pump; tau_upper params.tau_upper; sigma_ase params.sigma_ase; c 3e10; % 光速cm/s omega params.omega; dNdt zeros(2, 1); dNdt(1) Pump * (params.N_total - N(1)) - N(1) / tau_upper ... - sigma_ase * c * N(2) * N(1); dNdt(2) sigma_ase * c * N(2) * N(1) N(1) / tau_upper * params.beta ... - N(2) / params.tau_cavity; end这段代码里有两个容易忽略的细节。第一N(2)是光子数密度而非光强两者之间差一个光子能量hbar * omega。输出后处理时如果你直接把N(2)当光强用数值会差好几个数量级。第二方程里 N(1) / tau_upper * params.beta这一项是自发辐射耦合进 ASE 模式的比例beta通常取 1e-5 到 1e-3。这个参数决定了 ASE 的种子来源——没有它初始光子数为零时方程解恒为零永远算不出 ASE。我见过不少人第一次跑脚本得到全零结果就是忘了设置初始种子光子数。ode45 是自适应步长算法好处是精度高坏处是在快慢时间尺度差距大的刚性问题上可能很慢。如果tau_upper是纳秒量级而泵浦是飞秒量级建议换成 ode15s刚性求解器否则跑一次模拟要等很久。这个经验在我自己算氩气ASE时非常有用。4. 避坑指南ASE 模拟中五个我踩过的坑4.1 输出光强莫名其妙变成负值现象ode45 跑出来的N(1)出现负值导致后续增益计算全部出错。原因速率方程是刚性系统ode45 对这种快慢时间尺度相差悬殊的方程容易产生振荡步长控制不住时就会让粒子数密度越过零变成负值。解决换 ode15s或者在方程里加非负约束。MATLAB 里可以用odeset(NonNegative, [1, 2])强制非负。从那以后我所有粒子数密度模拟都默认加这个选项。4.2 增益系数算出来比实验大 100 倍现象小信号增益系数 g0 算出来是 10 cm⁻¹实验上实际只有 0.1 cm⁻¹。原因没有考虑碰撞展宽和 Doppler 展宽对受激辐射截面的影响。你用的sigma_ase可能是低温或低压下的值常温常压下谱线展宽会把截面拉低一到两个数量级。解决把sigma_ase换成有效截面sigma_eff sigma_0 * (gamma_line / (gamma_line gamma_coll))其中gamma_line是自然展宽gamma_coll是碰撞展宽。这个修正我在跑常压空气时必加不加必翻车。4.3 改变泵浦速率时输出不连续跳变现象扫 Pump 参数时ASE 强度曲线不光滑有突变。原因可能触发了数值不稳定性也可能是beta参数设置太小导致初始种子光子数不足数值噪声在低增益区被放大。解决检查beta是否在合理范围1e-5 到 1e-3另外把 ode45 的相对误差RelTol设到 1e-8 而不是默认的 1e-3。有时候“物理现象”其实是数值假象这个区别得靠收敛性检查判断。4.4 时间窗口太短ASE 还没发展到饱和就截断了现象模拟结果里 ASE 强度一直在涨没有平台期。原因t_end设置太短。ASE 从开始到饱和需要几个上能级寿命的时间如果你只跑到 10 ns 而寿命是 60 ns永远看不到饱和。解决把t_end设为10 * tau_upper然后检查最终时刻的dNdt(2)是否接近零。如果还在涨继续加长时间窗口。4.5 用同一个脚本算氩气结果完全不匹配现象把氮气参数换成氩气模拟结果和实验对不上。原因氩气离子 ASE 的机制和氮气不同氩气主要是 3p 到 4p 能级跃迁而且和泵浦脉冲的 carrier-envelope phase 有依赖。零维速率方程模型根本描述不了这个物理过程。解决认清这个模型的边界——它适合氮气、氢气和简单分子体系不适合强场电离机制复杂的惰性气体。遇到这类问题直接放弃这个脚本换更复杂的模型。5. 从模拟到实验ASE 参数标定和数据解读5.1 怎么用实验数据反过来定标参数模拟的最终目的是指导实验或者解释实验。我的标准流程是先用ASE.m扫参数看趋势然后在实验上测一个点比如某个气压下的 ASE 能量反推sigma_ase和tau_upper的实际有效值再把这个标定后的参数代回模拟算其他条件。具体做法是把脚本里的核心计算包成一个函数输入是气压、泵浦速率、长度输出是 ASE 能量密度。然后用 MATLAB 的lsqcurvefit或者fminsearch去拟合实验数据。这里贴一段简单的拟合框架%% 反演标定 sigma_ase % 实验数据p_vals 是气压P_exp 是测量的 ASE 脉冲能量 p_vals [0.1, 0.5, 1.0, 2.0]; % 单位 atm P_exp [1e-6, 5e-6, 8e-6, 1e-5]; % 单位 J % 定义拟合函数 fit_func (params, p) ase_predict(p, params.sigma_ase, params.tau_upper); % 初始猜测 params0 struct(sigma_ase, 1e-15, tau_upper, 60e-9); % 非线性最小二乘拟合 options optimset(Display, iter); [params_opt, resnorm] lsqcurvefit(fit_func, params0, p_vals, P_exp, ... [1e-17, 10e-9], [1e-13, 100e-9], options);拟合完要检查resnorm——如果残差大说明模型本身不足以描述你的实验这时候不要强行调参掩盖问题。我遇到过拟合结果很好但参数是负值的情况那就说明模型有 bug不是参数的问题。5.2 从 ASE 光谱宽度反推增益介质均匀性单个速率方程脚本只能输出强度随时间的变化不能直接给出光谱。但如果把不同时间的强度演化做傅里叶变换可以得到一个粗略的线宽估计。更实用的做法是用相位噪声模型后处理% 从 ASE 时间信号估计光谱宽度简化版 % I_t 是 ASE 强度时间序列dt 是采样间隔 I_fft fftshift(fft(I_t)); f_axis (-length(I_t)/2 : length(I_t)/2-1) / (length(I_t) * dt); % 取半高全宽作为估计线宽 [peak_val, peak_idx] max(abs(I_fft).^2); half_max peak_val / 2; % 找到左右两侧第一次低于半峰值的位置 left_idx find(abs(I_fft(1:peak_idx)).^2 half_max, 1, last); right_idx peak_idx find(abs(I_fft(peak_idx:end)).^2 half_max, 1, first); linewidth f_axis(right_idx) - f_axis(left_idx);这个方法的物理基础是ASE 的线宽和增益介质的均匀性、碰撞展宽有关如果实验上测到的线宽远大于模拟值说明增益介质里可能存在温度梯度或者湍流导致的折射率不均匀光在不同位置获得的增益不一致展宽了光谱。这个脚本没法直接模拟空间不均匀性但它能帮你建立一个“正常线宽”的预期——超出预期就是实验系统里有额外问题。5.3 判断 ASE 是否进入饱和区输出能量和增益长度的关系模拟里最有用的诊断是看输出能量随增益介质长度的变化关系。在小信号区ASE 能量随长度指数增长在深饱和区能量随长度线性增长。这个转折点就是实验中最佳介质长度的选择依据。在ASE.m里跑一组不同 L 的模拟记录稳态输出能量 P_out然后取对数画图——如果曲线明显弯曲说明在长介质里已经进入饱和区。实验中遇到这种情况加长介质不会带来收益反而增加损耗和成本。这时候应该优化泵浦的时空分布而不是无脑加长度。6. 一个进阶技巧把 ASE.m 改造成泵浦-探测参数扫描器脚本写到这里能算的东西已经够了但还有一个小技巧可以大幅提升它的实用价值把泵浦速率和时间延迟做成二维扫描模拟泵浦-探测实验中的时间分辨 ASE 信号。这个做法在飞秒光丝空气激光的 pump-probe 实验里很常用。改造思路是把Pump从常数改成时间函数。真实实验中泵浦脉冲是高斯形状探测脉冲的到达时间可以调节。代码改造如下%% 时间相关的泵浦速率泵浦-探测场景 function Pump_t pump_profile(t, t_delay, tau_pump, Pump_peak) % t_delay: 探测脉冲相对泵浦的延迟 % tau_pump: 泵浦脉冲宽度FWHM 近似 % 高斯泵浦脉冲峰值在 t_delay 时刻 Pump_t Pump_peak * exp(-2.77 * (t - t_delay).^2 / tau_pump^2); end然后用它替换原来固定的Pump扫不同t_delay值记录每个延迟下的 ASE 积分强度就能得到一张“ASE 强度 vs 泵浦-探测延迟”的曲线。这张曲线能告诉你增益介质里粒子数恢复的时间——如果实验上有相应的 pump-probe 数据正好可以和模拟对比校验你设定的能级寿命是否正确。我还习惯在这个扫描基础上加一个“泵浦-探测能级填充模型”的参考线。当延迟时间远大于上能级寿命时ASE 强度应该回到单脉冲泵浦的基线当延迟时间在零附近时两束光叠加泵浦ASE 强度会增强。如果模拟里出现负延迟时 ASE 反而比正延迟强说明你设置的泵浦速率时间分布有问题或者存在能级记忆效应。有一个细节务必检查时间延迟扫描的步长。上能级寿命如果是 60 ns扫描步长选 1 ns 会造成曲线过稀看不出单指数衰减的形状选 0.01 ns 又会让计算量暴增。我一般先粗扫步长 0.5 倍寿命定位特征时间再细化局部区域两种步长各跑一遍确认结果一致才放心。从那以后我每次改ASE.m的时间逻辑都会先做一个步长收敛性测试——这个习惯帮我避免了好几次“模拟结果随步长漂移”的乌龙。希望这个脚本和这些调参经验帮到你至少让你少走一些我当年走过的弯路。本文还有配套的精品资源点击获取