A direct numerical simulation of many interacting ions in a Penning trap with a rotating wall is presented. The ion dynamics is modeled classically. Both axial and planar Doppler laser cooling processes are modeled using stochastic momentum impulses based on two-level atomic scattering rates. The plasmas being modeled are ultracold two-dimensional crystals made up of hundreds of ions. We compare Doppler cooled results directly to a previous linear eigenmodes analysis. Agreement in both frequency and mode structure is obtained. Additionally, when Doppler laser cooling is applied, the laser cooled steady state plasma axial temperature agrees with the Doppler cooling limit. Numerical simulations using the approach described and benchmarked here will provide insights into the dynamics of large trapped-ion crystals, improving their performance as a platform for quantum simulation and sensing.

## I. INTRODUCTION

A Penning trap with a rotating wall potential and Doppler laser cooling can produce stable ultracold non-neutral ion plasma crystals with temperatures of a fraction of a millikelvin.^{1–4} Such ultracold ion crystals enable interesting research at the forefront of different areas of physics including atomic physics,^{5–9} quantum optics and metrology,^{10,11} quantum simulation,^{12–14} and basic plasma physics,^{15} including studies of strongly coupled plasmas that model dense astrophysical matter.^{16–18} For many applications, lower ion temperatures are desirable. For example, lower ion temperatures can improve the fidelity of quantum simulations with trapped ions^{19} and enable the engineering of stronger interactions by using shorter wavelength spin-dependent forces.^{13} Colder temperatures also improve the stability of the ion crystal, which advances the prospects for single site optical manipulation and detection in multidimensional crystals, akin to what has been achieved with linear ion crystal arrays^{20} and with an “atom microscope” for neutral atoms.^{21} Single site optical manipulation and detection will enable the preparation of complicated entangled states through the implementation of variational quantum simulation protocols^{22} on large trapped ion crystals.

The study of the thermodynamics of ultracold ions in Penning traps is not just a facilitator for certain experiments but also interesting basic plasma physics research in its own right. The ions in these systems form a strongly coupled plasma with complex collective modes of motion.^{15} The behavior of energy transfer by nonlinear mode coupling and other heating and cooling mechanisms is subtle and complex.^{23–25}

The primary means of cooling the ions in a Penning trap are various forms of laser cooling including Doppler cooling,^{26–29} sideband cooling,^{9} and electromagnetically induced transparency (EIT) cooling.^{1,30} The basic principles for laser cooling ultracold ions are similar to those used for laser-cooling neutral atoms in traps, but there are important differences that can make it more challenging to fully understand the cooling limits and to optimize the laser cooling geometry and parameters such as laser intensity and detuning.^{28} A significant difference from neutral atom experiments is that the ion motion is subject to the trapping electric and magnetic fields. In the axial magnetic field of the Penning trap, the ion plasma rotates, resulting in a change in the ion velocity relative to the cooling laser beam, which is stationary in the laboratory frame. The typical velocity of the ions gives rise to Doppler shifts that can be large relative to the natural atomic linewidth of the ions.

A second major difference is the strong Coulomb interaction between the ions. The Coulomb force couples the ions producing collective modes of motion. In contrast to the cooling dynamics of neutral atoms—which is largely a single particle phenomenon—it is necessary to take the collective nature of the ion motion into account to fully understand the cooling of the ion crystal. Even an analysis in terms of collective linear eigenmodes cannot fully account for the ion dynamics. Nonlinear coupling between modes gives rise to frequency shifts^{31,32} and broadening of resonances, resulting in nontrivial frequency spectra.

In this paper, we present results from a first principles classical model of the ion dynamics, analogous to a molecular dynamics simulation, including the trapping fields with a rotating wall potential, Coulomb interactions between ions, and Doppler laser cooling. We focus on the experimentally relevant case of a single-plane crystal generated with trapping parameters similar to those used in recent quantum simulation and sensing experiments.^{10–14} These simulations provide a more complete picture of the Doppler cooling dynamics and cooling limits in the complex Penning trap geometry and can help with optimizing the trap configuration and cooling parameters in current and future experiments. The simulation presented in this paper goes beyond previous modeling work^{12,27,28} in that it captures the combined effects of finite temperature, nonlinear coupling, and dynamics occurring at disparate time scales in a single unified simulation code. The rest of this article is organized as follows: Section II discusses the simulation model, including the Penning trap Hamiltonian, Doppler laser cooling, the time integration of the ion equations of motion, and numerical convergence with respect to time step. We then compare the simulation results with both the axial and the planar linear eigenmode analysis in Sec. III. In Sec. IV, we show that Doppler cooled simulations run to a steady state and compare with the results from equilibrium models.^{26–28} We conclude with a summary and discuss future extensions of this work.

## II. MODEL AND COMPUTATIONAL ALGORITHM

In this section, we describe our mathematical model and computational approach for the simulation of ultracold ions in a Penning trap with a rotating wall. For a sketch of a Penning trap and a detailed discussion of the confining fields, see Ref. 12.

### A. Trap forces and Hamiltonian

We treat the ions as classical point particles with velocity $vi$ and position $xi$. Excluding the cooling laser for now, the motion of the ions is governed by the Hamiltonian

where the Hamiltonian for the motion in the strong axial B-field is

including the vector potential **A** for the homogeneous axial magnetic field $B=(0,0,Bz)T$ in the Penning trap. To be noted, in the following, we use the terms axial, out-of-plane, and parallel as directions parallel to the magnetic field and planar, in-plane, and perpendicular as directions perpendicular to the magnetic field. We choose our coordinate system such that the Penning trap magnetic field $B=\u2207\xd7A$ is parallel to the *z* axis. With that choice of coordinate system, we may choose $A=\u2212yBzx\u0302$ with $x\u0302=x/x$ the unit vector along **x**. In Eq. (2), *N* is the number of ions in the trap and *m _{i}* and

*q*are the mass and charge of ion

_{i}*i*. The electrostatic potential

contains the potential $\phi trap$ for the electrodes in the Penning trap, the rotating wall potential $\phi wall$, and the Coulomb potential for the interaction between the ions. In the vicinity of the ion crystal, the trap potentials are well approximated by harmonic potentials. We parametrize $\phi trap(x)$ as

and the time-dependent rotating wall potential $\phi wall(x,t)$ as

where *δ* is the dimensionless parameter that characterizes the strength of the rotating wall potential to the trapping potential $\phi trap(x)$ and *θ* is the azimuthal angle in cylindrical coordinates. Next, we transform to the rotating frame,

Using the above transformation, Eq. (6), to the rotating frame, and combining Eqs. (4) and (5), we obtain

where

Note that in the rotating frame, the applied potential is time independent, Eq. (7). Due to the rotating wall potential, Eq. (5), the externally applied trapping fields produce a Hamiltonian that is time dependent in the laboratory frame. However, the applied potential is time independent in the rotating frame [Eq. (7)], resulting in energy conservation in this frame.^{15} This is a useful check for the validity of the model in the absence of laser cooling. The trap Hamiltonian in the rotating frame, *H _{R}*, is

where all “R” subscripts indicate the rotating frame and $rRi2=xRi2+yRi2$. The third term is new and is due to the Lorentz and centrifugal forces. The Coulomb interaction energy term only involves the difference between the coordinates of the ions. The transformation to the rotating frame conserves the length of vectors, so the expression is the same in the rotating or lab frame.^{15} After the system reaches equilibrium, all the particles nearly rigidly rotate as a crystal with the constant frequency *ω _{R}*. In Sec. IV, the measurement of temperature via kinetic energy is made in the rotating frame so as not to include the energy associated with the rigid rotation of the ion crystal.

### B. Doppler cooling

In addition to the conservative dynamics described by the Hamiltonian, the ions are subject to radiation pressure forces due to cooling lasers. Our model of laser cooling is based on resonance fluorescence for a driven two-level atom. We begin by discussing our approach for a single cooling laser. The photon scattering rate for a driven two-level atom is given by^{27}

where *γ*_{0} is the natural linewidth of the atom transition (in radians per second), $S(x)$ is the saturation parameter, and $\Delta (v)=\Delta 0\u2212k\xb7v$ is the detuning of the atomic transition from the laser frequency including the first order Doppler shift with $\Delta 0$ the detuning at rest. The vector **k** is the wave vector of the cooling laser. We assume that the atoms scatter photons with this rate with Poissonian number statistics. The saturation parameter is spatially dependent due to the intensity variations of the cooling laser beam. For a Gaussian beam, we have

where *ρ* is the distance of the atom from the axis of the beam and *W _{y}* is the 1/

*e*radius of the intensity of the beam.

To incorporate laser cooling into our numerical simulations, we proceed as follows: first, we compute the mean number of photons scattered by ion *i* in time interval Δ*t*,

The velocities and positions needed for computing $n\u0307i$ are evaluated at the center of the time step in accordance with the integration scheme discussed below in Sec. II C. We then compute the actual number of photons scattered by each ion, *n _{i}*, as a Poisson random number with mean $n\xafi$. Each ion receives a total momentum kick of

where $\Delta pi,absorb=ni\u210fk$ and $\Delta pi,emit$ is the recoil corresponding to *n _{i}* photons scattered in random directions with an isotropic probability distribution. To compute $\Delta pi,emit$, we generate

*n*vectors of length $\u210fk$ pointing in random directions. The recoil momentum $\Delta pi,emit$ is then obtained by adding up these vectors. The resulting velocity change in one time step is much smaller than the ion thermal velocity.

_{i}For simulating multiple cooling lasers, we find the momentum kick Eq. (13) for each laser individually, and then, we add up the results. This approach captures the salient features of laser cooling with the following approximations: first, in the case of strong saturation, $S\u22731$, the momentum kicks from multiple beams are not additive but are correlated in a more complicated way. We neglect these correlations. The other approximation in our model is that we assume that the ion motion is uniform during an excited state lifetime. Specifically, in treating the laser scattering, we neglect the change in an ion's velocity between absorption and re-emission. For the motion of ions in the Penning trap, the most rapid change of the ion velocity is due to the magnetic field. We can characterize the validity of the assumption of uniformity of motion by the dimensionless parameter $\eta =\omega B/2\pi \gamma 0$ being much small compared to 1. The parameter $\omega B=Bzq/m$ is the Larmor precession frequency. For the trap and ion parameters considered in this paper, we have $\eta \u223c0.05$. We also assume that the cooling laser intensity is approximately constant during an excited state lifetime. For the typical parameters considered in this publication, this approximation holds to a higher degree of accuracy than the assumption of a constant ion velocity during an excited state lifetime.

### C. Time integration of the equations of motion

The numerical integration of the ion motion is performed in the lab frame of reference. A time splitting technique is used to numerically integrate the equations of motion for the ions. This is analogous to the Buneman and Boris algorithms^{33–35} with an exact matrix rotation for the $v\xd7B$ force term. To advance the positions and velocities of the ions, ${xi,vi}$, from time *t* to *t* + Δ*t*, we apply the following time evolution operator:

where *U*(Δ*t*) is given by

$U0(\Delta t/2)$ is the time evolution operator corresponding to *H*_{0} [see Eq. (2)] which advances the state of the ions for a time interval of duration Δ*t*/2. Since *H*_{0} contains just the kinetic energy of the ions and the Lorentz force due to the axial magnetic field, the motion generated by *U*_{0} is given by the Larmor precession

where

The time evolution operator $Ukick(t+\Delta t/2;\Delta t)$ captures the forces due to the electrostatic potential [see Eq. (3)] and the radiation pressure forces due to laser cooling (Sec. II B). This operator is time dependent. We evaluate it at the midpoint of the time interval, *t* + Δ*t*/2. The operator *U*_{kick} changes the velocities of the particles only

where

### D. Convergence with respect to the time step

To verify the precision of our numerical integration procedure, we check that the simulations converge with the expected quadratic rate as the time step size is reduced. To evaluate the convergence, we initialize our simulation with an equilibrium configuration for 127 ions shown in Fig. 1. This is the lowest energy state equilibrium (zero temperature 2D crystal) in the rotating frame generated using the procedure previously discussed by Wang *et al.*^{12} The simulation is in the laboratory frame and so tracks the full dynamics, including the cyclotron and magnetron motion, and so, maintaining the 2D crystal equilibrium is a good test of the numerics. We study the numerical convergence without Doppler laser cooling. Here and in the rest of the paper, unless indicated otherwise, we use typical parameters for the Penning trap at the National Institute of Standards and Technology (NIST)^{10,13,36} with a strong homogeneous magnetic field of *B _{z}* = 4.4588 T, a trap rotation frequency of $\omega trap=2\pi \xd7$ 180 kHz, end cap voltages yielding a confining potential of

*k*= 9.21 MV/m

_{z}^{2}, and a 1 V rotating wall potential yielding

*δ*= 3.5 × 10

^{−4}. For

^{9}Be

^{+}ions, which is assumed throughout this paper, the cyclotron frequency is $qB/m=2\pi \xd77.596\u2009MHz$ and the axial confining frequency parallel to the magnetic field is 2

*π*× 1.58 MHz.

Starting from this initial 2D crystal configuration, we integrate the equations of motion for a total duration, *τ*, using different time step sizes, Δ*t*, to find the final positions $xi(\tau ;\Delta t)$. For each time step size, we compare the final solution with a reference solution, $xiref(\tau )=xi(\tau ,2\xd710\u221210\u2009s)$, computed using a time step size of Δ*t* = 2 × 10^{−10} s. We compute the error $\Delta x(\tau ;\Delta t)$ as the average Euclidean distance of the positions from the reference solution

Figure 2 shows the integration errors for an integration time of *τ* = 10 *μ*s. The error decreases quadratically for time step sizes below approximately $1.0\xd710\u22127\u2009s$. At a time step size of Δ*t* = 1 ns, the mean position error of the ions is $\Delta x(1\u2009ns)\u22481\u2009nm$. During this time interval, the ions move on the order of 1 mm, i.e., the ion motion is integrated with a relative error on the order of 1.0 × 10^{−6}. For time step sizes greater than $\Delta t\u22481.0\xd710\u22127\u2009s$, resonances occur where the integrator step size matches the period of one of the vibrational eigenmodes of the system. At these resonances, the time integration errors can grow large. The results presented in this paper were obtained with Δ*t* = 1 ns.

## III. COMPARISON WITH LINEAR THEORY NEAR ZERO TEMPERATURE

We now turn to an analysis of the vibrational eigenmodes of the 2D ion crystal. To find the eigenmodes, we initialize the ions in a low energy spatial configuration found by minimizing the potential energy in the rotating frame. We subject the ions to cooling lasers with a geometry typical for the NIST experiment.^{10,13,28,36} The details of the setup of the cooling lasers are described in Sec. IV. The cooling time for both in-plane and out-of-plane degrees of freedom is on the order of 1 ms. After 5 ms of laser cooling, the ions have settled into an equilibrium state. We then turn off the cooling lasers and record the trajectories of the ions $xi(tj)$ at discrete times $tj=j\Delta t$ for $j=0,1,\u2026,Nsample\u22121$.

From the trajectories, we compute power spectra of the ion motion by means of

In this equation, $z\u0303i$ is the Fourier transform of the *z* coordinate of ion *i*

with $\tau =Nsample\xb7\Delta t$ the total integration time of the simulation. In terms of the discretely sampled trajectories, the Fourier transform can be approximated by a discrete Fourier transform, which we evaluate with the aid of a Fast Fourier Transform numerically

Here, we have introduced the frequency resolution $\Delta \omega =2\pi /\tau $ and $l=0,1,\u2026,Nsample\u22121$. Frequencies exceeding the Nyquist frequency, $l\Delta \omega >\pi /\Delta t$, wrap to negative frequencies.

We compute spectra *P _{x}* and

*P*for the in-plane degrees of freedom similarly with a minor caveat: to eliminate the large coherent component due to the uniform rotation in the trap, we use the coordinates in the rotating frame,

_{y}*x*and

_{R}*y*.

_{R}### A. Out-of-plane eigenmodes

Figure 3(a) shows *P _{z}* for 7 ions with a sampling period of Δ

*t*= 0.25

*μ*s and a total integration time of

*τ*= 50

*μ*s corresponding to a frequency resolution of $\Delta \omega /(2\pi )=20\u2009Hz$. We superimpose the mode frequencies found by linearizing the equations of motion around the ground state as gray vertical lines. For the linear analysis (gray lines), we use a code developed by Wang

*et al.*

^{12}with minor modifications. We observed that the resonance frequencies obtained with our molecular dynamics simulations agree well with the modes from the linear theory. Several pairs of modes are nearly degenerate. For example, the degeneracy of the two tilt modes around 1.55 MHz is lifted only by the weak perturbation of the rotating wall potential. The out-of-plane eigenmodes are frequently referred to as drumhead modes because they resemble the vibrations of a membrane with open boundary conditions on the edge of the membrane.

Figure 3(b) shows *P _{z}* for 127 ions with otherwise identical parameters. Again, the frequencies of the highest frequency modes agree well with the zero temperature linear mode analysis. However, the frequencies of several of the lower frequency, shorter wavelength modes deviate from the linear frequencies. The reason for these shifts is that, at finite temperature, defects form in the ion crystal, especially for larger numbers of ions. Due to their shorter wavelength, the lower frequency modes more sensitively depend on the local crystal structure and therefore are perturbed by the crystal defects. In Sec. IV, we will discuss how the power spectrum in Fig. 3(b) can be used to estimate the temperature of the drumhead modes.

In addition to frequency spectra, the axial vibrational mode patterns can be analyzed.^{12} We obtain the eigenfunction shapes using a notch filter of width Δ*ω* about the mode frequency of interest. We then inverse the Fourier transform $z\u0303i(\omega )$ and plot the normalized displacement as a function of time. Figure 4 shows a snapshot of the first four (starting from the highest frequency) vibrational modes for 127 ions. Figure 4(i) shows the center-of-mass mode where all ions move together in the *z*-direction. Figures 4(ii) and 4(iii) show the next two highest frequency modes that are tilt modes that would be degenerate except for the weak structural anisotropy produced by the rotating wall potential. The fourth mode, shown in Fig. 4(iv), has two linear nodes. Figure 4 shows similar results as obtained in the eigenfunction plots presented in Ref. 12 (Fig. 11).

### B. In-plane spectrum

The spectrum of in-plane modes can be decomposed into two qualitatively different types of modes. There is a high frequency branch of modes near the cyclotron frequency and a low frequency branch near zero frequency. The high frequency modes are commonly referred to as cyclotron modes, and the low frequency modes are referred to as *E *×* B* modes. The spectra *P _{x}* and

*P*are very similar. Therefore, we consider just

_{y}*P*for simplicity.

_{x}Figure 5(a) shows the two narrow bands of modes for 7 ions together on the same plot. The mode frequencies predicted by the linear theory are again superimposed as gray vertical lines. Figures 5(b) and 5(c) show the cyclotron and *E *×* B* modes separately in more detail. The left most peak near zero frequency shown in Fig. 5(c), which is relatively higher than the other peaks, is associated with the slow periodic deformation of the elliptical ion crystal due to the confining potential of the rotating wall. This mode, sometimes called a rocking mode, is a zero frequency mode in the absence of a rotating wall. The presence of a wall potential leads to a small nonzero frequency for this mode. Figures 6(a) and 6(b) show the two branches of in-plane modes for the 127 ion crystal. Again, we observe good agreement between the zero-temperature linear analysis and the finite temperature direct numerical simulations.

## IV. DOPPLER COOLING GEOMETRY AND SIMULATION RESULTS

We now turn to an analysis of Doppler cooling of the ion motion. The numerical implementation of Doppler laser cooling was discussed in Sec. II B. Doppler laser cooling with a rotating wall is discussed in detail in Ref. 28. There, it was shown in a fluid model that assumed thermal equilibrium that low ion temperatures close to the Doppler laser cooling limit were possible over a range of experimental conditions. This is in contrast to a Doppler laser cooling model without a rotating wall,^{27} where the equilibrium temperature depended sensitively on the details of the setup. Here, we discuss in more detail the geometry of the laser beams and the achievement of a Doppler cooled ultracold steady-state 2D crystal of 127 ions using direct numerical simulation. Measurements of the average temperatures of the ion motion in the planar (perpendicular to the B-field) and axial (parallel to the B-field) directions are presented. In addition, mode resolved temperatures in the axial direction are investigated.

### A. Doppler cooling in a rotating wall Penning trap

We consider a cooling laser configuration similar to that employed in a typical NIST Penning trap experimental setup.^{10,13,28,36} We assume that the out-of-plane motion is cooled by two counter-propagating cooling beams along the positive and negative *z* directions, i.e., parallel and antiparallel to the trap magnetic field. These two beams have equal intensity, and they are detuned by $\Delta \u2225=\u2212\gamma 0/2$ relative to the atomic transition. The waist of the beams is assumed to be larger than the size of the ion crystal so that the spatial variation of the laser intensity can be neglected.

The in-plane motion is cooled by an additional laser directed along the *x* axis, i.e., perpendicular to the trap magnetic field. The arrangement of the planar cooling beam is illustrated in Fig. 7. This beam has a width *W _{y}* and is displaced by a distance

*d*from the center of the trap. To adjust for the rotational motion of the ions at the planar beam location, the beam is detuned by $\Delta \u22a5=\u2212\gamma 0/2+k\omega Rd$ from the atomic transition.

The in-plane and out-of-plane degrees of freedom are only weakly coupled to one another. Therefore, they are not typically in thermal equilibrium. To study the temperature anisotropy, we introduce parallel and perpendicular temperatures

To simulate the Doppler cooling process, we initialize the ions in a finite temperature state by setting their positions to the 2D crystal equilibrium, with random velocities drawn from a Maxwellian (Gaussian) distribution at a given starting temperature. We use a perpendicular cooling beam with a peak saturation intensity of *S*_{0} = 1 and a width of *W _{y}* = 5

*μ*m displaced by

*d*=

*5*

*μ*m from the trap center. The axial cooling beams have a uniform saturation intensity of

*S*

_{0}= 0.005. We deliberately choose a much smaller axial cooling beam intensity because cooling of the out-of-plane motion is very efficient since all ions interact continuously with the parallel beams. Larger axial cooling beam intensities would heat the in-plane degrees of freedom due to the recoil of the scattered photons, overpowering the planar cooling.

With all the cooling beams on, the system is left to evolve for 10 ms. The time histories of the in-plane and out-of-plane temperatures during the simulation are shown in Fig. 8. For these simulation parameters, both in-plane and out-of-plane temperatures relax to a steady state on a time scale on the order of 1 ms.

### B. Measured steady-state axial and planar temperatures

We now present steady-state temperature results. The axial temperature $T\u2225$ can be well diagnosed in the experiments.^{1} Also, the theoretical Doppler cooling limit^{27,28} gives an estimate for $T\u2225$, which agrees fairly well with temperatures achieved in the NIST experiment.^{10,13,36} The theoretical Doppler cooling limit for motion parallel to the magnetic field is

where $\gamma 0=2\pi \xb718\u2009MHz$ is the atomic transition linewidth for ^{9}Be^{+}. This is the cooling limit predicted for laser cooling along one axis (in this case parallel to the magnetic field) with minimal recoil heating from the perpendicular laser beam.^{26} This limit assumes isotropic scattering, which was also assumed in the numerical simulation of Doppler laser cooling. We estimate the steady-state $T\u2225$ by averaging the results with the lasers on continuously over the final millisecond, the time interval shown in the insets of Fig. 8. We find steady state temperatures in Fig. 8 of $T\u2225\u22480.37\u2009mK$ and $T\u22a5\u22480.92\u2009mK$. The lowest $T\u2225$ was obtained in a different simulation with *d *=* *0 (no planar beam offset) and was measured to be $0.28\u2009mK\xb10.04\u2009mK$, in good agreement with the Doppler cooling limit in Eq. (25). Also, it is important to note that the simulation shows fluctuations at the axial temperature to be approximately ±15%, which are produced by the interaction with the cooling laser beams and the dynamic exchange between potential and kinetic energies.

We can use the simulation to measure the thermal energy in each axial mode and investigate if an equipartition assumption is valid. Because the system is at very low temperature, nonlinear couplings are very weak. Additionally, the simulations are for a finite time. Therefore, one might not expect the system to be in equipartition. For an axial mode with frequency *ω _{m}*, we can obtain the mode temperature $T\u2225(\omega m)$ from the spectrum

where the integral is over a single resonance, and we evaluate the integral using the trapezoidal rule. Doppler cooling is turned off during the mode temperature recording time *τ*.

Table I shows the axial mode temperatures, relative to the theoretical parallel Doppler laser cooling limit, for the eight highest frequency axial modes for 127 ions. The axial mode temperatures shown in Table I are also obtained again with *d *=* *0, which gives the lowest parallel temperature steady state. The values are in the range of the Doppler cooling limit given the ±15% fluctuation, except for the two tilt modes (mode 1 and 2). More simulations and further analysis are required to determine the uncertainty of mode-resolved temperatures. Because this may be larger than the ±15% global temperature measurement uncertainty, it is premature to conclude that the tilt modes are colder than the Doppler laser cooling limit. However, Table I indicates the type of detailed information that should be possible with the reported simulations. Higher mode numbers are harder to diagnose using Eq. (26) due to the “grassy” spectrum in the range of 1.1 MHz to 1.5 MHz seen in Fig. 3(b). The spectral peaks become very close to each other with the increasing mode number, making them difficult to differentiate. Additionally, the lowest $T\u22a5$ is obtained by varying *W _{y}* and

*d*and found to be $T\u22a5\u2248$ 1 mK with $Wy=d=5\u2009\mu m$. Such small

*W*and

_{y}*d*, compared to the radius of the crystal

*R*= 70

_{c}*μ*m, allow the perpendicular beam to generate a modest torque, which is straight forwardly balanced by the rotating wall.

^{28}Additionally, the different $T\u22a5$ and $T\u2225$ demonstrate the temperature anisotropy and weak coupling between axial and planar directions.

m
. | ω (MHz)
. _{m} | $T\u2225(\omega m)/T\u2225,Doppler$ . |
---|---|---|

0 | 1.580532 | 1.05 |

1 | 1.553344 | 0.58 |

2 | 1.552781 | 0.31 |

3 | 1.532620 | 0.81 |

4 | 1.532467 | 0.81 |

5 | 1.515915 | 0.85 |

6 | 1.515464 | 0.94 |

7 | 1.511197 | 1.01 |

m
. | ω (MHz)
. _{m} | $T\u2225(\omega m)/T\u2225,Doppler$ . |
---|---|---|

0 | 1.580532 | 1.05 |

1 | 1.553344 | 0.58 |

2 | 1.552781 | 0.31 |

3 | 1.532620 | 0.81 |

4 | 1.532467 | 0.81 |

5 | 1.515915 | 0.85 |

6 | 1.515464 | 0.94 |

7 | 1.511197 | 1.01 |

## V. SUMMARY

We have described in detail a direct numerical simulation model of ultracold ions in a Penning trap with a rotating wall. The simulation includes a Doppler cooling model based on the microscopic physics of resonance fluorescence with a realistic laser beam geometry. Due to the low temperatures of interest, very good agreement is obtained with a linear eigenmode analysis, in both the in-plane and out-of-plane directions. The advantages of a direct numerical simulation are that weak nonlinear coupling between modes is accurately modeled and the laser cooling/heating processes. Additionally, it is difficult to diagnose in-plane motion experimentally, so simulations can help in better interpreting and understanding this physics.

Such a model is very useful for understanding the parametric dependences of the heating and cooling processes in the experiment and the relevant time scales. The Doppler cooling process was successfully incorporated in a molecular dynamics simulation for the first time, and the results were compared with the existing theory. Axial and planar temperatures were within the range of experimental results, and more work is needed to make detailed comparisons with the experiment. We varied the waist *W _{y}* and offset

*d*of the perpendicular laser beam to obtain the lowest planar temperature with the laser offset equal to the beam width $d=Wy=5\u2009\mu m\u226aRc\u224870\u2009\mu m$. The fluid equilibrium model

^{28}predicts low in-plane temperatures for small

*W*, with a very weak dependence on

_{y}*d*if $\Delta \u22a5=\u2212\gamma /2+k\omega Rd$. This prediction will be interesting to investigate because increasing

*d*produces a large laser torque which, when balanced with a similar large torque from the rotating wall, produces shear and potential instabilities in the crystal.

There are areas where the simulation model could be improved further. The Doppler cooling model treats photon-ion interactions instantaneously, neglecting the transition time. We also neglect infrequent interactions with background neutrals and associated ion loss, as well as impurities and trap error fields. These features of a more realistic experiment will be addressed in future work.

## ACKNOWLEDGMENTS

We thank Murray Holland and Athreya Shanker, JILA, University of Colorado, and Elena Jordan, NIST, for useful discussions. This work was supported in part by the U.S. Department of Energy under Award No. DE-FG02-08ER54954. This manuscript is a contribution of NIST and not subject to U.S. copyright.