4. 马赫-曾德尔干涉仪中相位估计精度的仿真分析
双模马赫-曾德尔干涉仪相位估计模型
本例在双模（双端口）马赫-曾德尔干涉仪（MZI）框架下，仿真相位扫描时的差分光电流及其噪声，并计算相位灵敏度，与散粒噪声极限（SNL）对比。MZI 由两个 50:50 分束器与两臂相位差 φ 构成：第一分束器把输入模 a、b 混合，两臂累积相对相位 φ（其中一臂含待估相位），第二分束器再混合，最后由平衡探测器取两输出端口光强之差，得到差分光电流算符
Î− = n̂c − n̂d = a†b e−iφ + b†a eiφ
对相干态输入，⟨Î−⟩ ∝ sin φ，方差为常数，相位灵敏度 δφ = ΔÎ− / |∂⟨Î−⟩/∂φ| 在 φ=π/2 附近达到 SNL：δφSNL = 1/√N（N 为总平均光子数）。当一臂注入压缩真空态时，噪声在特定相位被压缩，灵敏度可低于 SNL。下面用截断 Fock 空间数值实现该模型。
数值实现：截断 Fock 空间与算符构造
取每个模的截断维数 D（例如 D=40），在双模基 |na, nb⟩ 上构造湮灭算符 a、b 及分束器幺正变换。50:50 分束器写作 B = exp[i(π/4)(a†b + b†a)]，两臂相位算符 P(φ) = exp[iφ a†a]（相位累积在 a 臂）。差分光电流算符 Î− = a†b e−iφ + b†a eiφ 可直接由 a、b 与 φ 构造。输入态取相干态 |α⟩a ⊗ |0⟩b（相干态输入）或 |α⟩a ⊗ |ξ⟩b（一臂压缩真空）。
import numpy as np
from scipy.linalg import expm
D = 40 # 每模截断维数
N = D * D # 双模希尔伯特空间维数
a = np.zeros((D, D)); b = np.zeros((D, D))
for n in range(1, D):
a[n-1, n] = np.sqrt(n) # a|n> = sqrt(n)|n-1>
b[n-1, n] = np.sqrt(n)
I = np.eye(D)
# 双模算符（Kronecker 积，a 在前、b 在后）
A = np.kron(a, I); B = np.kron(I, b)
Ad = A.conj().T; Bd = B.conj().T
def bs_unitary(theta=np.pi/4):
# 50:50 分束器：B = exp[i*theta*(a†b + b†a)]
H = Ad @ B + Bd @ A
return expm(1j * theta * H)
def phase_op(phi):
# 相位累积在 a 臂：P = exp[i*phi*a†a]
return expm(1j * phi * (Ad @ A))
def coherent_state(alpha, D):
n = np.arange(D)
c = np.exp(-0.5*np.abs(alpha)**2) * alpha**n / np.sqrt(np.math.factorial(n))
return c
def squeezed_vacuum(r, D):
# 压缩真空态在 Fock 基的系数（偶数光子）
n = np.arange(0, D, 2)
c = np.zeros(D, dtype=complex)
c[n] = (1/np.cosh(r)) * (np.sqrt(np.math.factorial(n)) /
(2**(n/2) * np.math.factorial(n//2))) * (np.tanh(r))**(n/2)
return c
def mzi_output(state_in, phi):
# 输入 -> BS1 -> 相位 -> BS2
U = bs_unitary() @ phase_op(phi) @ bs_unitary()
return U @ state_in
def diff_current(state, phi):
# 差分光电流算符 Î_- = a†b e^{-iφ} + b†a e^{iφ}
Op = Ad @ B * np.exp(-1j*phi) + Bd @ A * np.exp(1j*phi)
mean = np.real(state.conj() @ Op @ state)
var = np.real(state.conj() @ (Op @ Op) @ state) - mean**2
return mean, var
# 相干态输入：|α>_a ⊗ |0>_b
alpha = 2.0
psi_coh = np.kron(coherent_state(alpha, D), np.array([1.0]+[0.0]*(D-1)))
# 压缩真空输入：|α>_a ⊗ |ξ>_b，r=0.8
r = 0.8
psi_sq = np.kron(coherent_state(alpha, D), squeezed_vacuum(r, D))
phis = np.linspace(0, 2*np.pi, 181)
res_coh = np.array([diff_current(mzi_output(psi_coh, p), p) for p in phis])
res_sq = np.array([diff_current(mzi_output(psi_sq, p), p) for p in phis])
# 相位灵敏度：δφ = ΔÎ_- / |∂⟨Î_-⟩/∂φ|
dmean = np.gradient(res_coh[:,0], phis)
sens_coh = np.sqrt(res_coh[:,1]) / np.abs(dmean)
dmean_sq = np.gradient(res_sq[:,0], phis)
sens_sq = np.sqrt(res_sq[:,1]) / np.abs(dmean_sq)
# 散粒噪声极限：δφ_SNL = 1/sqrt(N)，N = |α|^2
N_ph = np.abs(alpha)**2
snl = 1/np.sqrt(N_ph)
print("SNL =", snl)
print("相干态最小灵敏度 =", sens_coh.min())
print("压缩态最小灵敏度 =", sens_sq.min())相位扫描结果与灵敏度判读
运行上述代码（教学示例参数：α=2.0，N=4，r=0.8，D=40）可得到以下预期现象：
- 差分光电流均值：⟨Î−⟩ 随 φ 呈正弦振荡，在 φ=π/2 处斜率最大，是相位估计的最佳工作点。
- 噪声：相干态输入时 ΔÎ− 近似恒定；压缩真空输入时噪声随 φ 变化，在特定相位被压缩至低于相干态水平。
- 灵敏度：相干态输入的最小 δφ 应接近 SNL=1/√4=0.5；压缩态输入在最佳相位处 δφ 可低于 SNL，体现压缩态对相位估计的增强。
判读要点：若数值结果偏离理论值，优先检查截断维数 D 是否足够（D 过小会截断相干态尾部与压缩态高光子分量，导致方差偏小或均值偏差）；其次确认相位算符与差分电流算符的相位符号约定一致（本例相位累积在 a 臂，差分电流中 e−iφ 对应 a†b 项）。
结果验证与常见偏差
验证数值实现是否正确，可从三个层面检查：
- 算符一致性：在 φ=0 时，Î− 应等于 a†b + b†a；对相干态输入，⟨Î−⟩ 在 φ=0 处应为 0（两臂对称），在 φ=π/2 处取极值。
- SNL 校验：相干态输入下，最小灵敏度应接近 1/√N。若明显偏小，通常是截断维数 D 不足导致相干态被截断、方差被低估；增大 D 后结果应收敛。
- 压缩增强校验：压缩真空输入时，最佳相位处的 δφ 应低于相干态结果；若未出现，检查压缩态系数是否归一化、相位算符是否作用在正确模式上。