简介面向高次谐波产生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.17UpUp 是有质动力能所以截止光子能量写为[ 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) ]这里用到的是长度规范。对线偏振光体系绕激光偏振轴保持轴对称磁量子数 (m0) 守恒所以波函数可以按角动量 (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 关键参数波长、强度、脉宽、网格和气压所有物理量建议统一到原子单位制否则换算错误很难排查。下面给出一份我常用的示例配置和源码里的默认值基本一致参数示例值说明驱动波长800 nm换算成原子单位约 0.057峰值强度1.2e14 W/cm²1 a.u. 强度约 3.509e16 W/cm²脉宽30 fs功率半高宽约 1240 a.u.时间步0.05 a.u.约 1.2 as保证高频谐波不混叠径向格点150 r_max1024 点覆盖经典电子最大偏移量角动量截断l_max 8线偏振下已可收敛到高阶谐波气体靶长500 μm典型超声速喷嘴长度气压或密度1050 mbar决定电子密度和相位失配配置脚本的核心是单位换算。下面这一段就是我在每个新环境下都会替换的头部% 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.0572.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^2a^2}) 为例基态能量和 (a) 的取值强相关。如果你随便取 (a1)算出来的“基态”可能比实际目标气体的电离能高出几十个百分点。这样电子在势阱里被束缚得更紧电离率、返回轨迹和谐波相位全部偏离。解决源码里通常自带一个一维定态求解器用它扫 (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 全网格希望帮到你。本文还有配套的精品资源点击获取