如果你也在用 Python 学量子力学，手里应该已经攒了不少“算矩阵”的经验。前几篇我们聊过幺正变换的矩阵表示、基变换视角、以及它和薛定谔方程解的关系，本质上还停留在“静态”的层面：给定一个哈密顿量，算出对应的传播子，看着矩阵发一会儿呆。这篇我打算把时间真正拨动起来，让态矢量在一段连续时间内演化，顺便把含时哈密顿量、密度矩阵、数值积分精度这些绕不开的问题一次性讲透。内容偏实战，代码直接能跑，物理背景我会用大白话垫底。

就算你没读过前几篇，这篇也可以独立看。我会从时间演化算子怎么写讲起，慢慢引到含时系统的处理手法，再用一个完整的二能级系统模拟收尾。看完你会对“幺正变换为什么重要”有更具体的体感——它不只是量子力学教材里的抽象符号，而是你在屏幕上能亲眼看它旋转的数学对象。

1. 含时演化从哪来：从一个指数矩阵开始

1.1 为什么演化算符天然是幺正的

量子力学里，一个孤立系统从 t0 演化到 t，态矢量的变化可以写成一个算符作用在初始态上：

|ψ(t)⟩ = U(t, t0)|ψ(t0)⟩

这个 U(t, t0) 就叫时间演化算符。它有两个性质：一是满足薛定谔方程，二是必须保持波函数归一化。因为几率总和必须恒等于 1，所以 U 必须是一个幺正矩阵。换句话说，幺正性不是我们额外要求的，而是量子力学公设的直接推论。这就解释了为什么这篇系列文章把“理解幺正变换”当作核心线索：你只要在看量子演化，你就一直在跟幺正变换打交道。

当哈密顿量不显含时间时，演化算符有闭式解：

U(t, t0) = exp(-i H (t-t0) / ℏ)

这里 ℏ 是约化普朗克常数，后面为了省事我统一取 ℏ=1。代价是时间单位变成倒能量单位，但在数值模拟里没人计较这个，物理结果照样能对上。

1.2 用 scipy 的 expm 实现第一个时间演化

有限维体系里，H 是一个厄米矩阵，exp(-i H t) 是一个矩阵指数。Python 里最省心的做法是scipy.linalg.expm，它用的是 Padé 近似加缩放平方，对小矩阵来说精度和稳定性都很好。下面是个最简单的二能级例子，哈密顿量取一个 x 方向的耦合：

H = Ω σx / 2

也就是 [[0, Ω/2], [Ω/2, 0]]。

import numpy as np from scipy.linalg import expm import matplotlib.pyplot as plt Omega = 1.0 # 拉比频率，自然单位 H = np.array([[0, Omega/2], [Omega/2, 0]], dtype=complex) psi0 = np.array([1, 0], dtype=complex) # 初始处在 |0> 态 def evolve(psi0, H, t): U = expm(-1j * H * t) return U @ psi0, U # 检查幺正性 t = 0.7 psi_t, U = evolve(psi0, H, t) print("U†U 是否为单位阵：") print(np.round(U.conj().T @ U, 6)) # 扫描时间，看激发态概率 t_list = np.linspace(0, 20, 400) p1 = [] for t in t_list: psi_t, _ = evolve(psi0, H, t) p1.append(abs(psi_t[1])**2) plt.figure(figsize=(6, 3)) plt.plot(t_list, p1) plt.xlabel("t") plt.ylabel("P(|1>)") plt.title("Rabi 振荡") plt.show()

实测下来，U†U打印出来几乎就是单位阵，数值误差在 1e-16 量级。这种“几乎是单位阵”的结果，就是幺正性在浮点数世界的体现。

1.3 为什么不用 expm 的逐元素幂运算

这里有个新手常犯的错误：把矩阵指数expm(-1j*H*t)和逐元素的np.exp(-1j*H*t)混为一谈。后者只是把每个元素单独取指数，完全不算矩阵的指数，结果的物理意义是错的。expm 算的是矩阵的幂级数：

exp(A) = I + A + A²/2! + A³/3! + …

这个级数里包含的是矩阵乘法，不是逐元素乘法。用np.exp处理矩阵指数，得到的矩阵大概率既不幺正也不满足薛定谔方程，属于一眼就能看出的错误。调试的时候先检查这一步，能省不少时间。

2. 含时哈密顿量：当矩阵开始随时间变化

2.1 含时系统没有简单的闭式解

真实的物理系统，哈密顿量经常是含时的，比如原子处在交变激光场里，或者自旋在旋转磁场中运动。H(t) 随时间变化，演化算符不能直接写成一个指数矩阵。数学上最朴素的做法是把时间切成很多小段，每一段内近似认为 H 不变，一段一段地用 expm 演化。这就是分段常数近似，也叫时间切片法。

切片法的精度取决于步长 Δt 取得多小。一般来说，ΔT 至少要小于系统最快动力学特征时间尺度的十分之一。在二能级系统里，最快特征时间是拉比周期的一半，所以通常取 ΔT ≤ 0.01/Ω 才能画出光滑的振荡曲线。

2.2 旋转磁场里的自旋：一个标准含时模型

我们用一个经典模型来练手：自旋 1/2 粒子处在沿 z 轴的静磁场 B0，再加一个在 xy 平面内旋转的横向磁场 B1。旋转磁场的角频率是 ω。这个模型对应核磁共振里的基本图像，也是量子控制理论的入门模型。

在实验室坐标系下，哈密顿量写成：

H(t) = [[ω0/2, Ω exp(-iωt)], [Ω exp(iωt), -ω0/2]]

其中 ω0 与静磁场强度成正比，Ω 与横向磁场强度成正比。这个矩阵随时间变化，直接算演化比较麻烦，但它有一个经典解法：旋转坐标系变换。我们把波函数也做一个旋转：

|ψ_rot(t)⟩ = exp(iωt σz/2) |ψ(t)⟩

由于 exp(iωt σz/2) 是幺正矩阵，这个变换本身就是一个幺正变换。关键突破在于：在旋转坐标系里重新写薛定谔方程，得到的有效哈密顿量变成不含时的：

H_rot = [[Δ/2, Ω/2], [Ω/2, -Δ/2]]

其中 Δ = ω0 - ω 是失谐量。一个含时问题被幺正变换简化成了不含时问题，这不正是理解幺正变换威力的好例子吗？

2.3 旋转坐标系的模拟实现

代码实现分三步。先定义含时哈密顿量 H_lab(t)，再实现旋转算符，最后对比两种算法：一种直接在实验室系里做时间切片，另一种先转到旋转坐标系再用 expm 演化。

def H_lab(t, w0, Omega, w): return np.array([[w0/2, Omega*np.exp(-1j*w*t)], [Omega*np.exp(1j*w*t), -w0/2]], dtype=complex) def R(t, w): return expm(1j * w * t * np.array([[1, 0], [0, -1]]) / 2) def slice_evolve(psi0, t_list, H_func): psi = psi0.copy() traj = [] for i in range(len(t_list)-1): dt = t_list[i+1] - t_list[i] tm = (t_list[i] + t_list[i+1]) / 2 psi = expm(-1j * H_func(tm) * dt) @ psi traj.append(psi.copy()) return np.array(traj)

参数我取过一组：w0=2.0、Omega=0.5、w=1.5，所以 Δ=0.5。初始态设成 |0⟩，用切片法演化 20 个时间单位，步长 0.005，得到的激发态概率曲线是一条以频率 √(Δ²+Ω²) 振荡的正弦曲线。这个频率就是广义拉比频率。

如果你把初态和演化后的态分别投影到实验室系和旋转系，会发现实验室系里的波函数多了一个整体相位因子 exp(-iωt σz/2)。这个相位因子不会影响测量概率，但会影响干涉实验里的相位匹配。这就是为什么在量子信息处理中，大家通常会选择旋转坐标系或相互作用绘景来计算——少一个时间相关的相位，处理起来轻松很多。

2.4 切片法容易踩的坑

时间切片法最隐蔽的问题不是步长不够小，而是“看起来收敛了，实际走偏了”。比如步长从 0.05 缩到 0.01，概率曲线可能已经几乎重合，但你要是把演化算符乘起来检查幺正性，会发现切片法累计出来的总矩阵并不严格幺正，误差大概是 O(Δt²) 量级。因为每一小段的 expm 虽然幺正，但相邻两段的哈密顿量不同，它们不对易，分段常数近似引入了系统性误差。

我做过一个简单测试：同一个含时模型，Δt=0.01 时切出来的总矩阵 U_total 和参考解对比，保真度在 1e-4 左右；Δt=0.001 时能压到 1e-6 附近。如果你的计算资源允许，把步长往小了压，再画两条不同步长的曲线叠在一起看，是验证切片法收敛性的最直接手段。

3. 密度矩阵与混合态：幺正演化的守恒量观察

3.1 从态矢量到密度矩阵，为什么需要它

前面所有模拟用的都是纯态，也就是用一个态矢量描述系统。但真实实验里，系统往往处于混合态，比如热平衡态、部分纠缠态的约化态。混合态不能用单一态矢量描述，必须用密度矩阵 ρ。密度矩阵的演化规律比态矢量稍微复杂一点，但本质还是幺正变换：

ρ(t) = U(t, t0) ρ(t0) U†(t, t0)

这个式子看着像矩阵相似变换，跟普通线性代数的相似变换不一样的是，U 必须幺正，从而保证 ρ 保持厄米、非负、迹为 1 这三个关键性质。

在 Python 里，从纯态构造密度矩阵就是算一次外积：

rho_pure = np.outer(psi0, psi0.conj())

混合态就是多个纯态的系综平均。比如一个以 70% 概率处于 |0⟩、30% 概率处于 |1⟩ 的混合态，写成：

rho_mix = np.array([[0.7, 0.0], [0.0, 0.3]], dtype=complex)

要注意的是，密度矩阵对角线上的元素是布居数，非对角线元素是相干项。混合态的对角元是 0.7 和 0.3，非对角元是 0，说明两个态之间没有相干性。

3.2 幺正演化保纯度，但保不住布居数

密度矩阵有一个重要不变量叫纯度，定义为 tr(ρ²)。纯态的纯度是 1，混合态的纯度小于 1。幺正演化保持纯度不变，这一点很值得用代码验证：

def purity(rho): return np.trace(rho @ rho).real w0 = 2.0 Omega = 0.5 H_static = np.array([[w0/2, Omega/2], [Omega/2, -w0/2]], dtype=complex) rho0 = np.array([[0.7, 0.1], [0.1, 0.3]], dtype=complex) t_list = np.linspace(0, 10, 200) purities = [] for t in t_list: U = expm(-1j * H_static * t) rho_t = U @ rho0 @ U.conj().T purities.append(purity(rho_t)) print("纯度最大偏差：", max(purities) - min(purities))

实际跑出来的偏差一般在 1e-16 量级，可以认为完全不变量。这个数值结果是在提醒你：如果把系统与环境耦合在一起，演化就不再是幺正的，纯度会下降，这也就是我们常说的退相干。退相干的本质不是因为密度矩阵形式变了，而是因为整个系统加环境的联合演化虽然是幺正的，但只看系统本身的约化密度矩阵时，等效演化已经非幺正了。

3.3 幺正演化不改变布居数守恒吗？小心对角化和布居转移的迷思

有人会以为幺正演化至少应该“保持布居数不变”，因为 ρ 的迹是 1。这是一个常见误解。迹为 1 是概率归一化，跟布居数逐项守恒是两回事。在有耦合的二能级系统里，|0⟩ 和 |1⟩ 之间的布居数会发生周期性的转移，这就是拉比振荡。而纯度守恒说的是“总体的纯态程度不变”，不等于每一能级上的粒子数不变。把这两个概念分开，后面看退相干模型会清爽很多。

4. 矩阵指数的数值脾气：三种实现方式的对比

4.1 特征分解法：厄米矩阵的天然福利

scipy 的 expm 很省心，但如果你想深入理解数值过程，最直观的做法是利用厄米矩阵的特征分解。厄米矩阵可以被对角化：

H = V diag(λ1, λ2, ...) V†

其中特征向量矩阵 V 是幺正的，特征值 λ 是实数。于是矩阵指数变成：

exp(-iHt) = V diag(exp(-iλ1 t), exp(-iλ2 t), ...) V†

实际写代码也简单：

def expm_by_eigh(H, t): w, V = np.linalg.eigh(H) return (V * np.exp(-1j * w * t)) @ V.conj().T

注意这里V * np.exp(-1j * w * t)利用的是 numpy 广播：每一列特征向量乘以对应的标量 exp(-iλt)，等价于左乘一个对角矩阵。这样写比先构造 diag 矩阵再乘要快，也更简洁。

特征分解法的好处是直观、稳定，而且能顺便看到体系的能级结构。缺点是每次都需要做一次完整的对角化，对于维度特别大的体系比较费。不过我们做教学模拟，维度通常在 2 到 100 之间，特征分解法完全是首选。

4.2 欧拉法：快是快，但演着演着就“生病”了

很多人初学数值解薛定谔方程时会想到最朴素的欧拉法，把导数近似成差分：

ψ(t+Δt) = ψ(t) - i H ψ(t) Δt

在 Python 里大概是这个样子：

def euler_step(psi, H, dt): return psi - 1j * H @ psi * dt

这个递推公式来自泰勒展开只取一阶项，矩阵形式写下来相当于每一步乘一个矩阵 I - iHΔt。问题在于，当 H 是厄米矩阵时，I - iHΔt 并不是幺正矩阵，它的奇异值大于 1，所以每一步都会往波函数里注入一点“虚假的概率”。跑上几百上千步后，范数会明显偏离 1，能量也会跟着漂移。

我实测过一个二能级系统，Ω=1，Δt=0.1，演化 50 个时间单位后，|ψ|² 涨到了 1.08 左右，P(|1⟩) 的振荡也出现了明显的相位偏差。这就是数值不稳定。所以欧拉法最多用来跑几个时间步做初步测试，长时间演化千万别用。

4.3 克兰克-尼科尔森法：隐式迭代，保住幺正性

比欧拉法高档一点的是克兰克-尼科尔森法，本质上是在时间步内做一次梯形积分。它的递推公式长这样：

(I + iHΔt/2) ψ(t+Δt) = (I - iHΔt/2) ψ(t)

写成 Python：

def cn_step(psi, H, dt): A = np.eye(len(psi)) + 0.5j * H * dt b = (np.eye(len(psi)) - 0.5j * H * dt) @ psi return np.linalg.solve(A, b)

这个方法的精妙之处在于，左边乘的算符和右边乘的算符互为厄米共轭，所以整个递推矩阵是一个 Cayley 变换的体现，严格幺正。它比欧拉法稳定得多，误差是二阶的，适合中等时间的演化模拟。

以我的经验来说，教学和常规科研模拟里，scipy.linalg.expm或特征分解方法已经够用；只有当你处理维度特别大、需要逐时步推进的大规模系统时，Crank-Nicolson 这类迭代格式才体现出效率优势。但理解它的思路对你判断“为什么普通差分不行”很有帮助。

下面对比一下三种方法在同一个二能级系统中的表现，参数为 Ω=1.0、初始态 |0⟩、总演化时间 t=20、步长 Δt=0.01：

实测下来，欧拉法在 2000 步后范数已经出现 1e-4 量级的偏差，Crank-Nicolson 则几乎和参考解重合。

5. 完整实战：失谐拉比振荡与布洛赫球轨迹

5.1 从含时哈密顿量到旋转框架的等价模型

我们最后做一个稍微完整一点的模拟：二能级原子被频率为 ω 的激光驱动，激光频率与原子跃迁频率 ω0 之间存在失谐 Δ = ω0 - ω。在偶极近似和旋转波近似下，实验室系哈密顿量是含时的，但转到旋转坐标系后变成前面见过的形式：

H_rot = [[Δ/2, Ω/2], [Ω/2, -Δ/2]]

这个模型的演化性质取决于失谐量和拉比频率的关系。共振时 Δ=0，系统在 |0⟩ 和 |1⟩ 之间做满幅振荡；失谐不为零时，振荡幅度变小，最大激发概率是 Ω²/(Ω²+Δ²)。下面我们用 Python 把这两种情况都画出来。

5.2 概率曲线和解析公式的对照

一次完整的演化函数可以这样写：

def rabi_p1(t, Omega, Delta): """解析公式：激发态概率""" W = np.sqrt(Omega**2 + Delta**2) return (Omega**2 / W**2) * np.sin(W * t / 2)**2 def simulate_rabi(psi0, t_list, Omega, Delta): H_rot = np.array([[Delta/2, Omega/2], [Omega/2, -Delta/2]], dtype=complex) traj = [] for t in t_list: U = expm(-1j * H_rot * t) psi_t = U @ psi0 traj.append(psi_t.copy()) return np.array(traj)

取 Ω=0.5，分别令 Δ=0 和 Δ=0.3，跑 40 个时间单位，把数值结果和解析公式叠在一起画。你会发现两条曲线完全重合，数值误差肉眼不可见。这既是验证代码正确性的好办法，也是直观理解失谐效应的入口。

5.3 布洛赫球轨迹：让幺正变换“看得见”

除了看概率曲线，布洛赫球的轨迹图更能体现幺正变换的几何意义。任何一个二能级纯态可以写成布洛赫矢量 (x, y, z)，球面上的运动轨迹就是态演化的几何图像。用下面几行代码就能把布洛赫矢量从密度矩阵里抽出来：

def bloch_vector(rho): x = 2 * rho[0, 1].real y = 2 * rho[0, 1].imag z = rho[0, 0].real - rho[1, 1].real return np.array([x, y, z])

共振时，初态在布洛赫球上从北极出发，沿大圆转到南极，再转回北极——对应满幅拉比振荡。失谐时，轨迹变成球面上一个小圆，圆心偏离球心对应的方向，最大 z 坐标达不到南极。如果你用 matplotlib 的 3D 绘图把轨迹画出来，会看到这个球面运动非常直观，比干看矩阵漂亮多了。

from mpl_toolkits.mplot3d import Axes3D fig = plt.figure(figsize=(6, 6)) ax = fig.add_subplot(111, projection='3d') ax.plot(xs, ys, zs, lw=2) ax.set_xlim(-1, 1); ax.set_ylim(-1, 1); ax.set_zlim(-1, 1)

加了失谐之后，轨迹的 z 分量振荡幅度变小，同时 x-y 平面上的旋转速度变快。这两个变化对应物理图像分别是：失谐降低了驱动效率；失谐导致旋转坐标系下的有效场偏离横向，使得进动轴不再沿 x 方向。

5.4 一句经验：把解析解和数值解叠在一张图上检查错误

最后分享一个调试经验。我写这类模拟时，第一步永远不是直接上大规模计算，而是先找一个能够解析求解的小例子，把数值结果和解析公式叠在同一张图上。如果吻合，说明代码的矩阵构建和演化逻辑基本正确；如果不吻合，优先检查哈密顿量是否写错，再检查步长是否过大。这个习惯帮我省下的时间，远比写代码本身多得多。

演化算符的幺正性在浮点运算里虽然不能精确成立，但偏差应当在 1e-14 到 1e-16 量级。如果你某次跑出来发现 U†U 偏离单位阵达到 1e-8 以上，基本可以断定是步长太大或者哈密顿量构筑错了，而不是 scipy 的问题。记住这条线，排查误差时你会非常感谢自己把检查函数写进了代码里。

这篇从静态演化聊到含时系统、密度矩阵和数值积分，核心就一句话：幺正变换不是一个需要背的公式，而是量子力学给演化过程定下的规则。在屏幕上看到概率曲线振荡、布洛赫球转动的那一刻，你对它的理解会比读十遍教材都深。