开放量子系统与退相干建模：密度矩阵、Kraus 与 Lindblad 方程

引言

教科书里的量子力学是「封闭系统」的量子力学：一个理想化的态在酉演化下完美地旋转，永不丢失信息。真实的量子硬件不是这样。一个超导量子比特浸泡在充满热光子的电磁环境里，一个离子阱里的离子被激光相位噪声和电极电压涨落包围，一个光子比特在光纤里不断遭遇散射——没有任何量子系统是孤立的。

这种「与环境耦合」带来两个后果。第一是退相干（decoherence）：量子态的相干性（叠加）随时间衰减，相位信息扩散到环境里，宏观上表现为 $T_1$、$T_2$ 时间有限。第二是耗散（dissipation）：布居会自发地从高能级掉到低能级，最终趋向热平衡。

理解开放量子系统不是学术兴趣，而是量子工程的日常：你写的噪声模型决定了误差缓解能做到什么程度、纠错码的阈值在哪里、量子比特能跑多深的电路。本文系统讲解：先区分封闭与开放系统，再拆密度矩阵、量子信道与 Kraus 表示、Lindblad 主方程，然后讲三类标准噪声的建模、数值方法（Liouvillian、量子轨迹）与 QuTiP 实战，最后落到硬件上的表现与缓解。目标：能写出一段可信的噪声仿真代码，并理解每个参数对应硬件的哪个物理量。

前置：量子比特校准与噪声 、误差缓解 。

封闭系统与开放系统

封闭系统：酉演化

封闭系统的态是纯态 $|\psi\rangle$，演化由薛定谔方程给出：

iħ d|ψ>/dt = H|ψ>
解：|ψ(t)> = U(t)|ψ(0)>，U(t) = exp(-iHt/ħ)
关键性质：U 是酉的（U†U = I），所以内积守恒、纯度守恒、可逆。

酉演化有一个致命限制：它无法描述「信息丢失」。如果量子比特的相干性真的衰减了，那这个衰减过程不可能由某个酉变换描述——因为酉变换可逆，而退相干不可逆（信息跑进了你无法访问的环境自由度）。

开放系统：必须用密度矩阵

要描述开放系统，有两种等价视角：

1. 把环境也纳入：系统 + 环境构成大封闭系统，整体酉演化，然后**部分迹（partial trace）**掉环境自由度，得到系统的约化密度矩阵 $\rho_S = \mathrm{Tr}_E(|\Psi\rangle\langle\Psi|)$。
2. 只描述系统：用一个量子信道（quantum channel）——一个把输入密度矩阵映射到输出密度矩阵的完全正定保迹（CPTP）映射。

视角 1 是物理上的真相，视角 2 是工程上可用的抽象。两者的桥梁就是 Kraus 表示。

对比：
  封闭：|ψ>  →  U|ψ>
  开放：ρ   →  Σ_i K_i ρ K_i†      （Kraus 算符，Σ K_i†K_i = I）
  纯态 → 混合态（纯度 Tr(ρ²) < 1）

密度矩阵与部分迹

密度矩阵基础

密度矩阵 $\rho$ 是描述量子态最一般的工具：

性质：
  ρ = ρ†                （厄米）
  Tr(ρ) = 1             （归一）
  ρ >= 0                （半正定）
  纯态 ⇔ Tr(ρ²) = 1
  混合态 ⇔ Tr(ρ²) < 1
  期望值：<A> = Tr(ρ A)
  单比特 ρ = (I + x σx + y σy + z σz) / 2   （Bloch 向量表示）

Bloch 球是理解退相干最直观的图像：纯态在球面上，混合态在球内，完全退相干到球心（最大混合态 $I/2$）。

部分迹

部分迹的规则很直接：把环境自由度「求和掉」。

import numpy as np

def partial_trace_b(rho, dim_a, dim_b, keep="A"):
    """两体系统求部分迹。rho 为 (dim_a*dim_b) 维方阵。"""
    r = rho.reshape(dim_a, dim_b, dim_a, dim_b)
    if keep == "A":
        return np.einsum("ibjb->ij", r)     # 对 B 求迹
    return np.einsum("aiaj->ij", r)         # 对 A 求迹

# 例：Bell 态的部分迹是最大混合态
bell = np.array([1, 0, 0, 1]) / np.sqrt(2)
rho = np.outer(bell, bell)
print(np.round(partial_trace_b(rho, 2, 2, keep="A"), 3))   # ≈ [[0.5,0],[0,0.5]]

这个例子是纠缠与退相干联系的核心：对纠缠态求部分迹，得到的是混合态。反过来说，你的量子比特之所以退相干，正是因为它与环境的纠缠在增长——相干性没有消失，只是「漏」到了环境的自由度里。

量子信道与 Kraus 表示

信道必须完全正定

不是随便一个保迹映射都是合法的量子信道。合法信道必须完全正定（completely positive）——不仅它本身把正定矩阵映到正定矩阵，而且与任意辅助系统张量积后仍保持正定。转置映射是最著名的反例：它保迹、正定，但不完全正定（对纠缠态会给出负本征值），所以不是合法信道。

Kraus 定理

Kraus 定理：任何 CPTP 映射都可以写成

ρ' = Σ_i K_i ρ K_i†        （Kraus 表示）
约束：Σ_i K_i† K_i = I     （保迹条件）

Kraus 算符的个数 $i$ 最多是 $d^2$（$d$ 为系统维度）。Kraus 表示不唯一——同一信道可以有不同的 Kraus 集合（对应环境的不同「基底选择」）。

三类标准信道的 Kraus 表示

1. 振幅阻尼（amplitude damping，T1 过程）
   K0 = [[1, 0], [0, sqrt(1-γ)]]
   K1 = [[0, sqrt(γ)], [0, 0]]
   物理：|1> 以概率 γ 衰减到 |0>（自发辐射）

2. 相位阻尼（phase damping，纯退相位）
   K0 = sqrt(1-λ) I
   K1 = sqrt(λ) [[1,0],[0,0]]
   K2 = sqrt(λ) [[0,0],[0,1]]
   物理：丢失相位信息，不改变布居

3. 去极化（depolarizing）
   K0 = sqrt(1-p) I
   K1 = sqrt(p/3) X,  K2 = sqrt(p/3) Y,  K3 = sqrt(p/3) Z
   物理：以概率 p 被随机 Pauli 错误击中

# QuTiP: 用 Kraus 算符构造并作用信道
import qutip as qt
import numpy as np

gamma = 0.05
K0 = qt.Qobj([[1, 0], [0, np.sqrt(1 - gamma)]])
K1 = qt.Qobj([[0, np.sqrt(gamma)], [0, 0]])
rho0 = qt.basis(2, 1).proj()          # |1><1|
rho_out = K0 * rho0 * K0.dag() + K1 * rho0 * K1.dag()
print(rho_out)                        # 布居从 |1> 部分转到 |0>

主方程：Lindblad 形式

Kraus 表示描述「一段时间后的结果」，但不给出连续时间的演化方程。要描述连续演化，用主方程（master equation）。在玻恩-马尔可夫近似（环境无记忆、耦合弱）下，得到 Lindblad 方程：

dρ/dt = -i[H, ρ]/ħ  +  Σ_k γ_k ( L_k ρ L_k†  -  1/2 {L_k† L_k, ρ} )

第一项：酉演化（相干）
第二项：耗散项（非相干），L_k 为跳跃算符，γ_k 为速率
{} 为反对易子

Lindblad 形式是最一般的马尔可夫量子主方程：任何满足 CPTP 且马尔可夫的连续演化都能写成这个形式（Lindblad 1976 / Gorini-Kossakowski-Sudarshan 1976）。

三类噪声的 Lindblad 算符

注意热噪声要同时加吸收（$\sigma_+$，从环境吸能）和发射（$\sigma_-$）两个算符，比值由热占据数 $n_{th}$ 决定。超导量子比特在 10 mK 下 $n_{th} \approx 0.01$，接近零温，所以只用 $\sigma_-$ 就够；但离子阱或室温系统必须考虑热项。

Bloch 方程与 T1、T2

把 Lindblad 方程投影到 Bloch 向量上，就得到Bloch 方程——这是连接抽象主方程与硬件可测量量的关键：

dx/dt = -x / T2
dy/dt = -y / T2
dz/dt = -(z - z_eq) / T1

其中：
  1/T2 = 1/(2 T1) + 1/T_φ

物理含义：
  T1：纵向弛豫时间（布居衰减）
  T2：横向弛豫时间（相干性衰减）
  T_φ：纯退相位时间（与 T1 无关的相位噪声）
  z_eq：热平衡时的布居差

$1/T_2 = 1/(2T_1) + 1/T_\phi$ 这个关系式极其重要，它告诉你：即使没有任何纯相位噪声（$T_\phi \to \infty$），相干性也会因为 T1 衰减而消失，上限是 $T_2 = 2T_1$。 这就是前面基准测试里提到的 $T_2 \le 2T_1$ 硬约束的来源。

常见情形判读：
  T2 ≈ 2*T1   → 纯 T1 限制，相位噪声很小
  T2 << 2*T1  → 相位噪声主导，值得做动态解耦
  T2 > 2*T1   → 数据可疑（或测量方式特殊）

T2* 与 T2 的区别

$T_2^$ 是自由感应衰减（FID）时间，包含所有低频噪声；$T_2$ 是Hahn 回波测的，用 $\pi$ 脉冲把低频噪声抵消掉。所以 $T_2^ \le T_2$。差距越大，说明低频（准静态）噪声越强。

噪声谱与时间尺度：
  1/f 噪声（低频）  → 主导 T2*，可用回波/动态解耦压制
  白噪声（高频）    → 同时限制 T2 与 T2*，无法用回波消除
  测量方法：
    T2*：Ramsey 序列（π/2 - 延迟 - π/2）
    T2 ：Hahn echo（π/2 - 延迟 - π - 延迟 - π/2）
    噪声谱：Carr-Purcell 序列扫频率，反演功率谱密度 S(ω)

噪声谱学（noise spectroscopy） 是诊断的利器：用变间隔的 CPMG 序列扫描，拟合出环境的功率谱密度 $S(\omega)$。一旦知道 $S(\omega)$，就能预测不同电路的退相干行为，而不只是报一个 $T_2$ 数字。

数值方法：从 Liouvillian 到量子轨迹

向量化与 Liouvillian 超算子

Lindblad 方程对 $\rho$ 是线性的，所以可以「拉直」成一个线性常微分方程。做法是向量化：把 $d \times d$ 的密度矩阵拉成 $d^2$ 维的列向量，把超算子写成 $d^2 \times d^2$ 的矩阵。

向量化：ρ → vec(ρ)（按列堆叠）
超算子：dρ/dt = L ρ  →  d vec(ρ)/dt = L_super vec(ρ)

L_super = -i (I ⊗ H - H^T ⊗ I)/ħ
          + Σ_k γ_k [ conj(L_k) ⊗ L_k - 1/2 (I ⊗ L_k†L_k + (L_k†L_k)^T ⊗ I) ]

维度：d^2 × d^2
  1 比特：4 × 4       → 秒级
  10 比特：10^6 × 10^6 → 内存爆炸
  20 比特：10^12 × 10^12 → 不可行

这是开放系统仿真的根本瓶颈：Liouvillian 的维度是希尔伯特空间的平方。$n$ 个量子比特的密度矩阵有 $4^n$ 个元素，超算子有 $16^n$ 个。10 比特就已经是 $10^6$ 维，20 比特完全不可行。

Krylov 子空间方法

直接算矩阵指数 $\exp(L t)$ 对 $d^2 = 10^6$ 已经不可行（矩阵本身要存 $10^{12}$ 个数）。但 Lindblad 的超算子是稀疏的——每个跳跃算符只连接少数矩阵元。用 Krylov 子空间方法（如 scipy.sparse.linalg.expm_multiply）可以在不显式构造矩阵指数的情况下算 $\exp(Lt) \rho$，内存只需存稀疏矩阵。

import numpy as np
import scipy.sparse as sp
from scipy.sparse.linalg import expm_multiply

# 构造单比特 Lindblad 超算子（振幅阻尼）
sx = sp.csr_matrix(np.array([[0, 1], [0, 0]]))     # sigma_-
I2 = sp.identity(2, format='csr')
H = sp.csr_matrix(np.zeros((2, 2)))                # 无哈密顿量
gamma = 1.0 / 100e-6                                # 1/T1, T1 = 100 µs

L = sp.kron(I2, H) - sp.kron(H.T, I2)               # -i[H,·] 部分（H=0 略）
L = L - 1j * 0                                       # 占位
L = L + gamma * (sp.kron(sx.conj(), sx)
                 - 0.5 * (sp.kron(I2, sx.T @ sx) + sp.kron(sx.T @ sx, I2)))

rho0 = sp.csr_matrix(np.array([[0, 0], [0, 1]], dtype=complex))  # |1><1|
vec0 = rho0.reshape(-1)                              # 向量化
vec_t = expm_multiply(L, vec0, start=0, stop=1e-4, num=5)

量子轨迹（Quantum Trajectories）

对更大的系统，主流的办法是量子轨迹 / 蒙特卡洛波函数（MCWF）。核心思想：把主方程描述的「混合态系综演化」拆成大量「纯态的单次随机演化」。每条轨迹按薛定谔方程演化，但在随机时刻发生「跳跃」（对应环境探测到一个量子），跳跃算符就是 Lindblad 里的 $L_k$。

量子轨迹算法：
  1. 制备纯态 |ψ(0)>
  2. 用有效哈密顿量 H_eff = H - (i/2) Σ_k γ_k L_k† L_k 演化
     （非厄米，态范数随跳跃概率单调下降）
  3. 在每步 dt 内，以概率 dp_k = γ_k dt <ψ|L_k†L_k|ψ> 发生跳跃
       - 若跳跃：|ψ> → L_k|ψ> / ||L_k|ψ>||
       - 若无跳跃：|ψ> → 按 H_eff 演化并归一化
  4. 重复到目标时间，得一条轨迹
  5. 跑 N 条轨迹，用 ρ ≈ (1/N) Σ |ψ_i><ψ_i| 估计密度矩阵

优势：内存 O(d) 而非 O(d²)，轨迹天然并行
代价：统计误差随 1/sqrt(N) 收敛，需大量轨迹

# QuTiP: 主方程求解 vs 量子轨迹求解
import qutip as qt

H = qt.sigmaz() * 0.5 * 2 * np.pi * 1.0        # 1 MHz 拉比
c_ops = [np.sqrt(1/100e-6) * qt.sigmam()]      # T1 = 100 µs
rho0 = qt.basis(2, 1).proj()                   # |1><1>
tlist = np.linspace(0, 500e-6, 200)

# 方式一：主方程（小系统首选）
res_me = qt.mesolve(H, rho0, tlist, c_ops)

# 方式二：量子轨迹（大系统或需并行时）
res_mc = qt.mcsolve(H, rho0, tlist, c_ops, ntraj=500)

# 对比 |1> 布居的衰减，两者应一致（在统计误差内）
print(res_me.expect[0][-1], res_mc.expect[0][-1])

对低纠缠或 Clifford 主导的电路，还有更省的办法：稳定子/Pauli 传播类方法（如 Clifford 电路的 Pauli 框架、T 门计数的随机 Pauli 展开），能处理几十个比特，代价是精度受限。

需要强调的是：开放系统仿真的目标不是「复现硬件的全部细节」，而是抓住主导噪声项。 如果你的目标是评估一个哈密顿量模拟电路在退相干下的表现，那么把 $H$ 写成具体的自旋模型、再挂上 T1/T_φ 两个信道，往往比追求精细的噪声谱更实用。哈密顿量本身的构造与模拟方法，见 哈密顿量模拟 。

在硬件上的体现与缓解

各平台的退相干来源

超导平台的退相干来源最多样：介质损耗（TLS 二能级系统）限制 T1，通量噪声（1/f 谱）限制 T2*，准粒子隧穿造成偶发的大幅衰减。改进手段包括提高谐振腔质量因子、磁屏蔽、改进材料表面处理。

离子阱的问题集中在磁场漂移（可用磁不敏感态压制）和激光相位噪声（可用窄线宽激光与相位锁定）。因为 $T_2$ 可以做到秒级，离子阱在需要长相干时间的算法（如长序列的量子傅里叶变换）上有天然优势。

缓解手段与其对主方程的影响

动态解耦（如 XY-4、CPMG）的本质是在哈密顿量上叠加一串快速 $\pi$ 脉冲，把环境噪声在时域上「平均掉」。从频域看，它等价于把系统的有效噪声谱在低频段压出一个凹口——这正是 稀疏线性代数 里常见的「滤波函数」思路在量子控制中的对应物：设计一个滤波器响应 $F(\omega)$，让 $\int S(\omega) F(\omega) d\omega$ 尽可能小。

建模在误差缓解中的角色

误差缓解的效果高度依赖噪声模型的正确性。举例：零噪声外推（ZNE） 假设错误随噪声强度单调增长，如果实际噪声是相干的（可以相干增强也可以相干抵消），外推就会失效。概率误差抵消（PEC） 则需要精确的信道知识——你要知道每个门对应的 Kraus 算符，才能构造逆信道。

所以建模的质量直接决定缓解的上限：用一个错误标注的模型去设计缓解方案，比不缓解还糟。 这也是为什么噪声谱学、门集层析（GST）、循环基准这些诊断手段值得投入——它们给出的不只是「错误率是多少」，而是「错误是什么类型、怎么来的」。

小结

开放量子系统理论是量子工程的底层语言，它把「硬件为什么不好」翻译成可计算、可仿真的模型。

• 封闭系统用酉演化与纯态；开放系统必须用密度矩阵与 CPTP 信道，因为退相干不可逆。
• 约化密度矩阵通过部分迹得到；纠缠态的部分迹必然是混合态——退相干的本质是信息漏到环境。
• Kraus 表示是量子信道的一般形式，约束是 $\sum K_i^\dagger K_i = I$ 与完全正定性。
• Lindblad 主方程是马尔可夫开放系统的最一般形式；振幅阻尼、相位阻尼、去极化各有对应的跳跃算符。
• Bloch 方程给出 $1/T_2 = 1/(2T_1) + 1/T_\phi$，解释了 $T_2 \le 2T_1$ 的硬约束与 $T_2^*$ 和 $T_2$ 的差异来源。
• 数值上，Liouvillian 超算子维度是 $d^2$，10 比特以上必须转向量子轨迹或张量网络方法。
• 缓解手段（动态解耦、回波、纠错）都能在噪声模型里找到对应项，建模精度决定了缓解的上限。

实践建议：给任何量子硬件建噪声模型时，先测 $T_1$、$T_2^*$、$T_2$ 三个时间尺度，用 $1/T_2 = 1/(2T_1) + 1/T_\phi$ 反推纯退相位率，判断噪声是 T1 限制还是相位限制。然后做一次噪声谱学，看是 $1/f$ 还是白噪声主导——这决定了动态解耦能带来多少收益。最后用 QuTiP 的 mesolve 在小规模（≤6 比特）上验证模型，再用 mcsolve 或张量网络方法扩展到你的实际电路规模。模型不需要完美，但必须正确反映主导噪声的类型，否则误差缓解的设计会南辕北辙。