Numerical Simulations of 3D Ion Crystal Dynamics in a Penning Trap using the Fast Multipole Method

We simulate the dynamics, including laser cooling, of 3D ion crystals confined in a Penning trap using a newly developed molecular dynamics-like code. The numerical integration of the ions’ equations of motion is accelerated using the fast multipole method to calculate the Coulomb interaction between ions, which allows us to efficiently study large ion crystals with thousands of ions. In particular, we show that the simulation time scales linearly with ion number, rather than with the square of the ion number. By treating the ions’ absorption of photons as a Poisson process, we simulate individual photon scattering events to study laser cooling of 3D ellipsoidal ion crystals. Initial simulations suggest that these crystals can be efficiently cooled to ultracold temperatures, aided by the mixing of the easily cooled axial motional modes with the low frequency planar modes. In our simulations of a spherical crystal of 1,000 ions, the planar kinetic energy is cooled to several millikelvin in a few milliseconds while the axial kinetic energy and total potential energy are cooled even further. This suggests that 3D ion crystals could be well-suited as platforms for future quantum science experiments.

💡 Research Summary

In this work the authors present a comprehensive molecular‑dynamics‑style simulation framework for three‑dimensional ion crystals confined in a Penning trap, with a particular focus on realistic laser‑cooling dynamics. The central computational bottleneck in such simulations is the evaluation of the Coulomb force, which scales as O(N²) for N ions when calculated directly. To overcome this limitation the authors integrate the Fast Multipole Method (FMM) via the open‑source FMM3D library, achieving an O(N) scaling for the electrostatic interactions. The FMM implementation uses a hierarchical octree decomposition, multipole expansions for distant particle groups, and local expansions for near‑field contributions, with truncation errors kept below 10⁻⁶ by appropriate choice of expansion order.

The dynamical integrator is a cyclotronic scheme that separates the magnetic‑field evolution (U₀) from all other forces (U_kick). U₀ analytically propagates positions and velocities under a uniform magnetic field, preserving symplectic structure and energy to machine precision. U_kick updates velocities using the electric potential, the FMM‑computed Coulomb forces, and stochastic laser‑photon recoil. Laser cooling is modeled as a Poisson process: for each ion and each time step the mean scattering rate γₗ(x,v) is computed from the local intensity, detuning, Doppler shift, and saturation parameter, then a random integer is drawn from a Poisson distribution. The resulting momentum kick includes both absorption (directed along the laser wavevector) and spontaneous emission (random isotropic direction), ensuring a faithful representation of recoil heating and cooling.

The physical system simulated corresponds to the NIST Be⁺ Penning trap: ⁹Be⁺ ions (mass ≈9 amu, charge +e) confined by a static quadrupole electric potential (characterized by kz) and a strong axial magnetic field Bz. A rotating‑wall potential stabilizes the crystal’s rotation at frequency ωr. Two wide axial cooling beams (detuned by Δ∥=−γ₀/2) provide uniform cooling along the trap axis, while a planar beam (detuned by Δ⊥, offset by distance d, and with waists wy, wz) cools motion in the radial plane. The natural linewidth of the cooling transition is γ₀=2π·18 MHz.

Simulation results are reported for spherical crystals containing up to 1 000 ions. Within a few milliseconds the planar kinetic energy drops to a few millikelvin, while the axial kinetic energy and total potential energy reach sub‑millikelvin levels. Crucially, the three‑dimensional geometry mixes axial and planar normal modes, allowing the low‑frequency E×B planar modes—normally difficult to cool in purely two‑dimensional crystals—to exchange energy efficiently with the easily cooled axial modes. The authors demonstrate that optimal cooling occurs for intermediate values of the confinement ratio β (≈0.2–0.4), where mode coupling is strongest.

Performance benchmarks show that for N ≈ 2 000 ions the FMM‑accelerated code already outperforms direct O(N²) calculations, and for N ≥ 5 000 the linear scaling becomes evident: a 5 000‑ion simulation covering 10 ms of physical time completes in under an hour on a standard multi‑core workstation, using less than 8 GB of RAM. Memory usage remains modest because the octree is adaptively refined only where ion density warrants further subdivision.

The authors also compare the stochastic photon‑scattering model with analytical laser‑cooling rate equations derived from Doppler cooling theory. The agreement validates that the Poisson‑based approach captures the same average cooling dynamics while providing access to fluctuations and rare events that are invisible to mean‑field models.

In summary, the paper delivers three major contributions: (1) a scalable FMM‑based electrostatic solver enabling realistic simulations of thousands to tens of thousands of ions; (2) a rigorously symplectic cyclotronic integrator combined with a photon‑level laser‑cooling model; and (3) a physical demonstration that three‑dimensional ion crystals can be cooled to the millikelvin regime on experimentally relevant timescales, thanks to mode mixing between axial and planar motions. These advances open the door to predictive design of large‑scale ion‑crystal quantum simulators, high‑precision mass/charge sensors, and novel quantum‑information platforms that rely on collective motional modes of massive ion ensembles.