First principle simulation of ultra-cold ion crystals in a Penning trap with Doppler cooling and a rotating wall potential

TL;DR

This paper presents a classical simulation of ultra-cold ion crystals in a Penning trap, modeling Doppler cooling and rotating wall effects, and compares results with eigenmode analysis to validate the approach.

Contribution

It introduces a detailed numerical simulation method for large ion crystals, incorporating Doppler cooling and trap dynamics, advancing understanding of their behavior.

Findings

Simulation results agree with eigenmode analysis in frequency and mode structure.

Doppler cooling limits axial temperature to expected values.

The approach provides insights into large ion crystal dynamics for quantum applications.

Abstract

A direct numerical simulation of many interacting ions in a Penning trap with a rotating wall is presented. The ion dynamics is modelled classically. Both axial and planar Doppler laser cooling are modeled using stochastic momentum impulses based on two-level atomic scattering rates. The plasmas being modeled are ultra-cold two-dimensional crystals made up of 100's of ions. We compare Doppler cooled results directly to a previous linear eigenmodes analysis. Agreement in both frequency and mode structure are 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.

Figures

Tables

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

$m$	${\omega }_m$	$T_{\parallel }\left({\omega }_m\right)/T_{\parallel ,\mathrm{Doppler}}$
0	$\mathrm{1.580\hspace{0.17em}532}\mathrm{MHz}$	$1.05$
1	$\mathrm{1.553\hspace{0.17em}344}\mathrm{MHz}$	$0.58$
2	$\mathrm{1.552\hspace{0.17em}781}\mathrm{MHz}$	$0.31$
3	$\mathrm{1.532\hspace{0.17em}620}\mathrm{MHz}$	$0.81$
4	$\mathrm{1.532\hspace{0.17em}467}\mathrm{MHz}$	$0.81$
5	$\mathrm{1.515\hspace{0.17em}915}\mathrm{MHz}$	$0.85$
6	$\mathrm{1.515\hspace{0.17em}464}\mathrm{MHz}$	$0.94$
7	$\mathrm{1.511\hspace{0.17em}197}\mathrm{MHz}$	$1.01$

Equations

$$H=H_0+\sum _{i=1}^N{q}_i\phi ({\mathbf{x}}_i,t),$$

$$H_0=\sum _{i=1}^N\frac{1}{2m_i}{\left({\mathbf{p}}_i-q_i\mathbf{A}({\mathbf{x}}_i)\right)}^2,$$

$$\phi ({\mathbf{x}}_i,t)={\phi }_{\mathrm{trap}}({\mathbf{x}}_i)+{\phi }_{\mathrm{wall}}({\mathbf{x}}_i,t)+\sum _{j\ne i}\frac{1}{8\pi {\epsilon }_0}\frac{q_j}{\left|{\mathbf{x}}_i-{\mathbf{x}}_j\right|}$$

$${\phi }_{\mathrm{trap}}(\mathbf{x})=\frac{1}{4}{k}_z\left(2z^2-x^2-y^2\right),$$

$${\phi }_{\mathrm{wall}}(\mathbf{x},t)=\frac{1}{2}{k}_z\delta \left(x^2+y^2\right)\cos \left[2\left(\theta +{\omega }_{R}t\right)\right],$$

$$\left[\begin{array}[]{c}x_{R}\\ y_{R}\end{array}\right]=\left[\begin{array}[]{cc}\cos(\omega_{R}t)&-\sin(\omega_{R}t)\\ \sin(\omega_{R}t)&\cos(\omega_{R}t)\end{array}\right]\left[\begin{array}[]{c}x\\ y\end{array}\right]\;.$$

$${\phi }_{\mathrm{trap},\mathrm{R}}({\mathbf{x}}_R)+{\phi }_{\mathrm{wall},\mathrm{R}}({\mathbf{x}}_R)=\frac{1}{2}{k}_z{z}^2-\frac{1}{2}\left(k_x{x}_R^2+k_y{y}_R^2\right)$$

$$k_x=k_z(\frac{1}{2}-\delta ),k_y=k_z(\frac{1}{2}+\delta ).$$

$$H_R=$$

$$\sum _{i=1}^N\frac{1}{2}\left(q_i{B}_z{\omega }_R-m_i{\omega }_R^2\right)r_{Ri}^2+\sum _{i=1}^N{q}_i\sum _{j\ne i}\frac{1}{8\pi {\epsilon }_0}\frac{q_j}{\left|{\mathbf{x}}_{Ri}-{\mathbf{x}}_{Rj}\right|}$$

$$\dot{n}(\mathbf{x},\mathbf{v})=S(\mathbf{x}){\gamma }_0\frac{({\gamma }_0/2{)}^2}{({\gamma }_0/2{)}^2(1+2S(\mathbf{x}))+{\Delta }^2(\mathbf{v})},$$

$$S(\mathbf{x})=S_0{e}^{-{\rho }^2/W_y^2},$$

$${\bar{n}}_i={\dot{n}}_i\Delta t.$$

$$\Delta {\mathbf{p}}_i^{\mathrm{Laser}}=\Delta {\mathbf{p}}_{i,\mathrm{absorb}}+\Delta {\mathbf{p}}_{i,\mathrm{emit}},$$

$$\{{\mathbf{x}}_i(t+\Delta t),{\mathbf{v}}_i(t+\Delta t)\}=U(\Delta t)\{{\mathbf{x}}_i(t),{\mathbf{v}}_i(t)\},$$

$$U(\Delta t)=U_0(\Delta t/2)U_{\mathrm{kick}}(t+\Delta t/2;\Delta t)U_0(\Delta t/2).$$

$$\displaystyle\mathbf{x}_{i}(t+\Delta t/2)=\mathbf{x}_{i}(t)+\left[\begin{array}[]{ccc}s&c-1&0\\ -(c-1)&s&0\\ 0&0&\omega_{B}\Delta t/2\end{array}\right]\frac{\mathbf{v}_{i}(t)}{\omega_{B}}\;,$$

$$\displaystyle\mathbf{v}_{i}(t+\Delta t/2)=\left[\begin{array}[]{ccc}c&-s&0\\ s&c&0\\ 0&0&0\end{array}\right]\mathbf{v}_{i}(t)\;,$$

$$s=\sin ({\omega }_B\Delta t/2)andc=\cos ({\omega }_B\Delta t/2).$$

$${\mathbf{v}}_i(t+\Delta t)$$

$$\mathbf{E}(t+\Delta t/2,{\mathbf{x}}_i(t+\Delta t/2))=\nabla \phi ({\mathbf{x}}_i(t+\Delta t/2)).$$

$$\Delta x(\tau ;\Delta t)=\sqrt{N^{-1}\sum _{i=1}^N{\left({\mathbf{x}}_i(\tau ;\Delta t)-{\mathbf{x}}_i^{\mathrm{ref}}(\tau )\right)}^2}.$$

$$P_z(\omega )=\sum _{i=1}^N{\left|{\tilde{z}}_i(\omega )\right|}^2+\sum _{i=1}^N{\left|{\tilde{z}}_i(-\omega )\right|}^2.$$

$${\tilde{z}}_i(\omega )=\frac{1}{\tau }{\int }_0^{\tau }{e}^{-i\omega t}{z}_i(t)dt,$$

$$\tilde{z}(l\Delta \omega )$$

$$\approx$$

$$T_{\parallel }=\frac{m_{Be}}{N{k}_B}\sum _{i=1}^N{v}_{z,i}^2,$$

$$T_{\perp }=\frac{m_{Be}}{2N{k}_B}\sum _{i=1}^N\left(v_{x,i}^2+v_{y,i}^2\right)$$

$$T_{\parallel,Doppler}=\frac{1}{3}\frac{\hbar\gamma_{0}}{k_{B}}=$0.29\text{\,}\mathrm{mK}$$

$

$$T_{\parallel }({\omega }_m)=\frac{m{\omega }_m^2\tau }{\pi }{\int }_{res.}d\omega P_z(\omega ),$$

Full text

I Introduction

A Penning trap with a rotating wall potential and Doppler laser cooling can produce stable ultra-cold non-neutral ion plasma crystals with temperatures of a fraction of a millikelvin Jordan et al. (2019); Tan et al. (1995); Mavadia et al. (2013); Ball et al. (2018). Such ultra-cold ion crystals enable interesting research at the forefront of different areas of physics including atomic physics Gruber et al. (2001); Shiga et al. (2011); Andelkovic et al. (2013); Werth et al. (2018); Stutter et al. (2018), quantum optics and metrologyBohnet et al. (2016); Gilmore et al. (2017) , quantum simulation Wang et al. (2013); Britton et al. (2012); Gärttner et al. (2017), and basic plasma physicsDubin and O’neil (1999) , including studies of strongly coupled plasmas that model dense astrophysical matterVan Horn (1991); Ichimaru et al. (1987); Baiko (2009) . For many applications lower ion temperatures are desirable. For example, lower ion temperatures can improve the fidelity of quantum simulations with trapped ions Blatt and Roos (2012) and enable the engineering of stronger interactions by using shorter wavelength spin-dependent forcesBritton et al. (2012) . Colder temperatures also improve the stability of the ion crystal, which advances the prospects for single site optical manipulation and detection in multi-dimensional crystals, akin to what has been achieved with linear ion crystal arraysDebnath et al. (2016) and with an “atom microscope” for neutral atomsBakr et al. (2009). Single site optical manipulation and detection will enable the preparation of complicated entangled states through the implementation of variational quantum simulation protocolsKokail et al. (2018) on large trapped ion crystals.

The study of the thermodynamics of ultra-cold ions in Penning traps is not just a facilitator for certain experiments, it is 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 motionDubin and O’neil (1999). The behavior of energy transfer by nonlinear mode coupling, as well as other heating and cooling mechanisms is subtle and complexJensen et al. (2005); Dubin (2005); Anderegg et al. (2009).

The primary means of cooling the ions in a Penning trap are various forms of laser cooling including Doppler coolingItano and Wineland (1982); Itano et al. (1988); Torrisi et al. (2016); Asprusten et al. (2014), side band coolingStutter et al. (2018), and electromagnetically induced transparency (EIT) cooling Jordan et al. (2019); Shankar et al. (2019). The basic principles for laser cooling ultra-cold ions are similar to that 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 detuningTorrisi et al. (2016). 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 line width 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 shiftsMarquet et al. (2003); McAneny and Freericks (2014) and broadening of resonances resulting in non-trivial frequency spectra.

In this paper we present results from a first principle 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 that used in recent quantum simulation and sensing experimentsBohnet et al. (2016); Gilmore et al. (2017); Wang et al. (2013); Britton et al. (2012); Gärttner et al. (2017). 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 trap configuration and cooling parameters in current and future experiments. The simulation presented in this paper goes beyond prior modeling workWang et al. (2013); Itano et al. (1988); Torrisi et al. (2016) 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 timestep. We then compare the simulation results with both the axial and the planar linear eigenmode analysis in Sec. III. In Section IV, we show Doppler cooled simulations run to steady state and compare with results from equilibrium modelsItano and Wineland (1982); Itano et al. (1988); Torrisi et al. (2016). 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 ultra-cold ions in a Penning trap with rotating wall. For a sketch of a Penning trap and a detailed discussion of the confining fields, see Ref. Wang et al. (2013).

II.1 Trap forces and Hamiltonian

We treat the ions as classical point particles with velocity $\mathbf{v}_{i}$ and position $\mathbf{x}_{i}$. Excluding the cooling laser for now, the motion of the ions is governed by the Hamiltonian

[TABLE]

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

[TABLE]

including the vector potential $\mathbf{A}$ for the homogeneous axial magnetic field $\mathbf{B}=(0,0,B_{z})^{T}$ in the the Penning trap. To be noted, in the following we use the terms axial, out-of-plane and parallel as directions parallel to magnetic field, and planar, in-plane and perpendicular as directions perpendicular to magnetic field. We choose our coordinate system such that the Penning trap magnetic field $\mathbf{B}=\mathbf{\nabla}\times\mathbf{A}$ is parallel to the $z$ axis. With that choice of coordinate system we may choose $\mathbf{A}=-yB_{z}\mathbf{\hat{x}}$ with $\mathbf{\hat{x}}=\mathbf{x}/x$ the unit vector along $\mathbf{x}$. In Eq. (2), $N$ is the number of ions in the trap and $m_{i}$ and $q_{i}$ are the mass and charge of ion $i$. The electrostatic potential

[TABLE]

contains the potential $\varphi_{\rm{trap}}$ for the electrodes in the Penning trap, the rotating wall potential $\varphi_{\rm 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 $\varphi_{\rm{trap}}(\mathbf{x})$ as

[TABLE]

and the time-dependent rotating wall potential $\varphi_{\rm wall}(\mathbf{x},t)$ as

[TABLE]

where $\delta$ is the dimensionless parameter that characterizes the strength of the rotating wall potential to the trapping potential $\varphi_{\rm{trap}}(\mathbf{x})$ and $\theta$ is the azimuthal angle in cylindrical coordinates . Next we transform to the rotating frame

[TABLE]

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

[TABLE]

where

[TABLE]

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 frameDubin and O’neil (1999). 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

[TABLE]

where all ”R” subscripts indicate the rotating frame and $r^{2}_{Ri}=x^{2}_{Ri}+y^{2}_{Ri}$. 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 frameDubin and O’neil (1999). After the system reaches equilibrium, all the particles nearly rigidly rotate as a crystal with the constant frequency $\omega_{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.

II.2 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 byItano et al. (1988)

[TABLE]

where $\gamma_{0}$ is the natural line-width of the atom transition (in radians per second), $S(\mathbf{x})$ is the saturation parameter, and $\Delta(\mathbf{v})=\Delta_{0}-\mathbf{k}\cdot\mathbf{v}$ 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 $\mathbf{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

[TABLE]

where $\rho$ 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 $\Delta t$,

[TABLE]

The velocities and positions needed for computing $\dot{n}_{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, $n_{i}$, as a Poisson random number with mean $\bar{n}_{i}$. Each ion receives a total momentum kick of

[TABLE]

where $\Delta\mathbf{p}_{i,{\rm absorb}}=n_{i}\hbar\mathbf{k}$ and $\Delta\mathbf{p}_{i,{\rm emit}}$ is the recoil corresponding to $n_{i}$ photons scattered in random directions with an isotropic probability distribution. To compute $\Delta\mathbf{p}_{i,{\rm emit}}$, we generate $n_{i}$ vectors of length $\hbar k$ pointing in random directions. The recoil momentum $\Delta\mathbf{p}_{i,{\rm emit}}$ is then obtained by adding up these vectors.

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\gtrsim 1$, 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. For 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}/\gamma_{0}$ being small compared to $1$. The parameter $\omega_{B}=B_{z}q/m$ is the Larmor precession frequency. For the trap and ion parameters considered in this paper we have $\eta\sim 0.3$. 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.

II.3 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 Buneman (1967); Boris and Shanny (1972); Birdsall and Langdon (1985) with an exact matrix rotation for the $\mathbf{v}\times\mathbf{B}$ force term. To advance the positions and velocities of the ions, $\{\mathbf{x}_{i},\mathbf{v}_{i}\}$, from time $t$ to $t+\Delta t$ we apply the following time evolution operator

[TABLE]

where $U(\Delta t)$ is given by

[TABLE]

$U_{0}(\Delta t/2)$ is the time evolution operator corresponding to $H_{0}$ (see Eq. (2)) that advances the state of the ions for a time interval of duration $\Delta 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 Larmor precession,

[TABLE]

where

[TABLE]

The time evolution operator $U_{\rm kick}(t+\Delta t/2;\Delta t)$ captures the forces due to the electrostatic potential (see Eq. (3)) as well as the radiation pressure forces due to laser cooling (Sec. II.B). This operator is time dependent. We evaluate it at the mid-point of the time interval, $t+\Delta t/2$. The operator $U_{\rm kick}$ changes the velocities of the particles only,

[TABLE]

where

[TABLE]

II.4 Convergence with respect to timestep

To verify the precision of our numerical integration procedure, we check that the simulations converge with the expected quadratic rate as the timestep 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. Wang et al. (2013). 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 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)Bohnet et al. (2016); Britton et al. (2012); Sawyer et al. (2014) with a strong homogeneous magnetic field of $B_{z}=$4.4588\text{,}\mathrm{T}, a trap rotation frequency of $\omega_{\rm trap}=2\pi\times$180\text{\,}\mathrm{kHz}, end cap voltages yielding a confining potential of $k_{z}=$9.21\text{,}\mathrm{M}\mathrm{V}\mathrm{/}\mathrm{m}^{2}, and a $1\text{\,}\mathrm{V}$ rotating wall potential yielding $\delta=$3.5\text{\times}{10}^{-4}\text{\,}. For 9Be+ ions, which is assumed throughout this paper, the cyclotron frequency is $qB/m=2\pi\times$7.596\text{,}\mathrm{MHz} and the axial confining frequency parallel to the magnetic field is $2\pi\times$1.58\text{\,}\mathrm{MHz}.

Starting from this initial 2D crystal configuration we integrate the equations of motion for a total duration, $\tau$, using different time step sizes, $\Delta t$, to find the final positions $\mathbf{x}_{i}(\tau;\Delta t)$. For each time step size we compare the final solution with a reference solution, $\mathbf{x}_{i}^{\rm ref}(\tau)=\mathbf{x}_{i}(\tau,$2\text{\times}{10}^{-10}\text{,}\mathrm{s}$)$, computed using a time step size of $\Delta t=$2\text{\times}{10}^{-10}\text{,}\mathrm{s}$$. We compute the error $\Delta x(\tau;\Delta t)$ as the average Euclidean distance of the positions from the reference solution,

[TABLE]

Figure 2 shows the integration errors for an integration time of $\tau=$10\text{,}\mathrm{\SIUnitSymbolMicro s}. The error decreases quadratically for time step sizes below approximately $1.0\text{\times}{10}^{-7}\text{\,}\mathrm{s}$. At a time step size of $\Delta t=$1\text{\,}\mathrm{ns}, the mean position error of the ions is $\Delta x($1\text{,}\mathrm{ns}$)\approx$1\text{,}\mathrm{nm}. During this time interval the ions move on the order of $1\text{\,}\mathrm{mm}$, i.e. the ion motion is integrated with a relative error on the order of $1.0\text{\times}{10}^{-6}\text{\,}$. For time step sizes greater than $\Delta t\approx$1.0\text{\times}{10}^{-7}\text{\,}\mathrm{s} 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. Results presented in this paper were obtained with $\Delta t=$1\text{,}\mathrm{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 experimentBohnet et al. (2016); Britton et al. (2012); Sawyer et al. (2014); Torrisi et al. (2016). Details of the setup of the cooling lasers are described in section IV. The cooling time for both in-plane and out-of-plane degrees of freedom is on the order of $1\text{\,}\mathrm{ms}$. After $5\text{\,}\mathrm{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 $\mathbf{x}_{i}(t_{j})$ at discrete times $t_{j}=j\Delta t$ for $j=0,1,\ldots,N_{\rm sample}-1$.

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

[TABLE]

In this equation, $\tilde{z}_{i}$ is the Fourier transform of the $z$ coordinate of ion $i$,

[TABLE]

with $\tau=N_{\rm sample}\cdot\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,

[TABLE]

Here we have introduced the frequency resolution $\Delta\omega=2\pi/\tau$ and $l=0,1,\ldots,N_{\rm sample}-1$. Frequencies exceeding the Nyquist frequency, $l\Delta\omega>\pi/\Delta t$, wrap to negative frequencies.

We compute spectra $P_{x}$ and $P_{y}$ 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, $x_{R}$ and $y_{R}$.

III.1 Out-of-plane eigenmodes

Figure 3(a) shows $P_{z}$ for 7 ions with a sampling period of $\Delta t=$0.25\text{,}\mathrm{\SIUnitSymbolMicro s} and a total integration time of $\tau=$50\text{\,}\mathrm{ms} corresponding to a frequency resolution of $\Delta\omega/(2\pi)=$20\text{,}\mathrm{Hz}$$. We superimpose the mode frequencies found by linearizing the equations of motion around the ground state as grey vertical lines. For the linear analysis (grey lines) we use a code developed by Wang et al. Wang et al. (2013) 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\text{\,}\mathrm{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 analyzedWang et al. (2013). We obtain the eigenfunction shapes using a notch filter of width $\Delta\omega$ about the mode frequency of interest. We then inverse Fourier transform $\tilde{z}_{i}(\omega)$ and plot the normalized displacement as a function of time. Figure 4 shows a snaphot of the first four (starting from 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 (iii) show the next two highest frequency modes which are tilt modes that would be degenerate except for the weak structural anisotropy produced by the rotating wall potential. The fourth mode, showed in Fig. 4(iv), has two linear nodes. Figure 4 shows similar results as obtained in the eigenfunction plots presented in Ref. Wang et al. (2013), Fig.11.

III.2 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\times B$ modes. The spectra $P_{x}$ and $P_{y}$ are very similar. Therefore we consider just $P_{x}$ 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 grey vertical lines. Figures 5 (b) and  5 (c) show the cyclotron and $E\times 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 non-zero 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. Torrisi et al. (2016). There, it was shown in a fluid model that assumed thermal equilibrium that low ion temperatures close to the Doppler laser cooling limit was possible over a range of experimental conditions. This is in contrast with a Doppler laser cooling model without a rotating wall Itano et al. (1988), where the equilibrium temperature depended sensitively on the details of the set-up. Here, we discuss in more detail, the geometry of the laser beams and the achievement of a Doppler cooled ultra-cold 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 is investigated.

IV.1 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 setupBohnet et al. (2016); Britton et al. (2012); Sawyer et al. (2014); Torrisi et al. (2016). We assume the out-of-plane motion is cooled by two counter-propagating cooling beams along the positive and negative $z$ directions, i.e. parallel and anti-parallel to the trap magnetic field. These two beams have equal intensity and they are detuned by $\Delta_{\parallel}=-\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_{\perp}=-\gamma_{0}/2+k\omega_{R}d$ 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

[TABLE]

[TABLE]

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\text{,}\mathrm{\SIUnitSymbolMicro m} displaced by $d=$5\text{\,}\mathrm{\SIUnitSymbolMicro 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 recoil of the scattered photons, overpowering the planar cooling.

With all the cooling beams on, the system is left to evolve for $10\text{\,}\mathrm{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\text{\,}\mathrm{ms}$.

IV.2 Measured steady-state axial and planar temperatures

We now present steady-state temperature results. The axial temperature $T_{\parallel}$ can be well-diagnosed in the experiments Jordan et al. (2019). Also, the theoretical Doppler cooling limitTorrisi et al. (2016); Itano et al. (1988) gives an estimate for $T_{\parallel}$ that agrees fairly well with temperatures achieved in the NIST experimentBohnet et al. (2016); Britton et al. (2012); Sawyer et al. (2014). The theoretical Doppler cooling limit for motion parallel to the magnetic field is

[TABLE]

where $\gamma_{0}=2\pi\cdot$18\text{,}\mathrm{MHz} is the atomic transition line-width for ${}^{9}\textrm{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 Itano and Wineland ([1982](#bib.bib26)). This limit assumes isotropic scattering, which was also assumed in the numerical simulation of Doppler laser cooling. We estimate the steady-state $T_{\parallel}$ by averaging the results with the lasers on continuously over the final ms, the time interval shown in the inserts of Fig.  [8](#S4.F8). We find steady state temperatures in Fig. [8](#S4.F8) of $T_{\parallel}\approx$0.37\text{\,}\mathrm{mK} and $T_{\perp}\approx$0.92\text{,}\mathrm{mK}. The lowest $T_{\parallel}$ was obtained in a different simulation with $d=0$ (no planar beam offset) and was measured to be 0.28\text{,}\mathrm{mK}$\pm$0.04\text{,}\mathrm{mK}$$, in good agreement with the Doppler cooling limit in Eq. 32. Also, it is important to note that the simulation shows fluctuations in the axial temperature to be approximately $\pm 15\%$.

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 $\omega_{m}$, we can obtain the mode temperature $T_{\parallel}(\omega_{m})$ from the spectrum

[TABLE]

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 $\tau$.

Table 1 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 1 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 $\pm 15\%$ fluctuation, except for the two tilt modes (mode 1 and 2). Higher mode numbers are harder to diagnose using Eq. (26) due to the "grassy" spectrum in the range of $1.1\text{\,}\mathrm{MHz}$ to $1.5\text{\,}\mathrm{MHz}$ seen in Fig. 3.(b). The spectral peaks become very close to each other with increasing mode number making them difficult to differentiate. Additionally, the lowest $T_{\perp}$ is obtained by varying the $W_{y}$ and $d$ and found to be $T_{\perp}\approx$1\text{,}\mathrm{mK} with $W_{y}=d=$5\text{\,}\mathrm{\SIUnitSymbolMicro m}. Such small $W_{y}$ and $d$, compared to the radius of the crystal $R_{c}=$70\text{,}\mathrm{\SIUnitSymbolMicro m}$$, allow the perpendicular beam to generate a modest torque which is straight forwardly balanced by the rotating wallTorrisi et al. (2016). Additionally, the different $T_{\perp}$ and $T_{\parallel}$ demonstrate the temperature anisotropy and weak coupling between axial and planar directions.

V Summary

We have described in detail a direct numerical simulation model of ultra-cold 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 eigenmodes analysis, in both the in-plane and out-of-plane directions. The advantage of a direct numerical simulation is that weak nonlinear coupling between modes is accurately modeled, as well as 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 parametric dependences of the heating and cooling processes in the experiment, as well as 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 existing theory. Axial and planar temperatures were within the range of experimental results and more work is needed to make detailed comparisons with 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=W_{y}=$5\text{,}\mathrm{\SIUnitSymbolMicro m}$\ll R_{c}\approx$70\text{,}\mathrm{\SIUnitSymbolMicro m}$$. The fluid equilibrium model Torrisi et al. (2016) predicts low in-plane temperatures for small $W_{y}$, with a very weak dependence on $d$ if $\Delta_{\perp}=-\gamma/2+k\omega_{R}d$. 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.

VI Acknowledgment

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