4. 纯态与混态的判别与计算

在量子力学的计算实践中，拿到一个密度矩阵 ρ 之后，最先要回答的问题往往是：它描述的是一个纯态，还是一个混态？这个判断不是形式上的分类游戏，而是直接决定后续计算路径的分岔点——纯态可以用一个态矢量 |ψ⟩ 完整描述，混态则必须保留密度矩阵形式，且任何试图用单个态矢量替代它的做法都会丢失信息。本节围绕三个判别工具展开：纯度、冯·诺依曼熵和布洛赫球表示，并给出可直接运行的代码与可视化方法。

三个判别量的定义与判据

纯度定义为 γ = Tr(ρ²)。对于纯态，ρ² = ρ，因此 γ = 1；对于混态，γ < 1，且 γ 的最小值为 1/d，其中 d 是希尔伯特空间维数。单量子比特情形下 d = 2，纯度下限为 0.5，对应完全混态 ρ = I/2。纯度的优点是计算极简，只需一次矩阵平方和求迹，不涉及对数运算，数值上非常稳定。

冯·诺依曼熵定义为 S(ρ) = −Tr(ρ ln ρ)。若 ρ 的本征值为 {λᵢ}，则 S = −Σᵢ λᵢ ln λᵢ，约定 0 ln 0 = 0。纯态只有一个本征值为 1、其余为 0，故 S = 0；混态 S > 0，单量子比特的最大熵为 ln 2 ≈ 0.6931，对应完全混态。熵的物理含义是态的不确定度或混合程度，它比纯度更细致地刻画了本征值分布，但涉及对数，当本征值极小时需要做截断处理以避免数值噪声。

布洛赫球表示则提供了几何判据。任意单量子比特密度矩阵可写成 ρ = (I + r⃗·σ⃗)/2，其中 r⃗ = (rₓ, r_y, r_z) 是布洛赫矢量，σ⃗ 是泡利矩阵向量。纯态对应 |r⃗| = 1，即布洛赫球表面上的点；混态对应 |r⃗| < 1，位于球内部；完全混态对应 r⃗ = 0，即球心。三个判据在数学上等价，但各自适合不同场景：纯度适合批量快速筛查，熵适合信息论视角的分析，布洛赫球适合单比特的几何直觉与可视化。

用代码计算纯度与熵

下面用 Python 与 NumPy 实现三个判别量。示例选取三个密度矩阵：纯态 |0⟩、叠加态 |+⟩ 以及一个退相干后的混态 ρ = 0.75|0⟩⟨0| + 0.25|1⟩⟨1|。第三个矩阵的本征值为 0.75 和 0.25，预期纯度 γ = 0.75² + 0.25² = 0.625，熵 S = −0.75 ln 0.75 − 0.25 ln 0.25 ≈ 0.5623。以下数值均为按公式计算得到的示例结果，用于对照验证。

import numpy as np

def purity(rho):
    return np.real(np.trace(rho @ rho))

def von_neumann_entropy(rho):
    vals = np.linalg.eigvalsh(rho)
    vals = vals[vals > 1e-12]          # 截断数值噪声
    return -np.sum(vals * np.log(vals))

def bloch_vector(rho):
    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)
    return np.array([np.real(np.trace(rho @ s)) for s in (sx, sy, sz)])

ket0 = np.array([[1], [0]], dtype=complex)
ket_plus = np.array([[1], [1]], dtype=complex) / np.sqrt(2)

rho_pure = ket0 @ ket0.conj().T
rho_plus = ket_plus @ ket_plus.conj().T
rho_mix = np.array([[0.75, 0], [0, 0.25]], dtype=complex)

for name, rho in [("|0>", rho_pure), ("|+>", rho_plus), ("混态", rho_mix)]:
    r = bloch_vector(rho)
    print(f"{name:6s} 纯度={purity(rho):.4f} "
          f"熵={von_neumann_entropy(rho):.4f} "
          f"|r|={np.linalg.norm(r):.4f}")

运行后应得到类似输出：|0⟩ 与 |+⟩ 的纯度均为 1.0000、熵为 0.0000、|r⃗| 为 1.0000；混态纯度为 0.6250、熵为 0.5623、|r⃗| = |0.75 − 0.25| = 0.5000。三个判据给出的结论完全一致。这里有一个容易踩的坑：np.linalg.eigvalsh 要求输入是厄米矩阵，若密度矩阵因浮点误差出现极小的负本征值（如 −1e−17），直接取对数会产生 nan，因此代码中用 vals > 1e-12 做截断。另一个细节是 np.trace 返回复数类型，纯度取实部即可。

布洛赫球可视化

把布洛赫矢量画到球面上，能直观看出态的位置。用 Matplotlib 的 3D 坐标绘制单位球面，再把各态的 r⃗ 作为散点标出，并用从原点出发的箭头表示矢量方向。对于上面的混态，矢量沿 z 轴指向球内 0.5 处；对于 |+⟩，矢量沿 x 轴指向球面。可视化时建议把球面透明度调低（alpha 约 0.15），并画出三条坐标轴，避免散点被球面遮挡而误判位置。

import matplotlib.pyplot as plt

fig = plt.figure(figsize=(6, 6))
ax = fig.add_subplot(111, projection="3d")

u, v = np.mgrid[0:2*np.pi:40j, 0:np.pi:20j]
xs, ys, zs = np.cos(u)*np.sin(v), np.sin(u)*np.sin(v), np.cos(v)
ax.plot_surface(xs, ys, zs, color="#a9cdec", alpha=0.15, linewidth=0)

for name, rho, c in [("|0>", rho_pure, "red"),
                     ("|+>", rho_plus, "darkorange"),
                     ("混态", rho_mix, "purple")]:
    r = bloch_vector(rho)
    ax.quiver(0, 0, 0, r[0], r[1], r[2], color=c, arrow_length_ratio=0.12)
    ax.text(r[0], r[1], r[2], name, color=c)

ax.set_xlim(-1, 1); ax.set_ylim(-1, 1); ax.set_zlim(-1, 1)
ax.set_xlabel("x"); ax.set_ylabel("y"); ax.set_zlabel("z")
plt.show()

判别流程与常见误区

实际计算中，建议按“先纯度、后熵、再几何”的顺序处理：纯度只需矩阵乘法，适合对大量密度矩阵做初筛；当纯度接近 1 但又不精确等于 1 时，再计算熵确认是否为数值误差；对单量子比特，最后用布洛赫矢量模长做交叉验证，三者不一致通常意味着矩阵本身不合法（如迹不为 1 或非厄米）。需要特别注意的是，纯度小于 1 并不一定意味着“制备失败”，它可能来自与环境的纠缠导致的约化密度矩阵，此时子系统本身确实是混态，这是物理事实而非计算错误。另外，多量子比特情形下布洛赫球不再适用，应回到纯度与熵两个判据，且纯度的下限变为 1/d，判据阈值需随维数调整。