铯里德堡原子跃迁电子云演化 —— 物理与数值方法

本演示中所有形状都不是"画"出来的，而是逐点求解真实波动方程得到的。下面是完整的方程与算法链条。

1. 能量：铯的量子缺陷（Rydberg–Ritz）

碱金属价电子在离原子实较远处感受到近似 −1/r 的库仑势，但穿透/极化原子实会使能量偏离氢原子。用实验拟合的量子缺陷 δ 精确描述：

E_{nlj}/hc = − R_Cs / (n − δ_{lj})² , R_Cs = 109736.8626 cm⁻¹ δ_{lj}(n) = δ₀ + δ₂/(n−δ₀)² + δ₄/(n−δ₀)⁴ + δ₆/(n−δ₀)⁶ + δ₈/(n−δ₀)⁸

¹³³Cs 的 δ₀：nS_1/2 4.0493532，nP_1/2 3.5915871，nP_3/2 3.5590676，nD_3/2 2.475365，nD_5/2 2.4663144，nF 0.03339，nG 0.00704。有效主量子数 ν = n − δ 是所有尺度的核心：轨道半径 ~ 2ν² a₀，能级间距 ~ 1/ν³。

2. 径向波函数：Numerov 数值积分

u''(r) = [ l(l+1)/r² + 2V_l(r) − 2E ] u(r) , u(r) = r·R(r) （原子单位 ħ=m_e=e=1）

里德堡态的空间尺度横跨 5 个数量级（原子实 ~1 a₀，外壳 ~10⁴ a₀），均匀网格极不经济。做变量替换 x=√r、w = r^{3/4}R（消去一阶导数项）：

d²w/dx² = g(x)·w , g(x) = (4l(l+1) + 3/4)/x² − 8x²·[ E − V_l(x²) ]

从 x_max=√(2ν(ν+15)) 处的渐近解 w ∝ exp(−r/ν) 出发，用 Numerov 六阶递推

(1 − h²g₊/12)w₊ − (2 + 5h²g₀/6)w₀ + (1 − h²g₋/12)w₋ = 0

向内积分（积分过程中按需整体重标定以防指数溢出），再按 2∫x²w²dx = 1 归一化。反变换得 R(r)=w·r^{−3/4}，u(r)=w·r^{1/4}。左下图中可以直接数出径向节点。

关于内区截断：向内积分对"非本征解"（∝x^{−1/2}）天生不稳定，必须截断。本程序不是简单地砍在固定半径，而是自动判定：在 w 变量下经典允许区的振荡包络近似为常数（因 u 的振幅 ∝ r^{1/4}，而 w = u·r^{−1/4}），于是 (a) 若某处 |w| 超过外部已确立包络的 25 倍 → 判定已被污染； (b) 离心屏障内侧的禁戒区中物理解必须单调衰减，一旦包络回升即截断于极小点。此外使用模型势时以 α_c^{1/3} ≈ 2.5 a₀ 为下界（模型势在更内区不再可信，这也是 ARC 的惯例）。因此对 50S 这类穿透态，约 5 个位于原子实内部的节点会被截掉（右栏会明确标注"N 个在芯区"）—— 它们对里德堡态的半径、偶极矩与可视化的贡献均可忽略。用纯库仑势时程序会一直积到 ~10⁻³ a₀，此时氢原子的节点数与解析解逐一吻合。

3. 原子实：Marinescu 模型势

V_l(r) = − Z_l(r)/r − α_c/(2r⁴)·[1 − e^{−(r/r_c)⁶}] Z_l(r) = 1 + (Z−1)e^{−a₁r} − r(a₃ + a₄r)e^{−a₂r} (Z=55, α_c=15.6440 a₀³)

参数取自 Marinescu, Sadeghpour & Dalgarno, Phys. Rev. A 49, 982 (1994)，对 l=0,1,2,≥3 各一组。可在左栏关闭该势退化为纯库仑，对比高 l 态几乎不变、低 l（穿透）态内区差别明显。

4. 角向：严格球谐函数

Y_lm(θ,φ) = Θ_lm(cosθ)·e^{imφ} , Θ_lm = √[ (2l+1)(l−m)! / (4π(l+m)!) ] · P_l^m(cosθ)

P_l^m 用标准连带 Legendre 递推（含 Condon–Shortley 相位 (−1)^m）计算，非任何解析拟合。

5. 时间演化：两态时间依赖薛定谔方程

Ψ(r,θ,φ,t) = c_A(t)·R_A(r)Y_{l_A m_A}·e^{−iE_A t/ħ} + c_B(t)·R_B(r)Y_{l_B m_B}·e^{−iE_B t/ħ}

|Ψ|² = |c_A A_A|² + |c_B A_B|² + 2|c_A c_B| A_A A_B · cos( Δm·φ − ω_AB t + φ₀ ) ω_AB = (E_B − E_A)/ħ , A = R(r)Θ(cosθ)

第三项（干涉项）就是以跃迁频率振荡的振荡偶极，也正是原子辐射/吸收光子的经典对应。Δm=0（π 跃迁）时电子云沿 z 轴来回"晃动"；|Δm|=1（σ 跃迁）时干涉项含 e^{±iφ}，电子云绕 z 轴旋转——直接看到角动量与光子自旋的交换。

包络 c_A, c_B 由旋转坐标系中的二能级方程精确求解：

i·d/dt (c_A, c_B)ᵗ = ½ [ [−Δ, Ω], [Ω, Δ] ] (c_A, c_B)ᵗ c_A(t) = cos(Ω_R t/2) + i(Δ/Ω_R)sin(Ω_R t/2) , c_B(t) = −i(Ω/Ω_R)sin(Ω_R t/2) , Ω_R=√(Ω²+Δ²)

于是 |c_B|² = (Ω²/Ω_R²)sin²(Ω_R t/2)：共振时完全反转（π 脉冲），失谐时反转不完全——这就是右下角布居图上看到的行为。

6. 可观测的偶极矩

⟨z⟩(t) = 2·d_z·Re[ c_A* c_B e^{−iω_AB t} ] , d_z = ⟨A|z|B⟩ = (∫u_A u_B r dr)·⟨l_A m|cosθ|l_B m⟩ ⟨l+1,m|cosθ|l,m⟩ = √[ ((l+1)²−m²) / ((2l+1)(2l+3)) ]

径向积分用与波函数同一套网格上的 2∫x⁴w_Aw_B dx 计算。里德堡态的 d 可达 10³ a₀ e（比基态大 4~5 个数量级），这正是里德堡原子对微波场极端敏感、以及里德堡阻塞效应的根源。

7. 渲染如何保持"概率密度"的诚实性

点位置按混合密度 p = ½(|A_A|² + |A_B|²) 重要性抽样：径向从 x 网格上的精确 CDF 反变换抽样（无拒绝损失），角向按 |Θ_lm|² 拒绝抽样，φ 均匀。因此点的空间密度本身就是概率密度。
每帧在 GPU 顶点着色器中用该点预存的 (A_A, A_B, φ, 包络) 精确重算 |Ψ(t)|²，亮度 ∝ |Ψ|²/(2p) ∈ [0,1]。于是"点密度 × 亮度"始终 ∝ 真实 |Ψ(r,t)|²。
相位着色模式：色相 = arg Ψ（循环色环），可直接看到节点面处的 π 相位跳变与 σ 跃迁的相位绕数。

8. 时间尺度说明（唯一的"艺术处理"）

n≈50 的 S→P 微波跃迁频率约 10¹⁰–10¹¹ Hz（周期 ~10 ps），而实验 Rabi 频率约 10⁶–10⁷ Hz（周期 ~100 ns），两者相差约 10⁴。若用同一慢放倍率，一个 Rabi 周期需要几个小时。因此本演示对光学拍频与Rabi 包络使用两个独立的慢放倍率，二者的真实值与各自倍率都在右栏明确列出。其余一切（能量、波函数、节点、半径、偶极矩、失谐依赖）均为真实数值结果。

参考：Gallagher, Rydberg Atoms (CUP 1994)；Marinescu et al., PRA 49, 982 (1994)；Weber & Sansonetti, PRA 35, 4650 (1987)；Šibalić et al., ARC, Comput. Phys. Commun. 220, 319 (2017)。

¹³³Cs 里德堡跃迁 · 电子云演化 Numerov 数值求解真实径向薛定谔方程 + 球谐函数 + 时间依赖两态演化