2. 量子态与密度矩阵工程表示

在量子信息工程中，态矢 |ψ⟩ 只能描述“我们完全知道这个态”的情形。一旦系统与外界发生不可控耦合、制备过程存在经典随机性，或者我们只掌握某个子系统的信息，态矢就不够用了。密度矩阵 ρ 是更通用的表示：它既能描述纯态，也能描述混态，还能通过偏迹自然给出子系统的状态。本节围绕“如何用 ρ 做工程计算”展开，重点落在可编程实现上。

纯态、混态与密度矩阵的定义

纯态对应一个确定的态矢，其密度矩阵为 ρ = |ψ⟩⟨ψ|。混态则是若干纯态的经典概率混合：ρ = Σᵢ pᵢ |ψᵢ⟩⟨ψᵢ|，其中 pᵢ ≥ 0 且 Σpᵢ = 1。这里的“经典概率”是关键：混合来自制备信息缺失，而不是量子叠加。判断纯态与混态最直接的工程指标是纯度 Tr(ρ²)：纯态为 1，混态小于 1。另一个常用指标是冯·诺依曼熵 S(ρ) = −Tr(ρ log₂ρ)，纯态为 0，混态大于 0。

密度矩阵必须满足三条性质：厄米性 ρ = ρ†、半正定性 ρ ≥ 0（所有本征值非负）、迹为 1 Tr(ρ) = 1。工程实现中，这三条也是数值校验的入口：任何计算得到的 ρ 都应先检查厄米性误差、最小本征值和迹偏差，再进入后续熵计算。

约化密度矩阵：只看子系统

复合系统 AB 处于 ρ_AB 时，子系统 A 的状态由偏迹给出：ρ_A = Tr_B(ρ_AB)。偏迹的工程含义是“对 B 的自由度求和，丢弃 B 的信息”。以两比特 Bell 态 |Φ⁺⟩ = (|00⟩ + |11⟩)/√2 为例，其整体是纯态，但 ρ_A = I/2，是最大混态，熵为 1 bit。这正是纠缠的体现：整体信息完整，局部信息却完全丢失。约化密度矩阵是计算纠缠熵、互信息、条件熵的基础操作。

谱分解：把 ρ 变成可计算的形式

由于 ρ 是厄米半正定矩阵，它可以做谱分解：ρ = Σₖ λₖ |k⟩⟨k|，其中 λₖ 是本征值，|k⟩ 是对应本征矢。谱分解把“矩阵求对数”转化为“本征值求对数”：S(ρ) = −Σₖ λₖ log₂ λₖ。工程上，这意味着熵计算的核心步骤是：构造 ρ → 校验性质 → 求本征值 → 过滤零本征值 → 对正本征值求 −λ log₂ λ 并求和。注意 λ = 0 时 λ log λ 的极限为 0，数值实现中必须显式跳过，否则会出现 NaN。

工程实现：从 ρ 到熵的完整代码路径

下面给出一段可直接运行的 Python 示例，使用 NumPy 完成密度矩阵构造、性质校验、谱分解与熵计算。示例中的数值均为教学演示，用于说明计算流程与判读方式。

import numpy as np

def build_pure(psi):
    psi = np.asarray(psi, dtype=complex)
    psi = psi / np.linalg.norm(psi)
    return np.outer(psi, psi.conj())

def build_mixed(states, probs):
    rho = np.zeros((len(states[0]), len(states[0])), dtype=complex)
    for p, s in zip(probs, states):
        rho += p * build_pure(s)
    return rho

def check_rho(rho, tol=1e-9):
    herm = np.allclose(rho, rho.conj().T, atol=tol)
    evals = np.linalg.eigvalsh(rho)
    trace_ok = np.isclose(np.trace(rho).real, 1.0, atol=tol)
    psd = np.all(evals >= -tol)
    return herm, psd, trace_ok, evals

def vn_entropy(rho):
    evals = np.linalg.eigvalsh(rho)
    evals = evals[evals > 1e-12]          # 过滤零本征值
    return -np.sum(evals * np.log2(evals))

# 示例 1：纯态 |+⟩ = (|0⟩ + |1⟩)/√2
rho_pure = build_pure([1, 1])
print("纯态熵:", vn_entropy(rho_pure))   # 期望 0

# 示例 2：混态 ρ = 0.7|0⟩⟨0| + 0.3|1⟩⟨1|
rho_mix = build_mixed([[1, 0], [0, 1]], [0.7, 0.3])
print("混态熵:", vn_entropy(rho_mix))    # 期望约 0.881 bit

# 示例 3：Bell 态约化密度矩阵
bell = np.array([1, 0, 0, 1]) / np.sqrt(2)
rho_ab = build_pure(bell).reshape(2, 2, 2, 2)
rho_a = np.trace(rho_ab, axis1=1, axis2=3)   # 对 B 求偏迹
print("约化态熵:", vn_entropy(rho_a))    # 期望 1.0 bit

运行上述代码时，check_rho 返回的四个量应分别满足：厄米为 True、半正定为 True、迹为 True、本征值列表非负。若半正定校验失败，通常意味着构造过程中出现了负概率或数值舍入误差，需要检查输入概率是否归一化。若熵计算返回 NaN，几乎总是因为未过滤零本征值。

偏迹的工程实现要点

偏迹不是简单的“删掉 B 的行列”，而是对 B 的指标做求和。以两比特系统为例，ρ_AB 是 4×4 矩阵，按 (i_A, i_B) 重排为 2×2×2×2 张量后，对 B 的两个指标求迹，得到 2×2 的 ρ_A。NumPy 中可用 np.trace(rho_ab, axis1=1, axis2=3) 完成。工程上要特别注意指标顺序：如果系统由多个量子比特组成，偏迹前必须明确“哪个轴对应哪个子系统”，否则会得到错误的约化态。建议在代码中固定约定：张量轴按 (A, B, C, ...) 顺序排列，偏迹时显式写出 axis1/axis2。

数值稳定性与常见陷阱

• 零本征值处理：谱分解后本征值可能为 1e-16 量级的数值噪声，应设置阈值（如 1e-12）过滤，避免 log₂(0)。
• 厄米性校验：浮点运算可能使 ρ 轻微非厄米，建议用 np.allclose(rho, rho.conj().T) 而非严格相等。
• 概率归一化：构造混态前应检查 Σpᵢ = 1，否则 Tr(ρ) ≠ 1，后续熵值无物理意义。
• 本征值排序：np.linalg.eigvalsh 返回升序本征值，熵计算与顺序无关，但调试时按降序观察更直观。

从密度矩阵到纠缠度量的衔接

密度矩阵是后续所有纠缠度量的计算底座。约化密度矩阵的谱分解直接给出纠缠熵；两个子系统的联合谱分解可以计算互信息；部分转置后的本征值则用于判断 PPT 判据和计算负度。因此，本节的重点不是记住公式，而是建立一条稳定的工程路径：任何量子态先转成 ρ，校验三条性质，再做谱分解，最后按目标度量选择对应的本征值组合。只要这条路径在代码中固定下来，后续的熵计算与纠缠度量就只是替换求和表达式的问题。