简介：源自超短脉冲仿真需求，该资源面向光纤通信与光子学领域的科研人员及高年级学生，聚焦飞秒与皮秒激光高斯脉冲在光纤中的色散行为，兼顾啁啾效应与瞬态吸收光谱分析。压缩包仅含2个文件，包括用于无啁啾高斯脉冲传播仿真的Matlab脚本，以及系统阐释啁啾与色散概念的Word文档，整体包体约1.49MB。已有837人浏览学习，适合需要快速掌握脉冲展宽成因、正/负啁啾调节机制，或希望借助仿真工具验证理论结果的读者。通过运行脚本，可直观观察色散导致的脉冲展宽现象，并结合文档理解啁啾参数对脉冲压缩或扩大的影响；还可进一步延伸到自相位调制、交叉相位调制等非线性效应以及瞬态吸收光谱动力学分析的建模场景。整体上提供了从基础概念到仿真实验的完整闭环，有助于缩短光学超快现象的学习曲线。
1. 飞秒及皮秒激光进光纤前，先用MATLAB把色散和瞬态吸收仿真做对
飞秒脉冲进光纤后会发生什么，和你在示波器上看到的光谱完全是两回事：一个100 fs的高斯脉冲走完一米普通单模光纤，脉宽可能已经从0.1 ps展宽到几皮秒，峰值功率掉一个量级，而光谱几乎没变。这个“线性色散就能做到的事”，恰恰是飞秒瞬态吸收光谱实验里最难预估的背景项——泵浦光经过光纤压缩、展宽、加啁啾之后，激发密度的时间包络变了，ΔOD动力学也跟着变。本文用MATLAB把“高斯脉冲在光纤中的色散演化”和“飞秒瞬态吸收光谱的ΔOD仿真”串成一条可复现的链路，涉及分步傅里叶法、色散长度、泵浦-探测模型和参数取值。适合正要搭飞秒瞬态吸收系统、或想把光纤色散影响先仿清楚再动实验台的工程人员。
2. 高斯脉冲在光纤中的色散演化：从解析解到分步傅里叶法的选型
2.1 高斯脉冲的复包络与啁啾参数：先写对A(0,T)
所有仿真都从初始场写起。光纤里的脉冲用慢变复包络 A(z,T) 描述，其中 T 是随脉冲中心移动的本地时间，T = t - z/vg。无啁啾高斯脉冲的初始形式是：
A(0,T) = sqrt(P0) * exp(-T² / (2T0²))
这里 P0 是峰值功率，T0 是强度降到 1/e 处的半宽。注意很多人把 T0 直接当成半高全宽 FWHM，这是第一个单位坑：高斯脉冲的 FWHM ≈ 1.665 T0，所以“100 fs 高斯脉冲”实际对应的 T0 = 60 fs。写代码时要把这个换算关系写进注释，否则后面算色散长度会差 2.77 倍。
如果脉冲带有初始啁啾，写成：
A(0,T) = sqrt(P0) * exp(-(1+iC)/2 * (T/T0)²)
C 是无量纲啁啾参数。C > 0 表示上啁啾（瞬时频率随时间升高），C < 0 是下啁啾。飞秒振荡器直接输出的脉冲通常带几到几十的啁啾，经过压缩器补偿后接近 C = 0。仿真带啁啾脉冲时，最容易犯的错是把 C 写进功率项——啁啾只在相位指数里，不进幅度指数。记住这个，后面算展宽才不会翻车。
2.2 群速度色散项β2和色散长度LD：决定飞秒和皮秒的边界
高斯脉冲在光纤里传播，核心是群速度色散 GVD。频谱分量不同，走过的群速度不同，脉冲被拉伸。忽略损耗和非线性时，频域解是：
A(z,ω) = A(0,ω) * exp(iβ2 ω² z / 2)
时域里，无啁啾高斯脉冲在传播距离 z 后的脉宽是：
T1 = T0 * sqrt(1 + (z/LD)²)
其中 LD = T0² / |β2|，叫色散长度。这个公式是判断飞秒和皮秒仿真的第一个标尺：
- T0 = 100 fs，β2 = -20 ps²/km：LD = 0.01 ps² / 20 ps²/km = 0.0005 km = 0.5 m。走两米就展宽 4 倍以上。
- T0 = 1 ps：LD = 1 / 20 km = 50 m。走十米只展宽 2%。
这就是为什么标题里“飞秒及皮秒激光”必须分开对待：飞秒脉冲在实验室的几米光纤里就有明显色散，皮秒脉冲可能要几十米才看得出演化。做仿真时如果发现皮秒脉冲在短光纤里几乎不变，不一定是代码错了，先算 LD 再怀疑程序。
2.3 飞秒瞬态吸收光谱的信号来源：泵浦-探测与差分光学密度
飞秒瞬态吸收光谱（TAS）的基本图像是：泵浦脉冲先打到样品上，把一部分分子激发到激发态；探测脉冲延迟时间 τ 后到达，测量样品透过率变化。信号用差分光学密度表示：
ΔOD(λ,τ) = -log10( T_on(λ,τ) / T_off(λ,τ) )
其中 T_on 是有泵浦时的透过率，T_off 是没泵浦时的透过率。ΔOD 大于零对应激发态吸收增强或基态漂白，小于零对应受激发射或透过率增加。
把色散仿真和 TAS 仿真连起来的物理点在于：泵浦脉冲经过光纤和光学系统后，时域包络不再是最初的高斯，甚至变成带拖尾的啁啾脉冲。激发态布居数的建立速度、峰值密度、衰减恢复曲线都取决于泵浦光的时间包络。所以正确的工作流是：先用色散仿真算出到达样品处的泵浦脉冲形状，再用这个形状做激发密度的时间演化，最后换算成 ΔOD。这也是本文中间几章要逐步落地的内容。
3. MATLAB跑通色散演化：SSFM最小实现代码与时间窗、步长、采样点数
3.1 为什么选分步傅里叶法：比解析解更接近真实光纤
解析解只覆盖纯色散情形。一旦脉冲峰值功率高到非线性长度 L_NL = 1/(γP0) 和传播距离同量级，自相位调制（SPM）就会把光谱加宽，此时非线性薛定谔方程（NLSE）不再有简单解析解。工程上最常用的办法是分步傅里叶法（Split-Step Fourier Method, SSFM）：把一小段距离 dz 里的传播分成两步，先在线性算符里考虑色散，再在时域里乘一个非线性相位。好处是线性步在频域是精确的，非线性步在时域只涉及乘法，实现简单且对多数光纤问题足够稳定。
常见的做法是直接用对称 SSFM：一半线性步、一个非线性步、再一半线性步。对大多数实验室场景，这种两步对称格式的精度足够，且不用引入高阶算法。
3.2 完整代码：SSFM求解高斯脉冲在光纤中的色散
以下是可复现的最小 MATLAB 实现，只考虑 β2 色散和 Kerr 非线性，单位统一用 ps 和 km：
%% 参数设置：注意单位统一 T0 = 0.1; % 高斯脉冲T0，单位ps（1/e强度半宽） P0 = 0.1; % 峰值功率，单位W beta2 = -20; % 群速度色散，单位ps^2/km（1550nm单模光纤约-20） gamma = 1.3; % 非线性系数，单位W^-1 km^-1 L = 0.002; % 光纤长度，0.002 km = 2 m Nz = 200; % 纵向步数 %% 时间与频率网格 Ns = 2^12; % 采样点数，2的幂次 Tw = 50; % 时间窗半宽，单位ps dt = 2*Tw / Ns; % 时间步长 t = (-Ns/2 : Ns/2-1) * dt; dw = 2*pi / (Ns*dt); % 角频率间隔 w = (-Ns/2 : Ns/2-1) * dw; %% 初始高斯脉冲 A = sqrt(P0) * exp(-t.^2 / (2*T0^2)); % 无啁啾，幅度单位 sqrt(W) %% 对称分步傅里叶法 h = L / Nz; % 单步长度，单位km LinOp = exp(1i * beta2/2 * w.^2 * h/2); % 半线性步算子 Aout = A; for k = 1:Nz Aout = ifft( LinOp .* fft(Aout) ); % 线性半步（频域） Aout = Aout .* exp(1i * gamma * abs(Aout).^2 * h); % 非线性整步（时域） Aout = ifft( LinOp .* fft(Aout) ); % 线性半步（频域） end %% 输出：对比初始和末端的时域强度 I0 = abs(A).^2; Iout = abs(Aout).^2; figure; plot(t, I0/max(I0), 'DisplayName', 'z=0'); hold on; plot(t, Iout/max(Iout), 'DisplayName', ['z=', num2str(L*1000), 'm']); xlabel('t / ps'); ylabel('归一化强度'); legend show;为什么线性算子里频率项是1i * beta2/2 * w.^2？因为频域 NLSE 的线性项是∂A/∂z = -iβ2 ω² A / 2，积分一步得到exp(-iβ2 ω² h/2)，代码里用1i表示 i，β2 本身带符号，所以写出来就是这个形式。非线性项exp(iγ|A|²h)对应 SPM 相位，强度越高相位变化越快。这个代码中 P0 取值 0.1 W，γP0 约为 0.13 /km，在 2 m 长度上累积相位只有 0.00026 rad，SPM 可忽略，基本是纯色散演化。
3.3 参数怎么设：LD、时间窗Tw、采样点数Ns的配合
对于飞秒脉冲仿真，时间窗和采样点数必须一起定，不能只调一个。经验取值：
- 时间窗半宽 Tw 至少是 T0 的 20 倍以上。100 fs 脉冲用 Tw = 5 ps 够用，但如果走几米光纤后展宽到 1 ps，时间窗要放到 20～50 ps，否则脉冲尾部会被 FFT 周期性折回窗口另一边，看起来像出现了“鬼脉冲”。
- Ns 从 2^12 起步。飞秒脉冲频谱半宽约 0.44/T0（单位 THz），对 T0 = 0.1 ps 约 4.4 THz，对应角频率范围 ±27.6 rad/ps。时间步长 dt 必须满足奈奎斯特条件：dt < π / ω_max。Tw = 50 ps 时 Ns = 2^12 给出 dt = 0.0244 ps，对应 ω_max = 128 rad/ps，余量足够。
- 纵向步长 Nz 先取 200，然后在关键位置加密。判断标准是：步长减半后输出脉宽变化小于 0.1%。对纯色散，线性算子每个 z 步都是精确的，Nz 主要影响非线性项的积分精度；如果加了 SPM 之后出现振荡，把 Nz 提到 500 甚至 1000。
表：常见初始参数的推荐量级
| 脉冲类型 | T0 | Tw | Ns | Nz（2m光纤） |
|---|---|---|---|---|
| 飞秒 100 fs | 0.06 ps | 20 ps | 2^12 | 200 |
| 飞秒 300 fs | 0.18 ps | 50 ps | 2^12 | 200 |
| 皮秒 1 ps | 0.6 ps | 100 ps | 2^13 | 100 |
| 皮秒 10 ps | 6 ps | 400 ps | 2^14 | 50 |
皮秒脉冲的 LD 大，在实验室几米光纤里展宽很小，Nz 可以适当减少，但时间窗必须放大，否则频域分辨率不够，光谱细节会失真。
4. 飞秒瞬态吸收光谱仿真：泵浦-探测ΔOD的建模与MATLAB实现
4.1 泵浦-探测信号的计算：ΔOD怎么从脉冲里取
飞秒瞬态吸收光谱的常见做法分三步：先模拟泵浦脉冲到达样品处的时间包络 I_pump(t)，再算激发态布居数 N_e(t) 的建立与衰减，最后把布居数映射到探测波长上的 ΔOD。简化但物理正确的模型是速率方程：
dN_e/dt = W(t) * (N_g - N_e) - N_e / τ
其中 W(t) = σ_pump * I_pump(t) / (ℏω) 是激发速率，σ_pump 是泵浦波长处的吸收截面，τ 是激发态寿命。探测脉冲一般认为很弱，不会扰动布居数。这样 ΔOD 对延迟时间 τ_d 的依赖可以写成探测光与布居响应函数的重叠积分：
ΔOD(τ_d) = A * ∫ I_probe(t) * N_e(t - τ_d) dt
这里的 A 是比例系数，包含探测截面差和样品厚度等参数。真实实验里，探测脉冲也有一定宽度，当泵浦和探测脉宽都接近 100 fs 时，零延迟附近会出现相干尖峰（coherent artifact），这部分不是真实布居动力学，仿真时可以单独处理。
4.2 MATLAB代码：从色散输出到ΔOD时间扫描
下面的代码把第 3 章的色散结果直接用起来：把到达样品处的泵浦包络作为激发光源，用速率方程算出布居数，再卷积得到 ΔOD(τ_d) 曲线。探测脉冲取高斯形，脉宽设为 80 fs。
%% 输入：色散仿真得到的泵浦包络（来自上一节Aout） Ip = abs(Aout).^2 / max(abs(Aout).^2); % 归一化泵浦强度包络 tp = t; % 时间坐标，ps %% 速率方程参数 sigma_pump = 1e-17; % 泵浦吸收截面，cm^2 hbar_w = 1.55e-19; % 泵浦光子能量，J（对应800nm附近） P_peak = 1e-3; % 单位面积泵浦功率密度，W/cm^2（量级参考） tau_life = 5; % 激发态寿命，ps W0 = sigma_pump * P_peak / hbar_w; % 峰值激发速率，1/ps %% 数值求解布居数演化（欧拉法，步长取时间网格dt） Ne = zeros(size(tp)); Ng = ones(size(tp)); for k = 2:numel(tp) Wk = W0 * Ip(k-1); Ne(k) = Ne(k-1) + dt * (Wk*(1-Ne(k-1)) - Ne(k-1)/tau_life); end %% 探测脉冲与布居数卷积得到ΔOD(τ_d) T_probe = 0.05; % 探测脉冲T0，单位ps（FWHM约83fs） I_probe = exp(-tp.^2 / (2*T_probe^2)); tau_d = tp; % 延迟时间轴，单位ps delta_OD = zeros(size(tau_d)); for j = 1:numel(tau_d) % 卷积核：探测脉冲在延迟tau_d处采样布居数 delta_OD(j) = sum( I_probe .* Ne .* exp(-(tp - tau_d(j)).^2 / (2*T_probe^2)) ); end delta_OD = delta_OD / max(delta_OD); % 归一化，方便和实验对比 %% 画曲线 figure; plot(tau_d, delta_OD); xlabel('delay / ps'); ylabel('normalized \DeltaOD'); title('Simulated femtosecond transient absorption trace');这个代码里有两个参数要特别解释。第一，W0的绝对值很依赖样品浓度和泵浦能量，实际实验里很难先验确定，所以代码最后做了归一化；对比实验时应该把泵浦能量与实际测得的激发密度标定关联起来。第二，tau_d的扫描范围和步长要覆盖你关心的动力学尺度：寿命 5 ps，延迟扫描范围至少到 20 ps，步长要小于探测脉宽，这里用时间网格自身的 dt（约 0.024 ps），足够分辨 100 fs 量级的上升沿。
4.3 参数表：脉宽、延迟步长、色散参数与时间窗的配合
飞秒 TAS 仿真参数和实际实验系统的对应关系，可以用一张表说清楚。这里假设泵浦中心波长 800 nm，探测在可见区。
| 参数 | 典型值 | 仿真取值建议 | 说明 |
|---|---|---|---|
| 泵浦脉宽（FWHM） | 80～150 fs | T0 = 0.05～0.09 ps | 决定激发上升沿的极限分辨率 |
| 泵浦经过光纤后的脉宽 | 300 fs～2 ps | 由SSFM算出 | 直接决定ΔOD上升沿斜率 |
| 延迟扫描范围 | -5～50 ps | Tw ≥ 50 ps | 必须覆盖激发态寿命的3倍以上 |
| 延迟步长 | 实验延迟线精度 | dt ≤ 0.03 ps | 步长要能分辨脉冲上升沿 |
| β2 不确定度 | ±1 ps²/km | 跑两个边界值 | 对比上下界，评估对实验的影响 |
| 激发态寿命 τ | 1～100 ps | 分别扫1/5/20 ps | 看系统是寿命主导还是脉宽主导 |
皮秒泵浦的情形则相反：泵浦脉宽 1 ps 时，上升沿分辨率本来就只有皮秒量级，延迟步长取 0.1 ps 已经足够，时间窗拉到 200～500 ps 看长寿命衰减。飞秒和皮秒的区别不只是数字变了，时间窗、步长、是否要分辨相干尖峰这三个判断都要跟着变。
5. 仿真翻车排雷：色散和飞秒瞬态吸收的5个高频坑
5.1 现象：脉冲“穿墙”出现在时间窗另一侧
用 SSFM 计算飞秒脉冲传播时，如果设 Tw = 5 ps、T0 = 0.1 ps，在短距离内结果正常；一旦光纤长度增加到展宽超过时间窗一半，出射波形会在窗口左侧多出一个同样的脉冲，看起来像能量不守恒。
原因：FFT 隐含周期性边界条件，时域信号超出时间窗后从另一侧绕回。这不是物理现象，是数值周期性造成的混叠。
解决：把时间窗加大到预计展宽脉宽的 10 倍以上。工程做法是先粗算一遍：T1 = T0 * sqrt(1 + (z/LD)²)，然后取 Tw = max(10T1, 20T0)。对 100 fs 脉冲走 5 m 光纤，T1 可能到 1 ps，Tw = 20 ps 足够。
5.2 现象：皮秒脉冲在短光纤里“没有变化”，怀疑代码错了
T0 = 1 ps、β2 = -20 ps²/km 时 LD = 50 m。如果仿真光纤长度只有 2 m，展宽比例只有 sqrt(1 + 0.04) ≈ 1.02，几乎看不出来。这是物理，不是 bug。
原因：皮秒脉冲的色散长度远大于实验室常见光纤长度，色散效应很弱。
解决：不修改代码，先扫参数。把光纤长度从 1 m 扫到 100 m，看展宽随距离的变化是否符合 T0 * sqrt(1 + (z/LD)²)。如果符合，说明代码没问题；如果想在短光纤看到效果，要么减小 T0，要么改用正负 β2 交界附近的色散平坦光纤。
5.3 现象：加入SPM或更高阶色散项后波形振荡发散
在 3.2 代码里加入三阶色散 β3 或加大 P0 之后，出射波形出现高频振荡，甚至直接变成噪声。
原因：非线性步的相位γ|A|²h太大，单个 dz 内累积相位超过 π，分步傅里叶法的本地精度失效。对飞秒脉冲，频带很宽，时间步长不够细时频谱边缘会有混叠。
解决：增加 Nz，把 dz 缩小，直到输出不再变化；同时检查时间窗是否覆盖了整个光谱范围。经验法则是让 γP0h ≤ 0.01 rad。另一个手段是改用自适应步长：每一步比较半个步长和整步长的结果，误差超过阈值就把步长减半。
5.4 现象：ΔOD曲线在零延迟时间附近出现诡异的窄峰
泵浦-探测仿真中，当泵浦和探测脉宽都是亚百飞秒量级时，零延迟附近会出现一个很窄的正峰或负峰，宽度远小于激发态寿命。
原因：泵浦和探测脉冲在时域上直接重叠时，存在相干耦合信号（coherent artifact），包括交叉相位调制和频谱干涉，不是布居动力学响应。速率方程模型只描述布居数的建立，但如果代码里把探测卷积直接在零延迟处做了，相干项会被误算成信号。
解决：明确区分两个贡献。布居项用速率方程输出，相干项单独用一个高斯或双曲正割尖峰加在ΔOD曲线零延迟处。实验上也有办法区分：相干尖峰在泵浦和探测波长同时存在时才有，改变泵浦聚焦位置或交叉偏振可以抑制它。
5.5 现象：同样的代码在MATLAB里跑得很慢，循环内反复调用fft
如果直接用嵌套 for 循环逐点计算 4.2 里的 ΔOD 卷积，数值小还好，一旦时间窗到 100 ps、延迟步长到 0.02 ps，就是 5000 × 5000 的循环，几分钟起步。
原因：把卷积写成了逐元素求和，而 MATLAB 对实数矩阵运算的向量化能力远强于标量循环。
解决：把卷积改成矩阵乘法，一次算出所有延迟点。常见做法是构造 Toeplitz 矩阵或用conv函数代替循环，核心是把 I_probe 时间翻转后和 Ne 做标准卷积，再用tau_d取对应的中心部分。另外提前预分配delta_OD数组，避免循环里动态增长。
6. 用自相关曲线和二维扫参把仿真结果变成实验定标依据
仿真不是终点，能指导实验才有价值。这里给三个具体技巧，我每次做飞秒 TAS 实验规划时都会用到。
第一个技巧是用自相关曲线验证色散仿真。实验上很难直接测 100 fs 脉冲的时域波形，但可以用强度自相关仪测得宽度。仿真里造一个虚拟自相关器：对出射脉冲强度做自相关I_corr = xcorr(Iout)，再量自相关曲线半高宽，和实验值对比。自相关宽度和真实脉宽的比值对高斯脉冲是 sqrt(2)，如果实验和仿真在这个比值上一致，基本可以确认 β2、光纤长度这些输入是对的。
第二个技巧是把色散参数扫成一张二维图。固定 T0，以 z 为横轴、β2 为纵轴，画等脉宽线，能直接看出“这段光纤里色散要不要补偿”。例如 100 fs 脉冲，β2 = -20 ps²/km、z = 1 m 时展宽 2 倍；改用 β2 = -1 ps²/km 的色散位移光纤，展宽就降到可忽略。这个图我通常在实验方案里直接引用，比单条曲线有说服力得多。
第三个技巧是用仿真反推激发密度和 ΔOD 量级。做 TAS 之前不知道样品吸收截面时，先用 4.2 代码跑一个“最差情况”：泵浦经光纤后脉宽展宽 3 倍，峰值功率降 3 倍，激发密度不够时 ΔOD 会接近噪声基底。这时仿真会告诉你一个明确的方向：换色散更小的光纤、减短光纤长度、或者加大泵浦能量，而不是在实验台上盲试。这也是我为什么坚持把色散仿真和瞬态吸收仿真写在同一个流程里——分开算，每一个看起来都对，合起来才会发现泵浦到达样品处的脉宽已经毁掉了时间分辨率。
我个人的习惯是每次实验前花半小时跑一遍这套最小仿真，把三张图存下来：出射脉冲形状、ΔOD 预计曲线、展宽随光纤长度的等值线图。遇到实验数据和仿真趋势对不上的时候，先怀疑泵浦聚焦和色散补偿，再怀疑样品降解。仿真给不了你绝对正确的 ΔOD 数值，但能给你一个判断基线在哪里。希望帮到你。
本文还有配套的精品资源，点击获取
