1. 课程导入与全景概览
在理想化的量子力学教材里，一个封闭量子系统的演化由薛定谔方程完全决定：给定初态和哈密顿量，任意时刻的态矢量都是确定的。然而真实的量子系统几乎从不孤立——它会自发辐射光子、与热库交换能量、被环境测量、经历退相干。这类开放量子系统的演化不再是单一波函数的幺正演化，而是密度矩阵的完全正定保迹映射，其主方程（Lindblad 方程）为
dρ/dt = -i[H, ρ] + Σ_k ( L_k ρ L_k† - ½{L_k†L_k, ρ} )其中 H 是系统哈密顿量，L_k 是第 k 个耗散通道的跳跃算符（如 √γ σ₋ 描述自发辐射）。直接数值积分这个密度矩阵方程在维度上非常昂贵：一个 n 量子比特系统的 ρ 有 4ⁿ 个元素，而波函数只有 2ⁿ 个振幅。Monte Carlo 量子轨迹（Quantum Trajectories，也称随机波函数方法）正是为绕开这一维度瓶颈而生的核心算法。
核心思想：把"平均"拆成"随机样本"
量子轨迹方法的关键洞察是：密度矩阵可以写成大量纯态系综的平均，ρ(t) = E[ |ψ(t)⟩⟨ψ(t)| ]。于是我们不去演化 ρ，而是演化许多条随机纯态轨迹，每条轨迹在两次"跳跃"之间按有效非厄米哈密顿量连续演化：
H_eff = H - (i/2) Σ_k L_k† L_k在时间步 dt 内，系统以概率 δp = dt · Σ_k ⟨ψ|L_k†L_k|ψ⟩ 发生一次量子跳跃（如发射一个光子），跳跃后态被 L_k|ψ⟩ 重新归一化；若未跳跃，则用 H_eff 演化并归一化。对 N 条轨迹求平均，即可逼近 ρ(t)。误差随 N 按 1/√N 收敛——这正是"Monte Carlo"之名的由来。
在开放量子系统与量子光学中的定位
量子轨迹并非纯粹的数值技巧，它同时具有清晰的物理诠释。在量子光学中，它对应量子跳跃法（Quantum Jump Method）：把探测器记录到的光子到达时刻视为一次跳跃，未探测到光子的区间则对应条件演化。因此轨迹方法天然适合处理单次实验层面的问题——例如"给定一个持续监测的原子，某次运行中它在何时发射了光子、发射后处于什么态"。这与主方程只给出系综平均形成互补。
在方法谱系中，它位于几个经典工具之间：
| 方法 | 演化对象 | 适用场景 | 主要代价 |
|---|---|---|---|
| 薛定谔方程 | 纯态 |ψ⟩ | 封闭系统 | 无法描述耗散 |
| Lindblad 主方程 | 密度矩阵 ρ | 小维度开放系统 | 维度 4ⁿ 爆炸 |
| 量子轨迹 | N 条纯态 | 大维度、单次记录、光子计数 | 统计误差 1/√N |
典型应用包括：超导量子比特的弛豫与退相干模拟、腔 QED 中的光子统计、量子测量与反馈控制、以及量子纠错中的错误传播分析。当系统维度较大而耗散通道稀疏时，轨迹方法往往比直接积分主方程快几个数量级。
本课程的学习路径与工具准备
本课程共 30 章，整体按"原理—实现—应用"三层推进。前段建立数学与物理基础：Lindblad 形式、有效哈密顿量、跳跃概率的推导与归一化细节；中段进入实现层：离散时间步算法、跳跃判据的随机数处理、轨迹平均与误差估计、常见数值陷阱（如 dt 过大导致的概率越界）；后段面向应用：与量子光学可观测量对接、条件态与反馈、以及性能优化策略。
动手环境建议准备 Python 3.10+ 与 NumPy、SciPy、QuTiP。QuTiP 的 mcsolve 可直接作为参考实现，但课程会要求你从零写出核心循环，以便真正掌握跳跃判据与归一化的每一步。下面是一段最小可运行骨架，用于建立直觉（示例参数，非实测结果）：
import numpy as np
# 示例：二能级原子自发辐射，γ=1，dt=0.01
gamma, dt, Ntraj = 1.0, 0.01, 2000
sm = np.array([[0,1],[0,0]]) # σ₋
Heff = -0.5j*gamma*sm.conj().T@sm # H=0 时的有效哈密顿量
psi0 = np.array([0,1], dtype=complex) # 激发态
rho_acc = np.zeros((2,2), dtype=complex)
for _ in range(Ntraj):
psi = psi0.copy()
for step in range(500):
p_jump = dt*gamma*np.real(psi.conj()@sm.conj().T@sm@psi)
if np.random.rand() < p_jump:
psi = sm@psi
else:
psi = (np.eye(2) + Heff*dt)@psi
psi /= np.linalg.norm(psi)
rho_acc += np.outer(psi, psi.conj())
rho = rho_acc/Ntraj # 激发态布居应趋近 e^{-γt}运行后观察 rho[1,1]（激发态布居）随时间的衰减，应与解析解 e^{-γt} 在统计误差内吻合；若偏差明显偏大，优先检查 dt 是否过大、归一化是否遗漏、跳跃概率是否出现负值。这套"写循环—对解析解—查误差"的流程，将贯穿后续所有实操章节。