双光子吸收仿真精讲：方程、分布傅里叶法、参数与避坑全解

简介压缩包提供了一套与双光子吸收TPA计算和优化直接相关的MATLAB脚本集合主要面向光学、凝聚态物理、化学及生物光子学方向的研究生、科研人员与高年级本科生可用于TPA截面估算、吸收谱模拟、激发条件优化等任务。资源以zip格式打包共9个文件全部为.m脚本整体仅11KB体积非常精简适合直接按函数模块取用或集成到现有实验中。从内容结构看脚本按高斯贡献分解、策略一/策略二优化、截止功率计算、原始最大值定位、相位窗口分析、最优总和求解等功能划分构成一条从TPA谱形分解到参数寻优的完整处理链。已有130人学习尤其适合刚进入非线性光学模拟领域、需要可运行示例的读者可借助这套代码快速理解双光子吸收效应的数值表征方法并在此基础上结合具体材料模型进行二次开发。双光子吸收在生物成像、激光加工和光存储等领域应用广泛这份小而精的脚本包为相关研究提供了便捷的起点。1. 双光子吸收仿真打开 tpa.zip 前先想清楚的问题在代码资源站上tpa.zip 这类压缩包几乎和“双光子”三个字绑定出现解压后通常是一段非线性脉冲传播仿真的源码多数用 MATLAB 或 Python 写的一维分布傅里叶法。问题在于很多人把参数塞进去跑出来的透射率曲线却和实验对不上于是回头改步长、改分辨率、改边界折腾半天最后发现是双光子吸收项根本没被写进方程或者单位换算错了三个数量级。下面这套思路把过程讲透双光子吸收TPA在传播方程里为什么不是一个普通损耗项分布傅里叶法怎么把这项处理得又快又稳参数的物理单位怎么换算以及那些看起来在“修代码”、实则被数值坑卡住的常见流坑。适合正在做超短脉冲传输、半导体波导器件或非线性光谱测量的人。如果你只想在现有仿真包基础上把 TPA 补进去重点看第 3 章的写法和第 5 章的避坑清单要自己搭一套能拟合实验数据的模型从第 2 章开始读更稳。2. TPA 仿真先从方程立住非线性损耗项与求解器选型为什么先讲方程因为 tpa.zip 这类包里最容易出的问题就是“损耗项的单位和符号写错”。双光子吸收的物理并不复杂但放进传播方程时系数的来历、符号、归一化约定都会直接影响结果。这一章先把模型写清楚再解释分布傅里叶法为什么是默认解法。2.1 双光子吸收在传播方程里长什么样对于沿 z 方向传播的慢变包络 A(z, T)同时考虑色散、自相位调制SPM和双光子吸收时常见写法是∂A/∂z -α/2 A - iβ₂/2 ∂²A/∂T² iγ|A|²A - β_tpa/2 |A|²A四个项从左到右依次是线性吸收损耗、群速度色散GVD、自相位调制带来的非线性相位、双光子吸收带来的非线性强度损耗。其中 γ 是 Kerr 非线性系数β_tpa 是双光子吸收系数。这里最容易看错的是第四项SPM 项的系数 iγ 是纯虚数只改相位不耗能量TPA 项的系数 -β_tpa/2 是负实数对应“强度相关的损耗”它让脉冲在峰值处衰减得最快、两翼衰减得慢。一个高斯脉冲走一段距离后峰会被压低脉形从高斯慢慢变成平顶再塌陷。物理根源是两个光子同时被介质吸收吸收概率正比于光强的平方所以损耗功率密度写成 β_tpa I²。单位上β_tpa 常用 cm/GWI 用 GW/cm²乘起来恰好是无量纲速率。把强度 I 换成包络模平方 |A|² 并保持能量归一化约束时方程里就会出现 -β_tpa/2 |A|²A。写方程时最常踩的坑有两个。第一个是直接把 β_tpa 当线性吸收系数 α 用写成 -β_tpa A这样算出来的透射率跟强度无关完全失掉双光子吸收“越强越吸”的特征。第二个是系数 1/2 的来路方程描述的是场包络幅度而损耗发生在功率密度上功率正比于 |A|²所以幅度衰减率是功率衰减率的一半。这个 1/2 不能丢丢了最终透射率会差出一个平方关系。给出参考量级近红外波段典型半导体的 β_tpa 在 0.55 cm/GWKerr 系数 γ 在 1100 W⁻¹m⁻¹群速度色散 β₂ 在 -100100 ps²/km 之间。具体数值取决于材料与波长仿真时最好查材料库或实测值不要用某个包里的默认数据直接外推。很多 tpa.zip 里的默认参数来自某个特定实验波长和脉宽与手头场景往往不匹配。这里还有一个容易忽略的自检方法从方程可以推能量守恒。|A|² 对时间积分后对 z 求导SPM 和色散项贡献为零只有 α 项和 β_tpa 项贡献负值。如果你写的代码跑出来能量不降甚至上升那一定是符号或单位写错了。这个自检只需要在代码里加一行energy.append(np.sum(np.abs(A)**2)*dt)后面几章会一直用到它。2.2 分布傅里叶法为什么是 TPA 仿真的默认解法解决这个方程的标准做法是分布傅里叶法也叫分步傅里叶法。思路是把 z 方向切成很多小段 Δz在每一段内假定线性效应和非线性效应可以分开计算先忽略非线性只让色散和线性损耗在频域里作用再忽略色散只让 Kerr 和 TPA 在时域里作用。频域那一步用 FFT 完成时域那一步只是逐点乘指数因子所以单步开销是 O(N log N)比直接在时域做有限差分的网格推进快得多。为什么不用时域 FDTD超短脉冲在波导里传播几厘米到几十厘米空间网格要覆盖整个传播长度而时间网格要分辨光周期量级的变化网格数轻易到几百万甚至上亿一步时间推进还受 CFL 条件限制实际慢到没法用。分布傅里叶法把时间方向的偏导变成频域乘法色散项在频域是解析的对这类慢变包络方程几乎是标准答案。求解器选型可以简单对比如下求解方式单步复杂度色散处理适用场景主要限制分布傅里叶法O(N log N)频域解析精确慢变包络、长距离传播要求包络近似成立时域有限差分 FDTDO(N) 每网格点需要数值差分短距离、强非慢变网格量大CFL 限制显式 Runge-Kutta 配 FFTO(N log N)频域精确可用但步长偏保守每步多次 FFT较慢实际程序里我习惯用对称分布傅里叶法每一步先走 Δz/2 的线性段然后走 Δz 的非线性段再走 Δz/2 的线性段。对称安排的局部误差从 O(Δz) 降到 O(Δz²)同样的精度下步长可以放大好几倍。对脉冲传播来说步长 Δz 的经验上限是取色散长度 L_D T₀²/|β₂| 和非线性长度 L_NL 1/(γP₀) 里较小者的十分之一到五分之一。TPA 项加入后会额外消耗能量峰值功率沿 z 单调下降实际误差往往比纯 Kerr 情形小但保守起见仍按最小特征长度取步长。有一个细节容易被忽视非线性段里的 TPA 是实数损耗所以要乘的是 exp(-β_tpa/2 |A|² Δz)不是乘相位因子。很多从纯 SPM 仿真改过来的人在这里习惯性写了 exp(iγ|A|²Δz)把虚数单位带进去TPA 变成了一种“相位扰动”能量完全没掉曲线自然不对。判断有没有写对的方法是关掉色散跑一炮只有 TPA 和 SPM 时脉宽应该几乎不变总能量逐步下降频谱变宽来自 SPM而能量下降必须来自损耗项的实部。3. 把 TPA 仿真跑起来最小可运行代码与参数怎么设这一章直接给一个能运行的 Python 版本。用 Python 是因为 tpa.zip 里的 MATLAB 代码虽然多但 numpy/scipy 做 FFT 和后处理更顺手而且容易被集成到自动化扫描里比如后面第 6 章的参数反演。代码保持单文件没有第三方非线性依赖。3.1 先定介质与光源一张能直接抄的参数表仿真之前把下面这张参数表填好大多数“结果不对”都出在单位换算上。长度单位统一用 m功率单位用 W双光子吸收系数先换算成 m/W再结合有效模场面积 A_eff 变成一个组合系数避免在代码里到处乘 1e-4、1e-9。下面这套参数对应 1550 nm 波段、脉宽 100 fs、峰值功率 10 W 的脉冲在某半导体波导里传播 1 cm参数符号数值说明峰值功率P010 W决定 TPA 损耗率脉冲半高全宽T_FWHM100 fs换算 T0 T_FWHM/1.665群速度色散β₂-20 ps²/km反常色散换算为 -20e-27 s²/mKerr 系数γ20 W⁻¹m⁻¹高非线性波导典型量级双光子吸收系数β_tpa1 cm/GW换算为 1e-11 m/W有效模场面积A_eff1 μm²来自模式仿真或解析近似传播长度L1 cm时间窗口Twin10 ps约 100 个脉宽采样点数N2^14FFT 点数参数说明T0 是 1/e 强度的脉冲半宽对高斯脉冲 T0 T_FWHM/1.665。A 的归一化约定是 ∫|A|² dt 脉冲能量因此 |A|² 的单位是 W。组合 TPA 损耗系数 b2_eff β_tpa / A_eff单位是 (W·m)⁻¹。用这套参数算出来峰值处 TPA 等效损耗率 b2_eff × P0 10 × 10 100 m⁻¹传播 1 cm 总衰减 e⁻¹ ≈ 0.37效果很明显。如果峰值功率只有 0.1 W损耗率只有 1 m⁻¹1 cm 内基本看不到 TPA 痕迹。所以判断一个参数组合有没有意义先算 b2_eff × P0 × L这个量级决定透射率曲线变化幅度。有效模场面积 A_eff 不是随便填的。在波导里A_eff (∫|E|² dA)² / ∫|E|⁴ dA通常由模式求解器给出。如果手头没有模式仿真可以先按波导芯层几何尺寸估算但误差会直接换算进透射率绝对值。后面拟合 β_tpa 时A_eff 给错是最难排查的系统误差来源建议单独标定一次。3.2 分布傅里叶法主循环可直接改参数的脚本把下面代码存成 tpa_sim.py直接运行就会生成 tpa_out.npz里面有能量衰减曲线、输出脉冲和对应 z 坐标。import numpy as np from numpy.fft import fft, ifft, fftfreq # ---- 物理参数全部使用 SI 单位 ---- P0 10.0 # 峰值功率, W T_fwhm 100e-15 # 脉冲半高全宽, s T0 T_fwhm / 1.665 # 高斯脉冲 1/e 幅度半宽, s beta2 -20e-27 # 群速度色散, s^2/m-20 ps^2/km gamma 20.0 # Kerr 非线性系数, W^-1 m^-1 alpha 0.0 # 线性吸收, m^-1先关掉聚焦 TPA beta_tpa 1e-11 # 双光子吸收系数, m/W1 cm/GW A_eff 1e-12 # 有效模场面积, m^21 um^2 b2_eff beta_tpa / A_eff # 等效 TPA 损耗系数, (W·m)^-1 L 0.01 # 传播长度, m1 cm nz 500 # 纵向步数 N 2**14 # 时间采样点 Twin 10e-12 # 时间窗口, s约 100 个脉宽 # ---- 时间 / 频率网格 ---- t np.linspace(-Twin/2, Twin/2, N, endpointFalse) dt t[1] - t[0] f fftfreq(N, dt) omega 2 * np.pi * f # ---- 初始高斯包络峰值 |A|^2 P0 ---- A np.sqrt(P0) * np.exp(-0.5 * (t/T0)**2) # ---- 频域线性算子 ---- lin_op -alpha/2 - 1j * (beta2/2) * omega**2 dz L / nz energy [] for _ in range(nz): # 半个线性步频域乘色散线性损耗因子 A ifft(fft(A) * np.exp(lin_op * dz/2)) # 非线性步SPM 相位 TPA 实数损耗 A A * np.exp(1j * gamma * np.abs(A)**2 * dz - (b2_eff/2) * np.abs(A)**2 * dz) # 后半个线性步 A ifft(fft(A) * np.exp(lin_op * dz/2)) energy.append(np.sum(np.abs(A)**2) * dt) np.savez(tpa_out.npz, znp.arange(1, nz1) * dz, energynp.array(energy), A_finalA, tt)代码逻辑说明初始高斯包络用np.sqrt(P0)赋值保证峰值 |A|² P0这样 |A|² 的单位严格是 Wγ 和 b2_eff 直接按 SI 数值使用即可。线性算子lin_op里的 β₂ 已换算成 s²/m-20 ps²/km -20 × 1e-24 s² / 1e3 m -20e-27 s²/m。非线性步里第一项是纯虚相位第二项是实数损耗两项合并成一个 exp因为 exp(ab)exp(a)exp(b)顺序不重要。但要注意这里 Δz 用的是整个 dz而不是 dz/2这是对称分布傅里叶法的标准形式两个半线性步夹一个整非线性步。参数怎么改如果算到一半能量曲线出现锯齿或 NaN第一反应不是缩小 nz而是检查时间窗口 Twin 是否足够大。高斯脉冲在反常色散下会展宽窗口太小会让能量“流出”计算域能量积分越来越小看起来很像 TPA 损耗。一般把 Twin 设成初始脉宽的 200 倍以上并对非线性强的 case 每隔若干步检查一次 max(|A|²) 对应的半宽是否超过窗口的 1/10。还要记得频谱分辨率 Δf 1/Twin如果 Twin 太大而 N 不够频域采样间隔不足色散算子采样会失真。3.3 从仿真结果里读透射率与脉冲形变跑完脚本后先看能量曲线。把能量归一化 E(z)/E(0) 画出来调好参数后 1 cm 处大约是 0.37。如果同时开着 SPM峰值功率沿纵向不断降低TPA 损耗率跟着降低能量衰减会略微偏离单指数表现出“先快后慢”的特征。这是 TPA 区别于线性吸收的信号线性吸收的归一化能量衰减不随初始功率变化而 TPA 的衰减曲线随 P0 增大明显变陡。做一次功率扫描P0 分别取 1 W、3.16 W、10 W如果三条归一化能量曲线几乎重合说明 TPA 项没起作用如果间距明显拉开且功率越大掉得越快说明模型对了。第二个要看的是脉冲时域形状。由于 TPA 对峰值损耗更大脉冲峰会被削平高斯轮廓逐渐变成类似抛物线顶。如果 SPM 和 GVD 同时存在时域形状会和纯 TPA 不同但“峰值损耗最强”这个趋势不会变。第三个是频谱SPM 产生对称展宽TPA 因为损耗随强度变化频谱对称性会轻微破坏出现红蓝不对称这在实验中常作为 TPA 的诊断证据。仿真里如果看到频谱异常不对称且能量掉得很多先检查是不是窗口溢出再谈物理。4. 结果判读与边界条件双光子吸收与线性吸收、饱和效应的边界在哪跑通最小代码之后真正麻烦的是判读。仿真曲线“在动”不等于“物理对了”。这一章把 TPA 模型的有效边界讲清楚什么时候 TPA 会盖过线性损耗什么时候必须加三光子吸收和自由载流子吸收以及常见的边界条件误用。4.1 线性吸收与双光子吸收的分界临界功率判据用一个简单判据把两条损耗路径分开。线性吸收率是 α单位 m⁻¹TPA 的等效损耗率是 b2_eff × P(z,t)单位 m⁻¹。两者相等时对应一个临界功率P_crit α / b2_eff当瞬时功率远低于 P_crit 时总损耗主要由 α 决定能量曲线呈指数衰减功率超过 P_crit 后双光子吸收随功率增大而增强损耗不再满足指数衰减。实验里测透射率随入射功率的变化就能区分线性吸收区透射率基本不随 P0 变化进入 TPA 区后随 P0 增大透射率下降在 log-log 坐标里高功率段斜率和 β_tpa 有效长度相关。这个判据对实验设计很重要。在半导体波导中测 β_tpa入射光功率太低时信号被 α 淹没太高时又会出现三光子吸收和自由载流子吸收所以最佳测量窗口在 P_crit 附近一到两个量级内。仿真时同理如果 α 很大而 P0 很小你会看到漂亮的指数衰减但那是线性吸收跟双光子没有关系。还要区分峰值功率和平均功率。TPA 是瞬时强度平方效应只有峰值功率有意义。同一台激光器重频 1 MHz、脉宽 100 fs、平均功率 1 mW 时峰值功率是 10 W但如果重频提高到 100 MHz平均功率不变峰值功率还是 10 W——注意此时 TPA 损耗不变但热效应和载流子累积效应会显著增强。仿真里如果只给了平均功率必须先换算峰值功率P_peak P_avg / (f_rep × T_FWHM)对高斯脉冲还要乘一个约 1 的波形因子很多 tpa.zip 示例里直接拿平均功率当 P0曲线被高估或低估几个数量级都不自知。4.2 更高光强下的模型边界三光子吸收、自由载流子吸收与饱和双光子吸收系数 β_tpa 并不是常数。光强极高时会发生两种效应。第一是三光子吸收损耗项写成类似 -b3/3 |A|⁴Ab3 单位是 W⁻²m⁻¹强度再高时更高阶吸收也会依次出现。第二是自由载流子吸收FCATPA 产生的电子-空穴对被后续光子吸收形成一个与载流子密度成正比的附加损耗项。载流子寿命远大于脉冲宽度时这项会“记住”前面脉冲积累的能量表现为随脉冲重复频率升高而增强的损耗也就是载流子累积效应。如果实验测量显示透射率在低功率区符合 β_tpa 模型高功率区又有额外损耗最简单的扩展是在方程里增加 -σ_c N_c(z,t) A/2其中 N_c 由速率方程 dN_c/dt β_tpa I²/(2hν) 积分而来σ_c 是自由载流子吸收截面。加这项后能量曲线会多出一个与脉冲历史相关的损耗尾巴。更精细的做法是把 β_tpa 写成强度相关的饱和表达式β_tpa(I) β_tpa0 / (1 I/I_sat)这在宽带隙半导体里更符合实际。仿真阶段我的建议是先不加高阶项用最小模型把趋势跑出来再用实验数据判断偏差方向残差在低功率端大多半是 β_tpa 或 A_eff 给错了在高功率端大才需要引入 FCA 或三光子吸收。一开始就全部加上参数太多互相补偿反而不容易收敛拟合时更是灾难。边界条件也有坑。分布傅里叶法隐含周期性边界条件如果时间窗口不够宽脉冲展宽后碰到边界会“绕回”窗口另一端在能量曲线上表现为诡异的回升。检测方法是看最终时刻的 A_final 在窗口两侧是否趋近于 0不趋近就把 Twin 加倍重跑而不是靠加窗函数硬压。加窗会吃掉真实物理尤其是对 TPA 这种强度相关效应加窗后的峰值功率时间分布失真损耗被低估。5. 双光子吸收仿真避坑清单五个常见流坑的现象、原因与解决这一章集中写我实际跑 TPA 仿真时反复踩过的坑。每一条按现象、原因、解决三个层次写方便你在出问题时直接对照排查。5.1 能量曲线出现 NaN 或锯齿现象跑到中间某一步energy 突然变成 nan或者曲线像锯齿一样上下跳动。原因有二。一是非线性步的损耗因子 exp(-b2_eff × |A|² × dz/2) 在脉冲峰值处可能变得极小如果 dz 太大峰值样点一次衰减几十个数量级浮点数下溢二是时间窗口不足脉冲展宽后能量在 FFT 里混叠功率谱出现非物理振荡时域 |A|² 跟着乱跳。解决先把 nz 增大 10 倍看曲线是否平滑仍不平滑就把 Twin 增大到 100 倍脉宽以上N 提到 2^15。若问题依旧检查 β_tpa 换算1 cm/GW 1e-11 m/WA_eff 少写一个 1e-12损耗率会差出数量级。5.2 透射率只呈指数衰减双光子特征消失现象把 P0 从小调到大归一化能量曲线几乎重合或者只是整体按比例下移。原因把 TPA 项误写成 -b2_eff × A 的线性损耗形式或者在代码里把 β_tpa 当成 α 用另一个可能是 A_eff 设得过大峰值强度太低TPA 损耗率远小于 α。解决做一次功率扫描验证。P0 依次取 1、3.16、10如果归一化能量曲线明显不同且随 P0 增大变陡说明 TPA 项在起作用三条曲线完全重合就是写错了。这个测试应当写进回归用例以后改参数时随时跑一遍。5.3 仿真和实验的透射率绝对值对不上但趋势一致现象能量衰减的相对变化与实验一致绝对数值差 25 倍怎么调 β_tpa 都拉不平。原因脉冲的时间形状不是高斯。实验超短脉冲可能是双曲正割或带啁啾峰值功率相同但能量分布不同另外 β_tpa 实测值随波长变化明显用包里的默认值按另一个中心波长算自然对不上。解决先拟合实测脉冲自相关把初始包络 A(0,t) 换成 sech 或实测插值形状再反演 β_tpa而不是拿高斯结果硬套实验。啁啾也要注意如果输入脉冲有上啁啾前几厘米的峰值会暂时变低TPA 损耗比无啁啾情形小拟合时初始啁啾必须作为已知量带入。5.4 频谱异常不对称能量掉得反常快现象输出频谱左或右出现一条长尾巴能量损耗比理论估计大得多。原因时间窗口溢出导致 FFT 周期性混叠本质是频谱泄漏不是 TPA 物理。分布傅里叶法对色散和 SPM 的对称性比较敏感窗口边界反射会破坏频谱对称。解决Twin 至少取初始脉宽的 200 倍当脉冲被色散展宽 3 倍以上时要加大 Twin 后重跑。判断标准是输出脉冲两翼幅度低于峰值的 1e-6。如果确实需要长时间窗口可以考虑吸收边界层但嵌进分布傅里叶法会破坏周期性假设不如直接把窗口做大来得干净。5.5 步长找不到稳定收敛区间现象nz 从 500 加到 5000 曲线不变从 500 减到 100 就明显偏差中间找不到平滑收敛区间。原因对称分布傅里叶法的步长上限由最小特征长度决定。纯色散时 L_D 18 cm 很大纯 SPM 时 L_NL 0.5 cm但加上 TPA 后特征长度随功率实时变化不能用固定 L_NL 去算。解决把 dz 改成自适应。每走一步前先用试探 dz 算非线性相位增量 max(γ|A|²dz b2_eff|A|²dz/2)超过阈值就减半重算小于阈值的四分之一就加倍。这样既稳又省时间。阈值经验值单步非线性相位改变不超过 0.05 rad单步损耗因子不小于 0.9。6. 进阶用透射率曲线反演 β_tpa 的收敛流程仿真跑到稳定只是第一步真实需求往往是从测量数据里标定材料参数 β_tpa。常见做法是把第 3 章的传播函数封装成可调用形式输入峰值功率和 β_tpa输出透射率def solve_tpa(P0, b2_eff, gamma20.0, nz500): # 内部逻辑同第 3 章返回归一化透射率 Tr ... return Tr然后对一组实测入射功率和透射率做最小二乘拟合。反演里最容易翻车的不是优化器而是初始猜测和参数退化。我的习惯是先在 0.110 cm/GW 的 log 网格上粗扫找到残差最低的区间再用局部优化精扫同时把 A_eff 和 γ 固定只让 b2_eff 自由。第一次做这种反演时我只用了三个功率点结果三个点彼此矛盾后来才意识到实验里没扣除线性吸收 α。于是把 α 也作为自由参数发现 α 和 b2_eff 高度相关——低 α 高 b2 和高 α 低 b2 给出一模一样的曲线。解决办法是先测低功率段线性衰减独立定 α再在高功率段反演 b2否则拟合纯属看运气。反演时数据要在“TPA 效应占主导”的功率区间取每个功率点对应的 TPA 等效损耗至少是线性损耗的 2 倍以上。这个条件可以先用粗估参数预先扫描确认。还要保证测量模型和仿真模型用同一个能量归一化约定透射率才可比。示数里保持 nz 500 足够因为迭代中损耗率变化平缓粗拟合不需要高精度最终确定参数后再用自适应步长精跑一次做验证。这套流程我反复用过多次收敛稳定。如果你想更稳就把功率扫描点从 3 个加到 10 个并且每个点重复测量取平均拟合残差能小一个量级。希望帮到你。本文还有配套的精品资源点击获取