3. 自旋测量与关联函数
在 Bell 不等式实验中，真正被“测量”的不是一个预先存在的经典数值，而是自旋沿某一空间方向的投影结果。要编程实现这一过程，必须先把物理语言翻译成线性代数语言：测量方向是一个单位向量，测量基是该方向上的本征态，测量算符是投影算符，而关联函数则是两个投影结果在大量重复测量上的统计平均。本节围绕这条线索，给出可落地的算符构造与关联函数计算方法。
一、自旋投影算符与测量基
对自旋 1/2 粒子，沿单位方向 n = (nx, ny, nz) 的自旋算符可写成
σ·n = nx·σx + ny·σy + nz·σz其中 σx、σy、σz 为泡利矩阵。该算符有两个本征值 +1 和 −1，对应本征态 |n,+⟩ 与 |n,−⟩。测量“自旋沿 n 方向是 +1 还是 −1”，在数学上就是向这两个本征态投影。定义投影算符
P₊(n) = (I + σ·n) / 2
P₋(n) = (I − σ·n) / 2它们满足 P₊ + P₋ = I、P₊P₋ = 0、P₊² = P₊，正是正交投影算符的标准性质。于是“沿 n 方向测量自旋”这一可观测量可以写成
A(n) = (+1)·P₊(n) + (−1)·P₋(n) = σ·n这一步非常关键：测量算符就是 σ·n 本身，而投影算符给出了取 +1 或 −1 的概率分解。对任意态 |ψ⟩，取 +1 的概率为 ⟨ψ|P₊(n)|ψ⟩，取 −1 的概率为 ⟨ψ|P₋(n)|ψ⟩，两者之和为 1。
需要先约定参照系：本文取 z 轴为量子化轴，|0⟩ 对应 σz 的 +1 本征态，|1⟩ 对应 −1 本征态；方位角 φ 从 x 轴正向起算，极角 θ 从 z 轴正向起算。全文与图示均采用这一约定。于是
n = (sinθ·cosφ, sinθ·sinφ, cosθ)当 θ=0 时 n 指向 +z，测量基退化为计算基 {|0⟩, |1⟩}；当 θ=π/2、φ=0 时 n 指向 +x，测量基为 |±x⟩ = (|0⟩ ± |1⟩)/√2。编程时只需保存 (θ, φ) 两个参数，即可唯一确定一组测量基。
二、关联函数 E(a,b) 的定义
设 Alice 沿方向 a 测量粒子 1，Bob 沿方向 b 测量粒子 2，各自得到 ±1 的结果。定义单次联合结果的乘积为 A·B ∈ {+1, −1}，则关联函数为
E(a,b) = ⟨ψ| (σ·a) ⊗ (σ·b) |ψ⟩等价地，若用概率表示，设 P(i,j) 为 Alice 得 i、Bob 得 j 的联合概率（i, j ∈ {+1, −1}），则
E(a,b) = P(+,+) + P(−,−) − P(+,−) − P(−,+)两种写法在理想测量下完全等价：前者是算符期望值，后者是统计频率的极限。编程实现时，前者适合直接做矩阵运算，后者适合用采样数据估计，二者可以互相校验。对 Bell 态 |Φ⁺⟩ = (|00⟩ + |11⟩)/√2，解析结果为 E(a,b) = a·b = cosθab，其中 θab 是两测量方向夹角。这一结果将作为程序正确性的基准。
三、编程实现：任意方向测量算符
下面用 NumPy 给出最小可运行实现。核心只有三步：由 (θ, φ) 生成单位向量，用泡利矩阵线性组合出 σ·n，再与态做张量积并求迹。示例参数与结果均为教学演示，非实测数据。
import numpy as np
sx = np.array([[0, 1], [1, 0]], dtype=complex)
sy = np.array([[0, -1j], [1j, 0]], dtype=complex)
sz = np.array([[1, 0], [0, -1]], dtype=complex)
I2 = np.eye(2, dtype=complex)
def direction(theta, phi):
return np.array([np.sin(theta)*np.cos(phi),
np.sin(theta)*np.sin(phi),
np.cos(theta)])
def sigma_n(theta, phi):
nx, ny, nz = direction(theta, phi)
return nx*sx + ny*sy + nz*sz
def projectors(theta, phi):
sn = sigma_n(theta, phi)
return (I2 + sn)/2, (I2 - sn)/2
def E(rho, a, b):
A = sigma_n(*a)
B = sigma_n(*b)
return np.real(np.trace(rho @ np.kron(A, B)))
# Bell 态 |Φ⁺⟩
phi_plus = np.array([1, 0, 0, 1], dtype=complex) / np.sqrt(2)
rho = np.outer(phi_plus, phi_plus.conj())
a = (np.pi/2, 0.0) # 沿 +x
b = (np.pi/2, np.pi/2) # 沿 +y
print(E(rho, a, b)) # 演示输出约为 0.0
print(E(rho, a, a)) # 演示输出约为 1.0运行后可观察到：a 与 b 正交时 E≈0，同向时 E≈1，与解析式 a·b 一致。若把 b 改为 (π/2, π/4)，即两方向夹角 45°，则 E≈cos45°≈0.7071。这些数值可作为程序正确性的快速自检点。注意 np.trace 可能返回带微小虚部的复数，取实部是合理做法，因为厄米算符的期望值必为实数。
四、用采样估计关联函数
矩阵法给出的是理想期望值。若想模拟真实实验，需要按概率采样。对每个方向，先算四个联合概率 P(i,j)，再用多项式采样生成 N 次结果，最后用频率估计 E。示例代码如下：
def joint_probs(rho, a, b):
Pa_p, Pa_m = projectors(*a)
Pb_p, Pb_m = projectors(*b)
ops = {(1,1): np.kron(Pa_p, Pb_p), (1,-1): np.kron(Pa_p, Pb_m),
(-1,1): np.kron(Pa_m, Pb_p), (-1,-1): np.kron(Pa_m, Pb_m)}
return {k: np.real(np.trace(rho @ v)) for k, v in ops.items()}
def sample_E(rho, a, b, N=20000, seed=0):
rng = np.random.default_rng(seed)
probs = joint_probs(rho, a, b)
keys = list(probs.keys())
p = np.array([probs[k] for k in keys])
p = p / p.sum()
idx = rng.choice(len(keys), size=N, p=p)
vals = np.array([keys[i][0]*keys[i][1] for i in idx])
return vals.mean()以 N=20000、a=(π/2,0)、b=(π/2,π/4) 为例，采样估计通常落在 0.70 附近，与 0.7071 的偏差量级约为 1/√N≈0.007，属于正常统计涨落。判断实现是否正确，应看偏差是否随 N 增大而收缩，而不是要求单次结果与解析值完全相等。若偏差长期停留在 0.1 以上，应优先检查方向向量归一化、张量积顺序（Alice 在前、Bob 在后）以及投影算符是否用错本征态。