高次谐波数值模拟：MATLAB开源仿真如何耦合量子与宏观

简介：面向等离子体物理与超快光学方向研究者的MATLAB仿真源码包，围绕高次谐波产生（HHG）过程的数值模拟展开，涵盖量子力学、超音速流物理、相位匹配与吸收等核心物理机制。压缩包中共有451个文件，包括69个M脚本、371个TXT计算结果、3个PDF文档、2个FIG图形文件等，总体积约29.75MB，目录结构清晰，便于按需取用。已有714人学习下载。源码完整实现了从喷嘴出口超音速气流、等离子体动态离子、自由电子量子响应到相位匹配与吸收的整个模拟链路；图形用户界面（GUI）允许用户交互式调节输入参数，并即时查看电离、相位匹配和偶极响应的结果，所有计算数据自动保存为TXT文件，便于后续复现与提速。既可用于科研中的物理现象分析，也适合高等院校相关课程的实验教学与二次开发。

1. 高次谐波数值模拟：这份MATLAB开源的等离子体仿真是怎么把量子力学和超音速流耦合起来的

高次谐波产生（High-Harmonic Generation, HHG）是极紫外超快光源的核心机制，也是阿秒物理最经典的数值实验场景。很多人第一次接触它时被量子力学和超音速流物理两个词吓住，实际上拆开看就两件事：单原子层面求解含时薛定谔方程，宏观层面把气体喷流里成千上万个原子的辐射相干叠加。这套系统开源的MATLAB源码正好把两条链路都打通了，适合做超快光学、强场物理方向研究，也适合想快速验证理论模型的学生。它能帮你复现“平台区+截止频率”的典型谐波谱结构，也能让你看到密度、气压、喷流位置这些宏观参数怎么改变相位匹配，省掉从零写求解器的大量时间。

2. 物理模型与数值框架：三步模型、截止规则和SFA的取舍

2.1 为什么三步模型是理解这套代码的起点

高次谐波产生的物理图像可以用半经典的三步模型概括：强激光场把原子势垒压低，电子隧穿电离进入连续态；随后电子在激光场中加速，获得大量动能；当激光场反向时，一部分电子被拉回并回撞到母离子，复合过程中释放出截止能量为 Ip + 3.17Up 的高能光子。这里的 Ip 是电离势，Up 是 ponderomotive 能量，与激光强度 I 和波长 λ 的关系是 Up ≈ 9.33 × 10^-14 × I × λ^2（I 单位为 W/cm²，λ 单位为 μm）。

这套源码的第一步就是把这个三步模型变成可计算的参数输入。你会在初始化脚本里看到电离势、激光波长、峰值强度、脉宽和载波包络相位（CEP）这几个变量，它们直接决定隧穿电离率、电子回撞能量分布和谐波谱的截止位置。我一般会先把 Ip 和 Up 代入 3.17 规则算一遍理论截止阶次，再和代码输出比对，能快速判断参数是否设得合理。

三步模型的局限在于它没有包含量子干涉和色散效应。实际的谐波谱是大量量子轨道叠加的结果，短轨道贡献平台区，长轨道贡献截止附近的非定域结构。这就是为什么源码里除了半经典计算，还必须保留完整的量子力学求解路径。

2.2 TDSE与SFA两条路：怎么选、边界在哪

严格的高次谐波模拟需要求解含时薛定谔方程（TDSE），这是一条几乎不做近似但代价极高的路。以一维软核势模型为例：

% 软核势定义：Z是有效核电荷，a是软化参数，避免奇点 V = @(x) -Z ./ sqrt(x.^2 + a^2); % 初始波函数取基态，用虚时间传播求基态更稳 [Psi0, E0] = imaginary_time_propagation(V, Nx, dx);

这里的 Z 和 a 要按目标原子调节。比如模拟氩原子，Z 一般取 1 附近，a 取 1.5 左右能使基态能量接近 -0.58 a.u.（实验值）。虚时间传播是求基态的稳定办法，不依赖猜测初始波函数。我见过有人直接用高斯波包当基态，算出来的谐波谱容易出现虚假的宽峰，因为初始态不是本征态。

TDSE 的时间步进用分裂算符或 Crank-Nicolson 格式都行。分裂算符的好处是每个时间步里动能项和势能项分别处理，二阶精度且稳定：

% 动力学传播：二阶分裂算符 for t = 1:Nt % 1. 在动量空间传播 dt/2 Psi_k = fft(Psi_x); Psi_k = Psi_k .* exp(-1i * 0.5 * k.^2 * dt); Psi_x = ifft(Psi_k); % 2. 在坐标空间加激光场和势场传播 dt Psi_x = Psi_x .* exp(-1i * (V(x) + laser(t, E0, omega, CEP)) * dt); % 3. 再在动量空间传播 dt/2 收尾 Psi_k = fft(Psi_x); Psi_k = Psi_k .* exp(-1i * 0.5 * k.^2 * dt); Psi_x = ifft(Psi_k); end

这段代码里，k 是动量空间栅格，dt 是时间步长，laser() 函数返回当前时刻的电场值。注意激光项必须加在势能里，和库仑势一起作用，不能只加在动能项上，否则电子感受到的电场方向会出错。

强场近似（SFA）是另一条路。它假设电子一旦隧穿就忽略母离子势场对后续运动的影响，所有辐射偶极矩来自连续态波包的加速和回撞。这套源码里的 SFA 模块通常用 Lewenstein 积分实现，把每个时刻的电离振幅和传播相位写成积分形式，最后做傅里叶变换得到谐波谱。SFA 的优点是计算速度快一到两个数量级，适合扫参数找趋势；缺点是电离势较低或者激光强度过高的时候，多电子效应和库仑聚焦明显，SFA 结果会偏离 TDSE。

我的选择习惯是：单原子谱用于验证机制，用 SFA；要对接实验谱形、看干涉结构和长短轨道贡献，必须上 TDSE。这套源码把两条路都给了，中间数据格式统一，可以直接交叉验证。

3. 代码落地：源码包的结构和核心参数怎么读

3.1 源码包里都有什么

拿到这套系统开源的仿真资源后，第一件事不是打开 main.m 就跑，而是先摸清文件结构。它与常见的学术代码包一样，按模块划分了初始化、核心求解器、后处理三层。我建议你按下面的顺序阅读：

. ├── init_params.m # 所有物理参数集中管理：原子、激光、网格 ├── HHG_main.m # 主脚本，控制整个流程 ├── tdsse_solver.m # 含时薛定谔方程求解器 ├── sfa_spectrum.m # 强场近似谐波谱计算 ├── step3_semiclass.m # 三步模型经典轨迹计算 ├── macro_propagation.m # 宏观介质传播与相位匹配 ├── utils │ ├── laser_pulse.m # 余弦或sine载波包络脉冲生成 │ ├── atomic_potential.m # 软核库仑势 │ └── spectral_filter.m # 时间窗函数和谱平滑 └── output ├── HHG_spectrum.dat # 谐波谱输出 └── dipole_acc.dat # 偶极加速度时域数据

文件数不多，每个脚本职责单一。init_params.m 是全局参数的唯一入口，改参数不要直接去求解器里改。这种情况在科研代码里非常常见——主脚本调用子函数时从全局结构体读参数，保持跨模块一致性。

3.2 核心脚本逐段拆：从激光脉冲到谐波谱

先看激光脉冲生成函数。高次谐波对载波包络相位极其敏感，代码里通常会提供 cosine 和 sine 两种载波选择：

function E = laser_pulse(t, params) % params: E0峰值场强, lambda波长, tau脉宽, CEP载波包络相位 omega = 2*pi*299.792458 / params.lambda; % 波长单位micron env = exp(-2*log(2)*(t/params.tau).^2); % 高斯包络 E = params.E0 * env .* cos(omega*t + params.CEP); end

E0 的单位是原子单位（a.u.），1 a.u. 对应约 5.14×10^9 V/m。实际实验中 800nm 激光强度达到 1×10^14 W/cm² 时，E0 大约在 0.05 a.u. 附近。这里的 omega 乘以 299.792458 是为了把波长微米换算成角频率，1 微米波长对应的角频率约为 1.55 a.u.，换算过程要盯紧单位，这个坑我后面会专门讲。

再看谐波谱计算。单原子偶极加速度的傅里叶变换是标准做法：

% 偶极加速度: a(t) = <Psi(t)| -dV/dx + E(t) |Psi(t)> acc = -grad_potential(x).*abs(Psi).^2 + E_field.*abs(Psi).^2; Acc_w = fftshift(fft(acc)); freq = ( (0:Nt-1) - Nt/2 ) / (Nt*dt); % 频率轴，单位a.u. harm_order = freq ./ omega_laser; % 谐波阶次

这里 acc 是每时间步的偶极加速度，对它的傅里叶变换得到频谱。用偶极加速度而不是偶极矩，是因为加速度算符里自然包含激光场项，能压制低频噪声，得到更干净的谐波谱。谐波阶次用频率除以基频得到，整数阶次处出现峰值。平台区和截止位置在这个频谱里直接可读。

3.3 参数怎么调：网格步长、时间步、脉宽和强度

参数设置是仿真能否真实反映物理的关键。以下是我在实际调试中总结的参数范围表，直接对应 init_params.m 里的变量：

需要特别注意的是 dx 和时间步长 dt 的耦合关系。数值稳定性条件大致要求 dt < dx²/2，这是 Crank-Nicolson 和分裂算符格式的共同约束。如果 dx = 0.2 a.u.，dt 必须小于 0.02 a.u.。有人为了省时间把 dt 调到 0.1，结果谐波谱高频段全是噪声。这时候我会先缩小 dx 再调 dt，保证两者的比值满足稳定性条件。

强度参数直接影响 Up 和谐波截止。用 800nm 激光，强度从 1×10^14 升到 3×10^14 W/cm²，截止阶次从 35 阶左右提升到 70 阶以上。调参数时要同步观察电离率，如果中性原子耗尽太快，谐波效率会下降，这时候要降低强度或改用更长波长。

4. 跑通与调参：从阿秒时间尺度的解算到谐波截止验证

4.1 运行环境和主脚本入口

这套 MATLAB 源码不依赖额外工具箱，基础版本的 fft、ode45 和矩阵运算就能跑。主脚本 HHG_main.m 的执行流程是：读参数 → 求基态 → 时间循环求解偶极加速度 → 傅里叶变换得到谐波谱 → 调用宏观模块叠加。直接运行前先在 MATLAB 里设置当前目录为源码包根目录，然后执行：

% 运行整个高次谐波仿真流程 HHG_main;

如果一切正常，会看到输出目录下生成 HHG_spectrum.dat 文件。第一次运行建议用氩原子参数、800nm 波长、2×10^14 W/cm² 强度，这套参数被反复验证过，不容易出幺蛾子。

4.2 输出结果怎么解读：平台区、截止频率和干涉结构

谐波谱输出是双列数据，第一列是谐波阶次，第二列是强度相对值。画出来后你会看到几个典型特征：低阶谐波快速衰减，中间阶次形成一段强度起伏较小的“平台区”，在某个阶次之后强度急剧下降，这个转折点就是截止频率位置。

解读时有个细节要盯紧——平台区的强度并不完全平坦，而是叠加了准周期振荡。这个振荡来自长短两条量子轨道之间的干涉，短轨道产生平台区主要贡献，长轨道在截止附近与短轨道干涉形成调制结构。如果谱图里振荡周期非常规则，说明数值网格和时间步长设置合理；如果振荡杂乱无章，先查边界吸收条件是否生效。

截止频率验证是仿真是否可信的第一道关。把谱图里最后一个明显峰值的阶次记录下来，和理论值 Ip + 3.17Up 对应的阶次对比。偏差在 5% 以内算合理，超过 10% 就要回头查强度单位或波长单位是否换算错。

4.3 宏观传播的初步叠加：为什么单原子谱不够用

单原子谐波谱只反映微观响应，实验探测到的是整个气体喷流中辐射的相干叠加。宏观传播模块在这个资源里做了一件关键的事：把每个空间格点上的单原子偶极辐射按相位匹配条件叠加起来。

% 宏观叠加核心思想：每个格点的辐射相位 = 激光相位 + 几何相位 + 色散相位 for iz = 2:Nz % 累积相位差：泵浦激光和谐波在介质中传输速度不同 phase_mismatch = phase_mismatch + (k_harmonic - q*k_fundamental)*dz; % 相干叠加，带动吸收 E_harmonic(iz) = E_harmonic(iz-1) + dipole(iz) * exp(1i*phase_mismatch); end

宏观模块里最重要的参数是气体密度、介质长度和焦点位置。密度太高，基波在介质中快速耗尽；密度太低，信号强度不足。常见做法是扫描密度和介质长度，找到相位匹配曲线上的峰值位置。这个扫描用 SFA 快速谱预先估算趋势，再用 TDSE 在峰值附近做精细计算，可以节省大量时间。

5. 避坑：高次谐波仿真中反复出现的五个翻车点

5.1 边界反射污染谱线

现象：谐波谱平台区出现周期性尖刺，低阶谐波被淹没在周期性伪峰里。 原因：网格边界没有吸收边界条件，波函数传播到边界后被反射回来，和物理信号形成干涉。这种伪信号在时间上表现为边界往返周期的规律回波。 解决：在网格两端各加 10%~15% 宽度的余弦吸收层。每步传播后把边界区域的波函数乘一个衰减因子。我一般用带肩宽的 cos² 函数，平滑衰减不产生新的反射。加了吸收层后，短时间窗口内的反射污染会明显消失。

5.2 谐波谱高频段的数值色散噪声

现象：截止频率以上本该是噪声的地方出现规则的准周期振荡，看起来像真实的物理信号。 原因：时间步长 dt 相对空间步长 dx 偏大，分裂算符的高阶误差项被放大，高频分量出现非物理色散。这种问题在 400nm 短波长做高次谐波时尤其明显，因为激光频率高，要求更小的时间步。 解决：先固定 dx，把 dt 缩小一半重跑，看谐波谱截止位置附近的噪声是否成倍下降。如果噪声减半，说明是时间步收敛问题；如果不变，则要检查边界反射。收敛性验证是数值仿真最基本的习惯，我每次换激光波长都会强制检查一遍 dt 的收敛窗口。

5.3 SFA 里相位积分的参考时刻选错

现象：SFA 算出的谐波谱截止位置正确，但平台区的干涉结构完全错乱，周期忽密忽疏。 原因：Lewenstein 积分里的 saddle point 时间用的是电离时刻，而不是回撞时刻。这两者之差正好是电子在连续态的传播时间，通常是半个光周期量级，错用会引入一个额外的载波相位误差。 解决：核对代码里 saddle point 方程的求解变量含义。正确的 saddle point 条件是对动作积分的时间端点求极值，得到的是电离时刻 ti 和回撞时刻 tr 两个独立变量。如果实现的是一变量版本，只能捕获半经典收集角附近的结果，平台区干涉细节必然丢失。

5.4 宏观叠加里忽略了电离损耗

现象：密度增加但谐波产出反而下降，谱形还出现奇怪的凹陷。 原因：气体密度高时，前几层的强场电离会产生大量自由电子，这些电子引入负色散，破坏相位匹配条件。更严重的是，电离后的介质对谐波有共振吸收，能量被内层原子吸收后转成热。 解决：宏观传播模块里需要引入电离率模型，每传播一层就更新该层的电离程度和自由电子密度，然后修改波矢失配量。这套源码里用的是 ADK 隧穿电离率公式，和实际实验吻合度较高。密度扫描时要把电离损耗放在模型里一起算，不然扫描曲线和实验结果对不上。

5.5 双原子分子模拟里固定核间距导致谱结构偏移

现象：用双原子分子做高次谐波时，同位素替换后谐波平台变窄，截止往低阶移动。 原因：分子高次谐波里核间距在电离-回撞过程中会增大，这个动态过程对回撞电子能量有显著调制。固定核间距近似忽略了核波包的运动和振动波包的重叠积分，导致长轨道贡献被严重低估。 解决：在 TDSE 求解器里引入核运动自由度，至少用多振动能级展开。这部分计算量会增加一个数量级，但分子系统本来就不是单原子可比的。如果只是验证机制，可以先用固定核间距把电子动力学弄对，再逐步引入核运动。

6. 验证技巧：三步确认结果可信，再试着合成阿秒脉冲

仿真跑完不等于结果能用。我拿到一套新源码或新参数时，会强制走三个验证步骤，每一步都能筛掉大部分隐藏错误。

第一步是对齐截止频率。用理论截止阶次 N_cutoff = (Ip + 3.17Up)/ω_laser 和仿真谱图对比，误差限设在 5% 以内。这个验证最关键的是单位换算——从实验的微米、瓦特每平方厘米、飞秒换算到原子单位的流程，任何一处出错都会让截止位置大幅偏离。很多看起来“物理不对”的仿真结果，最后查到都是单位系数错了一位。

第二步是检查平台区的强度比值。仿真的单原子 HHG 平台区谐波强度从前 5 阶到截止位置应该只下降一到两个数量级，如果中间出现超过三个数量级的跌落，多半是数值色散或边界反射在作怪。实验上测到的平台会有介质的吸收和色散调制，但单原子仿真的平台比较平坦。

第三步是改变激光强度后看截止随 Up 的线性变化趋势。强度从 I 提到 2I，Up 翻倍，截止阶次应该向右移动约 3.17 倍 Up 对应的阶次。这一步能检验 SFA 或 TDSE 实现里激光强度项的标度关系是否正确。用三个不同的强度点做线性拟合，斜率在 2.8~3.5 之间都算合格，超出这个范围要回头查电场峰值 E0 的计算。

三步都过了之后，进阶玩法是合成阿秒脉冲。把高通滤波后的谐波谱做逆傅里叶变换，就能在时间域得到一串阿秒脉冲序列。谐波谱里相位被压缩的频段越宽，阿秒脉冲越窄。你可以尝试改变 CEP 参数在 0 到 π 之间扫描，观察阿秒脉冲是单脉冲还是脉冲串，这直接对应实验里用载波包络相位锁定来控制孤立阿秒脉冲的原理。

% 合成阿秒脉冲：对截取的高次谐波谱做逆傅里叶变换 harm_filtered = spectrum .* (freq > freq_cutoff); atto_pulse = abs(ifft(harm_filtered)).^2; % 时间轴换算: 1 a.u. = 24.2 as t_as = ( (0:Nt-1) - Nt/2 ) * dt * 24.2;

这里 freq_cutoff 要取截止位置向内 5~10 阶的位置，因为截止处的谐波相位不平滑，强行包含会引入时间旁瓣。做完这些验证和脉冲合成，就能确认这套源码在你自己的参数范围内是可复现、可扩展的。从那以后我每次拿到一套新的高次谐波仿真资源，都会先跑这三个验证步骤再碰宏观模型，省掉的排查时间远比多花的几分钟运行时间值。希望这些经验能帮你的仿真之路少一点无谓的调试，多一点扎实的结果。

本文还有配套的精品资源，点击获取