简介一份用于解析硅微环调制器中非线性效应的时域非线性模型 Matlab 代码面向光通信、电子信息、计算机及数学专业学生和研究人员。代码采用参数化编程参数可灵活调整注释明细完整给出 odesystem、ring_params、run_sims 等核心脚本支持 Matlab 2014a 至 2024a 多版本运行。资源包共 11 个文件以 8 个 .m 源文件为主另含 README.md 说明文档、数据文件及运行辅助文件压缩包整体仅 16KB轻量易部署。配套案例数据可直接运行便于使用者模拟强度依赖折射率变化、热效应、载流子等离子体色散等非线性过程观察不同输入光强、调制频率及环参数下的调制器响应。代码结构清晰、注释丰富适用于课程设计、期末大作业及毕业设计中的器件建模与性能分析已有 73 人参与学习可快速上手开展硅光调制器非线性效应研究。 硅光子入门的人多半早晚会碰到硅微环调制器Silicon Microring Modulator这个名字。数据中心光互联做WDM的场合它几乎是绕不开的选型尺寸小、功耗低、工艺和CMOS兼容动辄几十吉波特率的调制也不在话下。不过一旦你真正开始仿真它的动态行为就会发现一个麻烦——硅这个材料在高光强下“不太老实”热光效应、双光子吸收、自由载流子色散这些非线性效应会一起冒出来把谱线、损耗、谐振点搅得一团糟。我最近在处理一份用MATLAB实现硅微环调制器时域非线性模型的代码包就是专门用来描述这种“非线性效应”的动态行为。把代码跑通、又对照文献把各模块过了一遍之后我觉得里面的建模思路和调试经验很值得写出来。这篇博文就顺着“物理机理—数学模型—MATLAB实现—调试避坑”这条线把整个代码包的逻辑拆开讲透。不管你是刚接手微环仿真任务的研究生还是需要在系统级仿真里加一个非线性调制器模型的光电工程师照着这篇文章应该都能明白怎么让这套代码跑起来、又能怎么改。1. 为什么要给硅微环调制器建时域非线性模型1.1 硅微环里那些“看不见”的非线性效应很多刚接触微环的人第一反应是微环不就是个谐振器吗我算一下谐振波长、Q值、透过率不就行了这话在小信号、低光功率下没错但硅微环调制器实际工作时腔内会建立起很高的功率密度几种非线性效应就会轮番登场。第一种是热光效应。硅的折射率随温度变化很明显热光系数大约在1.8e-4 /K左右。光在波导里传播总有吸收哪怕是本征吸收很小加上耦合损耗和侧壁散射损耗也会让一小部分光功率变成热量。热量在180 ns量级的热时间常数内累积起来把波导局部温度抬高几十甚至上百开尔文谐振波长就往长波方向漂移。这就像水槽的壁被晒热以后变形了水面高度自然就不准了。第二种是双光子吸收TPA加自由载流子色散FCD。当腔内光强足够高硅材料会同时吸收两个光子产生电子空穴对。自由载流子一方面通过等离子体色散效应改变折射率——在近红外波段表现为折射率下降也就是FCD另一方面自身会带来额外的自由载流子吸收FCA增大波导损耗。载流子寿命一般在几纳秒到几十纳秒量级比热效应快得多但又比光场周期慢得多所以它带来的谐振漂移和损耗扰动会跟着调制信号来回变化。第三种是克尔效应。硅的三阶非线性系数不算小在极高峰值功率下克尔效应会瞬间改变折射率。它在皮秒量级响应跟光场几乎是同步的主要影响高速脉冲的峰功率位置。不过在大多数调制器工作条件下克尔效应相对热光和FCD来说贡献较小很多时域模型会先忽略它把它留作高阶修正项。这些效应的时间尺度差异非常关键光场周期是飞秒到皮秒量级载流子动态是纳秒量级热效应是百纳秒到微秒量级。要用一个统一的模型把三种快慢不同的物理过程都装进去就只有在时域里做耦合求解这正好回答了“为什么要写时域模型”这个核心问题。1.2 稳态仿真搞不定的场景只能交给时域模型如果你只关心偏置点固定、功率很低、调制速率很慢的情况用频域稳态模型算一次透过率曲线就够了。但硅微环调制器的真实工作场景恰恰是稳态模型最尴尬的几个地方。首先是高速调制动态。调制信号驱动PN结改变载流子浓度谐振点跟着信号来回摆动。这时候热效应也在慢悠悠地漂移谐振点两者叠加的结果是同一个驱动电压在不同时刻看到的透过率是完全不同的。稳态模型给不出这种动态响应只有把光场方程、热方程、载流子方程一起在时域里推才能看到输出光波形怎么被非线性“拧”变形。其次是高输入光功率。现代光互联为了压低系统噪声常常把激光器功率往上提几毫瓦甚至几十毫瓦的片上光功率很常见。腔内功率密度一上来双光子吸收和FCD立刻变得不可忽略调制器插损、消光比、啁啾全部偏离设计值。你拿稳态线性模型去算系统预算算出来跟实测对不上根因就在这。还有一个场景是突发模式或开关瞬态。热效应的时间常数很长意味着当通道从无光切换到大功率光谐振点要经过几十微秒才能稳定下来。这段“热建立时间”内调制器工作点是漂移的对判断系统能不能快速恢复、能不能做突发接收至关重要。时域非线性模型能够自然捕捉这个过程而任何频域打包模型都做不到。2. 时域非线性模型的数学框架从物理到方程2.1 耦合模方程用复振幅还是腔内能量时域模型最常见的出发点是耦合模理论CMT。核心思想是把微环看成一个谐振腔用腔内复振幅 a(t) 来描述光场输入光场记为 s_in(t)。复振幅的模平方 |a|² 对应腔内能量它跟输入功率之间通过耦合系数 κ 联系起来。用复振幅比用纯能量更有优势因为相位信息被保留了下来。调制器在PN结驱动下不仅会有幅度变化还会伴随相位调制啁啾如果你想仿真啁啾对光纤传输链路的影响就必须用复振幅而不是实时功率。方程写出来大体是这个形式[ \frac{da(t)}{dt} \left[ j 2\pi (\Delta f_{0} \Delta f_{th} \Delta f_{FCD}) - \frac{\gamma_{l} \gamma_{TPA} \gamma_{FCA}}{2} \right] a(t) \kappa s_{in}(t) ]这里的 j 是虚数单位Δf0 是静态失谐频率Δf_th 是热效应引起的谐振频率偏移Δf_FCD 是载流子色散引起的偏移。γ_l 对应线性损耗速率γ_TPA 和 γ_FCA 是双光子吸收和自由载流子吸收带来的额外损耗速率。整体上看这个方程就是一个带复频移的阻尼谐振子方程输入激励是 κ s_in(t)。在MATLAB里实现时a(t) 要拆成实部和虚部两个实数状态量方便ode系列求解器处理。不要图省事把复数状态直接硬编码进去后面接射频驱动信号、算输出光功率时会很被动。2.2 温度、载流子怎么“塞”进谐振方程光场的方程看似简单真正的难点在于把温度T和载流子浓度N这两个“慢变量”也变成求解器里的状态量。温度方程可以写成[ \frac{dT(t)}{dt} \frac{P_{abs}(t) \cdot R_{th} - (T(t) - T_{amb})}{\tau_{th}} ]P_abs 是腔内吸收光功率R_th 是热阻τ_th 是热时间常数T_amb 是环境温度。载流子方程类似[ \frac{dN(t)}{dt} -\frac{N(t)}{\tau_c} \frac{\beta_{TPA} |a(t)|^4}{2 \hbar \omega V_{eff}} ]右边第一项是载流子复合τ_c 是载流子寿命第二项是双光子吸收对载流子的产生速率β_TPA是TPA系数V_eff是有效作用体积。然后把这两个状态量映射回光场方程的频移项和损耗项。热光效应让谐振频率正比于温度偏移[ \Delta f_{th} - f_0 \cdot \frac{1}{n_g} \frac{dn}{dT} \cdot (T - T_{amb}) ]FCD效应让折射率下降对应频率偏移[ \Delta f_{FCD} - f_0 \cdot \frac{8.8 \times 10^{-22} N}{n_g} ]其中N的单位是cm-3。自由载流子吸收可以用额外损耗速率的一半来折算[ \gamma_{FCA} \frac{c}{n_g} \cdot 8.5 \times 10^{-18} \cdot N ]注意单位换算。这里N乘以前面系数得到的是吸收系数单位是1/cm再换成损耗速率就需要乘光速除以群折射率。很多复现代码跑出奇怪结果九成是这一串单位换算出了问题。2.3 一整套可求解的常微分方程组把2.1和2.2合在一起就得到了一个包含五个实数状态量的常微分方程组a的实部、a的虚部、温度T、载流子浓度N再加上最后一个对外部驱动频率或驱动相位的跟踪如果需要模拟chirp可以再补一个状态。写成矩阵形式在MATLAB里反而不好维护最好的做法是把它们打包在一个函数里比如function dstate ring_nl_ode(t, state, p) % state: [ar, ai, T, N] a complex(state(1), state(2)); T state(3); N state(4); dT (p.P_abs_func(a) * p.R_th - (T - p.T_amb)) / p.tau_th; dN -N / p.tau_c p.beta_tpa * abs(a)^4 / (2 * p.hbar * p.omega * p.Veff); df_th -p.f0 / p.ng * p.dndT * (T - p.T_amb); df_fcd -p.f0 / p.ng * 8.8e-22 * N; gamma_fca (p.c / p.ng) * 8.5e-18 * N; omega_shift 2 * pi * (p.df0 df_th df_fcd); gamma_total p.gamma_l p.gamma_tpa * abs(a)^2 gamma_fca; da (1j * omega_shift - gamma_total / 2) * a p.kappa * p.s_in(t); dstate [real(da); imag(da); dT; dN]; end这段代码的核心是把光场、热、载流子三个物理过程放在同一个时钟下推进。状态变量的维度不要想着用复数和实数混着存MATLAB的ode求解器对纯实数向量更友好。3. MATLAB代码实现从方程到可运行脚本3.1 核心参数定义与单位换算拿到代码包后我建议先把参数表完整筛查一遍因为它直接决定仿真结果对不对。常见的参数表长这样参数符号典型值单位备注环半径R5μm半径决定FSR和谐振位置群折射率n_g4.2-硅波导典型值热光系数dn/dT1.8e-41/K硅材料常数热时间常数τ_th200ns由衬底散热决定载流子寿命τ_c3ns受注入和复合机制影响TPA系数β_TPA0.5cm/GW波长相关的经验值线性损耗速率γ_l1e101/s由Q值反推耦合系数κ2e111/s由耦合区结构决定有效体积V_eff1e-12cm³粗略估算即可几个特别注意的单位坑。FCD公式里的载流子浓度N代码里通常以cm-3为单位但MATLAB变量里大家习惯直接写科学计数法一不小心就把1e17误写成1e15结果差了100倍。热阻R_th单位是K/W热方程里 P_abs 要用 W而不是mW。耦合系数κ到底取多少不能只看结构通常要跟实验测得的带宽和消光比做联立标定。如果代码包里参数文件写得比较随意我建议把参数结构体 p 单独写成一个脚本所有数值集中管理并用中文注释标注数据来源实验测试、文献、估算。这样后续调试才能快速追踪每一个可疑参数是从哪来的。3.2 主运行脚本怎么组织一个规范的主运行脚本大致分四步第一步定义参数结构体 p第二步定义输入信号 s_in通常是脉冲序列或正弦调制的光载波第三步设置初始状态 state0第四步调用ode求解器并做后处理。输入光场是调制器仿真的关键。实际代码包里通常会把激光器连续光写成 s_in(t) sqrt(P_in) * exp(1j * 2 * pi * f_opt * t)但如果你把光载波频率 f_opt 真实地放到方程里仿真步长会被迫到飞秒量级谁都跑不动。正确做法是只追踪“相对光频”的慢变包络把光载波频率从方程里提出来只保留相对谐振点的失谐量。这也是CMT模型的最大优势——你不需要真的去解光频振动只要解包络就行。光场包络的时间步长因此可以放到皮秒甚至纳秒量级极大地加速仿真。主函数代码可以这样写% 主脚本 ring_sim_main.m p ring_params(); % 加载参数 tspan [0 20e-6]; % 仿真20微秒观察热效应建立过程 s_in (t) sqrt(p.P_in) * exp(1j * 2 * pi * p.df_drive * t); state0 [0; 0; p.T_amb; 0]; % 初始无光、常温、无载流子 [t, state] ode15s((t, s) ring_nl_ode(t, s, p, s_in), ... tspan, state0, odeset(RelTol, 1e-6));这里的 s_in 写成匿名函数好处是可以随时换调制信号形式比如改成归零码脉冲、高斯脉冲或者直接从外部文件读入的任意波形。3.3 运行结果和后处理到底要看哪些图代码包跑通后第一件事不是去调参而是确认输出波形是否符合物理直觉。最基本的三个后处理图是时域输出光功率、谐振频率偏移随时间变化、透过率曲线。输出光功率从腔内场算出带耦合输出系数的平方即可。谐振频率偏移要把 Δf_th 和 Δf_FCD 分别画出来这样能直观看到热效应和载流子效应谁主导。如果输入功率只有1 mW你会看到频率偏移主要是热效应如果把功率拉到20 mWFCD的贡献会明显抬头而且它有更快的上升沿。透过率曲线可以做两种一是固定频率扫描观察到的是稳态谱线被非线性“压宽”或“劈裂”二是在时间轴上取快照看到的是瞬态谱线。后者才是时域模型独有的优势稳态仿真永远给不出来。调试阶段我还习惯在代码里加一个能量守恒检查输入的能量应该等于输出能量加腔内损耗加吸收转换的热量。如果仿真结果里能量差超过百分之几说明时间步长不够或者方程少了损耗项。4. 常见报错与调试经验4.1 数值发散与刚性系统怎么处理第一次跑这个模型最常见的问题就是数值爆炸输出功率变成NaN或者无穷大。原因基本出在两个地方。第一个是时间尺度差异导致的刚性。光场包络的响应时间在几十皮秒量级热效应却是几百纳秒载流子又是几纳秒三者可能差四五个数量级。你用普通的ode45跑它会为了满足精度把步长拖得非常小然后出现步长收敛失败或者干脆发散。解决办法是换用ode15s或者ode23t这类刚性求解器它们对慢变量和快变量耦合的处理好得多。第二个是初始腔内能量过高。如果初始状态给了一个很大的腔内场但输入光功率很小方程会先把腔内能量耗散掉这个瞬态过程在某些显式求解器里会震荡。保险做法是初始状态直接用0让光从零开始建立。还有一个很多人忽略的问题是输入信号本身有突变。如果你给的是理想方波脉冲光场包络在上升沿和下降沿会有非常陡的梯度数值上等效于高频分量容易激发振荡。解决办法是把脉冲边沿用高斯或余弦滚降函数修一下哪怕只做1 ps的rise time数值稳定性都会有明显改善。提示遇到NaN时先检查是不是时间步长太大导致载流子浓度出现负数。载流子浓度在物理上不可能小于0但在数值上可能越过零点继续下降最终把折射率扰动项推成一个荒谬的正数系统就炸了。处理方法是在状态方程函数里对负数载流子浓度做钳位N max(N, 0)。4.2 参数标定从文献和实验里拿到靠谱数值代码包给了一组默认参数但你的实际器件参数大概率跟它不同。我见过不少复现者直接拿默认参数去对实验曲线对不上就开始怀疑程序写错了。其实问题往往出在参数标定环节。Q值到损耗速率的换算是第一道坎。微环的加载Q值 Q_loaded 跟总损耗速率的关系是[ Q \frac{\omega_0}{\gamma_{l} \kappa_c} ]这里的 γ_l 是本征损耗速率κ_c 是耦合损耗速率。实验上测到的透过率谱线深度和3 dB带宽能够反推耦合条件和本征损耗但要注意区分过耦合、欠耦合和临界耦合三种状态。不同状态下同样的谱线深度对应完全不同的耦合系数。热时间常数最好不要直接抄文献数值。它主要由衬底和埋氧层厚度决定不同工艺平台可以差到一个数量级以上。自己用泵浦-探测实验测一下热建立时间或者查同一家Foundry的PDK报告得到的数据会比文献平均值可靠很多。载流子寿命同样和工艺强相关。PN结调制器的耗尽区、PIN结构的注人区载流子寿命都不一样。如果代码包的默认值是3 ns而你的器件是10 ns仿真得到的寄生损耗和频移会明显偏大。4.3 仿真结果跟实验对不上怎么办这可能是全流程最磨人的一环。仿真和实验对不上通常先检查坐标轴和归一化条件。很多实验测的是相对透过率而仿真输出的是绝对功率。如果没有把耦合损耗、光纤到芯片的耦合效率算进去两条曲线就差一个固定比例。其次是驱动信号的定义。实验里调制器是电压驱动PN结的载流子变化是通过耗尽区宽度改变实现的时域模型里你用的是折射率扰动或者失谐频率变化。这两个量之间的换算关系不是线性的尤其在大信号调制下由PN结C-V曲线带来的非线性会让调制器的啁啾和消光比都发生形变。如果代码包只给了一个线性系数你在做大信号验证时必然会失配。还有一个隐蔽问题是波长参考点。仿真里你设定的是绝对波长或绝对频率但实验里激光器锁定的是某个相对位置比如谐振峰的-3 dB点。如果两者对不上初始失谐量就会差很多输出的动态特性自然天差地别。排查时把你的初始失谐量跟实验校准点的谱线位置换算到同一个基准上再比较。调试时我习惯先做三个“固定测试”小功率连续光扫描验证线性谱线、低速率方波调制验证热效应和载流子效应的基本时间响应、大功率连续光下的光谱畸变验证非线性强度。这三个测试能过再进高速调制和系统级仿真能省掉大量摸索时间。5. 一点个人的建模心得这个模型跑通以后大多数人会想把它接到更高层级的系统仿真里比如做眼图、算BER、看非线性串扰。我的建议是先不要让这个时域模型直接跟DSP算法跑在一起那样计算量会很大而且容易把问题搅浑。正确做法是先用这个模型生成一组关键工况下的响应数据做成查找表或者简化等效模型再嵌入到链路仿真里。具体来说可以在仿真结束后把不同输入功率、不同失谐频率下的透过率包络存成三维数组然后在上层仿真里用插值方式调用。实测下来这种做法能把计算速度提升两个数量级以上同时保留非线性模型的关键动态特性。如果后面想深入研究载流子效应和热效应的相互作用还可以在代码包里加上功放效应或两光子吸收饱和的修正项扩展性很好。最后再分享一个小技巧用ode15s跑刚性系统时如果愿意提供状态方程对状态变量的解析雅可比矩阵求解速度还能再快一倍。手动推导这些偏导确实有点烦但跑长仿真或者做参数扫描的时候省下来的时间绝对值得。本文还有配套的精品资源点击获取