如果你也在用 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]], dtypecomplex) psi0 np.array([1, 0], dtypecomplex) # 初始处在 |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]], dtypecomplex) 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[i1] - t_list[i] tm (t_list[i] t_list[i1]) / 2 psi expm(-1j * H_func(tm) * dt) psi traj.append(psi.copy()) return np.array(traj)参数我取过一组w02.0、Omega0.5、w1.5所以 Δ0.5。初始态设成 |0⟩用切片法演化 20 个时间单位步长 0.005得到的激发态概率曲线是一条以频率 √(Δ²Ω²) 振荡的正弦曲线。这个频率就是广义拉比频率。如果你把初态和演化后的态分别投影到实验室系和旋转系会发现实验室系里的波函数多了一个整体相位因子 exp(-iωt σz/2)。这个相位因子不会影响测量概率但会影响干涉实验里的相位匹配。这就是为什么在量子信息处理中大家通常会选择旋转坐标系或相互作用绘景来计算——少一个时间相关的相位处理起来轻松很多。2.4 切片法容易踩的坑时间切片法最隐蔽的问题不是步长不够小而是“看起来收敛了实际走偏了”。比如步长从 0.05 缩到 0.01概率曲线可能已经几乎重合但你要是把演化算符乘起来检查幺正性会发现切片法累计出来的总矩阵并不严格幺正误差大概是 O(Δt²) 量级。因为每一小段的 expm 虽然幺正但相邻两段的哈密顿量不同它们不对易分段常数近似引入了系统性误差。我做过一个简单测试同一个含时模型Δt0.01 时切出来的总矩阵 U_total 和参考解对比保真度在 1e-4 左右Δt0.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]], dtypecomplex)要注意的是密度矩阵对角线上的元素是布居数非对角线元素是相干项。混合态的对角元是 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]], dtypecomplex) rho0 np.array([[0.7, 0.1], [0.1, 0.3]], dtypecomplex) 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Δt0.1演化 50 个时间单位后|ψ|² 涨到了 1.08 左右P(|1⟩) 的振荡也出现了明显的相位偏差。这就是数值不稳定。所以欧拉法最多用来跑几个时间步做初步测试长时间演化千万别用。4.3 克兰克-尼科尔森法隐式迭代保住幺正性比欧拉法高档一点的是克兰克-尼科尔森法本质上是在时间步内做一次梯形积分。它的递推公式长这样(I iHΔt/2) ψ(tΔt) (I - iHΔt/2) ψ(t)写成 Pythondef 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⟩、总演化时间 t20、步长 Δt0.01方法最终范数最终激发态概率相对误差特征分解参考1.0000000.1439600expm1.0000000.1439600欧拉法1.0000670.1442181.8e-4Crank-Nicolson1.0000000.143961约 1e-6实测下来欧拉法在 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]], dtypecomplex) 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, projection3d) ax.plot(xs, ys, zs, lw2) 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 的问题。记住这条线排查误差时你会非常感谢自己把检查函数写进了代码里。这篇从静态演化聊到含时系统、密度矩阵和数值积分核心就一句话幺正变换不是一个需要背的公式而是量子力学给演化过程定下的规则。在屏幕上看到概率曲线振荡、布洛赫球转动的那一刻你对它的理解会比读十遍教材都深。