# 离子晶体与射频微运动数值模型（v2）

本文件对应 `dist/dynamics.js` 与 `tests/dynamics.test.mjs`。这是具有明确方程、单位、参考解与收敛检查的经典数值模型。界面中的真实物理时间与动画播放速度必须分别标注；20 MHz 射频不能被当作数秒周期的实际离子运动。

## 1. 单位、常数与参考依据

模块内部统一使用 SI；接口位置使用 µm、设备距离使用 mm、RF 使用 MHz、模态使用 Hz、电势使用 V、电场使用 V/m、能量使用 eV。`rfVpeak` 是四极势定义中的零到峰值系数，并非 Vpp 或示波器在 50 Ω 端口测得的电压；实际电极连接方式必须按相同定义标定。

默认粒子是单电荷的 ⁴⁰Ca⁺，质量 `39.962042283 u`，`u=1.66053906660e−27 kg`。这与原平台物理模块一致。质量可通过 `massKg` 覆盖；目前不提供混合质量/混合电荷的晶体求解。

理论依据和可追溯原始文献：

- D. F. V. James, *Quantum dynamics of cold trapped ions, with application to quantum computation*, Applied Physics B **66**, 181–190 (1998), [DOI](https://doi.org/10.1007/s003400050373)，[作者预印本](https://arxiv.org/abs/quant-ph/9702053)。依据其无量纲 Coulomb 长度、轴向平衡与正则模定义，扩展到完整三维 Hessian。
- D. Leibfried, R. Blatt, C. Monroe, D. Wineland, *Quantum dynamics of single trapped ions*, Reviews of Modern Physics **75**, 281–324 (2003), [期刊原文](https://journals.aps.org/rmp/abstract/10.1103/RevModPhys.75.281)，[NIST 全文](https://www.nist.gov/system/files/documents/2017/05/09/rmp1.pdf)。依据线性 Paul 阱、Mathieu 方程、赝势和微运动的推导。
- D. J. Berkeland et al., *Minimization of ion micromotion in a Paul trap*, Journal of Applied Physics **83**, 5025–5033 (1998), [DOI](https://doi.org/10.1063/1.367318)，[NIST 作者全文](https://tf.nist.gov/general/pdf/1226.pdf)。依据残余直流电场引起的位移和超额微运动的理论。补偿电极的场转换系数仍是平台教学标定值，不能从该文献当成实际装置实测常数。

## 2. 三维 Coulomb 平衡与正则模

接口：

```js
solveCrystal({ count, radialXHz, radialYHz, axialHz,
  massKg?, temperatureK? })
```

`count` 必须为 0…10 整数；`count=0` 显式返回空位置与空模态，`structure:'empty'`。非空晶体要求三个频率和质量为有限正数。数值较大并不代表给定RF设备一定能实现该赝势；实际频率由主物理模块的 Mathieu/Floquet 计算后传入。

势能是

\[
U=\frac{m}{2}\sum_i(\omega_x^2x_i^2+\omega_y^2y_i^2+\omega_z^2z_i^2)
 +\frac{e^2}{4\pi\epsilon_0}\sum_{i<j}\frac{1}{|\mathbf r_i-\mathbf r_j|}.
\]

令 `l³=e²/(4πε₀mωz²)`、`u=r/l`，以 `mωz²l²` 为能量单位。全部优化在无量纲空间进行，不用肉眼调整的离子间距。数值梯度和 Hessian 使用 Coulomb 项的解析导数。

求解顺序为：

1. 对每个离子数缓存轴向线链 Newton 解。
2. 在该线链计算完整三维 Hessian，直接检查最软横向特征值。临界零模不判为严格稳定。
3. 线链失稳时沿软模、两个确定性三维扰动初始化，使用正定移位 Newton 及下降线搜索最小化势能；若仍有负曲率，则再沿负模扰动求解。
4. 取候选中的最低能量解，核查完整三维残余力与 Hessian。仅当无量纲残余力小于 `1e−7`，且没有低于相对 `1e−7` 阈值的负曲率时，`equilibriumConverged=true`。

这提供可核查的局部稳定平衡；确定性多初值搜索不是对任意参数下全局基态的数学证明。近结构转变可能存在亚稳构型；教学中应比较能量和稳定性，而不能把每次返回值都称为唯一的晶体结构。该模型是理想谐赝势模型，不包括 RF 周期内的多离子位置重排、非谐电极场、镜像电荷、有限温度相变或量子零点位移。

返回字段：

- `positionsUm:[{x,y,z},…]`：零温三维平衡位置，µm。
- `normalModes:[{frequencyHz,eigenvalue,omegaSquaredSI,axis,vector,unstable},…]`：按 Hessian 特征值排序。`eigenvalue=ωmode²/ωz²`，`omegaSquaredSI` 单位为 s⁻²。`vector` 长度为 `3N`，按 `x₁,y₁,z₁,x₂,y₂,z₂,…` 排列且单位范数。`axis` 为占比超过 98% 的 x/y/z，否则为 mixed。出现实质负特征值时 `unstable=true`，显示频率为零，不能把零频当作普通稳定模态。
- `forceResidualNorm`：全部 `3N` 分量合成的残余力二范数，单位 **N**；另返回 `forceResidualDimensionless`。小力并不替代 Hessian 稳定性核查。
- `linearChainStable`：沿指定 z 轴的线链是否严格稳定。若晶体旋转成沿 x 的线列，此字段仍为 false。
- `softModeHz`：实际返回平衡的最软模频率；`linearSoftModeHz`/`linearSoftEigenvalue` 专门描述原 z 线链最软横向模。
- `structure`：中心化位置协方差的秩判定，秩≤1 为 linear，秩2 为 zigzag，秩3 为 3d。在强各向异性外，zigzag 标签泛指平面构型，并不保证交替锯齿拓扑。
- `hasRotationalZeroMode`：存在临界/旋转近零模；不能据此宣称所有模严格稳定。
- `harmonicThermalRmsUm`：由经典谐模方差 `kBT/(mωmode²)` 相加得到的各离子三轴 RMS。`temperatureK` 不改变零温平衡；这不是有限温度分子动力学。零模方差被省略，此时正温度的 `thermalApproximationValid=false`，不得把省略后的有限数值当成该方向的真实热宽度。接近软化、大热振幅或非晶体态时也不宜使用此近似。
- `energyEV`、`coulombLengthUm`、`iterations`、`equilibriumConverged`。

缓存按离子数与两个频率比建立；绝对频率/质量仅做精确尺度变换。缓存最多160种频率比。常规线链缓存调用适合实时读数；10离子强三维重排首次计算约 0.1–0.2 s，应避免把连续拖动时大量独立优化阻塞在绘图循环。

### 独立解析基准

两个离子的轴向坐标为 `z/l=±2^(−2/3)`，轴向 COM 频率为 `fz`，伸缩模为 `√3 fz`，x 横向摇摆模为 `√(fx²−fz²)`。所以两离子的 z 线链在 `fx=fz` 具有真正的零模。三离子轴向位置为 `0, ±(5/4)^(1/3)l`，轴向三个频率是 `fz,√3 fz,√(29/5)fz`。测试还独立重算三维 Coulomb 力、检查 COM 频率、特征向量正交性以及尺度律 `l∝m^(−1/3)fz^(−2/3)`。

## 3. 瞬时电势、电场与平均势

接口 `potentialAt({xUm,yUm,zUm,rfPhase,...trapParameters})`；RF 相位为弧度。采用一致的局部理想四极势约定：

\[
\Phi=\frac{V_{RF}\cos(\Omega t)+U_r}{2r_0^2}(x^2-y^2)
 +\frac{\kappa U_z}{z_0^2}\left[z^2-\frac{x^2+y^2}{2}\right]-E_{res}x.
\]

`Eres = strayFieldVm − compensationV × 0.05/r0`，`r0` 为米。0.05 是显式的教学几何场系数，需要用真实电极场求解/测量替换才能用于具体硬件设计。

`potentialV`/`electricPotentialV` 单位V；`electricFieldVm:[Ex,Ey,Ez]=−∇Φ` 单位V/m。返回的 `rfPseudopotentialEV` 仅RF平均势，`pseudopotentialEV` 包含静态电势和RF平均势：

\[
U_{ps}=e\Phi_{DC}+\frac{e^2|\mathbf E_{RF,amp}|^2}{4m\Omega^2}.
\]

能量转换为 eV 后输出。`pseudopotentialEV` 包含端帽的径向反束缚和杂散场线性项，允许为负；不要为了画图裁成处处正值。

`q=2eV_RF/(mr0²Ω²)`，正RF势约定下 `qx=−q`、`qy=+q`。小q赝势预测的三个频率由以下 ω² 计算，而不是替代主模块的数值 Floquet 结果：

\[
\omega_{x,y}^2=\frac{e^2V_{RF}^2}{2m^2\Omega^2r_0^4}
\pm\frac{eU_r}{mr_0^2}-\frac{e\kappa U_z}{mz_0^2},\qquad
\omega_z^2=\frac{2e\kappa U_z}{mz_0^2}.
\]

此模型满足真空中的 Laplace 方程；完整边界电极形状、实际机械开口及制造误差尚未求解。热图显示理想局部解析场，不能称作 Blender 几何的有限元场计算。

## 4. 真实 RF 时间轨迹与收敛

```js
integrateRF({rfMHz,rfVpeak,r0mm,z0mm,endcapV,kappa,
  radialDcV,strayFieldVm,compensationV,
  initialPositionUm:[x,y,z],initialVelocityMS:[vx,vy,vz],
  cycles:160,stepsPerCycle:96,massKg?})
```

RK4 直接积分 `r̈=eE(r,t)/m`，没有用赝势谐运动冒充射频轨迹。一个 20 MHz 周期为0.05 µs，160周期就是8 µs；绘图可放慢播放，但横坐标必须仍为µs。单离子轨迹不加未声明的阻尼、激光冷却、随机碰撞或其它离子的 Coulomb 力。

允许1…2000周期、每周期16…1024步且请求总步数≤1,000,000。额外以 `h/2` 再积分，返回匹配物理终点时的 `convergence.positionErrorUm`、`velocityErrorMS` 与 `relativePositionError`。相对位置误差以初始位置/最终位置的较大范数（最小1e−6 µm）归一化，低于1e−4标为 passed。单个终点差是数值收敛诊断，不是严格全程误差界，也不是实验参数的绝对准确度。粗细轨迹因越界在不同时间停止时，差值为 null。

返回采样至多3000点：`timesUs`、`positionsUm:[[x,y,z],…]`、`velocitiesMS`、`kineticEnergyEV`、`instantaneousEnergyEV`、`secularEnergyEV`。后者仅为“瞬时速度动能+平均势”的诊断；它仍含微运动，**不是守恒的久期能量**。真实RF驱动持续做功，瞬时总能量同样不守恒；不以此虚构能量漂移误差。纯轴向且无RF径向运动的测试单独检验谐振子能量守恒。

轨迹达到 `sqrt(x²+y²)>0.7r0` 或 `|z|>0.7z0` 时终止，`termination='left-local-model-boundary'`。这是离开局部理想场有效区，**不是计算出的电极碰撞位置或真实逃逸鞍点**。非有限运算同样终止且保留最后有限状态。返回 boundary 字段明确该解释。

`secularFrequencyEstimateHz` 从每RF周期一次的等相位位置取正向过零周期中位数估计，剔除了快速RF摆动。每轴若振幅不足或不足三次正向过零，则为 null；`frequencyResolutionHz=1/Twindow` 提示有限观测窗口，不意味着达到该统计精度。当前单离子模型在线性稳定区可和独立 Floquet 特征指数比较；接近零模、RF Nyquist边界或短窗口时应增加周期数，不能把 null 写成0 Hz。

小q、静态残余场主导时，`x0≈eEres/(mωx²)`，超额微运动幅度 `Amm≈|q x0|/2`，在本电势相位约定下 `x≈x0[1+(q/2)cosΩt]`。接口同时返回解析近似 `displacementApproxUm`、`analyticExcessMicromotionUm`；完整轨迹与该近似的偏差包含有限q修正和初始条件引起的久期运动。实际过剩微运动对冷却和光谱的影响由独立光学/冷却模块处理。

## 5. 实际验证

运行 `node --test tests/dynamics.test.mjs` 可验证16项可重复基准。参考参数 `20 MHz、200 Vpeak、r0=0.55 mm、Uz=8 V、κ=0.12、z0=2.5 mm`：

- `q=0.20217549`，小q径向近似约1.426307 MHz；完整积分估计约1.437895 MHz；主模块独立 Floquet 约1.437943 MHz，轨迹过零估计相对差约3.3e−5。
- 400RF周期（20µs）、96步/周期时，h与h/2末端位置差约 `3.66e−5 µm`，相对约 `9.8e−6`（该初始条件为 `[3,2,1] µm`、零速度，无残余场）。
- 24→48→96步/周期的末端误差近似按四阶下降，而不是把界面帧率当作积分精度。
- 两离子解析模基准相对误差为浮点数舍入量级；独立三维 Coulomb 力和势梯度有限差分测试用于发现单位/符号错误。

这些数字验证的是所声明方程的求解可靠性；几何系数、真实激光参数、电极边界和仪器噪声仍须实测标定后才能声称与某台真实装置定量一致。