3. 分束器对相干态与压缩态干涉的仿真实现

任务拆解与物理约定

本案例要完成三件事：用光子数基截断构造 50:50 分束器的幺正矩阵 U_BS；分别把相干态 |α⟩ 与压缩真空态 S(r)|0⟩ 送入分束器的一个输入端口（另一端口为真空）；在输出端统计两模光子数分布，并对每个输出模做正交分量方差计算，观察干涉对量子噪声的压缩/反压缩效应。

约定：分束器两输入模记为 a（信号）、b（真空），输出模记为 c、d。采用对称分束器变换

c = (a + b)/√2
d = (a − b)/√2

对应幺正矩阵（作用于列向量 [a, b]ᵀ）：

U_BS = (1/√2) * [[1, 1],
                  [1, -1]]

正交分量定义为 X = (a + a†)/√2、P = (a − a†)/(i√2)，真空方差为 1/2（取 ħ=1）。压缩真空态满足 Var(X) = e^(−2r)/2、Var(P) = e^(2r)/2。

数值实现：截断维数、算符与态矢量

取截断维数 N = 40（单模 0…39 个光子）。湮灭算符用次对角线 √(n+1) 构造。相干态系数为 c_n = e^(−|α|²/2) α^n / √(n!)；压缩真空态系数为 c_(2m) = (tanh r)^m √((2m)!)/(2^m m! √(cosh r))，奇数项为 0。两模态用 Kronecker 积展开为 1600 维矢量，再作用 U_BS ⊗ U_BS 的共轭转置完成分束器演化。

import numpy as np
from math import factorial, sqrt, tanh, cosh

N = 40
a = np.diag(np.sqrt(np.arange(1, N)), 1).astype(complex)   # 湮灭算符
I = np.eye(N, dtype=complex)

def coherent(alpha):
    n = np.arange(N)
    c = np.exp(-abs(alpha)**2/2) * alpha**n / np.sqrt([factorial(k) for k in n])
    return c / np.linalg.norm(c)

def squeezed_vacuum(r):
    c = np.zeros(N, dtype=complex)
    for m in range(N//2):
        c[2*m] = tanh(r)**m * sqrt(factorial(2*m)) / (2**m * factorial(m) * sqrt(cosh(r)))
    return c / np.linalg.norm(c)

def bs_unitary():
    return np.array([[1, 1], [1, -1]], dtype=complex) / sqrt(2)

def two_mode_op(U, op_a, op_b):
    # 把 2x2 模变换 U 作用到 (op_a, op_b) 上，返回两模算符
    return np.kron(U[0,0]*op_a + U[0,1]*op_b, I) + np.kron(I, U[1,0]*op_a + U[1,1]*op_b)

# 两模初态：a 模为输入态，b 模为真空
psi_in = np.kron(coherent(1.2), np.array([1]+[0]*(N-1), dtype=complex))
U = bs_unitary()
# 分束器演化：输出 = U† 作用于输入（等价于对模做 U 变换）
psi_out = two_mode_op(U.conj().T, a, a).conj().T @ psi_in  # 见下方说明

import numpy as np
from math import factorial, sqrt, tanh, cosh

N = 40
a = np.diag(np.sqrt(np.arange(1, N)), 1).astype(complex)
I = np.eye(N, dtype=complex)

def coherent(alpha):
    n = np.arange(N)
    c = np.exp(-abs(alpha)**2/2) * alpha**n / np.sqrt([factorial(k) for k in n])
    return c / np.linalg.norm(c)

def squeezed_vacuum(r):
    c = np.zeros(N, dtype=complex)
    for m in range(N//2):
        c[2*m] = tanh(r)**m * sqrt(factorial(2*m)) / (2**m * factorial(m) * sqrt(cosh(r)))
    return c / np.linalg.norm(c)

def bs_matrix():
    return np.array([[1, 1], [1, -1]], dtype=complex) / sqrt(2)

def evolve_two_mode(psi_in, U):
    # U 为 2x2 模变换矩阵；两模幺正为 kron(U, U)
    U2 = np.kron(U, U)
    return U2 @ psi_in

def photon_distribution(psi, N):
    # 把两模态 reshape 成 (N, N)，分别求各模的边缘分布
    M = psi.reshape(N, N)
    p_a = np.sum(np.abs(M)**2, axis=1)   # 模 a 的光子数分布
    p_b = np.sum(np.abs(M)**2, axis=0)   # 模 b 的光子数分布
    return p_a, p_b

def quadrature_variance(psi, N, mode='a', phase=0.0):
    # 计算单模正交分量 X_θ = (a e^{-iθ} + a† e^{iθ})/√2 的方差
    M = psi.reshape(N, N)
    if mode == 'a':
        rho = M @ M.conj().T
    else:
        rho = M.T @ M.conj()
    X = (a * np.exp(-1j*phase) + a.conj().T * np.exp(1j*phase)) / sqrt(2)
    mean = np.trace(rho @ X).real
    var = np.trace(rho @ X @ X).real - mean**2
    return var

# 主流程
alpha, r = 1.2, 0.8
psi_a = coherent(alpha)
psi_b = np.array([1]+[0]*(N-1), dtype=complex)
psi_in = np.kron(psi_a, psi_b)

U = bs_matrix()
psi_out = evolve_two_mode(psi_in, U)

p_c, p_d = photon_distribution(psi_out, N)
var_X_c = quadrature_variance(psi_out, N, 'a', 0.0)
var_P_c = quadrature_variance(psi_out, N, 'a', np.pi/2)
var_X_d = quadrature_variance(psi_out, N, 'b', 0.0)
var_P_d = quadrature_variance(psi_out, N, 'b', np.pi/2)

print("输出模 c 平均光子数:", np.sum(np.arange(N)*p_c))
print("输出模 d 平均光子数:", np.sum(np.arange(N)*p_d))
print("c 模 Var(X):", var_X_c, " Var(P):", var_P_c)
print("d 模 Var(X):", var_X_d, " Var(P):", var_P_d)

结果判读：干涉如何改变噪声

把相干态（α=1.2）与真空态送入 50:50 分束器时，输出两模仍为相干态，平均光子数各约 |α|²/2 ≈ 0.72，光子数分布为泊松型，正交分量方差保持真空水平 0.5，即分束器对相干态不产生额外噪声，也不压缩。

把压缩真空态（r=0.8）与真空态送入分束器时，输出两模是纠缠的压缩态（双模压缩真空），单模边缘分布不再是纯压缩态：Var(X_c) = Var(X_d) = (e^(−2r) + e^(2r))/4，代入 r=0.8 得约 0.44，比真空 0.5 略低；而 Var(P_c) = Var(P_d) = (e^(−2r) + e^(2r))/4 同样约 0.44。也就是说，单模正交分量方差被“平均化”，压缩被两模间的关联“藏”了起来——这正是后续马赫-曾德尔干涉与平衡零差探测能利用的关联资源。

验证要点：① 输出总光子数守恒，⟨n_c⟩+⟨n_d⟩ = ⟨n_a⟩+⟨n_b⟩；② 相干态输入时输出方差应回到 0.5（截断 N=40 下误差 < 1e-3）；③ 压缩态输入时两模交叉关联 ⟨X_c X_d⟩ 不为零，单模方差下降但两模联合方差可低于真空极限。

参数扫描与验证脚本

为了确认“压缩被关联隐藏”这一结论，可对 r 做扫描，同时计算单模方差与两模联合正交分量 X_c − X_d 的方差。后者在理想情况下应等于 e^(−2r)/2，即压缩度被完整保留在差分信号中。

def joint_variance(psi, N, phase=0.0):
    M = psi.reshape(N, N)
    X = (a * np.exp(-1j*phase) + a.conj().T * np.exp(1j*phase)) / sqrt(2)
    Xc = np.kron(X, I)
    Xd = np.kron(I, X)
    Xdiff = Xc - Xd
    mean = np.vdot(psi, Xdiff @ psi).real
    var = np.vdot(psi, Xdiff @ Xdiff @ psi).real - mean**2
    return var

for r in [0.0, 0.4, 0.8, 1.2]:
    psi_in = np.kron(squeezed_vacuum(r), np.array([1]+[0]*(N-1), dtype=complex))
    psi_out = evolve_two_mode(psi_in, bs_matrix())
    v_single = quadrature_variance(psi_out, N, 'a', 0.0)
    v_joint  = joint_variance(psi_out, N, 0.0)
    print(f"r={r:.1f}  单模Var(X)={v_single:.4f}  差分Var(Xc−Xd)={v_joint:.4f}  理论e^(−2r)/2={np.exp(-2*r)/2:.4f}")

预期输出（示例，实际数值随截断 N 略有偏差）：r=0 时单模与差分方差均为 0.5；r=0.8 时单模约 0.44、差分约 0.10；r=1.2 时单模约 0.55（已超过真空，单模出现反压缩）、差分约 0.045。这说明分束器把压缩态的单模压缩“转移”为两模关联，单模方差不再直接反映压缩度，必须用差分或联合测量才能恢复。