2. 二能级系统与哈密顿量构建

单量子比特是量子计算中最小的信息载体，而它的物理实现几乎总可以抽象为一个二能级系统。要把“受驱动演化”变成可运行的仿真，第一步不是写演化算符，而是把物理对象翻译成数学对象：用态矢量描述量子态，用泡利矩阵搭建可观测量与哈密顿量，再区分静态项与含时驱动项。本节围绕这条主线，给出可直接落地的构建方法与代码实现。

态矢量与布洛赫球约定

在计算基下，一个量子比特的纯态写成二维复矢量：

|ψ⟩ = α|0⟩ + β|1⟩ = [α, β]^T,   |α|² + |β|² = 1

其中 |0⟩ = [1,0]^T，|1⟩ = [0,1]^T。为避免符号混乱，本节统一约定：泡利 Z 矩阵的本征态就是计算基，Z|0⟩ = +|0⟩，Z|1⟩ = −|1⟩；因此 |0⟩ 对应布洛赫球北极，|1⟩ 对应南极。这个约定会直接影响后面驱动项中 σx 的相位与旋转方向，全文与图示保持一致。

任意纯态还可用布洛赫矢量表示：⟨σx⟩、⟨σy⟩、⟨σz⟩ 三个期望值构成单位球面上的一个点。对 |ψ⟩ = cos(θ/2)|0⟩ + e^{iφ}sin(θ/2)|1⟩，有 ⟨σz⟩ = cosθ，⟨σx⟩ = sinθcosφ，⟨σy⟩ = sinθsinφ。这一表示的价值在于：二能级系统的演化在布洛赫球上就是刚体旋转，后面验证仿真结果时，可以直接看矢量是否保持单位长度、旋转轴是否正确。

泡利矩阵与哈密顿量的物理含义

三个泡利矩阵为：

σx = [[0, 1], [1, 0]]
σy = [[0, -i], [i, 0]]
σz = [[1, 0], [0, -1]]

它们与单位矩阵 I 一起构成 2×2 厄米矩阵的一组完备基。任何二能级哈密顿量都可写成 H = h₀I + hxσx + hyσy + hzσz，其中 h₀ 只贡献整体相位、不改变可观测的能级差，真正决定动力学的是 (hx, hy, hz) 这个“等效磁场”矢量。于是“构建哈密顿量”本质上就是确定这个三维矢量如何随时间变化。

以驱动型超导量子比特或核磁共振中的自旋为例，典型形式为：

H(t) = -(ħωq/2)σz + ħΩ(t)/2 · [cos(ωd t)σx + sin(ωd t)σy]

第一项是静态项，ωq 为比特频率（能隙/ħ），决定 |0⟩ 与 |1⟩ 的能量差；第二项是含时驱动项，Ω(t) 为拉比频率（驱动幅度），ωd 为驱动频率。取 ħ = 1 是仿真中的常见简化，本节代码也采用该约定，此时能量与角频率同量纲。注意符号：静态项前的负号使 |0⟩ 能量高于 |1⟩，若你的物理系统约定相反，只需整体改号，但必须在全文统一。

从物理参数到矩阵：构建流程

实际编码时，建议把“参数—矩阵”的映射拆成三层，避免把物理量硬编码进演化函数。第一层是参数层：ωq、ωd、Ω、相位 φ 与总时长 T；第二层是矩阵层：由泡利矩阵组合出 H_static 与 H_drive(t)；第三层是求解层：对本征值问题或含时演化做数值处理。下图展示了这一构建关系与两类哈密顿量的分支。

静态哈密顿量的本征值求解

先处理不含驱动的静态情形，这是理解能级结构的基础。取示例参数 ωq = 2π×5.0 GHz，ħ = 1，则 H₀ = −(ωq/2)σz。由于 σz 已是对角矩阵，本征值可直接读出：E₀ = −ωq/2（对应 |0⟩），E₁ = +ωq/2（对应 |1⟩），能隙 ΔE = ωq。用代码求解时，应走通用本征值路径，以便后续加入失谐或耦合项后仍能复用：

import numpy as np

sx = np.array([[0, 1], [1, 0]], dtype=complex)
sy = np.array([[0, -1j], [1j, 0]], dtype=complex)
sz = np.array([[1, 0], [0, -1]], dtype=complex)
I2 = np.eye(2, dtype=complex)

wq = 2 * np.pi * 5.0e9      # 比特角频率，示例值
H0 = -0.5 * wq * sz

evals, evecs = np.linalg.eigh(H0)   # 厄米矩阵用 eigh
print(evals)                        # [-1.5708e10, 1.5708e10]
print(evecs)                        # 列向量为本征态

这里有两个工程细节值得强调。其一，必须用 eigh 而非 eig，前者利用厄米性保证本征值为实数且本征向量正交归一；其二，eigh 返回的本征值按升序排列，因此索引 0 对应低能级。若你的物理约定是 |0⟩ 为低能级，就需要在绘图时把标签与索引对应清楚，否则能级图会“上下颠倒”。

绘制能级图时，横轴可放参数（如失谐 δ = ωd − ωq），纵轴放能量。当 H = −(δ/2)σz + (Ω/2)σx 时，本征值为 ±√(δ² + Ω²)/2，能级在 δ = 0 处出现最小间隔 Ω（避免交叉），这正是驱动下能级排斥的直观体现。用一段循环扫描 δ 并记录 evals，即可得到能级随失谐变化的曲线，比单点数值更能说明问题。

含时哈密顿量的离散化实现

含时项无法一次性对角化，标准做法是把总时长 T 切成 N 步，每步 Δt = T/N，在每步内把 H(t) 近似为常数（一阶分段常数近似），再用矩阵指数推进态矢量：

from scipy.linalg import expm

def H_drive(t, Om, wd, phi):
    return 0.5 * Om * (np.cos(wd*t + phi) * sx + np.sin(wd*t + phi) * sy)

def evolve(psi0, T, N, wq, Om, wd, phi=0.0):
    dt = T / N
    psi = psi0.astype(complex).copy()
    for k in range(N):
        t = (k + 0.5) * dt          # 中点采样，提高一阶精度
        H = -0.5 * wq * sz + H_drive(t, Om, wd, phi)
        psi = expm(-1j * H * dt) @ psi
    return psi

几个参数选择直接影响结果可信度。第一，采样点取步中点而非步起点，可把分段常数近似从一阶提升到二阶精度，代价几乎为零。第二，Δt 必须远小于驱动周期 2π/ωd 和拉比周期 2π/Ω，经验上每周期至少取 50～100 步；若只关心布洛赫矢量慢变包络，可先做旋转坐标系变换去掉 ωd 的快振荡，再对等效哈密顿量积分，这样步长可放大一个量级以上。第三，expm 对 2×2 矩阵开销很小，但对长时演化仍建议预分配数组、避免在循环内反复构造矩阵。

结果判读与常见陷阱

运行上述代码后，建议做三项检查。其一，态矢量模长应始终为 1，若出现明显漂移，说明步长过大或矩阵非厄米。其二，在共振条件 ωd = ωq、初态 |0⟩ 下，⟨σz⟩ 应呈余弦振荡，振荡频率约等于 Ω（拉比频率），这是判断驱动项符号与幅度是否正确的最直接依据。其三，关闭驱动（Ω = 0）时，⟨σz⟩ 应保持为常数，仅整体相位随时间变化，若出现布洛赫矢量进动，说明静态项符号或频率设置有误。

把以上环节串起来，就得到一条清晰的构建链：物理参数 → 泡利基展开 → 静态/含时分支 → 本征值或矩阵指数求解 → 能级图与布洛赫轨迹判读。后续加入脉冲包络、失谐扫描或噪声项时，只需替换参数层与含时项，求解层与判读逻辑可以原样复用。