2. 自由粒子高斯波包扩散的数值模拟

问题拆解：自由波包的两种等价视角

自由粒子（V(x)=0）的高斯波包有一个非常干净的解析结论：位置空间中的概率密度始终是高斯型，其宽度按

σ(t) = σ₀ · √(1 + (ℏt / (2mσ₀²))²)

增长。这个式子说明扩散完全由初始宽度 σ₀ 和质量 m 决定，与波包中心动量 p₀ 无关（p₀ 只让整个包平移）。本案例要做的是：用两种数值路径分别演化同一个初始波包，再把测得的宽度曲线与上式逐点比较。

两条路径的数学依据不同，正好互相校验：

• 动量空间路径：自由粒子动量本征态是平面波，演化只是给每个动量分量乘相位 e^{−iEₖt/ℏ}，Eₖ=ℏ²k²/(2m)。这一步是精确的，误差只来自离散化与 FFT。
• 位置空间路径：用分裂算符法，把传播子 e^−iHΔt/ℏ 近似为 e^{−iVΔt/2ℏ}e^−iTΔt/ℏe^{−iVΔt/2ℏ}。自由情形 V=0，势能因子退化为 1，只剩动能因子，因此每步都是“FFT→乘相位→IFFT”。

下面这张图给出两条路径共享同一初始态、各自演化、最后汇合到宽度比较环节的数据流结构，并标出关键参数与判据。

数值实现：两条路径的代码

下面给出一个自洽的脚本，包含初始态构造、两条演化路径、宽度测量与误差输出。所有参数在脚本内定义，可直接运行。

import numpy as np

# ---------- 参数 ----------
hbar, m = 1.0, 1.0
sigma0, x0, p0 = 1.0, 0.0, 2.0
L, N = 200.0, 2048                 # 空间范围与网格点数
x = np.linspace(-L/2, L/2, N, endpoint=False)
dx = x[1] - x[0]
k = 2*np.pi*np.fft.fftfreq(N, d=dx)  # 对应 FFT 顺序的波数

# ---------- 初始高斯波包 ----------
psi0 = (1/(2*np.pi*sigma0**2))**0.25 \
       * np.exp(-(x-x0)**2/(4*sigma0**2)) * np.exp(1j*p0*x/hbar)
psi0 /= np.sqrt(np.sum(np.abs(psi0)**2)*dx)   # 归一化

def width(psi):
    rho = np.abs(psi)**2
    norm = np.sum(rho)*dx
    xm = np.sum(x*rho)*dx/norm
    x2 = np.sum(x**2*rho)*dx/norm
    return np.sqrt(x2 - xm**2)

def sigma_theory(t):
    return sigma0*np.sqrt(1 + (hbar*t/(2*m*sigma0**2))**2)

# ---------- 路径 A：动量空间一次演化 ----------
def evolve_momentum(t):
    phi0 = np.fft.fft(psi0)
    phi_t = phi0 * np.exp(-1j*hbar*k**2*t/(2*m))
    return np.fft.ifft(phi_t)

# ---------- 路径 B：位置空间分裂算符逐步演化 ----------
def evolve_split(dt, nsteps):
    psi = psi0.copy()
    phase = np.exp(-1j*hbar*k**2*dt/(2*m))   # V=0，势能因子为 1
    for _ in range(nsteps):
        psi = np.fft.ifft(np.fft.fft(psi)*phase)
    return psi

# ---------- 比较 ----------
dt, nsteps = 0.01, 2000
t_final = dt*nsteps
psiA = evolve_momentum(t_final)
psiB = evolve_split(dt, nsteps)

sA, sB, sT = width(psiA), width(psiB), sigma_theory(t_final)
print(f"t={t_final:.2f}  σ_A={sA:.6f}  σ_B={sB:.6f}  σ_th={sT:.6f}")
print(f"相对误差 A: {abs(sA-sT)/sT:.2e}   B: {abs(sB-sT)/sT:.2e}")

运行后应看到 σ_A 与 σ_B 在 1e−3 量级内一致，且都贴近 σ_th。若把 dt 改成 0.5 而保持总时间不变，路径 B 的误差会明显上升——这是分裂算符在自由情形下唯一会引入的误差来源：相位因子 exp(−iℏk²Δt/2m) 在 k 较大时随 Δt 快速振荡，采样不足就会混叠。

扩散规律的定量观察

为了看清“扩散”本身，而不是只看终点，建议在多个时刻采样并记录宽度。下表是示例参数下由解析式直接算出的预期值（σ₀=1, m=1, ℏ=1），可作为数值结果的对照基准：

两个特征时间值得记住：短时 σ(t)≈σ₀[1+(ℏt/2mσ₀²)²/2]，增长是二次的；当 t 超过 2mσ₀²/ℏ（本例为 2）后，σ 近似线性增长，斜率 ℏ/(2mσ₀)。也就是说，初始越窄的波包扩散越快——这是色散的直接体现，也是为什么动量空间路径里高频分量相位转得更快。

常见偏差与排查方向

• 宽度偏小且随时间恶化：网格范围 L 不够。扩散后波包尾部被截断，⟨x²⟩ 被低估。本例 t=10 时 σ≈5，波包有效范围约 ±4σ≈±20，L=200 足够；若把 L 缩到 40 就会出问题。
• 路径 B 与 A 不一致但各自都接近解析式：检查两条路径是否用了同一个归一化初始态，以及 FFT 的 k 排列是否与相位因子一致（用 np.fft.fftfreq 而非 linspace）。
• 中心位置漂移异常：自由波包中心应以 p₀/m 匀速平移，本例 p₀=2、m=1，t=10 时中心应在 x=20 附近。若明显不符，检查初始相位 exp(ip₀x/ℏ) 的符号约定。

把上述检查做成一个循环，对每个采样时刻同时输出 σ_A、σ_B、σ_th 与波包中心，就能在一张表里同时确认“扩散规律正确”和“两条路径等价”两件事。