光学频率梳这个词这几年在光子学里几乎是绕不开的。无论是光通信里的多波长光源、微波光子学里的低噪声振荡器还是精密光谱测量大家都会提到它。而微环谐振腔作为片上实现光学频率梳的核心平台之一更是圈内人默认的高起点。我最初上手这个方向时用的就是Matlab做全链路仿真从单频泵浦光在微环里一点点演化成几十根等间隔梳齿整个过程非常“上瘾”但踩坑也不少。这篇内容就把我自己的实现路径、理论基础、关键参数和排查方法完整捋一遍写给正在或者准备做微环频率梳仿真的朋友参考。1. 为什么是微环谐振腔先弄懂它在光学频率梳里的位置微环谐振腔本质上是集成光子学里最常见的结构之一一根直波导加上一个闭合的环形波导光通过倏逝波耦合从直波导进入环内在环里绕圈共振。它的名字听起来很学术但结构上可以类比成一个“光学蓄水池”——光在环内不断循环积累能量只要满足谐振条件腔内功率会远远高于输入功率。这种功率增强效应是后面产生非线性效应的根基。1.1 微环的基本结构与核心参数标准的微环结构有两种一种是全通型all-pass一根直波导加一个环只有一个耦合区另一种是上下载型add-drop两根直波导加一个环两个耦合区。频率梳仿真里常用的通常是全通型结构因为它结构简单、损耗路径清晰便于分析腔内功率积累过程。和微环相关的几个核心参数决定了整个仿真的走向自由光谱范围FSR环内相邻谐振峰之间的频率间隔由环周长和群折射率决定公式是 FSR c / (ng · L)单位是Hz。半径越大FSR越小。这个参数直接决定了梳齿之间的间距。品质因子Q值衡量谐振腔储能能力的指标Q值越高腔内光子寿命越长功率增强越明显。Q值受耦合损耗、波导传输损耗等因素约束。功率增强因子腔内功率和输入功率的比值近似正比于Q值。这是微环能作为非线性平台的核心原因。耦合系数κ直波导和环之间能量交换的强度决定了腔的加载Q值和外耦合率。仿真中很多人上来就写方程结果参数一团乱最后出了结果也不知道对不对。我的建议是先固定物理结构参数——选一个常见氮化硅Si₃N₄或者硅Si波导平台把有效折射率、群折射率、非线性折射率这些材料参数确定下来再推导耦合方程。否则Matlab脚本里每个变量都悬空后面调参就是灾难。1.2 为什么不能只用“线性腔模”解释频率梳线性状态下微环的频率响应就是一系列等间隔的洛伦兹峰泵浦光耦合进某个谐振峰后输出仍是单频光。频率梳的产生需要“新频率”出现而这必须依赖非线性过程。微环频率梳的核心物理机制是四波混频FWM这是由材料的三阶非线性极化率χ⁽³⁾引起的。两个泵浦光子湮灭、产生一红一蓝两个新光子信号光与闲频光只要满足能量守恒和相位匹配约束能量就会不断从泵浦光转移到新的频率分量上。随着腔内光强增大新产生的频率又作为新的“泵浦”继续产生更多频率成分最终形成跨越数个谐振峰的光学频率梳。这也是为什么仿真频率梳一定要用非线性模型线性传输方程只能给你梳状滤波器的透射谱永远得不到梳状光源的光谱。你在Matlab里既要处理色散、损耗、耦合这些线性效应也要处理克尔非线性带来的相位调制和频率搬移这两类效应在腔内往返过程中是同时发生的。2. 从物理机制到Matlab模型核心方程与仿真原理频率梳的时域仿真圈内最常用的模型是Lugiato-Lefever方程LLE。它本质上是一个带边界条件的非线性薛定谔方程专门用来描述被动微环谐振腔内连续波CW泵浦下的非线性动力学。可以说掌握了LLE就掌握了微环频率梳仿真的“主心骨”。2.1 Lugiato-Lefever方程的含义LLE的归一化形式大概长这个样子dE(t,τ)/dt -(1 i·θ)E i·|E|²E i·(∂²E/∂τ²) F这行式子看起来复杂但拆开看每一项物理意义非常清晰左边 dE/dt 是腔内慢变幅度随时间慢时间域的演化t代表慢时间也就是在腔内的往返次数累计。右边第一项 -(1i·θ)E 中的实部1是腔内损耗项虚部θ是泵浦光与谐振峰的失谐量。第二项 i·|E|²E 是克尔非线性项对应四波混频和自相位调制效应这是产生新频率的根源。第三项 i·(∂²E/∂τ²) 代表群速度色散GVDτ是快时间坐标对应一个腔周期内的波形细节。这个色散项决定了梳齿包络的形状和孤子形成的可能性。最后一项 F 是连续波泵浦项。在实际仿真中你不需要从第一性原理推导这个方程但要能写出每一项的物理量纲理解它在“腔往返”这个物理图像里的位置。很多人拿着别人代码跑通了但不知道每个系数对应实验中哪个参数一遇到异常结果就抓瞎。2.2 为什么用分步傅里叶法求解LLE是非线性偏微分方程线性项损耗、色散、失谐和非线性项克尔效应在物理上是耦合的。直接做时域有限差分可以求解但数值稳定性和效率都差一些。工程里最流行的方案是分步傅里叶法Split-Step Fourier MethodSSFM。核心思想很朴素在一个极小的时间步长内先假设非线性效应单独作用再用傅里叶变换把线性效应放到频域里单独处理最后再把它们拼起来。因为色散项在频域里只是一个乘法运算效率远高于时域差分。用生活化类比来说这就像你做一道菜——先腌制非线性步再加热线性步但加热过程中腌制的效果也在同步进行只要每一步时间足够短顺序误差就可以忽略。在Matlab里实现SSFM主要的循环结构是初始化腔内场E通常从一个小噪声或弱泵浦场开始。非线性步E E · exp(i·γ·|E|²·dz)直接在时域完成。FFT变换到频域。线性步E E · exp(线性算子·dz)包含损耗和色散。IFFT回时域检查收敛或粒子数守恒。这里要注意Matlab的fft和ifft在处理光场包络时频率轴的顺序是f(1)对应直流然后把负频率放在后半段。很多人第一次画光谱时发现左右不对称就是没做fftshift。这是入门阶段最典型的低级错误但不仔细看真的很容易忽略。2.3 归一化与无量纲化的作用LLE里最关键的设计决策是归一化。物理参数波导损耗α、非线性系数γ、色散系数β₂、环周长L数量级差异极大非线性系数可能是W⁻¹·m⁻¹量级而色散可能是ps²/m量级。直接代入计算会让Matlab里的数值巨大或极小极易导致误差和溢出。归一化的思路是把时间、长度、场强都缩放到无量纲单位。常见的做法是引入一个参考时间尺度T₀通常取快时间窗口内的某个特征值把快时间τ归一化为τ/T₀把腔传播距离归一化为环周长L场强E则用非线性系数和损耗的比值来归一化。这样方程里只剩下几个无量纲系数失谐θ、归一化色散、归一化泵浦幅度F。这一步做好的好处是同一套代码可以扫描任意物理参数组合只需要改归一化系数不需要改求解核心。我自己的习惯是先在脚本开头写一个“归一化计算”段落把所有物理参数转换成无量纲参数再用这些参数去跑SSFM。调试时如果输出异常优先检查归一化是否正确而不是怀疑求解算法本身。3. Matlab代码实现从空白脚本到能出梳状谱仿真环境的搭建虽然不难但很多细节会影响效率甚至影响正确性。我在做这个项目时用的是较新版本的MatlabWindows平台上跑中途也经历过换机器、重装环境、工具箱缺失的问题。这里把环境准备、参数设置和核心代码骨架分开说。3.1 环境准备版本、工具箱和安装要点Matlab本身的分发包很大早期版本装起来比较折腾。现在的安装流程已经友好很多但有几个点还是值得提醒版本选择如果你手上是新电脑可以直接用较新的Matlab版本比如R2025b或R2026b新版本对GPU加速、并行工具箱的支持更好。老版本跑SSFM也可以只是大规模参数扫描时会慢得让人崩溃。工具箱要求核心仿真只需要基础Matlab和Signal Processing Toolbox用于滤波和信号分析。如果要画漂亮的光谱图和动态演化图还需要DSP和基本绘图功能——这些在标准安装包里都有。Optimization Toolbox可用于自动拟合参数但不是必须。许可证和激活安装时经常遇到许可证问题比如License Manager报错或不识别License文件通常需要检查Host ID是否匹配、许可证路径是否设置正确。网上很多教程提到下载和安装步骤其实关键就两条安装路径不要带中文激活时需要管理员权限。遇到license.lic的Host ID对不上多半是网卡顺序变化或虚拟机环境导致的重新生成License文件就行。并行计算如果CPU核心多建议启用Parallel Computing Toolbox的parfor做参数扫描。频率梳仿真本质上是一个大参数循环扫描泵浦功率、失谐、色散等参数每个参数点都要跑几千个腔往返串行计算很浪费时间。在Matlab里跑的脚本结构我习惯分成四个文件参数定义脚本params.m、归一化计算脚本normalize.m、主仿真循环main_sim.m、绘图与数据分析脚本plot_results.m。这样改参数时不用反复翻代码出图也更快。千万不要把所有代码塞在一个脚本里后期调参时你会后悔的。3.2 参数设定从一个“能出梳”的典型值开始对于仿真新手最大的坑是参数设置不合理导致永远无法产生频率梳。我的建议是先找到文献里已经被验证过的参数范围然后在附近微调。以下是一组典型的氮化硅微环参数基于常见文献和商用SOI/SiN平台指标参数数值说明微环半径 R40 μm决定FSR半径越小FSR越大环波导截面0.5 μm × 0.8 μm常见脊波导尺寸有效折射率 n_eff1.98与波长相关需用模式求解器算群折射率 n_g2.1决定FSR仿真中常近似非线性折射率 n₂2.4×10⁻¹⁹ m²/W氮化硅材料典型值波导损耗 α0.5 dB/cm代表损耗水平功率增强倍数10~50由Q值决定估算用泵浦功率 P_in100 mW ~ 1 W低于阈值则无法起振泵浦波长1550 nmC-band常见泵浦实际上在仿真里你并不会直接输入n₂和α而是换算成γ非线性系数和腔损耗系数。如果嫌换算麻烦可以用近似公式γ ≈ 2π·n₂/(λ·A_eff)A_eff为有效模面积损耗系数则需要把dB/cm换算成m⁻¹。这里用到了Matlab的数据处理能力写个小函数批量换算体力活很值得一提频率梳仿真大多数时间不是耗在求解上而是耗在参数单位换算上。3.3 核心代码骨架分步傅里叶法的Matlab实现下面是经过简化但仍可直接复现的Matlab代码框架我用它跑出了初步的CI梳状光谱。代码里保留了关键步骤注释方便你对照上面的LLE解释。% 主仿真脚本流程简化版 % 参数初始化 lambda0 1550e-9; % 泵浦波长 c 3e8; % 光速 R 40e-6; % 环半径 L 2*pi*R; % 环周长 n_g 2.1; % 群折射率 FSR c/(n_g*L); % 自由光谱范围 (Hz) % 非线性/损耗参数换算后 alpha 0.5; % dB/cm alpha_m alpha/8.686/0.01; % 换算为 1/m gamma 2.0; % W^-1 m^-1需结合材料与模场面积估算 P_pump 0.5; % 泵浦功率 W % 快时间窗口与网格 T_window 20/FSR; % 时间窗口覆盖20个FSR N 4096; % 网格点数 tau linspace(-T_window/2, T_window/2, N); d_tau tau(2)-tau(1); omega 2*pi*fftshift(-N/2:N/2-1)/(N*d_tau); % 频域角频率轴 % 色散参数以二阶色散为主单位换算后代入 beta2 -20e-24; % s^2/m择优异常色散负号 disp_operator -1i*0.5*beta2*(omega.^2); % 频域色散算子 % 内腔场初始化小噪声 E sqrt(2e-3)*exp(-(tau/T0).^2); % 弱高斯种子 噪声 % 实际中常用白噪声初始化更好避免包络偏差 % 腔往返迭代 N_roundtrip 2000; % 往返次数 for k 1:N_roundtrip % 非线性步时域相位调制 E E .* exp(1i*gamma*abs(E).^2*L/2); % 线性步频域色散和损耗 E_hat fftshift(fft(E)); E_hat E_hat .* exp(disp_operator * L/2) .* exp(-alpha_m*L/2); E ifft(ifftshift(E_hat)); % 每一步加入泵浦项对应LLE的外部注入 E E sqrt(P_pump/P_round) * exp(-1i*delta*...); end % 输出光谱取输出场做FFT S abs(fftshift(fft(E))).^2; plot((omega-omega_pump)/(2*pi), 10*log10(S/max(S)));这里必须说明上面的代码是“半伪代码”因为完整的泵浦注入项、失谐项和边界条件的处理需要对照你自己的归一化方式完成。但从结构上它已经覆盖了SSFM的核心流程非线性一步、线性一步、泵浦注入往复迭代。与NLSE求解的一个区别普通的光纤脉冲传输是单向推进而微环频率梳是“绕圈循环”每一步都必须把场绕回起点接着走所以时间窗口要设计成恰好覆盖整数个FSR否则会因为边界不匹配产生人为调制。这也是很多人仿真结果里出现不明纹波的原因之一——窗口没对齐。3.4 怎么判断仿真“成功”了跑完仿真你得到的是一个二维矩阵快时间变量×往返次数或者只取稳态时域输出。判断是否形成频率梳最简单的方法是看光谱出现了多根等间距的峰值间距恰好等于FSR或FSR的整数倍这是频率梳的标志性特征。峰值底部有一定的连续光谱或旁瓣说明存在孤子或混沌态具体形态取决于参数。如果只有泵浦附近一两根峰说明能量没有有效转移到其他模式非线性作用不足需要提高功率或降低损耗。如果你还能进一步观察到稳定的高信噪比等间距梳齿并且每根梳齿的线宽远小于泵浦线宽那恭喜这个仿真已经接近真实微环频率梳的行为。4. 参数扫描与物理洞察让仿真为你解释背后的规律仿真最大的价值不只是“跑出一个梳状谱”而是通过扫描关键参数理解频率梳的生成边界和演化特征。这部分我建议你养成做参数扫描的习惯固定其他参数单独扫描泵浦功率、失谐量和色散值观察阈值、梳齿包络和混沌区的变化。4.1 泵浦功率的阈值行为微环频率梳的产生有明确阈值。低于阈值时腔内非线性相移不足以补偿色散失配四波混频被抑制输出几乎只有泵浦峰。超过阈值后调制不稳定性MI开始起作用泵浦附近率先出现对称的边带这就是频率梳的“第一对梳齿”。我在仿真时观察到泵浦功率从0.2 W升到0.6 W的过程中光谱是“阶跃”式变化的一开始只有边带然后边带附近出现更多并峰最后横跨整个费米窗口形成整齐梳齿。这里有个很容易被忽略的点阈值功率不是固定的它和腔内Q值强烈相关。Q值每提高一倍有效腔内功率可能提升一个量级阈值功率会显著下降。这也是为什么实际实验中大家都在拼高Q值和低损耗波导。4.2 色散符号异常色散是梳子的“润滑剂”LLE里的色散项是i·(∂²E/∂τ²)符号决定了频率梳的形态。多数稳定梳状谱生成都发生在异常色散区域β₂ 0因为此时色散和克尔非线性可以平衡支持孤子和类梳结构。如果你把色散设成正常色散β₂ 0你大概率得到的是弱的、不稳定的、包络呈“暗孤子”形式的梳状结构甚至直接得不到梳齿。在Matlab里发现问题就这么简单把beta2符号改一下跑完光谱肉眼立刻能看出现差异。这也是验证你的仿真代码是否真正捕捉到非线性效应的一个“仪器级别”指标——如果正常色散和异常色散的结果几乎一样说明代码里的色散项根本没生效检查频域算子的符号和单位吧。色散类型预期光谱特征物理原因异常色散β₂0明显的等间距梳齿能量带宽宽支持亮孤子和调制不稳定性增益正常色散β₂0弱而碎的边带梳齿不明显缺少调制不稳定性增益相位平衡困难4.3 失谐量与孤子形成泵浦失谐θ是另一个决定性参数。这里可以观察到一个非常有趣的“孤子步”现象随着失谐从零增大腔内场从小噪声逐步演化当失谐超过某个临界值系统会跳变到孤子态光谱呈现双曲正割包络梳齿整齐排列。但失谐过大也不行否则泵浦项太弱腔内场衰减频率梳像蜡烛一样被吹灭。在Matlab里扫θ把每次迭代的腔内能量画出来你会看到明显的台阶结构——这一小段阶跃正是孤子形成的证据。这部分内容我在刚开始仿真时完全不理解觉得“不就是把泵浦波长调偏一点吗”后来扫了上百个失谐点才真正明白失谐不仅仅影响谐振耦合效率也是频率梳稳定性的“开关”。扫参时建议用热图展示横轴是失谐纵轴是往返次数颜色是腔内能量能非常直观地看到阈值和孤子台阶的出现。5. 常见问题、排查思路与实操心得仿真做到后面遇到的大量问题都不是物理问题而是数值和工程问题。这里把我踩过的坑按“症状—原因—解法”列成速查表方便大家对照。5.1 症状一仿真结果发散场强变成NaN或Inf可能原因解决方案时间窗口太小边界反射干扰增加窗口长度或改用吸收边界网格点数不够色散算子误差大增大N如4096→8192单次步长过大减小分步长度或增加往返迭代次数初始场包含高频分量用平滑高斯种子或小幅噪声非线性系数γ设置过大检查单位换算可能是量纲错误新手最容易犯的错误是直接用物理参数比如10⁻²⁴量级的β₂去计算算子结果数值太小在双精度浮点下直接损失精度。正确的做法是先在无量纲化框架下设定系数再换算回物理参数显示。5.2 症状二光谱里只有泵浦峰没有梳齿这是最常见的情况尤其是第一次仿真的新手。先不要怀疑代码按顺序检查物理参数泵浦功率是否超过阈值如果只有泵浦峰先用刚才提的扫功率方法看看是否有边带从噪声中冒出来。Q值是否足够高腔内功率增强不够时非线性效应根本起不来。失谐量是否过大失谐过大时泵浦无法有效耦合进环内。色散是否恰好是合适符号如果你设置的是正常色散大概率等不到梳齿。初始场是否过于“干净”四波混频需要种子“噪声”来打破对称性实际实验里噪声来源于量子涨落仿真里如果初始场毫无噪声可能因对称性而无法起振。特别想强调第5点物理上频率梳的产生需要从微弱噪声中自发建立如果你用一个完全平滑的连续波作为初始条件系统可能因为对称性守恒而维持在无梳状态。解决方法是加一个极低幅度的高斯随机噪声幅度通常是泵浦幅度的10⁻⁶到10⁻⁴倍。这个细节在文献里很少有人提但在仿真里绝对关键。5.3 症状三梳齿间隔不对或出现额外包络调制出现这种问题多半是“快时间窗口”设置有问题。频率梳仿真要求快时间窗口恰好对应整数个FSR周期否则在窗口边缘会出现不连续导致频谱泄露和额外峰谷。解决办法是精确计算T_window m / FSR其中m为正整数。同时在Matlab里用fftshift之前检查频率轴是否按顺序定义不然画出来的光谱左右颠倒、梳齿间隔错乱。5.4 我踩过的“数值坑”和提速技巧仿真调参过程中我自己摸索出了几个提升效率和稳定性的技巧保存中间态每次扫描参数时把腔内场演化的中间态存成.mat文件。这样出问题以后不用从头跑直接加载某个中间态继续迭代能节省大量时间。我通常在每50个往返周期存一帧配合MAT格式的压缩存储文件不大但非常有价值。GPU加速分步傅里叶法里占计算量最大的是FFT。如果你的机器有NVIDIA显卡可以考虑用gpuArray把场E放到GPU上FFT计算可以快5到10倍。注意GPU上的ifft/fft结果要记得传回CPU再画图否则内存管理混乱。并行参数扫描用parfor循环扫描泵浦功率或失谐参数每个参数点启动一个worker并行计算。这需要Parallel Computing Toolbox但在多核机器上的提速非常明显尤其在扫上百个点时。不要过早优化先跑通单点仿真再考虑加速。很多新手一上来就写并行循环结果参数还没摸清代码已经在并发崩溃。先确保物理模型正确再优化性能。5.5 个人经验从“跑通代码”到“理解物理”做这个项目最深的体会是仿真代码跑通只是第一步真正的门槛是把光学语义和Matlab数值操作一一对应。你没有必要背下LLE的推导但你必须知道每个矩阵乘法在物理上意味着什么。比如非线性步的 exp(iγ|E|²L) 是在模拟同一点上的自相位调制频域的色散算子是在模拟不同频率分量以不同速度传播——一旦脑子里有了这幅“物理图像”排查问题就快得多。另一个实用建议是刚开始不要把目标定成“复现一篇Nature论文里的完美频率梳”而是先复现最基本的色散效应和四波混频现象。比如先跑一个不含非线性项的版本确认色散导致的脉冲展宽与解析解一致再加非线性项确认会出现自陡峭和频谱展宽最后加入微环边界条件才可能看到真正的频率梳。这种“分阶段验证法”帮你把每一步的误差隔离在一个可控范围避免最后锅底下的问题根本不知道出在哪根柱子。另外Matlab本身的脚本化特性非常适合这种“假设—仿真—验证”的循环。我习惯在主脚本里加一堆assert比如检查输入功率是否为正数、时间窗口是否覆盖整数个FSR、色散算子的最大值是否在安全范围内稍微花点时间但能避免深夜debug的绝望感。真实项目里跑一遍1000个腔往返大概只要几分钟但如果你在参数初始化时就埋了坑后面所有结果都是垃圾。数据与图形分析脚本分离也很重要我用的plot脚本里只管读取.mat文件生成频谱演化热图和稳态光谱图减少重复劳动。最后再分享一个小技巧想确认频率梳是否真的“梳”起来不要只看稳态光谱要把每个往返周期的光谱叠在一张图上看看稳定性。真实的微环频率梳在孤子态下是周期稳定的如果有人让你看“动态演化图”你会看到等间距梳齿像琴弦一样从头到尾稳定在那里。如果梳齿在整个迭代过程中位置漂移或者幅度抖动说明系统还在瞬态或者锁模不稳定输出光谱看着漂亮但落到实验里是没法工作的。这个判断方法在论文里很难找到是我自己反复跑了上百组参数后总结出来的。 
长河君 · 2026/10/3 15:17:21 · 阅读 1.2k 微环谐振腔光学频率梳Matlab仿真：从LLE方程到分步傅里叶法全解析
光学频率梳这个词这几年在光子学里几乎是绕不开的。无论是光通信里的多波长光源、微波光子学里的低噪声振荡器还是精密光谱测量大家都会提到它。而微环谐振腔作为片上实现光学频率梳的核心平台之一更是圈内人默认的高起点。我最初上手这个方向时用的就是Matlab做全链路仿真从单频泵浦光在微环里一点点演化成几十根等间隔梳齿整个过程非常“上瘾”但踩坑也不少。这篇内容就把我自己的实现路径、理论基础、关键参数和排查方法完整捋一遍写给正在或者准备做微环频率梳仿真的朋友参考。1. 为什么是微环谐振腔先弄懂它在光学频率梳里的位置微环谐振腔本质上是集成光子学里最常见的结构之一一根直波导加上一个闭合的环形波导光通过倏逝波耦合从直波导进入环内在环里绕圈共振。它的名字听起来很学术但结构上可以类比成一个“光学蓄水池”——光在环内不断循环积累能量只要满足谐振条件腔内功率会远远高于输入功率。这种功率增强效应是后面产生非线性效应的根基。1.1 微环的基本结构与核心参数标准的微环结构有两种一种是全通型all-pass一根直波导加一个环只有一个耦合区另一种是上下载型add-drop两根直波导加一个环两个耦合区。频率梳仿真里常用的通常是全通型结构因为它结构简单、损耗路径清晰便于分析腔内功率积累过程。和微环相关的几个核心参数决定了整个仿真的走向自由光谱范围FSR环内相邻谐振峰之间的频率间隔由环周长和群折射率决定公式是 FSR c / (ng · L)单位是Hz。半径越大FSR越小。这个参数直接决定了梳齿之间的间距。品质因子Q值衡量谐振腔储能能力的指标Q值越高腔内光子寿命越长功率增强越明显。Q值受耦合损耗、波导传输损耗等因素约束。功率增强因子腔内功率和输入功率的比值近似正比于Q值。这是微环能作为非线性平台的核心原因。耦合系数κ直波导和环之间能量交换的强度决定了腔的加载Q值和外耦合率。仿真中很多人上来就写方程结果参数一团乱最后出了结果也不知道对不对。我的建议是先固定物理结构参数——选一个常见氮化硅Si₃N₄或者硅Si波导平台把有效折射率、群折射率、非线性折射率这些材料参数确定下来再推导耦合方程。否则Matlab脚本里每个变量都悬空后面调参就是灾难。1.2 为什么不能只用“线性腔模”解释频率梳线性状态下微环的频率响应就是一系列等间隔的洛伦兹峰泵浦光耦合进某个谐振峰后输出仍是单频光。频率梳的产生需要“新频率”出现而这必须依赖非线性过程。微环频率梳的核心物理机制是四波混频FWM这是由材料的三阶非线性极化率χ⁽³⁾引起的。两个泵浦光子湮灭、产生一红一蓝两个新光子信号光与闲频光只要满足能量守恒和相位匹配约束能量就会不断从泵浦光转移到新的频率分量上。随着腔内光强增大新产生的频率又作为新的“泵浦”继续产生更多频率成分最终形成跨越数个谐振峰的光学频率梳。这也是为什么仿真频率梳一定要用非线性模型线性传输方程只能给你梳状滤波器的透射谱永远得不到梳状光源的光谱。你在Matlab里既要处理色散、损耗、耦合这些线性效应也要处理克尔非线性带来的相位调制和频率搬移这两类效应在腔内往返过程中是同时发生的。2. 从物理机制到Matlab模型核心方程与仿真原理频率梳的时域仿真圈内最常用的模型是Lugiato-Lefever方程LLE。它本质上是一个带边界条件的非线性薛定谔方程专门用来描述被动微环谐振腔内连续波CW泵浦下的非线性动力学。可以说掌握了LLE就掌握了微环频率梳仿真的“主心骨”。2.1 Lugiato-Lefever方程的含义LLE的归一化形式大概长这个样子dE(t,τ)/dt -(1 i·θ)E i·|E|²E i·(∂²E/∂τ²) F这行式子看起来复杂但拆开看每一项物理意义非常清晰左边 dE/dt 是腔内慢变幅度随时间慢时间域的演化t代表慢时间也就是在腔内的往返次数累计。右边第一项 -(1i·θ)E 中的实部1是腔内损耗项虚部θ是泵浦光与谐振峰的失谐量。第二项 i·|E|²E 是克尔非线性项对应四波混频和自相位调制效应这是产生新频率的根源。第三项 i·(∂²E/∂τ²) 代表群速度色散GVDτ是快时间坐标对应一个腔周期内的波形细节。这个色散项决定了梳齿包络的形状和孤子形成的可能性。最后一项 F 是连续波泵浦项。在实际仿真中你不需要从第一性原理推导这个方程但要能写出每一项的物理量纲理解它在“腔往返”这个物理图像里的位置。很多人拿着别人代码跑通了但不知道每个系数对应实验中哪个参数一遇到异常结果就抓瞎。2.2 为什么用分步傅里叶法求解LLE是非线性偏微分方程线性项损耗、色散、失谐和非线性项克尔效应在物理上是耦合的。直接做时域有限差分可以求解但数值稳定性和效率都差一些。工程里最流行的方案是分步傅里叶法Split-Step Fourier MethodSSFM。核心思想很朴素在一个极小的时间步长内先假设非线性效应单独作用再用傅里叶变换把线性效应放到频域里单独处理最后再把它们拼起来。因为色散项在频域里只是一个乘法运算效率远高于时域差分。用生活化类比来说这就像你做一道菜——先腌制非线性步再加热线性步但加热过程中腌制的效果也在同步进行只要每一步时间足够短顺序误差就可以忽略。在Matlab里实现SSFM主要的循环结构是初始化腔内场E通常从一个小噪声或弱泵浦场开始。非线性步E E · exp(i·γ·|E|²·dz)直接在时域完成。FFT变换到频域。线性步E E · exp(线性算子·dz)包含损耗和色散。IFFT回时域检查收敛或粒子数守恒。这里要注意Matlab的fft和ifft在处理光场包络时频率轴的顺序是f(1)对应直流然后把负频率放在后半段。很多人第一次画光谱时发现左右不对称就是没做fftshift。这是入门阶段最典型的低级错误但不仔细看真的很容易忽略。2.3 归一化与无量纲化的作用LLE里最关键的设计决策是归一化。物理参数波导损耗α、非线性系数γ、色散系数β₂、环周长L数量级差异极大非线性系数可能是W⁻¹·m⁻¹量级而色散可能是ps²/m量级。直接代入计算会让Matlab里的数值巨大或极小极易导致误差和溢出。归一化的思路是把时间、长度、场强都缩放到无量纲单位。常见的做法是引入一个参考时间尺度T₀通常取快时间窗口内的某个特征值把快时间τ归一化为τ/T₀把腔传播距离归一化为环周长L场强E则用非线性系数和损耗的比值来归一化。这样方程里只剩下几个无量纲系数失谐θ、归一化色散、归一化泵浦幅度F。这一步做好的好处是同一套代码可以扫描任意物理参数组合只需要改归一化系数不需要改求解核心。我自己的习惯是先在脚本开头写一个“归一化计算”段落把所有物理参数转换成无量纲参数再用这些参数去跑SSFM。调试时如果输出异常优先检查归一化是否正确而不是怀疑求解算法本身。3. Matlab代码实现从空白脚本到能出梳状谱仿真环境的搭建虽然不难但很多细节会影响效率甚至影响正确性。我在做这个项目时用的是较新版本的MatlabWindows平台上跑中途也经历过换机器、重装环境、工具箱缺失的问题。这里把环境准备、参数设置和核心代码骨架分开说。3.1 环境准备版本、工具箱和安装要点Matlab本身的分发包很大早期版本装起来比较折腾。现在的安装流程已经友好很多但有几个点还是值得提醒版本选择如果你手上是新电脑可以直接用较新的Matlab版本比如R2025b或R2026b新版本对GPU加速、并行工具箱的支持更好。老版本跑SSFM也可以只是大规模参数扫描时会慢得让人崩溃。工具箱要求核心仿真只需要基础Matlab和Signal Processing Toolbox用于滤波和信号分析。如果要画漂亮的光谱图和动态演化图还需要DSP和基本绘图功能——这些在标准安装包里都有。Optimization Toolbox可用于自动拟合参数但不是必须。许可证和激活安装时经常遇到许可证问题比如License Manager报错或不识别License文件通常需要检查Host ID是否匹配、许可证路径是否设置正确。网上很多教程提到下载和安装步骤其实关键就两条安装路径不要带中文激活时需要管理员权限。遇到license.lic的Host ID对不上多半是网卡顺序变化或虚拟机环境导致的重新生成License文件就行。并行计算如果CPU核心多建议启用Parallel Computing Toolbox的parfor做参数扫描。频率梳仿真本质上是一个大参数循环扫描泵浦功率、失谐、色散等参数每个参数点都要跑几千个腔往返串行计算很浪费时间。在Matlab里跑的脚本结构我习惯分成四个文件参数定义脚本params.m、归一化计算脚本normalize.m、主仿真循环main_sim.m、绘图与数据分析脚本plot_results.m。这样改参数时不用反复翻代码出图也更快。千万不要把所有代码塞在一个脚本里后期调参时你会后悔的。3.2 参数设定从一个“能出梳”的典型值开始对于仿真新手最大的坑是参数设置不合理导致永远无法产生频率梳。我的建议是先找到文献里已经被验证过的参数范围然后在附近微调。以下是一组典型的氮化硅微环参数基于常见文献和商用SOI/SiN平台指标参数数值说明微环半径 R40 μm决定FSR半径越小FSR越大环波导截面0.5 μm × 0.8 μm常见脊波导尺寸有效折射率 n_eff1.98与波长相关需用模式求解器算群折射率 n_g2.1决定FSR仿真中常近似非线性折射率 n₂2.4×10⁻¹⁹ m²/W氮化硅材料典型值波导损耗 α0.5 dB/cm代表损耗水平功率增强倍数10~50由Q值决定估算用泵浦功率 P_in100 mW ~ 1 W低于阈值则无法起振泵浦波长1550 nmC-band常见泵浦实际上在仿真里你并不会直接输入n₂和α而是换算成γ非线性系数和腔损耗系数。如果嫌换算麻烦可以用近似公式γ ≈ 2π·n₂/(λ·A_eff)A_eff为有效模面积损耗系数则需要把dB/cm换算成m⁻¹。这里用到了Matlab的数据处理能力写个小函数批量换算体力活很值得一提频率梳仿真大多数时间不是耗在求解上而是耗在参数单位换算上。3.3 核心代码骨架分步傅里叶法的Matlab实现下面是经过简化但仍可直接复现的Matlab代码框架我用它跑出了初步的CI梳状光谱。代码里保留了关键步骤注释方便你对照上面的LLE解释。% 主仿真脚本流程简化版 % 参数初始化 lambda0 1550e-9; % 泵浦波长 c 3e8; % 光速 R 40e-6; % 环半径 L 2*pi*R; % 环周长 n_g 2.1; % 群折射率 FSR c/(n_g*L); % 自由光谱范围 (Hz) % 非线性/损耗参数换算后 alpha 0.5; % dB/cm alpha_m alpha/8.686/0.01; % 换算为 1/m gamma 2.0; % W^-1 m^-1需结合材料与模场面积估算 P_pump 0.5; % 泵浦功率 W % 快时间窗口与网格 T_window 20/FSR; % 时间窗口覆盖20个FSR N 4096; % 网格点数 tau linspace(-T_window/2, T_window/2, N); d_tau tau(2)-tau(1); omega 2*pi*fftshift(-N/2:N/2-1)/(N*d_tau); % 频域角频率轴 % 色散参数以二阶色散为主单位换算后代入 beta2 -20e-24; % s^2/m择优异常色散负号 disp_operator -1i*0.5*beta2*(omega.^2); % 频域色散算子 % 内腔场初始化小噪声 E sqrt(2e-3)*exp(-(tau/T0).^2); % 弱高斯种子 噪声 % 实际中常用白噪声初始化更好避免包络偏差 % 腔往返迭代 N_roundtrip 2000; % 往返次数 for k 1:N_roundtrip % 非线性步时域相位调制 E E .* exp(1i*gamma*abs(E).^2*L/2); % 线性步频域色散和损耗 E_hat fftshift(fft(E)); E_hat E_hat .* exp(disp_operator * L/2) .* exp(-alpha_m*L/2); E ifft(ifftshift(E_hat)); % 每一步加入泵浦项对应LLE的外部注入 E E sqrt(P_pump/P_round) * exp(-1i*delta*...); end % 输出光谱取输出场做FFT S abs(fftshift(fft(E))).^2; plot((omega-omega_pump)/(2*pi), 10*log10(S/max(S)));这里必须说明上面的代码是“半伪代码”因为完整的泵浦注入项、失谐项和边界条件的处理需要对照你自己的归一化方式完成。但从结构上它已经覆盖了SSFM的核心流程非线性一步、线性一步、泵浦注入往复迭代。与NLSE求解的一个区别普通的光纤脉冲传输是单向推进而微环频率梳是“绕圈循环”每一步都必须把场绕回起点接着走所以时间窗口要设计成恰好覆盖整数个FSR否则会因为边界不匹配产生人为调制。这也是很多人仿真结果里出现不明纹波的原因之一——窗口没对齐。3.4 怎么判断仿真“成功”了跑完仿真你得到的是一个二维矩阵快时间变量×往返次数或者只取稳态时域输出。判断是否形成频率梳最简单的方法是看光谱出现了多根等间距的峰值间距恰好等于FSR或FSR的整数倍这是频率梳的标志性特征。峰值底部有一定的连续光谱或旁瓣说明存在孤子或混沌态具体形态取决于参数。如果只有泵浦附近一两根峰说明能量没有有效转移到其他模式非线性作用不足需要提高功率或降低损耗。如果你还能进一步观察到稳定的高信噪比等间距梳齿并且每根梳齿的线宽远小于泵浦线宽那恭喜这个仿真已经接近真实微环频率梳的行为。4. 参数扫描与物理洞察让仿真为你解释背后的规律仿真最大的价值不只是“跑出一个梳状谱”而是通过扫描关键参数理解频率梳的生成边界和演化特征。这部分我建议你养成做参数扫描的习惯固定其他参数单独扫描泵浦功率、失谐量和色散值观察阈值、梳齿包络和混沌区的变化。4.1 泵浦功率的阈值行为微环频率梳的产生有明确阈值。低于阈值时腔内非线性相移不足以补偿色散失配四波混频被抑制输出几乎只有泵浦峰。超过阈值后调制不稳定性MI开始起作用泵浦附近率先出现对称的边带这就是频率梳的“第一对梳齿”。我在仿真时观察到泵浦功率从0.2 W升到0.6 W的过程中光谱是“阶跃”式变化的一开始只有边带然后边带附近出现更多并峰最后横跨整个费米窗口形成整齐梳齿。这里有个很容易被忽略的点阈值功率不是固定的它和腔内Q值强烈相关。Q值每提高一倍有效腔内功率可能提升一个量级阈值功率会显著下降。这也是为什么实际实验中大家都在拼高Q值和低损耗波导。4.2 色散符号异常色散是梳子的“润滑剂”LLE里的色散项是i·(∂²E/∂τ²)符号决定了频率梳的形态。多数稳定梳状谱生成都发生在异常色散区域β₂ 0因为此时色散和克尔非线性可以平衡支持孤子和类梳结构。如果你把色散设成正常色散β₂ 0你大概率得到的是弱的、不稳定的、包络呈“暗孤子”形式的梳状结构甚至直接得不到梳齿。在Matlab里发现问题就这么简单把beta2符号改一下跑完光谱肉眼立刻能看出现差异。这也是验证你的仿真代码是否真正捕捉到非线性效应的一个“仪器级别”指标——如果正常色散和异常色散的结果几乎一样说明代码里的色散项根本没生效检查频域算子的符号和单位吧。色散类型预期光谱特征物理原因异常色散β₂0明显的等间距梳齿能量带宽宽支持亮孤子和调制不稳定性增益正常色散β₂0弱而碎的边带梳齿不明显缺少调制不稳定性增益相位平衡困难4.3 失谐量与孤子形成泵浦失谐θ是另一个决定性参数。这里可以观察到一个非常有趣的“孤子步”现象随着失谐从零增大腔内场从小噪声逐步演化当失谐超过某个临界值系统会跳变到孤子态光谱呈现双曲正割包络梳齿整齐排列。但失谐过大也不行否则泵浦项太弱腔内场衰减频率梳像蜡烛一样被吹灭。在Matlab里扫θ把每次迭代的腔内能量画出来你会看到明显的台阶结构——这一小段阶跃正是孤子形成的证据。这部分内容我在刚开始仿真时完全不理解觉得“不就是把泵浦波长调偏一点吗”后来扫了上百个失谐点才真正明白失谐不仅仅影响谐振耦合效率也是频率梳稳定性的“开关”。扫参时建议用热图展示横轴是失谐纵轴是往返次数颜色是腔内能量能非常直观地看到阈值和孤子台阶的出现。5. 常见问题、排查思路与实操心得仿真做到后面遇到的大量问题都不是物理问题而是数值和工程问题。这里把我踩过的坑按“症状—原因—解法”列成速查表方便大家对照。5.1 症状一仿真结果发散场强变成NaN或Inf可能原因解决方案时间窗口太小边界反射干扰增加窗口长度或改用吸收边界网格点数不够色散算子误差大增大N如4096→8192单次步长过大减小分步长度或增加往返迭代次数初始场包含高频分量用平滑高斯种子或小幅噪声非线性系数γ设置过大检查单位换算可能是量纲错误新手最容易犯的错误是直接用物理参数比如10⁻²⁴量级的β₂去计算算子结果数值太小在双精度浮点下直接损失精度。正确的做法是先在无量纲化框架下设定系数再换算回物理参数显示。5.2 症状二光谱里只有泵浦峰没有梳齿这是最常见的情况尤其是第一次仿真的新手。先不要怀疑代码按顺序检查物理参数泵浦功率是否超过阈值如果只有泵浦峰先用刚才提的扫功率方法看看是否有边带从噪声中冒出来。Q值是否足够高腔内功率增强不够时非线性效应根本起不来。失谐量是否过大失谐过大时泵浦无法有效耦合进环内。色散是否恰好是合适符号如果你设置的是正常色散大概率等不到梳齿。初始场是否过于“干净”四波混频需要种子“噪声”来打破对称性实际实验里噪声来源于量子涨落仿真里如果初始场毫无噪声可能因对称性而无法起振。特别想强调第5点物理上频率梳的产生需要从微弱噪声中自发建立如果你用一个完全平滑的连续波作为初始条件系统可能因为对称性守恒而维持在无梳状态。解决方法是加一个极低幅度的高斯随机噪声幅度通常是泵浦幅度的10⁻⁶到10⁻⁴倍。这个细节在文献里很少有人提但在仿真里绝对关键。5.3 症状三梳齿间隔不对或出现额外包络调制出现这种问题多半是“快时间窗口”设置有问题。频率梳仿真要求快时间窗口恰好对应整数个FSR周期否则在窗口边缘会出现不连续导致频谱泄露和额外峰谷。解决办法是精确计算T_window m / FSR其中m为正整数。同时在Matlab里用fftshift之前检查频率轴是否按顺序定义不然画出来的光谱左右颠倒、梳齿间隔错乱。5.4 我踩过的“数值坑”和提速技巧仿真调参过程中我自己摸索出了几个提升效率和稳定性的技巧保存中间态每次扫描参数时把腔内场演化的中间态存成.mat文件。这样出问题以后不用从头跑直接加载某个中间态继续迭代能节省大量时间。我通常在每50个往返周期存一帧配合MAT格式的压缩存储文件不大但非常有价值。GPU加速分步傅里叶法里占计算量最大的是FFT。如果你的机器有NVIDIA显卡可以考虑用gpuArray把场E放到GPU上FFT计算可以快5到10倍。注意GPU上的ifft/fft结果要记得传回CPU再画图否则内存管理混乱。并行参数扫描用parfor循环扫描泵浦功率或失谐参数每个参数点启动一个worker并行计算。这需要Parallel Computing Toolbox但在多核机器上的提速非常明显尤其在扫上百个点时。不要过早优化先跑通单点仿真再考虑加速。很多新手一上来就写并行循环结果参数还没摸清代码已经在并发崩溃。先确保物理模型正确再优化性能。5.5 个人经验从“跑通代码”到“理解物理”做这个项目最深的体会是仿真代码跑通只是第一步真正的门槛是把光学语义和Matlab数值操作一一对应。你没有必要背下LLE的推导但你必须知道每个矩阵乘法在物理上意味着什么。比如非线性步的 exp(iγ|E|²L) 是在模拟同一点上的自相位调制频域的色散算子是在模拟不同频率分量以不同速度传播——一旦脑子里有了这幅“物理图像”排查问题就快得多。另一个实用建议是刚开始不要把目标定成“复现一篇Nature论文里的完美频率梳”而是先复现最基本的色散效应和四波混频现象。比如先跑一个不含非线性项的版本确认色散导致的脉冲展宽与解析解一致再加非线性项确认会出现自陡峭和频谱展宽最后加入微环边界条件才可能看到真正的频率梳。这种“分阶段验证法”帮你把每一步的误差隔离在一个可控范围避免最后锅底下的问题根本不知道出在哪根柱子。另外Matlab本身的脚本化特性非常适合这种“假设—仿真—验证”的循环。我习惯在主脚本里加一堆assert比如检查输入功率是否为正数、时间窗口是否覆盖整数个FSR、色散算子的最大值是否在安全范围内稍微花点时间但能避免深夜debug的绝望感。真实项目里跑一遍1000个腔往返大概只要几分钟但如果你在参数初始化时就埋了坑后面所有结果都是垃圾。数据与图形分析脚本分离也很重要我用的plot脚本里只管读取.mat文件生成频谱演化热图和稳态光谱图减少重复劳动。最后再分享一个小技巧想确认频率梳是否真的“梳”起来不要只看稳态光谱要把每个往返周期的光谱叠在一张图上看看稳定性。真实的微环频率梳在孤子态下是周期稳定的如果有人让你看“动态演化图”你会看到等间距梳齿像琴弦一样从头到尾稳定在那里。如果梳齿在整个迭代过程中位置漂移或者幅度抖动说明系统还在瞬态或者锁模不稳定输出光谱看着漂亮但落到实验里是没法工作的。这个判断方法在论文里很难找到是我自己反复跑了上百组参数后总结出来的。