2. 含退相干二能级系统的布洛赫方程仿真

物理模型与仿真思路

自由感应衰减（FID）描述的是：二能级系统被激励到相干叠加态后，撤去驱动，横向相干随退相干速率衰减、纵向布居向热平衡弛豫的过程。用布洛赫方程描述时，引入两个唯象时间常数：纵向弛豫时间 T1（布居恢复）与横向退相干时间 T2（相干衰减）。在旋转坐标系下，无驱动（Ω=0）且共振（Δ=0）时，方程退化为三个解耦的一阶线性方程：

du/dt = -u / T2
dv/dt = -v / T2
dw/dt = -(w - w_eq) / T1

其中 u、v 为横向布洛赫分量，w 为纵向布居差（归一化后热平衡值 w_eq = -1 或按约定取 0，本例取 w_eq = -1 表示基态布居占优）。初始条件取激励后 u(0)=0, v(0)=1, w(0)=0，即完全横向相干态。解析解为 v(t)=exp(-t/T2)、w(t)=-1+exp(-t/T1)。仿真任务就是用数值积分复现这两条曲线，并从衰减曲线中反提取 T2。

数值实现（Python）

采用四阶龙格-库塔法积分布洛赫方程，扫描三组退相干速率，同时用解析指数模型做对照。关键点：时间步长需远小于最短的 T2，本例取 dt = 0.01（最短 T2=2，约 200 步/衰减常数，足够）。

import numpy as np

def bloch_fid(T1, T2, t_max=20.0, dt=0.01, w_eq=-1.0):
    """无驱动共振条件下的自由感应衰减，RK4 积分布洛赫方程。"""
    n = int(t_max / dt)
    t = np.linspace(0, t_max, n + 1)
    u = np.zeros(n + 1); v = np.zeros(n + 1); w = np.zeros(n + 1)
    u[0], v[0], w[0] = 0.0, 1.0, 0.0          # 激励后的初始相干态

    def deriv(u, v, w):
        return (-u / T2, -v / T2, -(w - w_eq) / T1)

    for i in range(n):
        k1 = deriv(u[i], v[i], w[i])
        k2 = deriv(u[i] + 0.5*dt*k1[0], v[i] + 0.5*dt*k1[1], w[i] + 0.5*dt*k1[2])
        k3 = deriv(u[i] + 0.5*dt*k2[0], v[i] + 0.5*dt*k2[1], w[i] + 0.5*dt*k2[2])
        k4 = deriv(u[i] + dt*k3[0],     v[i] + dt*k3[1],     w[i] + dt*k3[2])
        u[i+1] = u[i] + dt/6*(k1[0] + 2*k2[0] + 2*k3[0] + k4[0])
        v[i+1] = v[i] + dt/6*(k1[1] + 2*k2[1] + 2*k3[1] + k4[1])
        w[i+1] = w[i] + dt/6*(k1[2] + 2*k2[2] + 2*k3[2] + k4[2])
    return t, u, v, w

def fit_T2(t, v):
    """对 v(t) 取对数做线性拟合，斜率即 -1/T2。"""
    mask = v > 1e-3
    slope, _ = np.polyfit(t[mask], np.log(v[mask]), 1)
    return -1.0 / slope

# 三组退相干速率：T2 从 2 到 8，T1 固定为 10
for T2 in (2.0, 4.0, 8.0):
    t, u, v, w = bloch_fid(T1=10.0, T2=T2)
    T2_fit = fit_T2(t, v)
    print(f"设定 T2={T2:4.1f}  拟合 T2={T2_fit:6.3f}  "
          f"v(0)={v[0]:.3f}  v(5)={v[500]:.4f}  w(20)={w[-1]:.4f}")

预期输出（示例运行结果，非实测）：设定 T2=2.0 时拟合值约 2.00，T2=4.0 时约 4.00，T2=8.0 时约 8.00；v(5) 分别约为 0.082、0.287、0.535；w(20) 均趋近 -1（因 T1=10，20 时间单位约 2 个 T1，恢复约 86%）。拟合值与设定值一致，说明数值积分与指数衰减模型吻合。

衰减曲线与弛豫时间提取

下图把三组 T2 下的横向分量 v(t) 衰减曲线、纵向分量 w(t) 恢复曲线，以及从 ln v(t) 斜率反提取 T2 的判读逻辑放在一起。注意 T2 越小，曲线越陡，ln v 斜率绝对值越大；w(t) 的恢复只由 T1 决定，与 T2 无关，这正是区分两种弛豫的判据。

结果判读与验证要点

• 曲线形状：v(t) 应为单调指数衰减，半对数坐标下呈直线；若出现振荡，说明驱动项未置零或失谐未清零，方程未真正退化为 FID 形式。
• 拟合一致性：拟合 T2 与设定值相对误差应小于 1%（本例步长下约 0.1% 量级）；误差偏大时优先检查 dt 是否过大或拟合区间是否包含 v 接近数值噪声的尾部。
• T1 与 T2 解耦：改变 T2 时 w(t) 曲线应完全重合，这是验证纵向与横向弛豫项未相互串扰的直接判据。
• 物理约束：通常 T2 ≤ 2T1；若设定 T2 > 2T1，虽数值上可积分，但已超出纯退相位与弛豫的物理可实现范围，仅作数学练习。