2. 量子态与密度矩阵的仿真表示

在量子噪声信道仿真中，我们最终要处理的对象几乎都不是“一个漂亮的态矢量”，而是密度矩阵。原因很直接：噪声会把纯态变成混合态，而混合态无法用单个 ket 描述；此外，多体系统中的部分迹、纠缠、退相干，也都必须借助密度矩阵才能自然表达。因此，本节的任务不是重复线性代数定义，而是把“量子态”落到可计算的矩阵对象上：如何用数值工具构造纯态、混合态和多体态，并验证迹为 1、半正定这些基本性质。

先明确约定。本节所有密度矩阵均采用计算基 \(|0\rangle, |1\rangle, \dots\) 的列向量表示，即 \(|0\rangle=(1,0)^T\)、\(|1\rangle=(0,1)^T\)；密度矩阵定义为 \(\rho=\sum_i p_i|\psi_i\rangle\langle\psi_i|\)，满足 \(\mathrm{Tr}(\rho)=1\)、\(\rho=\rho^\dagger\)、\(\rho\succeq 0\)。后文提到的“迹”“半正定”“本征值”都基于这一约定，图示中的矩阵元位置也与此一致。

纯态：从态矢量到秩一密度矩阵

单量子比特纯态可写成 \(|\psi\rangle=\cos(\theta/2)|0\rangle+e^{i\varphi}\sin(\theta/2)|1\rangle\)。在仿真中，我们通常先构造态矢量，再做外积。以 \(\theta=\pi/2,\ \varphi=0\) 为例，得到 \(|+\rangle=(|0\rangle+|1\rangle)/\sqrt{2}\)，其密度矩阵为

import numpy as np

# 计算基
ket0 = np.array([[1.0], [0.0]], dtype=complex)
ket1 = np.array([[0.0], [1.0]], dtype=complex)

# 构造 |+> = (|0> + |1>)/sqrt(2)
plus = (ket0 + ket1) / np.sqrt(2)
rho_plus = plus @ plus.conj().T

print(rho_plus)
# [[0.5+0.j 0.5+0.j]
#  [0.5+0.j 0.5+0.j]]
print(np.trace(rho_plus))          # (1+0j)
print(np.linalg.eigvalsh(rho_plus))# [0. 1.]
print(np.allclose(rho_plus, rho_plus.conj().T))  # True

这里有两个值得注意的观察点。第一，纯态的密度矩阵是秩一的，本征值只有一个 1 和其余 0，因此“纯态”在数值上等价于 \(\mathrm{Tr}(\rho^2)=1\)。第二，非对角元 \(0.5\) 携带的是相位相干信息，它正是后续噪声信道要削弱的对象。若把 \(\varphi\) 改为 \(\pi/2\)，对角元不变而非对角元变为 \(\pm 0.5i\)，这说明仅看布居数（对角元）无法区分 \(|+\rangle\) 与 \(|+i\rangle\)，必须保留复矩阵表示。

混合态：概率加权与“最大混合”基准

混合态最典型的构造方式是概率凸组合。例如以 \(p=0.7\) 的概率制备 \(|0\rangle\)、以 \(0.3\) 的概率制备 \(|1\rangle\)，得到 \(\rho=0.7|0\rangle\langle0|+0.3|1\rangle\langle1|=\mathrm{diag}(0.7,0.3)\)。它的迹为 1、本征值为 0.7 和 0.3，均非负，因此半正定；但 \(\mathrm{Tr}(\rho^2)=0.58<1\)，说明它是混合态。另一个极端是最大混合态 \(I/2=\mathrm{diag}(0.5,0.5)\)，本征值全为 0.5，纯度 \(\mathrm{Tr}(\rho^2)=0.5\)，可作为衡量退相干程度的参照基准。

实际仿真中，混合态往往不是手工写对角元，而是由一组带权纯态生成。下面这段代码同时构造“经典混合”和“相干叠加再混合”两种情形，用于对比：

def density_from_ensemble(states, probs):
    rho = np.zeros((states[0].shape[0], states[0].shape[0]), dtype=complex)
    for psi, p in zip(states, probs):
        rho += p * (psi @ psi.conj().T)
    return rho

# 情形A：经典混合 diag(0.7, 0.3)
rhoA = density_from_ensemble([ket0, ket1], [0.7, 0.3])

# 情形B：|+> 与 |-> 各半，注意非对角元相互抵消
minus = (ket0 - ket1) / np.sqrt(2)
rhoB = density_from_ensemble([plus, minus], [0.5, 0.5])
# rhoB ≈ diag(0.5, 0.5)，非对角元为 0

情形 B 很能说明问题：两个纯态叠加概率各半，结果却是完全无相干的最大混合态。这正是“混合”与“叠加”的本质区别——叠加发生在概率幅层面，混合发生在概率层面。噪声信道中的退相位过程，本质上就是把非对角元按某个因子衰减，使 \(\rho\) 从情形 A 的相干形态逐步滑向情形 B 的对角形态。

多体量子态：张量积、部分迹与纠缠判据

两比特系统的态空间是 4 维，基矢按 \(|00\rangle, |01\rangle, |10\rangle, |11\rangle\) 排列。构造多体态的关键操作是张量积：单比特算符 \(A\)、\(B\) 作用在各自子系统上时，整体算符为 \(A\otimes B\)。以 Bell 态 \(|\Phi^+\rangle=(|00\rangle+|11\rangle)/\sqrt{2}\) 为例，其密度矩阵在 4×4 矩阵中只有四个非零元：\(\rho_{00,00}=\rho_{00,11}=\rho_{11,00}=\rho_{11,11}=0.5\)。

def kron(*ops):
    out = ops[0]
    for op in ops[1:]:
        out = np.kron(out, op)
    return out

# Bell 态 |Phi+>
phi = (kron(ket0, ket0) + kron(ket1, ket1)) / np.sqrt(2)
rho_phi = phi @ phi.conj().T

# 对第二比特做部分迹，得到第一比特的约化密度矩阵
def partial_trace_b(rho, dim=2):
    r = rho.reshape(dim, dim, dim, dim)
    return np.trace(r, axis1=1, axis2=3)

rho1 = partial_trace_b(rho_phi)
print(rho1)                 # [[0.5, 0], [0, 0.5]]
print(np.linalg.eigvalsh(rho1))  # [0.5 0.5]

约化密度矩阵是最大混合态，说明单看任一比特都得不到确定信息，这正是纠缠的典型特征。工程上常用的快速判据是纯度：若整体为纯态而子系统约化态纯度小于 1，则两体态纠缠。对 Bell 态，\(\mathrm{Tr}(\rho_1^2)=0.5<1\)，判据成立。需要提醒的是，该判据只对“整体纯、子系统混合”的两体情形直接有效，对整体本身为混合态的情况不能简单套用，此时应改用部分转置本征值（PPT 判据）等更一般的方法。

性质验证：把“检查”写成可复用的函数

构造完密度矩阵后，必须做三项基本验证，否则后续信道仿真会建立在错误对象上。第一是迹为 1，容差建议取 \(10^{-10}\) 量级；第二是厄米性，即 \(\rho=\rho^\dagger\)；第三是半正定，最稳妥的做法是求本征值并检查最小值是否不小于 \(-\varepsilon\)，而不是只看对角元非负——对角元非负只是必要条件，存在对角元全正但本征值为负的反例。下面给出一个可直接复用的校验函数：

def validate_density_matrix(rho, tol=1e-10):
    checks = {}
    checks["shape"] = rho.shape[0] == rho.shape[1]
    checks["trace"] = abs(np.trace(rho) - 1.0) < tol
    checks["hermitian"] = np.allclose(rho, rho.conj().T, atol=tol)
    eigs = np.linalg.eigvalsh((rho + rho.conj().T) / 2)
    checks["positive_semidefinite"] = eigs.min() >= -tol
    checks["purity"] = float(np.real(np.trace(rho @ rho)))
    return checks, eigs

for name, r in [("|+>", rho_plus), ("diag(0.7,0.3)", rhoA),
                ("I/2", np.eye(2)/2), ("Bell", rho_phi)]:
    checks, eigs = validate_density_matrix(r)
    print(name, checks, "min_eig=", eigs.min())

预期输出中，四个对象的 trace、hermitian、positive_semidefinite 均为 True；纯度分别为 1.0、0.58、0.5、1.0。若某个构造结果出现 trace 偏离 1，通常源于归一化遗漏；若出现最小本征值为负且超出容差，则多半是概率权重为负或态矢量未归一化。把这三项检查固化为函数并在每次构造后调用，是避免“错误态悄悄进入信道”的最低成本手段。