简介：面向高次谐波产生（HHG）数值模拟的开放源码 MATLAB 仿真包，聚焦强激光场与气体靶相互作用中真空紫外/极紫外辐射的产生机制，适合等离子体物理、超音速流动力学及量子原子响应方向的研究者与学生使用。代码库提供图形界面，可交互调整输入参数并观察电离、相位匹配和偶极响应等关键物理量；所有计算过程会自动保存为 txt 文件，再次运行同类参数时速度更快。资源共 451 个文件，其中 69 个 m 脚本为核心仿真代码，371 个 txt 文件为中间/结果数据，另有少量 PDF、PNG、FIG 辅助说明，整体大小约 29.75MB。仿真覆盖喷嘴出口超音速气流、等离子体动态离子、自由电子量子响应、相位匹配与吸收等完整物理链条，是研究 HHG 全过程的较系统实现。已有 716 人学习/下载，适合需要深入理解高次谐波产生机理或在这一方向扩展 MATLAB 算法的使用者。

1. 高次谐波数值模拟：一套把 TDSE、SFA 和宏观传播串起来的 Matlab 源码

高次谐波产生（HHG）的数值模拟，是阿秒物理里最绕不开的一道坎：经典三步模型能解释截止能量，但真实实验里那条平台区的谐波谱，需要同时解量子力学含时薛定谔方程（TDSE）、做强场近似（SFA）对照，再把超声速气体靶的密度剖面代入宏观传播方程。这套 Matlab 源码正好把这几段串成一条完整链路。它解决的是“单原子响应算对了，一加相位匹配就翻车”的痛点。适合刚入门强场物理、需要在论文里放一张可信谐波谱的硕博生，也适合做气体靶实验、想快速估算谐波产率的工程师。你不需要重写求解器，主要工作是改对参数、边界条件和网格尺寸。

2. 物理模型与选型：高次谐波为什么不能只靠经典三步模型

2.1 三步模型只给截止能量，给不了平台区频谱

很多第一次跑 HHG 仿真的人会问：三步模型也能画出一条漂亮的截止曲线，为什么源码里还要放一个 TDSE 求解器？我的建议是，三步模型只适合做量级估算。它把“电离—加速—复合”当作一条经典轨迹：电子在激光场峰值附近隧穿电离，之后被电场加速，等电场反向时以一定概率回到母离子复合，放出一个高能光子。最大回撞动能是 3.17Up，Up 是有质动力能，所以截止光子能量写为：

[ E_{\text{cutoff}} = I_p + 3.17U_p ]

其中 (I_p) 是电离能，(U_p = \frac{e^2 E_0^2}{4m_e \omega^2})。这个公式能帮你在改波长、改强度时快速估算谐波能推到多少阶，实测误差通常能控制在百分之十以内。但它解释不了两件事：第一，为什么谐波谱上会出现清晰的奇次阶峰，而不是连续谱；第二，为什么平台区的谐波强度并不单调下降。这两个问题都来自连续态波包在激光场中的量子干涉，经典轨迹模型没有相位信息，自然也复现不了平台区。所以，只要你的目标是发论文、对比实验产率，终归要回到量子模型。

在这套源码里，三步模型被放在预处理模块里，专门用来快速扫描截止边界。真正干活的是 TDSE 和 SFA。我的习惯是先用三步模型过一遍参数，确定强度范围有没有可能把谐波推到目标阶次，再开 TDSE，避免一上来就在不该花的网格上烧时间。

2.2 含时薛定谔方程：长度规范、软核势与分裂算符

单原子响应的核心求解器，是在单电子近似下解三维含时薛定谔方程：

[ i\frac{\partial}{\partial t}\psi(\mathbf r, t) = \left[-\frac12 \nabla^2 + V(r) - \mathbf E(t)\cdot \mathbf r\right]\psi(\mathbf r, t) ]

这里用到的是长度规范。对线偏振光，体系绕激光偏振轴保持轴对称，磁量子数 (m=0) 守恒，所以波函数可以按角动量 (l) 展开。源码里常见的做法是把径向格点取成均匀网格，每个 (l) 分量单独推进，再用分裂算符处理动能项和势能项。分裂算符的典型形式是：

[ \psi(t+\Delta t) \approx e^{-i \hat T \Delta t/2} e^{-i V \Delta t} e^{-i \hat T \Delta t/2} \psi(t) ]

动能项在动量空间用 FFT 完成，势能项在坐标空间做乘法，这样做的好处是每步误差在 (\Delta t^3) 量级，比直接用 Crank-Nicolson 快很多。要注意的是，势函数 (V(r)) 不能直接用氢原子库仑势，因为惰性气体原子的基态并不简单地是 (-1/r) 形式。源码里用的是软核势：

[ V(r) = -\frac{1}{\sqrt{r^2 + a^2}} ]

其中参数 (a) 要由目标气体的实际电离能反推。这一步如果省略，出来的截止位置会偏。我一般会先用一个一维定态求解器扫描 (a)，把基态能量调到和实验 (I_p) 一致后再做全维 TDSE。另外，长度规范里的 (-E(t)r) 项在大 (r) 处会迅速累积相位，所以网格边界必须加吸收层，否则反射波包会污染谐波谱。这个点后面避坑章会细说。

TDSE 的优点是物理完备：它天然包含束缚态、连续态、库仑修正和所有量子干涉，只要网格够细、时间步够小，结果就能当作理论基准。缺点也明显，计算量上天。对三维轴对称问题，径向格点 1000 到 2000，角动量 (l) 取到 8 左右，一个 30 飞秒的 800 纳米脉冲，推进步数大约几万步，单核跑需要几个小时，这也是源码里同时放 SFA 的原因。

2.3 强场近似 Lewenstein 模型：快但需要校准

SFA 的做法是只保留基态和连续态，忽略中间激发态，并且把连续态当成平面波处理。Lewenstein 模型给出的偶极矩是一个含时积分，本质上是电子在连续态里的传播子，可以写成：

[ d(t) = i \int_0^\infty d\tau \left(\frac{\pi}{\epsilon + i\tau/2}\right)^{3/2} E(t-\tau) d_c^* E(t) d_g \exp[-i I_p \tau] \exp[-i S(t, \tau)] ]

其中 (S(t,\tau)) 是作用量，(d_g)、(d_c) 是偶极矩阵元。这个公式跑起来非常快，一个波长扫描几分钟就能出结果。SFA 在长波长、高强度的条件下非常准，因为库仑势对返回电子轨迹的修正相对变小。但在短波长、低强度或低电离能情况下，SFA 和 TDSE 的差距会拉开，尤其是低阶谐波的相位，SFA 经常给出明显偏差，这会影响后续宏观传播的相干叠加。

在这套源码里，SFA 和 TDSE 的输出格式统一，都是时间域的偶极矩，后续 FFT 和宏观传播模块不需要区分数据来源。我会把它当“快照”工具：先跑 SFA 扫描强度、波长参数空间，把可疑区域标出来，再在 TDSE 里精算几个关键点。如果你的课题只是做平台区产率趋势，不纠结绝对相位，SFA 基本够用。

2.4 宏观传播：超声速气体靶密度剖面如何进入谐波相位匹配

单原子偶极矩乘上气体密度，再沿着整个介质累积场强和相位，才得到实验能测到的谐波信号。这个过程用麦克斯韦方程在慢变包络近似下逐阶推进：

[ \frac{\partial A_q}{\partial z} = -i\frac{\omega_q}{2 c} \left[ \frac{\omega_q}{c} \chi_q^{(1)} A_q + \frac{\omega_q}{c} \chi_q^{(3)} A_1^3 + \cdots \right] ]

简化处理时，每个谐波阶次 (q) 单独传播，最重要的物理量是相位失配：

[ \Delta k_q = k_q - q k_1 = \Delta k_{\text{neutral}} + \Delta k_{\text{plasma}} + \Delta k_{\text{geom}} ]

中性原子色散让谐波和驱动光走得不一致，电离出来的自由电子（也就是等离子体）贡献一个负色散项，几何项则来自聚焦光的高斯相移。当这三项之和接近零时，谐波产率出现相干增长；否则谐波在介质里一边产生一边抵消，最后探测到的强度大打折扣。

这时候超声速流物理就进来了。实验上为了降低气体密度梯度，常用超声速喷嘴形成准均匀靶。但超声速喷流的密度剖面并不是简单矩形，它有边界层、激波和轴向密度平台。源码里预留了一个密度剖面输入接口，你既可以填解析表达式，也可以直接读流体模拟导出的离散数据。我的经验是，密度剖面取平均密度会出大问题，相位失配对轴向梯度极其敏感，尤其是电子密度非线性增长的高强度区，一个格点内 (\Delta k L) 就可能跳过相干峰，后面避坑章我会再展开。所以宏观传播网格要和密度剖面的变化尺度匹配，不能随手用二三十个格点敷衍。

3. 把源码跑起来：配置参数、主流程和结果验证

3.1 功能模块怎么拆

拿到源码包之后，不要急着在 Matlab 里双击运行。先把模块分清楚。这套源码通常由四块组成：

• 单原子响应模块：TDSE 求解器、SFA 计算器、软核势参数标定。
• 宏观传播模块：沿 z 方向逐层推进驱动场和谐波场，支持自定义气体密度剖面。
• 后处理模块：傅里叶变换、谐波谱画图、截止位置提取、时频分析。
• 配置脚本：集中放波长、强度、脉宽、网格大小、气体靶长度等参数。

我习惯把它们当成四个独立的“实验台”：改配置脚本，跑单原子模块，存下偶极矩后再跑宏观模块，最后统一后处理。这样每次调参数不需要从头算，尤其是 TDSE 最烧时间的一步。源码里模块之间用 .mat 文件传递数据，单原子响应算完一次，可以反复用在不同的宏观条件里，这个设计很实用。

3.2 关键参数：波长、强度、脉宽、网格和气压

所有物理量建议统一到原子单位制，否则换算错误很难排查。下面给出一份我常用的示例配置，和源码里的默认值基本一致：

配置脚本的核心是单位换算。下面这一段就是我在每个新环境下都会替换的头部：

% config_hhg.m % 驱动波长 lambda_nm = 800; lambda = lambda_nm * 1e-9; omega = 2*pi*2.998e8 / lambda; % 角频率，SI omega_au = omega / 4.134e16; % 角频率，原子单位 % 峰值强度与电场幅度 I_peak = 1.2e14; % W/cm^2 E0 = sqrt(I_peak / 3.509e16); % 峰值电场幅度，原子单位 % 脉宽与传播时间 tau_fs = 30; tau_au = tau_fs * 1e-15 / 2.418e-17; % 时间，原子单位 t_total = 4 * tau_au; % 传播总时长 dt = 0.05; % 时间步长 Nt = ceil(t_total / dt); % 径向网格 r_max = 150; Nr = 1024; l_max = 8;

这里要注意omega_au的计算：原子单位下时间单位是 2.418e-17 秒，角频率要除以这个数。很多新手直接用 SI 角频率去做量子力学计算，结果所有参数都飘了。tau_au * 4是因为脉冲前后有很长的零场区域，需要让波包充分演化，避免人为截断导致频谱泄漏。如果你只算 1~2 个周期的极短脉冲，可以改成 2 倍脉宽，但要多周期脉冲，4 倍比较稳。

3.3 一次完整复现：从驱动场到谐波谱

核心脚本大致是下面这个流程。先建立一个时间轴和驱动电场包络，然后调用单原子模块得到偶极矩，再加窗函数做 FFT 得到谐波谱。

% run_hhg_single_atom.m t = (0:Nt-1) * dt; mid = t(end)/2; % 高斯包络，峰值在 mid env = exp(-((t - mid) / (tau_au/1.177)) .^ 2); E_drive = E0 * env .* cos(omega_au * (t - mid)); % 调用 TDSE 求解器 dipole = tdse_solve(E_drive, Ip_au, r_max, Nr, l_max, dt); % 加汉宁窗后算谐波谱 Nfft = Nt; dip_window = dipole .* hann(Nfft).'; spec = fft(dip_window); freq = (0:Nfft-1) / (Nfft * dt); q_axis = freq / omega_au; % 画图，去掉直流分量 semilogy(q_axis(2:Nfft/2), abs(spec(2:Nfft/2)).^2); xlabel('Harmonic order q'); ylabel('Intensity (a.u.)');

这段脚本里的tdse_solve是源码包里的核心求解函数，名字可能略有差异，但参数列表基本就是电场序列、电离能、径向网格参数和时间步长。加汉宁窗这一步非常关键，不窗的话低频连续背景会拖得很长，谐波峰会被泄漏信号淹没。q_axis是谐波阶次，它等于频率除以驱动光频率。画图时我习惯看强度对数坐标，平台区会直观地呈现为一段平缓下降的斜坡。

如果你更想先用 SFA 试跑，只需要替换调用行：

dipole = sfa_dipole(E_drive, Ip_au, lambda_nm);

SFA 没有网格参数，速度至少快两个数量级，代码包里的接口是一样的。唯一要注意的是，SFA 对短脉冲的波形比较敏感，需要保证电场序列里没有直流偏置。如果E_drive的包络在时间轴上不对称，你会看到一个虚假的偶次谐波平台，那不是物理，是数值问题。

3.4 验证结果：先看奇偶性和截止位置

跑完第一张谐波谱，先别急着调参数进宏观模块，先做三个静态检查。

第一个检查是奇偶性。对于中心对称气体介质，在长脉冲、低强度下，谐波谱应当只有奇次峰。如果偶次谐波明显存在，通常是驱动电场的时间包络不对称，或者时间轴上直流分量没清零。偶次谐波在实验里也可能通过其他机制出现，但在单原子仿真里，它大概率是数值对称性破缺。这时候我会检查sum(E_drive)是不是接近零，而不是去调吸收边界。

第二个检查是截止位置。用公式 (E_{\text{cutoff}} = I_p + 3.17U_p) 估算一下，再在谱图上找平台结束的位置。如果偏差超过百分之十，优先怀疑单位换算。特别是omega_au或E0这两个量，一个符号错整个谱就废了。

第三个检查是网格收敛。把 (N_t) 翻倍、(r_{\max}) 增加 20%，重算一次。谐波峰的高度变化应该在小数点后两位以内，如果明显上升，说明之前网格没有收敛。这一步很多人会跳过，等到宏观传播跑完才发现产率异常，回头排查成本极高。我现在每换一组强度或波长，都强制自己先跑收敛测试，再用固定网格做后续扫描，算是吃够了“网格不收敛就发图”的亏。

4. 高次谐波仿真避坑：四个最容易翻车的数值细节

4.1 边界反射假信号

现象：TDSE 跑完后，谐波谱低阶部分出现一条很宽的连续背景，甚至能看到等间隔的小振荡，把平台区起点都遮住了。更典型的是，连续背景随着传播时间增加而变强，让人误以为物理上产生了连续辐射。

原因：波函数在激光场中电离出的电子向大 (r) 方向运动，碰到r_max硬边界后被反射回来，反射波包与仍在束缚态区域的波函数发生干涉。这个干涉体现在偶极矩上就是非物理的低频振荡，FFT 之后变成连续背景。

解决：在网格外部加吸收区。常用做法是在 (r > 0.8 r_{\max}) 的区域乘一个掩膜 (M(r))，形状取 (M(r)=\cos^{1/8}(\pi (r-r_{\text{abs}})/(2(r_{\max}-r_{\text{abs}}))))，每个分裂算符子步结束都乘一次。掩膜强度不能太陡，否则反射变成散射；也不能太平，否则吸收不完。我的调法是先设吸收区占比 20%，看反射背景是否消失，如果还在，就把占比提到 30% 并减小掩膜指数。

4.2 时间步长导致频率混叠

现象：在谐波谱的高频端，截止峰附近突然出现一个看起来像“镜像”的上升结构，而且峰位和理论截止完全不匹配。把谱图横轴拉长后会发现，这个结构本质上是真实谱以奈奎斯特频率为轴的镜像。

原因：时间步长 (\Delta t) 太大，最高可分辨频率是 (f_{\max}=1/(2\Delta t))。如果谐波阶次对应的频率超过这个上限，频谱就会折叠混叠。800 纳米驱动的第 50 阶谐波，频率约为 50×0.057=2.85 a.u.，那 (\Delta t) 至少要小于 (1/(2\times 2.85)\approx 0.175) a.u.。实际中包络上升沿还有更高瞬态频率，0.05 a.u. 才保险。

解决：把时间步长减小一半，重算谐波谱，看截止峰位置是否稳定。如果稳定，说明之前只是临界；如果还在动，继续减半。要注意的是，单纯把Nt翻倍但保持dt不变是没用的，因为你并没有提高最高采样频率，只是把总时长拉长。改参数时先固定时间窗，再改采样率，才是有效收敛测试。

4.3 软核势基态能量对不上实验电离能

现象：SFA 和 TDSE 算出的截止阶次不一致，差出一两个 eV。更隐蔽的是，TDSE 的平台区斜率明显偏陡，和实验光谱趋势对不上。

原因：软核势参数 (a) 没标定。以 (V(r) = -1/\sqrt{r^2+a^2}) 为例，基态能量和 (a) 的取值强相关。如果你随便取 (a=1)，算出来的“基态”可能比实际目标气体的电离能高出几十个百分点。这样电子在势阱里被束缚得更紧，电离率、返回轨迹和谐波相位全部偏离。

解决：源码里通常自带一个一维定态求解器，用它扫 (a)，让最低本征能量等于目标气体的实验电离能 (I_p)（原子单位）。比如想模拟某个惰性气体，就找到对应 (I_p) 的 (a) 值，再把这个 (a) 代回到三维 TDSE。标定完成的标准是：基态能量误差小于 0.01 a.u.，约 0.27 eV。如果源码没带定态求解器，可以在主脚本里加一段 Lanczos 迭代，十几行就能算出来，不值得在这上面省时间。

4.4 宏观传播 z 步长与密度剖面的不匹配

现象：宏观模块算出来的谐波产率随气压变化不是平滑曲线，而是出现无规律抖动；或者实验上某段气压区间产率上升，模拟里却下降。

原因：传播步长 (\Delta z) 太大，相位失配在一个步长里就累积了超过 (\pi) 的相位，数值上把相干增长平均掉了。超声速气体靶的密度剖面存在平台与边界层的快速变化区，平均密度近似会抹掉相位匹配的敏感点。尤其是电离率随 z 变化时，等离子体色散项逐层不同，粗网格完全捕捉不到。

解决：宏观传播模块里强制用自适应 z 步长，每步都检查 (\Delta k(z)\Delta z < 0.2)。如果你的代码不支持自适应，就把 (z) 方向网格设到几千层，再配合线性插值读取密度剖面。还有一个容易被忽略的点：背压升高时，自由电子密度非线性增长，(\Delta k_{\text{plasma}}) 会迅速超过中性色散项，相位失配从正变负。这个拐点想抓准，密度剖面必须按实际超声速喷嘴的参数给定，不能为了省事用矩形靶。

5. 进阶：把谐波频谱还原成阿秒脉冲序列

5.1 短时傅里叶重建

谐波谱本身是频域信息，但做阿秒实验的人最终关心的是时域脉冲序列。高频谐波之间相位锁定，叠加起来就是一串阿秒脉冲。这个变换我用短时傅里叶变换来做：

% stft_reconstruct.m Nwin = 256; Nhop = 32; Nfft = 1024; [S, f_axis, t_axis] = spectrogram(dipole, hann(Nwin), Nhop, Nfft, 1/dt);

spectrogram需要信号处理工具箱，如果没有，手写一个滑动窗 FFT 也可以。重点不是工具，而是物理：首先在频谱上圈出要重建的谐波范围，例如第 15 到第 29 阶，然后把范围外的频谱直接置零，再做逆傅里叶变换：

q_min = 15; q_max = 29; mask = zeros(size(q_axis)); mask(q_axis >= q_min & q_axis <= q_max) = 1; E_att = ifft(fft(dipole) .* mask); t_fs = t * 2.418e-2; % 时间从 a.u. 转成飞秒 plot(t_fs, real(E_att));

检查重建结果时，我会重点看两个值：第一，主脉冲间隔是不是接近驱动光半周期，800 纳米驱动对应约 1.33 飞秒；第二，脉冲宽度是不是远小于驱动光周期。如果脉冲间隔对不上，说明你选的谐波范围内相位关系没有锁定；如果重建出连续振荡，说明频谱滤波后的窗口里包含过多连续背景。我第一次跑这步时，滤波窗口取得太宽，把 14 阶到 31 阶全拉进去，结果混进一段低阶连续谱，脉冲序列完全看不出结构。从那以后，每次重建前我都会先看一眼谐波谱的相图，只选相位曲线接近线性的那一段。

5.2 “SFA 快扫，TDSE 精算”的日常节奏

这套源码用久了，我总结出最省时间的工作流：凡是需要扫参数空间的，先用 SFA 跑一遍，把强度、波长、脉宽画成产率色图，标出有希望的参数区；再在选区里挑两三个点，开 TDSE 精算，用来校准宏观传播的绝对产率。SFA 和 TDSE 完全对齐的时候很少，但趋势一致就够用了。宏观传播模块里密度剖面固定住之后，批量扫描气压和靶长，不用反复动单原子模块。这样每次下完源码，一个下午就能把一篇论文用的参数图画得差不多。从那以后，我每次改驱动波长，都强制先跑一遍 5 分钟的 SFA 快扫，再决定是否开动 TDSE 全网格，希望帮到你。

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