4. 二能级系统与拉比振荡实战

二能级系统是量子动力学中最简单、也最有教学价值的模型。一个量子比特、一个自旋、一个原子中的两个能级，在忽略其它自由度后都可以抽象成二维希尔伯特空间。本节我们不再停留在公式推导，而是从零构建一个含时哈密顿量，用 scipy.integrate.solve_ivp 求解概率幅随时间的演化，并绘制布居数振荡曲线，最后通过改变驱动频率观察共振与失谐带来的差异。

一、物理模型与约定

取二能级基矢为 |0⟩ 与 |1⟩，对应能量 E₀ 与 E₁。为方便，令 E₀ = 0，能级差为 ℏω₀，即 ω₀ = (E₁ − E₀)/ℏ。我们采用一个经典驱动场（半经典近似）作用于系统，驱动频率为 ω，耦合强度（拉比频率）为 Ω。在旋转波近似下，系统在相互作用绘景中的有效哈密顿量为：

H(t) = (ℏ/2) * [ [Δ, Ω],
                  [Ω, -Δ] ]

其中 Δ = ω − ω₀ 为失谐量。当 Δ = 0 时称为共振，此时布居数在 |0⟩ 与 |1⟩ 之间完全翻转；当 Δ ≠ 0 时，振荡幅度被压缩，且振荡频率变为广义拉比频率 Ω_R = √(Δ² + Ω²)。这里我们统一约定：初始时刻 t = 0 系统处于 |0⟩，即概率幅向量为 [1, 0]ᵀ；驱动从 t = 0 开始持续作用。全文与图示均采用这一约定。

二、构建含时哈密顿量与求解器

在数值实现中，我们把概率幅写成复数向量 psi = [c0, c1]，薛定谔方程为 iℏ dψ/dt = H(t)ψ。取 ℏ = 1 简化单位，则 dψ/dt = −i H(t) ψ。由于 H(t) 依赖时间，solve_ivp 的右端函数必须接收 (t, y) 两个参数，并在内部实时构造哈密顿量。下面给出一个可直接运行的完整示例：

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# 参数（示例值，单位取 ℏ = 1）
omega0 = 5.0        # 能级差对应频率
Omega  = 1.0        # 拉比频率（耦合强度）
omega  = 5.0        # 驱动频率，此处先取共振
Delta  = omega - omega0

def hamiltonian(t, Delta, Omega):
    return 0.5 * np.array([[Delta, Omega],
                           [Omega, -Delta]], dtype=complex)

def rhs(t, y, Delta, Omega):
    psi = y[0:2] + 1j * y[2:4]      # 实部 + i*虚部
    H = hamiltonian(t, Delta, Omega)
    dpsi = -1j * (H @ psi)
    return np.concatenate([dpsi.real, dpsi.imag])

# 初始态 |0>
psi0 = np.array([1.0, 0.0], dtype=complex)
y0 = np.concatenate([psi0.real, psi0.imag])

t_span = (0.0, 20.0)
t_eval = np.linspace(0.0, 20.0, 2000)

sol = solve_ivp(rhs, t_span, y0, t_eval=t_eval,
                args=(Delta, Omega), rtol=1e-9, atol=1e-11)

psi_t = sol.y[0:2] + 1j * sol.y[2:4]
p0 = np.abs(psi_t[0])**2
p1 = np.abs(psi_t[1])**2

这里有一个关键技巧：solve_ivp 默认处理实数向量，而概率幅是复数。我们通过把实部与虚部拼接成四维实向量 y = [Re(c0), Re(c1), Im(c0), Im(c1)] 来绕过这一限制，在右端函数内部再还原成复数。这是量子仿真中非常常用的写法。此外，rtol 与 atol 设得较严，是因为布居数振荡对相位误差敏感，容差过松会导致振幅漂移。

三、布居数振荡曲线与结果判读

运行上述代码后，绘制 p0 与 p1 随时间的变化，会看到两条互补的余弦曲线：p0 从 1 开始下降，p1 从 0 开始上升，在 t = π/Ω 附近完成第一次完全翻转，随后周期性返回。共振时理论解为 p1(t) = sin²(Ωt/2)，周期 T = 2π/Ω。以 Ω = 1 为例，第一次完全翻转发生在 t ≈ π ≈ 3.14，完整周期约 6.28。数值曲线应与这一解析结果高度吻合，这本身就是一次有效的正确性验证。

判断仿真是否可信，可以检查三点：其一，任意时刻 p0 + p1 应恒等于 1（概率守恒），若明显偏离说明积分误差过大；其二，共振时最低点 p0 应接近 0，若始终无法归零，可能是失谐未设为零或耦合项符号有误；其三，振荡周期应与 2π/Ω_R 一致。这三点构成了本节最实用的结果判读依据。

四、共振与失谐的对比分析

把驱动频率从 ω = ω₀ 改为 ω = ω₀ + 0.5，即 Δ = 0.5，保持 Ω = 1 不变，重新运行同一段代码，只需修改 omega 一个参数。此时广义拉比频率 Ω_R = √(0.25 + 1) ≈ 1.118，振荡周期缩短到约 5.62；同时 p1 的最大值不再是 1，而是 Ω²/Ω_R² = 1/1.25 = 0.8，即最多只有 80% 的布居能转移到 |1⟩。这正是失谐的物理含义：驱动频率偏离能级差，能量交换不再完全，系统在 |0⟩ 与 |1⟩ 之间做“不彻底”的往复。

从工程角度看，这一现象直接对应量子比特操控中的频率校准问题。若实际驱动频率与比特频率存在偏差，门操作的保真度会下降，因为目标态无法被完全到达。因此在实际实验中，往往需要先扫描驱动频率，找到 p1 峰值最高的那个频率点，才认为实现了共振。数值仿真在这里的作用，是让我们在动手之前就清楚“失谐会以什么方式、多大程度地破坏翻转”，从而为参数扫描设定合理范围。

五、参数扫描的实用做法

与其逐个手改参数，更高效的方式是把失谐量做成数组，循环求解并记录每个 Δ 下 p1 的最大值。示例中可取 Δ 从 −2 到 2 共 41 个点，对每个点调用一次 solve_ivp，记录 np.max(p1)。结果会呈现一条以 Δ = 0 为中心的洛伦兹型峰线：Δ = 0 时峰值 1，随 |Δ| 增大而单调下降。这条曲线就是最朴素的“共振峰”，它把抽象的失谐效应转化成了可直接观察、可量化比较的数据。

至此，我们完成了一个从哈密顿量构建、含时求解、曲线绘制到共振/失谐对比的完整闭环。这套“构造 H(t) → 复数转实向量 → solve_ivp → 提取概率幅 → 计算布居数 → 参数扫描”的流程，是后续处理更复杂量子系统时的通用骨架，值得反复练习到能独立默写。