1. 为什么马赫-曾德干涉仪值得单独建模仿真我真正下决心把马赫-曾德干涉仪Mach-Zehnder InterferometerMZI放进波动光学仿真环境里完整跑一遍是在一次折射率传感实测翻车之后。当时测出来的所有特征曲线都指向同一个结论——理论上输出光强应该随溶液浓度近似线性变化但实验数据里明显叠了一层周期性的起伏。排查了很久才意识到问题出在合束端两条高斯光束并没有严格同轴夹角虽然只有毫弧度级别却直接决定了干涉条纹的空间频率而条纹一旦落到探测器靶面上就成了那个多余的周期项。这种问题靠解析公式手算非常痛苦换成复振幅场仿真之后几分钟就能把整个因果链条理清楚。这篇文章我会从波动光学的角度把马赫-曾德干涉仪的建模思路、数学基础、Python仿真实现和结果分析完整拆开讲一遍。内容适合三类人一是刚接触干涉测量、想知道条纹到底怎么形成的初学者二是在实验室里调过MZI、但对仿真手段不熟的实验物理或光学工程师三是做集成光子学设计、需要快速评估结构参数的芯片设计者。无论你是哪一种看完之后应该都能自己动手搭出一套可复用的MZI仿真流程。1.1 双光路结构为什么偏偏选MZI马赫-曾德干涉仪的基本结构很直观一束光先经过第一个分束器BS1被分成两束分别走两条独立的光路然后由第二个分束器BS2重新合束在两个输出端口形成干涉光场。与迈克尔逊干涉仪不同MZI的参考光和信号光是在空间上分开的两条臂而且都是透射式结构不依赖反射镜的回程对准。这个特点让它在工程上有几个不可替代的优势两条臂可以分别独立改造比如在一条臂里放样品池、加电光晶体或者引入折射率变化另一条臂保持参考状态测量时互不干扰没有返回光避免了对光源的反馈扰动对激光器稳定性友好两个输出端口天然带有互补特性一路增强时另一路减弱配合平衡探测器就能做共模抑制把光源噪声压下去。正因为两条臂物理上分开了仿真模型也必须把两条独立传播路径这个结构如实还原出来。不能像处理法布里-珀罗腔那样只在一个方向上来回叠而是要分别计算两束光各自的振幅、相位和空间分布最后再做场的叠加。这个差别听起来不大但它决定了整个仿真代码的组织方式——两条臂就是两个独立传播分支。1.2 波动光学与几何光学的分界线在哪里说到仿真干涉仪会有人问直接用几何光学光线追迹不就行了答案是不行。干涉本身是一种波动现象它的数学本质是复振幅的叠加而光线追迹只跟踪能量传播方向完全不携带相位信息。换句话说几何光学可以把光走到哪算清楚却算不出两束光叠在一起是变亮还是变暗。具体到MZI场景两条臂的光程差哪怕只差一个波长量级输出光强就会在极大和极小之间来回摆动。这个周期性行为完全由相位差驱动属于波动的特征。只有两种情况下几何光学够用一是特征尺寸远大于波长且不需要关心光场分布细节二是系统中没有相干叠加环节只看能量分配。MZI显然两个条件都不满足所以从一开始就应该用波动光学的场仿真来做。1.3 仿真到底能回答哪些工程问题我搭这套仿真最想回答的问题不是理想情况下条纹长什么样——那个用平面波公式几分钟就能算出来而是几个更接地气的问题当入射光是高斯光束而不是理想平面波时输出干涉图样与教科书公式差多少分束比如果偏离50:50条纹可见度会掉到多少合束端存在微小夹角时条纹间距和方向如何变化一条臂里加载了折射率变化后输出光强随相位的变化曲线是什么形状哪个参数对可见度最敏感设计容限在哪里这些问题里有些可以用解析公式估算但一旦涉及二维光场分布、高斯包络、非理想对准解析解就会变得非常繁琐甚至没法闭式表达。仿真好在能把所有实际因素放进同一个框架里改一个参数重跑一遍五分钟之内就能看到趋势。这就是我觉得MZI值得专门建模仿真的核心原因——它的结构足够简单物理图像清晰但工程细节足够丰富是练习波动光学仿真的好题目。2. 从复振幅到干涉条纹MZI的数学建模要点仿真做之前先把数学模型理清楚。MZI的一切行为都可以归结为一句话输出场是两束复振幅的相干叠加。理解这句话后面所有代码都只是它的翻译。2.1 干涉项从哪来假设两束光到达合束端时的复振幅分别是E₁和E₂输出光场就是E_out (E₁ E₂) / √2这里除以√2来自50:50合束器的场振幅分配。光强是振幅模平方I |E_out|² I₁ I₂ 2√(I₁I₂)·cos(Δφ)其中Δφ是两束光在合束端的相位差。最关键的是第三项——交叉项它不依赖于任何光束重叠之外的物理过程单靠波的线性叠加就自然出现。如果中间没有这个余弦项两束光合起来只是能量相加那就退化成非相干光的行为也就没有干涉可言了。所以做仿真时一定要保证在复振幅层面做加法而不是在光强层面做加法。这也是很多初学者第一版代码跑不出来条纹的原因——他们分别算了I₁和I₂再相加交叉项永远出不来。2.2 相位差的三种来源与条纹形态在MZI里相位差Δφ通常有三个来源它们对应的干涉图样特征完全不同第一是臂长差。两条臂的物理长度差ΔL直接换算成相位Δφ (2π/λ)·ΔL如果两束光完全同轴合束这个相位差只会让输出端整体变亮或变暗不产生空间条纹。这也是MZI做传感的基础——只要一条臂的有效光程发生变化输出光强就在相长和相消之间切换。第二是折射率差。一条臂里介质折射率从n₁变成n₂等效光程变化是ΔL·(n₂ − n₁)。生物传感、气体检测用的都是这个机制。折射率的变化往往很小10⁻⁴量级所以系统需要很高的相位灵敏度这又反过来要求干涉条纹对比度足够高。第三是合束端的夹角。两束光不是严格平行而是存在一个夹角θ时它们在探测器平面上的波前就有了一个线性的相对相位梯度。假设夹角在x方向干涉光强分布变成I(x) I₁ I₂ 2√(I₁I₂)·cos(kx·sinθ Δφ₀)此时你会看到沿x方向的等间距直线条纹条纹间距d λ / sinθ这个公式非常实用。θ1 mrad、λ632.8 nm时d约0.63 mm——在毫米量级的探测器靶面上就能看到好几条条纹。实验里MZI出条纹通常不是因为故意制造夹角而是因为合束镜没调好残留了一个毫弧度级别的角度。仿真里可以精确控制这个量直接看出它对图样的影响。2.3 分束比决定可见度干涉条纹的质量通常用可见度V描述V (I_max − I_min) / (I_max I_min)假设两束光场振幅比例由第一个分束器决定光强分束比是t₁:t₂t₁t₂1那么可见度可以推导成V 2√(t₁t₂) / (t₁ t₂)下面这个表格是直接算出来的分束比 t₁:t₂理论可见度 V50:501.000060:400.979870:300.916580:200.800090:100.6000注意一个反直觉的点分束比偏到60:40时可见度仍然有0.98看上去没太大损失但到80:20就掉到0.8到90:10只剩0.6。如果你的应用对对比度要求很高比如平衡探测系统分束器的一致性反而是比臂长差更需要盯紧的参数。2.4 高斯光束传播选择角谱法的理由实际光源输出的是高斯光束不是平面波。高斯光束在自由空间的传播可以用解析公式描述但一旦涉及二维光场、任意相位分布、非理想波前解析公式就不够通用了。仿真里我选择角谱法Angular Spectrum MethodASM来处理传播问题。角谱法的思路是把任意光场分解成一组平面波的叠加每个平面波有一个传播方向对应的空间频率为(fx, fy)。自由空间传播只是给每个平面波乘一个相位因子H(fx, fy) exp(ikz·√(1 − λ²fx² − λ²fy²))然后逆变换回空间域。整个过程就是两次傅里叶变换各乘一个因子一次调用就完成任意距离的传播。相比菲涅尔衍射积分角谱法对传播距离没有近轴近似强加的限制只要采样满足要求它既能处理近距离也能处理远距离非常适合MZI这种两臂分别传播、最后合成的链路型仿真。3. 用Python搭建波动光学仿真环境理论模型有了接下来就是把数学翻译成代码。我用的工具栈非常朴素Python NumPy Matplotlib。没有选商用光学软件因为这个场景本质上是二维复振幅数组的运算NumPy天然就是干这个的。3.1 技术选型为什么不复杂有人会觉得仿真干涉仪应该上FDTD或者COMSOL。其实要看仿真对象是什么层级。FDTD是在麦克斯韦方程层面求解适合处理亚波长结构、散射、模式耦合这类问题但计算代价大一个MZI结构里如果大部分是厘米量级的自由空间光路FDTD会跑到天荒地老。MZI的经典玩法——自由空间传播加相干叠加——在标量衍射理论框架下就能很准确地描述Python原生的FFT操作完全够用。如果你的场景换成了集成光子学芯片上的波导型MZI那是另一个话题那确实需要FDTD或者波束传播法BPM后面我会单独提到。但对于空间光路MZINumPy是最快出成果的路径。3.2 网格、波长与采样率的匹配仿真前三个核心参数必须定下来波长λ、仿真窗口尺寸L、网格点数N。我常用的配置是这样的参数取值说明λ632.8 nm氦氖激光波长L8 mm窗口覆盖光束直径的3~4倍N1024每维网格数dx7.8125 μm空间采样间隔网格间距dx要能分辨你关心的最小结构。对MZI来说最小结构通常是干涉条纹条纹间距dλ/sinθ。当θ1 mrad时d0.633 mm远大于dx7.8 μm采样绰绰有余。但如果你做的是高倍显微镜成像或者亚波长结构dx就得往λ/2量级走计算量会急剧上升。还有一个容易被忽略的约束来自角谱法的传递函数采样。菲涅耳近似下角谱传递函数是空间频率的二次相位项要正确采样它传播距离z需要满足z ≤ N·dx² / λ代入上面参数N·dx²/λ 1024×(7.8125e-6)²/632.8e-9 ≈ 0.098 m。所以我的传播距离取0.05 m是安全的如果非要传播0.5 m而不改网格结果会出现明显的数值畸变。这个约束很多人不知道是仿真结果莫名其妙发散的经典原因。3.3 角谱传播函数与FFT坐标陷阱角谱传播的实现只有三句话但里面有一个人人都踩过的坑FFT的坐标约定。NumPy的fft2默认把数组的第0个索引当作原点而我们构造的高斯光束中心是在数组正中间索引N//2。如果直接把数组喂给fft2空间中心就对不上频域原点结果就是光场整体偏移或者传播后出现错位。正确的做法是用ifftshift把中心挪到原点做完变换再用fftshift把结果挪回中心。下面这段是我实际在用的传播函数import numpy as np def propagate(E, z, wavelength, dx): 角谱法自由空间传播 E: 空间域复振幅数组 z: 传播距离米 wavelength: 波长米 dx: 空间采样间隔米 N E.shape[0] k 2 * np.pi / wavelength # 频域坐标FFT ordering fx np.fft.fftfreq(N, ddx) FX, FY np.meshgrid(fx, fx) F2 FX**2 FY**2 # 角谱传递函数 tmp 1 - (wavelength**2) * F2 H np.zeros_like(tmp, dtypecomplex) mask tmp 0 H[mask] np.exp(1j * k * z * np.sqrt(tmp[mask])) # tmp 0 的区域是倏逝波传播中衰减直接置零 # 域切换ifftshift 让光场中心对准 FFT 原点 E_f np.fft.fft2(np.fft.ifftshift(E)) Ez np.fft.ifft2(E_f * H) return np.fft.fftshift(Ez)这里对倏逝波的处理我直接置零了。对自由空间毫米级光路来说倏逝波只出现在近场波长尺度内传播几十毫米完全可以忽略。如果仿真对象里有亚波长缝隙结构那就得另说。高斯光束初始场的构造同样要注意坐标原点# 空间坐标网格 x (np.arange(N) - N // 2) * dx X, Y np.meshgrid(x, x) # 初始高斯光束束腰半径 w0 w0 1.2e-3 E_in np.exp(-(X**2 Y**2) / w0**2)这段代码生成的E_in最高点在数组正中央物理坐标上对应x0、y0。配合前面的ifftshift坐标约定就完全统一了。我最早在这上面浪费过整整一天——光场传播后不是应该在中心继续发亮吗结果出来一个往角落跑的亮斑排查半天就是少了ifftshift。4. 完整仿真从分束到合束的MZI实现数学和工具都准备好了现在把整个MZI链条串起来。这一步最爽的地方在于真实实验里每个环节都是物理器件仿真里每个环节都只是一行NumPy操作。4.1 分束器和合束器如何抽象理想分束器的作用是把入射场的振幅按比例拆分并在其中一个出口引入一个π/2相位细节取决于器件实现。但在标量模型中我们通常只关心强度分配所以分束器可以抽象成一个简单的系数乘法# 第一个分束器光强按 t1:t2 分配 E_arm1 np.sqrt(t1) * E_in E_arm2 np.sqrt(t2) * E_in如果t1t20.5就是标准的50:50分束。为什么要开根号因为系数作用在振幅上而分束比通常按光强定义所以要用√t。这个细节搞反了后面算可见度会对不上。合束器稍微复杂一点因为它有两个输出端口# 第二个分束器合束两个输出端口互为反相 E_out1 (E_arm1 E_arm2) / np.sqrt(2) E_out2 (E_arm1 - E_arm2) / np.sqrt(2)端口1对应同相叠加端口2对应反相叠加。注意当两条臂相位差为0时端口1是相长干涉亮端口2完全相消暗能量守恒就在这两个端口之间流动。4.2 两条臂的传播与相位加载分束完成后两束光分别走各自的臂。对于自由空间MZI每条臂就是一截传播距离对于集成光学MZI每条臂是一段波导。仿真时两类情况都能处理——自由空间调传播距离波导就调相位因子。我把臂长差、折射率变化和合束夹角都做了进去代码分三步# 臂1传播 z1不做额外处理 z1 0.05 E_arm1 propagate(E_arm1, z1, wavelength, dx) # 臂2传播 z2与臂1有细微长度差再叠加折射率相位和合束夹角 z2 0.05 E_arm2 propagate(E_arm2, z2, wavelength, dx) # 折射率变化导致的全局相位比如相位改变 0.8 rad E_arm2 E_arm2 * np.exp(1j * 0.8) # 合束端夹角 theta_x给一个线性相位斜坡模拟波前倾斜 theta_x 1e-3 # 1 mrad E_arm2 E_arm2 * np.exp(1j * k * X * np.sin(theta_x))这里的逻辑是折射率变化让整束光整体带一个常数相位合束端夹角让波前带上一个空间线性相位。它们分别对应着条纹整体明暗移动和条纹空间分布两个截然不同的现象可以独立开关、单独观察。4.3 输出端口分析与能量守恒校验第一次跑通之后别急着去看条纹先做一步校验能量守恒。两条臂的能量加起来应该等于输入能量两个输出端口的能量之和也应该等于输入能量。在复振幅模型里这是个很好的自检手段。I_total_in np.sum(np.abs(E_in)**2) * dx**2 I_total_out (np.sum(np.abs(E_out1)**2) np.sum(np.abs(E_out2)**2)) * dx**2 print(f输入能量: {I_total_in:.6e}, 输出能量: {I_total_out:.6e})如果这两者对不上排除数值误差后最常见的原因就是前面的ifftshift/fftshift坐标问题导致光场在传播时发生了人为的边界截断或者相位畸变。能量守恒校验通过后再看条纹才有意义。完整代码我整理成了这样一个主流程# MZI 完整仿真主流程 wavelength 632.8e-9 L 8e-3 N 1024 dx L / N x (np.arange(N) - N // 2) * dx X, Y np.meshgrid(x, x) w0 1.2e-3 E_in np.exp(-(X**2 Y**2) / w0**2) # 50:50 分束 t1, t2 0.5, 0.5 E_arm1 np.sqrt(t1) * E_in E_arm2 np.sqrt(t2) * E_in # 两臂传播 E_arm1 propagate(E_arm1, 0.05, wavelength, dx) E_arm2 propagate(E_arm2, 0.05, wavelength, dx) # 信号臂加载相位0.8 rad 夹角1 mrad E_arm2 E_arm2 * np.exp(1j * 0.8) E_arm2 E_arm2 * np.exp(1j * k * X * np.sin(1e-3)) # 合束 E_out1 (E_arm1 E_arm2) / np.sqrt(2) E_out2 (E_arm1 - E_arm2) / np.sqrt(2) I1 np.abs(E_out1)**2 I2 np.abs(E_out2)**2跑完这十几行你就拿到了一张包含全部干涉信息的二维光强图。接下来的问题就是怎么看它、怎么量化它。5. 结果解读与仿真中的常见陷阱仿真跑通只是开始真正的功夫在结果解读上。这一节我把几个典型场景的结果特征和量化方法讲清楚顺便把我实际踩过的坑都倒出来。5.1 等倾条纹的频率验证第一个必做的验证是仿真出来的条纹间距和理论值dλ/sinθ对不对得上。取输出面中心行做一维强度切片就能看到周期性的峰和谷。我上面的参数θ1 mradd理论值约0.633 mm在8 mm窗口内应该出现大约12到13个周期。用这个验证最大的好处是能暴露整套仿真的系统误差。如果d偏大或偏小说明频域坐标或者传递函数写错了如果条纹方向不对说明线性相位斜坡的坐标轴搞反了如果条纹数量对但对比度异常说明合束系数或者分束比写错了。一次对不上通常就能顺着d的偏差方向猜出是哪个环节的问题。我见过一个很隐蔽的错误有人把相位斜坡写成np.exp(1j * k * X * theta)但X的单位是米、theta的单位是度没转换结果条纹间距莫名其妙大了57倍还以为是程序没错、理论错了。参数单位一致性永远是第一排查对象。5.2 分束比与可见度的定量对照第二个必做的验证是把仿真可见度和2.3节的公式对照。我把分束比从50:50一路改到90:10记录中心区域的可见度和理论曲线叠在一起画基本是重合的。这验证了分束器抽象和合束系数是正确的。这里要提醒一句可见度的计算窗口很讲究。如果在整幅图像上取全局max和min高斯光束的边缘暗区会把I_min拉得极低导致可见度虚高。正确做法是在光强包络相对平坦的中心区域取局部的峰谷值或者沿中心行在单个条纹周期内做滑动窗口统计。这个细节在真实实验里同样存在——探测器如果覆盖了整个光斑边缘测出来的对比度往往偏乐观。5.3 高斯包络对可见度测量的干扰高斯光束和理想平面波最大的区别是光强分布不均匀。用平面波公式算出来的可见度是纯常数而高斯光束的可见度在空间上是包络调制的。光束中心附近条纹峰值高、边缘条纹被高斯包络压低如果你只取离中心远的一行算可见度会明显低于理论值。做仿真分析时要么明确说自己测量的是中心区域局部可见度要么用二维条纹图样做空间滤波把载波频率以外的包络成分分离再在载波信号上算对比度。实验室里如果真的要用MZI定量测相位最省事的办法其实是把高斯光束先经过空间滤波扩束成近似均匀的平顶光再进干涉仪避免包络影响。5.4 我踩过的几个坑这一节是血泪总结按坑的典型程度排序。第一个坑一维数组当二维用。NumPy的meshgrid用不对X和Y的shape不匹配所有乘法都会广播成奇怪的矩阵。我那次的症状是条纹只在对角线方向出现怎么调都调不直。检查才发现Y坐标忘记转置了。第二个坑倏逝波区域的sqrt。角谱传递函数里当λ²(fx²fy²)1时根号下面是负数np.sqrt会返回nan或者触发warning然后整张图报废。解法就是我前面代码里写的用mask把倏逝波区域直接置零。这个问题在传播距离大、网格细的时候几乎必然出现。第三个坑传播距离超限。3.2节那个z ≤ N·dx²/λ的约束不是摆设。有一次我为了模拟长臂差直接设z0.5 m忘了改网格出来的结果条纹数量对但每个条纹边缘都带锯齿状的数值噪声。后来把N从1024提到4096问题立刻消失。记住改传播距离时先检查这个不等式。第四个坑把能量加在光强上而不是振幅上。分束器的√t忘记开根号直接np.sqrt(t)是对的但有人会写成t*E_in然后可见度公式怎么都对应不上因为光强分配变成立方关系了。每个系数都要想清楚它作用在振幅还是光强。6. 从单元仿真走向系统应用MZI的仿真模型一旦跑通你会发现它像一个乐高积木一样可以往各个方向拓展。最后这部分聊聊我在实际项目中怎么用这套仿真以及往更复杂系统走的时候会遇到什么。6.1 折射率传感把相位变化变成测量信号我最初做MZI仿真就是为折射率传感。思路很简单样品臂放入待测溶液折射率变化Δn让这一臂的有效光程变化等效相位变化Δφ (2π/λ)·Δn·L_sample。在仿真里这条臂就多乘一个exp(iΔφ)。关键是要搞清楚工作点的问题。MZI输出光强对相位是余弦响应灵敏度最高的点出现在Δφπ/2附近余弦的拐点但那里响应也是非线性的。为了稳定工作通常要人为把工作点偏置到正交位置。仿真里你能很直观地看到如果只扫Δφ输出是一条余弦曲线但如果加上合束夹角的条纹你在某个探测器位置采到的强度随Δφ的曲线还叠加了探测器落在条纹暗区带来的衰减。这种条纹相位和全局相位耦合的效应不仿真很容易忽略。6.2 电光调制器MZM怎么仿真在通信领域MZI最常见的身份是马赫-曾德调制器。铌酸锂或者硅波导里一条臂上加了电极电场改变折射率从而改变相位输出就随调制电压在亮暗之间切换。这个场景里光路不再是自由空间而是波导。我的做法是把上一节的空间传播换成波导模式的有效折射率修正传播相位用neff·k·L调制臂的相位加载用电压产生的Δneff换算。空间自由度和波导自由度在这个层面可以解耦——先算横向模式再算轴向传播。这套半解析方法在器件设计初期比FDTD快几个数量级精度对大多数工程评估完全够用。6.3 集成光子学场景下的扩展方向如果你真的在开发硅光芯片上的MZI那自由空间角谱法就到头了。硅波导的尺寸在几百纳米量级和非线性材料、波导侧壁粗糙度等因素耦合需要用FDTD或BPM来求解。这时候自由度从二维标量场变成三维矢量场计算资源需求会跳好几个量级。我的建议是分层次走先用本文这种快速模型把系统级参数分束比、臂长差、相位加载扫一遍缩小参数空间再用FDTD对少数关键结构做精细验证。这样的组合拳在工程上既快又不失精度。另外现在很多集成光子学设计还会配合逆向设计算法自动优化分束器结构这时候仿真模型的速度直接决定了优化能不能在可接受时间内收敛。仿真做到最后我的感受是MZI是个好老师它结构简单到你能完全掌握每一个参数的物理含义又复杂到能覆盖波动光学里几乎所有核心概念——复振幅、相位、干涉、衍射、相干性。把这套仿真吃透以后遇到更复杂的干涉系统建模思路基本是相通的。至少对我来说搭完这套流程之后再回头看那些实测曲线里的奇怪周期项一眼就能认出它们的真面目了。