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.

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 ions19 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 arrays20 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 protocols22 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 shifts31,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 work12,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.

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.

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

H=H0+i=1Nqiφ(xi,t),
(1)

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

H0=i=1N12mi(piqiA(xi))2,
(2)

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=×A is parallel to the z axis. With that choice of coordinate system, we may choose A=yBzx̂ with x̂=x/x the unit vector along x. In Eq. (2), N is the number of ions in the trap and mi and qi are the mass and charge of ion i. The electrostatic potential

φ(xi,t)=φtrap(xi)+φwall(xi,t)+ji18πε0qj|xixj|
(3)

contains the potential φtrap for the electrodes in the Penning trap, the rotating wall potential φ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 φtrap(x) as

φtrap(x)=14kz(2z2x2y2)
(4)

and the time-dependent rotating wall potential φwall(x,t) as

φwall(x,t)=12kzδ(x2+y2)cos[2(θ+ωRt)],
(5)

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

[xRyR]=[cos(ωRt)sin(ωRt)sin(ωRt)cos(ωRt)][xy].
(6)

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

φtrap,R(xR)+φwall,R(xR)=12kzz212(kxxR2+kyyR2),
(7)

where

kx=kz(12δ),ky=kz(12+δ).
(8)

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, HR, is

HR=i=1N12mi(pRiqiA(xRi))2+i=1Nqi(φR(xRi)+ϕwall,R(xRi))+i=1N12(qiBzωRmiωR2)rRi2+i=1Nqiji18πε0qj|xRixRj|,
(9)

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.

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 by27 

ṅ(x,v)=S(x)γ0(γ0/2)2(γ0/2)2(1+2S(x))+Δ2(v),
(10)

where γ0 is the natural linewidth of the atom transition (in radians per second), S(x) is the saturation parameter, and Δ(v)=Δ0k·v is the detuning of the atomic transition from the laser frequency including the first order Doppler shift with Δ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

S(x)=S0eρ2/Wy2,
(11)

where ρ is the distance of the atom from the axis of the beam and Wy 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,

n¯i=ṅiΔt.
(12)

The velocities and positions needed for computing ṅi 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, ni, as a Poisson random number with mean n¯i. Each ion receives a total momentum kick of

ΔpiLaser=Δpi,absorb+Δpi,emit,
(13)

where Δpi,absorb=nik and Δpi,emit is the recoil corresponding to ni photons scattered in random directions with an isotropic probability distribution. To compute Δpi,emit, we generate ni vectors of length k pointing in random directions. The recoil momentum Δ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.

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, S1, 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 η=ωB/2πγ0 being much small compared to 1. The parameter ωB=Bzq/m is the Larmor precession frequency. For the trap and ion parameters considered in this paper, we have η0.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.

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 algorithms33–35 with an exact matrix rotation for the v×B 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:

{xi(t+Δt),vi(t+Δt)}=U(Δt){xi(t),vi(t)},
(14)

where Ut) is given by

U(Δt)=U0(Δt/2)Ukick(t+Δt/2;Δt)U0(Δt/2).
(15)

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

xi(t+Δt/2)=xi(t)+[sc10(c1)s000ωBΔt/2]vi(t)ωB,vi(t+Δt/2)=[cs0sc0000]vi(t),
(16)

where

s=sin(ωBΔt/2)andc=cos(ωBΔt/2).
(17)

The time evolution operator Ukick(t+Δt/2;Δ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 Ukick changes the velocities of the particles only

vi(t+Δt)=vi(t)+Δt(qi/mi)E(t+Δt/2,xi(t+Δt/2))+ΔpiLaser/mi,

where

E(t+Δt/2,xi(t+Δt/2))=φ(xi(t+Δt/2)).

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 Bz = 4.4588 T, a trap rotation frequency of ωtrap=2π× 180 kHz, end cap voltages yielding a confining potential of kz = 9.21 MV/m2, and a 1 V rotating wall potential yielding δ = 3.5 × 10−4. For 9Be+ ions, which is assumed throughout this paper, the cyclotron frequency is qB/m=2π×7.596MHz and the axial confining frequency parallel to the magnetic field is 2π × 1.58 MHz.

FIG. 1.

Top view of the steady state configuration of 127 ions in a Penning trap used for the convergence study.

FIG. 1.

Top view of the steady state configuration of 127 ions in a Penning trap used for the convergence study.

Close modal

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(τ;Δt). For each time step size, we compare the final solution with a reference solution, xiref(τ)=xi(τ,2×1010s), computed using a time step size of Δt = 2 × 10−10 s. We compute the error Δx(τ;Δt) as the average Euclidean distance of the positions from the reference solution

Δx(τ;Δt)=N1i=1N(xi(τ;Δt)xiref(τ))2.
(18)

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×107s. At a time step size of Δt = 1 ns, the mean position error of the ions is Δx(1ns)1nm. 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 Δt1.0×107s, 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.

FIG. 2.

Time integration errors as a function of time step size for a total integration time of τ = 10 μs for a typical Penning trap simulation. The dashed orange line is quadratic for the orientation. See text for details.

FIG. 2.

Time integration errors as a function of time step size for a total integration time of τ = 10 μs for a typical Penning trap simulation. The dashed orange line is quadratic for the orientation. See text for details.

Close modal

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Δt for j=0,1,,Nsample1.

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

Pz(ω)=i=1N|z̃i(ω)|2+i=1N|z̃i(ω)|2.
(19)

In this equation, z̃i is the Fourier transform of the z coordinate of ion i

z̃i(ω)=1τ0τeiωtzi(t)dt,
(20)

with τ=Nsample·Δ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

z̃(lΔω)=1τ0τeilΔωtz(t)dt,
(21)
1NsnNsei2πln/Nsz(nΔt).
(22)

Here, we have introduced the frequency resolution Δω=2π/τ and l=0,1,,Nsample1. Frequencies exceeding the Nyquist frequency, lΔω>π/Δt, wrap to negative frequencies.

We compute spectra Px and Py 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, xR and yR.

Figure 3(a) shows Pz 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 Δω/(2π)=20Hz. 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.

FIG. 3.

Spectra Pz(ω) of out-of-plane motion for 7 ions (a) and 127 ions (b). Some modes are nearly degenerate, and individual peaks are not discernible.

FIG. 3.

Spectra Pz(ω) of out-of-plane motion for 7 ions (a) and 127 ions (b). Some modes are nearly degenerate, and individual peaks are not discernible.

Close modal

Figure 3(b) shows Pz 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̃i(ω) 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).

FIG. 4.

Vibrational mode patterns (eigenfunctions) for the first 4 out-of-plane modes using 127 ions. Each dot represents the ion position, zi. The color represents the scale value of the normalized out-of-plane motion of the particular ion, as described in Sec. III A.

FIG. 4.

Vibrational mode patterns (eigenfunctions) for the first 4 out-of-plane modes using 127 ions. Each dot represents the ion position, zi. The color represents the scale value of the normalized out-of-plane motion of the particular ion, as described in Sec. III A.

Close modal

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 Px and Py are very similar. Therefore, we consider just Px for simplicity.

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.

FIG. 5.

Spectra of the in-plane motion of a 7 ion crystal, Px(ω). The top figure (a) shows the entire frequency range with the E × B modes near zero frequency and the cyclotron modes near ωB=eB/m=(2π)×7.596 MHz. (b) and (c) Closeups of the cyclotron (upper branch) and E × B (lower branch) modes, respectively. For this 7 ion case, we use a different trap axial frequency, to better display the cyclotron modes.

FIG. 5.

Spectra of the in-plane motion of a 7 ion crystal, Px(ω). The top figure (a) shows the entire frequency range with the E × B modes near zero frequency and the cyclotron modes near ωB=eB/m=(2π)×7.596 MHz. (b) and (c) Closeups of the cyclotron (upper branch) and E × B (lower branch) modes, respectively. For this 7 ion case, we use a different trap axial frequency, to better display the cyclotron modes.

Close modal
FIG. 6.

Two branches of the in-plane motion spectrum for a 127 ion crystal. (a) and (b) Closeups of the E × B (lower branch) and cyclotron (upper branch) modes, respectively.

FIG. 6.

Two branches of the in-plane motion spectrum for a 127 ion crystal. (a) and (b) Closeups of the E × B (lower branch) and cyclotron (upper branch) modes, respectively.

Close modal

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.

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 Δ=γ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 Wy 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 Δ=γ0/2+kωRd from the atomic transition.

FIG. 7.

Schematic of the perpendicular laser cooling geometry. The black circle represents the approximate area occupied by 127 ions, which have a collective rotating motion at frequency ωR. Two parallel beams (not displayed here) are perpendicular to the X-Y plane. The perpendicular laser shown in red has a Gaussian intensity profile with a width of Wy and a displacement of d from the trap center, and it is directed along the x axis.

FIG. 7.

Schematic of the perpendicular laser cooling geometry. The black circle represents the approximate area occupied by 127 ions, which have a collective rotating motion at frequency ωR. Two parallel beams (not displayed here) are perpendicular to the X-Y plane. The perpendicular laser shown in red has a Gaussian intensity profile with a width of Wy and a displacement of d from the trap center, and it is directed along the x axis.

Close modal

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

T=mBeNkBi=1Nvz,i2,
(23)
T=mBe2NkBi=1N(vx,i2+vy,i2).
(24)

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 S0 = 1 and a width of Wy = 5 μm displaced by d =5 μm from the trap center. The axial cooling beams have a uniform saturation intensity of S0 = 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.

FIG. 8.

A simulation of Doppler cooling starting with an initial temperature T =20 mK. The blue line represents the axial motion, and the green line represents the planar motion. The ion crystal achieves a steady-state temperature in qualitative agreement with the experiment. The two insets show a zoom-in of the axial and planar temperatures during the last millisecond of cooling. For these parameters (see the text), the steady state axial temperature is T=0.37mK, and the planar temperature is T=0.92mK.

FIG. 8.

A simulation of Doppler cooling starting with an initial temperature T =20 mK. The blue line represents the axial motion, and the green line represents the planar motion. The ion crystal achieves a steady-state temperature in qualitative agreement with the experiment. The two insets show a zoom-in of the axial and planar temperatures during the last millisecond of cooling. For these parameters (see the text), the steady state axial temperature is T=0.37mK, and the planar temperature is T=0.92mK.

Close modal

We now present steady-state temperature results. The axial temperature T can be well diagnosed in the experiments.1 Also, the theoretical Doppler cooling limit27,28 gives an estimate for T, 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

T,Doppler=13γ0kB=0.29mK,
(25)

where γ0=2π·18MHz is the atomic transition linewidth for 9Be+. 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 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 T0.37mK and T0.92mK. The lowest T was obtained in a different simulation with d =0 (no planar beam offset) and was measured to be 0.28mK±0.04mK, 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(ωm) from the spectrum

T(ωm)=mωm2τπres.dωPz(ω),
(26)

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 is obtained by varying Wy and d and found to be T 1 mK with Wy=d=5μm. Such small Wy and d, compared to the radius of the crystal Rc = 70 μm, allow the perpendicular beam to generate a modest torque, which is straight forwardly balanced by the rotating wall.28 Additionally, the different T and T demonstrate the temperature anisotropy and weak coupling between axial and planar directions.

TABLE I.

Mode temperatures of the first eight (starting with m =0) axial eigenmodes.

mωm (MHz)T(ωm)/T,Doppler
1.580532 1.05 
1.553344 0.58 
1.552781 0.31 
1.532620 0.81 
1.532467 0.81 
1.515915 0.85 
1.515464 0.94 
1.511197 1.01 
mωm (MHz)T(ωm)/T,Doppler
1.580532 1.05 
1.553344 0.58 
1.552781 0.31 
1.532620 0.81 
1.532467 0.81 
1.515915 0.85 
1.515464 0.94 
1.511197 1.01 

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 Wy 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μmRc70μm. The fluid equilibrium model28 predicts low in-plane temperatures for small Wy, with a very weak dependence on d if Δ=γ/2+kω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.

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.

1.
E.
Jordan
,
K. A.
Gilmore
,
A.
Shankar
,
A.
Safavi-Naini
,
J. G.
Bohnet
,
M. J.
Holland
, and
J. J.
Bollinger
,
Phys. Rev. Lett.
122
,
053603
(
2019
).
2.
J. N.
Tan
,
J. J.
Bollinger
,
B.
Jelenkovic
, and
D.
Wineland
,
Phys. Rev. Lett.
75
,
4198
(
1995
).
3.
S.
Mavadia
,
J. F.
Goodwin
,
G.
Stutter
,
S.
Bharadia
,
D. R.
Crick
,
D. M.
Segal
, and
R. C.
Thompson
,
Nat. Commun.
4
,
2571
(
2013
).
4.
H.
Ball
,
C. D.
Marciniak
,
R. N.
Wolf
,
A. T.-H.
Hung
,
K.
Pyka
, and
M. J.
Biercuk
,
Rev. Sci. Instr.
90
,
053103
(
2019
).
5.
L.
Gruber
,
J.
Holder
,
J.
Steiger
,
B.
Beck
,
H.
DeWitt
,
J.
Glassman
,
J.
McDonald
,
D.
Church
, and
D.
Schneider
,
Phys. Rev. Lett.
86
,
636
(
2001
).
6.
N.
Shiga
,
W. M.
Itano
, and
J. J.
Bollinger
,
Phys. Rev. A
84
,
012510
(
2011
).
7.
Z.
Andelkovic
,
R.
Cazan
,
W.
Nörtershäuser
,
S.
Bharadia
,
D.
Segal
,
R.
Thompson
,
R.
Jöhren
,
J.
Vollbrecht
,
V.
Hannen
, and
M.
Vogel
,
Phys. Rev. A
87
,
033423
(
2013
).
8.
G.
Werth
,
S.
Sturm
, and
K.
Blaum
, in
Advances in Atomic, Molecular, and Optical Physics
(
Elsevier
,
2018
), Vol.
67
, pp.
257
296
.
9.
G.
Stutter
,
P.
Hrmo
,
V.
Jarlaud
,
M.
Joshi
,
J.
Goodwin
, and
R.
Thompson
,
J. Mod. Opt.
65
,
549
(
2018
).
10.
J. G.
Bohnet
,
B. C.
Sawyer
,
J. W.
Britton
,
M. L.
Wall
,
A. M.
Rey
,
M.
Foss-Feig
, and
J. J.
Bollinger
,
Sci.
352
,
1297
(
2016
).
11.
K. A.
Gilmore
,
J. G.
Bohnet
,
B. C.
Sawyer
,
J. W.
Britton
, and
J. J.
Bollinger
,
Phys. Rev. Lett.
118
,
263602
(
2017
).
12.
C.-C. J.
Wang
,
A. C.
Keith
, and
J.
Freericks
,
Phys. Rev. A
87
,
013422
(
2013
).
13.
J. W.
Britton
,
B. C.
Sawyer
,
A. C.
Keith
,
C.-C. J.
Wang
,
J. K.
Freericks
,
H.
Uys
,
M. J.
Biercuk
, and
J. J.
Bollinger
,
Nature
484
,
489
(
2012
).
14.
M.
Gärttner
,
J. G.
Bohnet
,
A.
Safavi-Naini
,
M. L.
Wall
,
J. J.
Bollinger
, and
A. M.
Rey
,
Nat. Phys.
13
,
781
(
2017
).
15.
D. H.
Dubin
and
T.
O'neil
,
Rev. Mod. Phys.
71
,
87
(
1999
).
16.
17.
S.
Ichimaru
,
H.
Iyetomi
, and
S.
Tanaka
,
Phys. Rep.
149
,
91
(
1987
).
18.
D.
Baiko
,
Phys. Rev. E
80
,
046405
(
2009
).
19.
R.
Blatt
and
C. F.
Roos
,
Nat. Phys.
8
,
277
(
2012
).
20.
S.
Debnath
,
N. M.
Linke
,
C.
Figgatt
,
K. A.
Landsman
,
K.
Wright
, and
C.
Monroe
,
Nature
536
,
63
(
2016
).
21.
W. S.
Bakr
,
J. I.
Gillen
,
A.
Peng
,
S.
Fölling
, and
M.
Greiner
,
Nature
462
,
74
(
2009
).
22.
C.
Kokail
,
C.
Maier
,
R.
van Bijnen
,
T.
Brydges
,
M. K.
Joshi
,
P.
Jurcevic
,
C. A.
Muschik
,
P.
Silvi
,
R.
Blatt
,
C. F.
Roos
 et al,
Nature
569
,
355
360
(
2019
).
23.
M.
Jensen
,
T.
Hasegawa
,
J. J.
Bollinger
, and
D.
Dubin
,
Phys. Rev. Lett.
94
,
025001
(
2005
).
24.
D. H.
Dubin
,
Phys. Rev. Lett.
94
,
025002
(
2005
).
25.
F.
Anderegg
,
D.
Dubin
,
T.
O'Neil
, and
C.
Driscoll
,
Phys. Rev. Lett.
102
,
185001
(
2009
).
26.
W. M.
Itano
and
D.
Wineland
,
Phys. Rev. A
25
,
35
(
1982
).
27.
W. M.
Itano
,
L.
Brewer
,
D.
Larson
, and
D. J.
Wineland
,
Phys. Rev. A
38
,
5698
(
1988
).
28.
S. B.
Torrisi
,
J. W.
Britton
,
J. G.
Bohnet
, and
J. J.
Bollinger
,
Phys. Rev. A
93
,
043421
(
2016
).
29.
M.
Asprusten
,
S.
Worthington
, and
R. C.
Thompson
,
Appl. Phys. B
114
,
157
(
2014
).
30.
A.
Shankar
,
E.
Jordan
,
K. A.
Gilmore
,
A.
Safavi-Naini
,
J. J.
Bollinger
, and
M. J.
Holland
,
Phys. Rev. A
99
,
023409
(
2019
).
31.
C.
Marquet
,
F.
Schmidt-Kaler
, and
D. F.
James
,
Appl. Phys. B
76
,
199
(
2003
).
32.
M.
McAneny
and
J.
Freericks
,
Phys. Rev. A
90
,
053405
(
2014
).
33.
O.
Buneman
,
J. Comput. Phys.
1
,
517
(
1967
).
34.
J. P.
Boris
and
R. A.
Shanny
, in
Proceedings: Fourth Conference on Numerical Simulation of Plasmas
, November 2, 3, 1970 (
Naval Research Laboratory
,
1972
), pp.
3
67
.
35.
C.
Birdsall
and
A.
Langdon
,
Plasma Physics via Computer Simulation
(
McGraw-Hill
,
New York
,
1985
), p.
62
.
36.
B. C.
Sawyer
,
J. W.
Britton
, and
J. J.
Bollinger
,
Phys. Rev. A
89
,
033408
(
2014
).