Creation of two-dimensional coulomb crystals of ions in oblate Paul traps for quantum simulations

We develop the theory to describe the equilibrium ion positions and phonon modes for a trapped ion quantum simulator in an oblate Paul trap that creates two-dimensional Coulomb crystals in a triangular lattice. By coupling the internal states of the ions to laser beams propagating along the symmetry axis, we study the effective Ising spin-spin interactions that are mediated via the axial phonons and are less sensitive to ion micromotion. We find that the axial mode frequencies permit the programming of Ising interactions with inverse power law spin-spin couplings that can be tuned from uniform to $r^{-3}$ with DC voltages. Such a trap could allow for interesting new geometrical configurations for quantum simulations on moderately sized systems including frustrated magnetism on triangular lattices or Aharonov-Bohm effects on ion tunneling. The trap also incorporates periodic boundary conditions around loops which could be employed to examine time crystals.


Introduction
Using a digital computer to predict the ground state of complex many-body quantum systems, such as frustrated magnets, becomes an intractable problem when the number of spins becomes too large. The constraints on the system's size become even more severe if one is interested in the (nonequilibrium) quantum dynamics of the system. Feynman proposed the use of a quantum-mechanical simulator to efficiently solve these types of quantum problems [1]. One successful platform for simulating lattice spin systems is the trapped ion quantum simulator, which have already been used to simulate a variety of scenarios [2,3,4,5,6,7,8,9]. In one realization [10], ions are cooled in a trap to form a regular array known as a Coulomb crystal and the quantum state of each simulated electron spin can be encoded in the internal states of each trapped ion. Laser illumination of the entire crystal then can be used to program the simulation (spin-spin interactions, magnetic fields, etc.) via coupling to phonon modes, and readout of the internal ion states at the end of the simulation corresponds to a projective measurement of each simulated spin on the measurement basis.
To date, the largest number of spins simulated in this type of device is about 300 ions trapped in an a rotating approximately-triangular lattice in a Penning trap [11]. In that experiment, a spin-dependent optical dipole force was employed to realize an Ising-type spin-spin coupling with a tunable power law behavior. However, the Penning trap simulator was not able to perform certain desirable tasks such as the arXiv:1406.5545v1 [quant-ph] 20 Jun 2014 adiabatic state preparation of the ground state of a frustrated magnet because it did not include a time-dependent transverse magnetic field. The complexity of the Penning trap apparatus also creates a barrier to adoption and therefore does not seem to be as widespread as radio-frequency (RF) Paul traps.
Experiments in linear Paul traps have already performed a wide range of different quantum simulations within a one-dimensional linear crystal [2,3,4,5,6,7,8,9]. The linear Paul trap is a mature architecture for quantum information processing, and quantum simulations in linear chains of ions have benefited from a vast toolbox of techniques that have been developed over the years. Initially, the basic protocol [10] was illustrated in a two-ion trap [2], which was followed by a study of the effects of frustration in a three-ion trap [3]. These experiments were scaled up to larger systems first for the ferromagnetic case [4] and then for the antiferromagnetic case [5]. Effective spin Hamiltonians were also investigated using a Trotter-like stroboscopic approach [6]. As it became clear that the adiabatic state preparation protocol was difficult to achieve in these experiments, ion trap simulators turned to spectroscopic measurements of excited states [7] and global quench experiments to examine Lieb-Robinson bounds and how they change with long-range interactions (in both Ising and XY models) [8,9].
It is therefore desirable to be able to use the demonstrated power of the Paul trap systems to extend access to the 2D physics that is native to the Penning trap systems. However, extension of this technology to higher dimensions is hampered by the fact that most ions in 2D and 3D Coulomb crystals no longer sit on the RF null. This leads to significant micromotion at the RF frequency and can couple to the control lasers through Doppler shifts if the micromotion is not perpendicular to the laser-illumination direction, leading to heating and the congestion of the spectrum by micromotion sidebands.
In an effort to utilize the desirable features of the Paul trap system to study the 2D physics, arrays of tiny Paul traps are being pursued [12,13,14,15,16,17]. It has also been shown that effective higher-dimensional models may be programmed into a simulator with a linear crystal if sufficient control of the laser fields can be achieved [18]. In this article, we study an alternative approach to applying Paul traps to the simulation of frustrated 2D spin lattices. We consider a Paul trap with axial symmetry that forms an oblate potential, squeezing the ions into a 2D crystal. The micromotion in this case is purely radial due to symmetry, and lasers that propagate perpendicular to the crystal plane will therefore not be sensitive to Doppler shifts from micromotion. We study the parameter space of a particular model trap geometry to find the crystal structures, normal modes, and programmable spin-spin interactions for 2D triangular crystals in this trap. We find a wide parameter space for making such crystals, and an Ising spin-spin interaction with widely-tunable range, suitable for studying spin frustration and dynamics on triangular lattices with tens of ions.
In the future, we expect the simulation of larger systems to be made possible within either the Penning trap, or in linear Paul traps that can stably trap large numbers of ions. It is likely that spectroscopy of energy levels will continue to be examined, including new theoretical protocols [19]. Designing adiabatic fast-passage experiments along the lines of what needs to be done for the nearest-neighbor transverse field Ising model [20] might improve the ability to create complex quantum ground states. Motional effects of the ions are also interesting, such as tunneling studies for motion in the quantum regime [21]. It is also possible that novel ideas like time crystals [22,23] could be tested (although designing such experiments might be extremely difficult).
The organization of the remainder of the paper is as follows: In Sec. 2, we introduce and model the potential energy and effective trapping pseudopotential for the oblate Paul trap. In Sec. 3, we determine the equilibrium positions and the normal modes of the trapped ions, with a focus on small systems and how the geometry of the system changes as ions are added in. In Sec. 4, we show representative numerical results for the equilibrium positions and the normal modes. We then show numerical results for the effective spin-spin interactions that can be generated by a spin-dependent optical dipole force. In Sec. 5, we provide our conclusions and outlook.

Oblate Paul Trap
The quantum simulator architecture we study here is based on a Paul trap with an azimuthally symmetric trapping potential that has significantly stronger axial confinement than radial confinement, which we call an "oblate Paul trap" since the resulting effective potential resembles an oblate ellipsoid. As we show below, this can create a Coulomb crystal of ions that is a (finite) two-dimensional triangular lattice. 2D Coulomb crystals in oblate Paul traps have been studied since the 1980's and were used, for instance, by the NIST Ion Storage Group [24] to study spectroscopy [25,26], quantum jumps [27], laser absorption [28,29] and cooling processes [30]. 2D crystals in oblate Paul traps have also been studied by other groups both experimentally [31] and theoretically [32,33,34,35,36].
The particular oblate Paul trap we study has DC "end cap" electrodes above and below a central radio-frequency (RF) ring, as depicted in Fig. 1. The trap we propose uses modern microfabrication and lithography technology (manufactured by Translume, Ann Arbor, MI) to realize the DC end cap electrodes as surface features on a monolithic fused silica substrate, providing native mechanical indexing and easier optical access to the ions than discrete end cap traps. Similar to Penning traps, oblate Paul traps can be used to study frustration effects when the lattice of ions has multiple rings. However in a Penning trap, the lattice of ions is rotating and the ions are in a strong magnetic field, which can add significant complications. In an oblate Paul trap, the lattice of ions is stationary except for the micromotion of each ion about its equilibrium positions (which is confined to the plane of the crystal by symmetry) and the qubits can be held in nearly zero magnetic field, permitting the use of the m = 0 "clock state" used in linear Paul trap quantum simulators [3]. For trapped ions in a crystal that is a single polygon (N = 3, 4, or 5), we can study periodic boundary conditions applied to a linear chain of trapped ions in the linear Paul quantum simulators. Oblate Paul traps can also potentially be used to perform experiments that are similar to those recently exploring the Aharonov-Bohm effect [21] with more ions.
It is well known that Maxwell's equations forbid the possibility to use a static electric field to trap charged particles in free space through Earnshaw's theorem. However, a static electric field can create a saddle-point, which confines the charged particles in some directions and deconfines them in the other directions. A static electric field that provides a saddle-point is where A is a constant andê i are the perpendicular unit vectors with i = 1, 2, 3. Using a static electric field with a saddle-point, both Penning and radio-frequency (RF) Paul traps have successfully trapped charged particles in free space by applying an additional field. In the Penning trap, one applies a strong uniform magnetic field, such that the charged particles are confined to a circular orbit via the Lorentz force, qv × B/c. The RF Paul trap applies a time-varying voltage to its electrodes, which produce a saddle potential that oscillates sinusoidally as a function of time. This rapid change of sign allows for certain ions to be trapped because for particular charge to mass ratios, the effective focusing force is stronger than the defocussing force.
If the ions remain close to the nulls of the potential, then the micromotion of the ions is small, and it is a good approximation to describe the system via a static pseudopotential that approximates the trapping effect of the time-varying potential. We use the numerical modeling software Comsol to simulate this effective pseudopotential that arises from applying a time-varying voltage to the RF ring and additional DC voltages on the other electrodes. The effective total potential energy, V (x 1 ,x 2 ,x 3 ), of an ion in the oblate Paul trap can be approximated bỹ where ψ(x 1 ,x 2 ,x 3 ) is the effective pseudopotential due to the RF fields and φ(x 1 ,x 2 ,x 3 ) is the additional potential due to the DC voltage applied on the top and bottom electrodes and the RF ring. The resulting pseudopotential at a certain point in space will depend upon the RF frequency, Ω RF , and the RF electric field amplitude, E o,RF (x 1 ,x 2 ,x 3 ), at that point [37]and is given by which depends on the charge, q, and the mass, m, of the particular ion being trapped. After simulating the field using Comsol, we find that the electric field amplitude from the RF field near the trap center can be approximated by where V o,RF is the amplitude of the RF voltage. Plugging this into Eq. (3) yields where r o = 512 µm is a fitting parameter, that is determined by grounding the top and bottom electrodes and numerically modeling the square of the RF electric field amplitude, as shown in Fig. 2. We calculate the DC electric field as having 3 contributions: one from the voltage applied to the RF ring, φ r (x 1 ,x 2 ,x 3 ), one from the voltage applied to the top electrodes, φ t (x 1 ,x 2 ,x 3 ) and one from the bottom where V r is the DC voltage on the ring. We numerically model the electrostatic potential due to the DC voltage applied to the either the top or bottom electrodes, as shown in Fig. 3. We find that near the trap center, the numerical results for the electrostatic potential produced by a voltage of V t,b the top or bottom electrodes is reasonably modeled by the polynomial We can use the results from Eqs. (5-7) in Eq. (2) to yield the final effective potential energy of an ion in this oblate Paul trap Since the effective potential energy is just a function ofx 2 1 +x 2 2 , it is rotationally symmetric around theê 3 -axis and we would expect there to be a zero frequency rotational mode in the phonon eigenvectors. That mode can be lifted from zero by breaking the symmetry, which can occur by adding additional fields that do not respect the cylindrical symmetry, and probably occur naturally due to imperfections in the trap, the optical access ports, stray fields, etc.

Equilibrium structure and normal modes
Using Eq. (8) (the calculated pseudopotential), we solve for the equilibrium structure of the crystal in the standard way. We first construct an initial trial configuration for the ions and then minimize the total potential energy of the oblate Paul trap (including the trap potential and the Coulomb repulsion between ions), as summarized in Eq. (9); MatLab is used for the nonlinear minimization with a multidimensional Newton's method. We rewrite the total potential energy of the oblate Paul trap in a conventional form (up to a constant) viã where k e = 1/4π o . Herex in is the i th component of the n th ion's location and The frequencies in Eq. (9) are defined via We will express all distances in terms of a characteristic length, l o , which satisfies and we will work with dimensionless coordinates x =x/l o when calculating the equilibrium positions. Furthermore, we measure all frequencies relative to ω ψ,3 . The normalized frequencies are for i = 1, 2, β r,3 = ω r,3 /ω ψ,3 , β t,3 = ω t,3 /ω ψ,3 , and β b,3 = ω b,3 /ω ψ,3 . The dimensionless total potential energy becomes where we have defined x o,t = a 2 /(2l o b t ) and x o,b = a 2 /(2l o b r ). To find the equilibrium positions, we use the gradient of the total potential energy and numerically minimize the total potential energy using a multidimensional Newton's method. The gradient of Eq. (9) is The force on ion m in theî direction will be − ∇V ·ê im . We seek the solution in which all ions lie in a plane parallel to theê 1 −ê 2 plane, such that x 3m = x 3 for all m ∈ [1, N ]. As a result of this condition, x 3n = x 3m , and there is no x 3 contribution to the Coulomb potential term. The value of x 3 is determined by setting theê 3m term equal to zero in Eq. (17) and is given by the condition Using x 3 , the ion equilibrium positions are numerically obtained when all 3N components of the force on each ion are zero, which is given by ∇V | equilib. = 0. After the equilibrium positions {x in , n = 1, . . . , N, i = 1, 2, 3} are found, we expand the total potential about the equilibrium positions up to quadratic order The nonzero terms of the expansion are the zeroth order and the quadratic terms; the first order term is zero because the equilibrium position is defined to be where the gradient of the total potential is zero, however the zeroth term is also neglected since it is a constant. We calculate the Lagrangian of the trapped ions using the quadratic term of the Taylor expanded total potential, with q in being the dimensionless displacement from the equilibrium position for the n th ion in the i th direction. The expanded Lagrangian becomes where K ij mn represents the elements of the effective spring constant matrices which are given by   3 and r mn = (x 1m − x 1n ) 2 + (x 2m − x 2n ) 2 is the planar interparticle distance between ions n and m. Note that motion in the 3-direction (axial direction) is decoupled from motion in theê 1 −ê 2 plane.
After applying the Euler-Lagrange equation to Eq. (20) and substituting the eigenvector solution q im = Re(b α im e iωαt ), we are left to solve the standard eigenvalue equation There are two sets of normal modes: eigenvectors of the N ×N matrix K 33 yield the "axial" modes (those corresponding to motion perpendicular to the crystal plane) and eigenvectors of the 2N × 2N matrix K ij , i, j ∈ [1, 2], yield the "planar" modes (those corresponding to ion motion in the crystal plane).

Results
Now that we have constructed the formal structure to determine the equilibrium positions, phonon eigenvectors, and phonon frequencies, and we have determined the total pseudopotential of the trap, we are ready to solve these systems of equations to determine the expected behavior of the trapped ions. We present several numerical examples to illustrate the equilibrium structure, eigenvalues of the normal modes, and the effective spin-spin interaction J mn for the axial modes with a detuning of the spin-dependent optical dipole force above the axial center-of-mass phonon frequency, ω CM . We use an ytterbium ion with mass m = 171u, where u is the atomic mass unit, and a positive charge q = e. For the frequency of the RF voltage, we use Ω RF = 2π × 35MHz and the amplitude of the potential applied to the RF ring is V o,RF ≈ 500V. The DC voltage applied to the RF ring and to the top and bottom electrodes will be |V DC | ≤ 100V. We work in a region where the trapped ion configurations are stable. The ion crystal is stable only when both β 1 and β 2 are real and nonzero, as defined in Eq. (14) and this region is shown in Fig. 4, which depends on the voltages applied to the RF ring and to the top and bottom electrodes. We work with ion crystals that contain up to 20 ions.

Equilibrium configurations
A single ion will sit in the center of the trap. As more ions are added, because the ions repel each other, a "hard core" like structure will form, starting with rings of ions until it becomes energetically more stable to have one ion in the center of the ring, and then additional shells surrounding it, and so on. We find that as we increase the number of ions, the single ring is stable for N = 3, 4, and 5. Increasing N further creates more complex structures. We show the common equilibrium configurations for N = 5, 10, 15, 20 with DC voltages of V r = 46.3V and V t = V b = 50V, in Fig. 5. As mentioned above, N = 5 is the last configuration that is comprised of a single ring of ions, as depicted in Fig. 5(a). The N = 5 configuration is ideal to use in order to study when periodic boundary conditions are applied to the linear chain, this is due to the configuration being in a single ring. For configurations with N = 6 through 8, the additional ions are added to an outer ring. When N = 9 the additional ions are added to the center. In Fig. 5(b), N = 10 is the first configuration the forms a ring in the center, with three ions. The equilibrium configuration of N = 15 is the maximum number to have two rings, as shown in Fig. 5(c). Ion configurations with N > 15 have a single ion at the center, as an example of this, we show N = 20 in Fig. 5(d). The common configurations for N > 5 are nearly formed from triangular lattices (up to nearest neighbor) and this could be used to study frustration in the effective spin models (except, of course, that due to the finite number of ions there are many cases where the coordination number of an interior ion is not equal to 6, as seen in Fig. 5). The shape of all of these clusters for small N agree with those found in Ref. [33], except for N = 10, 12, and 14, which have small differences due to the different potential that describes the oblate Paul trap from the potential used in [33]. We next show the dependence of the equilibrium positions on the DC voltages applied to the RF ring and independently applied to the top and bottom electrodes. We fix N = 5. As each DC voltage is independently varied, the shape of the equilibrium configuration for N = 5 remains the same and only the distances between ions change, as shown in the four cases in Fig. 6.

Normal modes
After determining the equilibrium positions, we can find the spring constants and then solve the eigenvalue problem to find the normal modes. Note that due to rotational symmetry, there always is a zero frequency planar mode corresponding to the free rotation of the crystal. In an actual experiment, however, we expect that the rotational symmetry of the trap will be broken by stray fields, the radial optical access tunnels, imperfections in the electrodes, etc., so that mode will be lifted from zero.
We show the eigenvalues of the normal modes for N = 5 in Fig. 7. The axial phonon frequencies decrease as the DC voltage on the RF ring increases and the planar phonon frequencies decrease as the DC voltage on the RF ring decreases. For the majority of the combinations of V t and V b the axial phonon frequencies lie in a narrow band that is separated from the planar mode frequencies, which also lie in a narrow band. As V r increases the axial band broadens and eventually overlaps the planar band, which is also broadening. When the axial band has an eigenvalue that goes soft, the system is no longer stable within one plane (which is the equivalent of the zig-zag transition in the linear Paul trap). When V t = V b = 0 the initial clustering of the axial modes and planar modes is not present, as shown in Fig. 7(a).

Ising spin-spin interaction
The ions in our trap have two hyperfine states that are separated by a frequency ω 0 . Three laser beams with two beatnotes at frequencies ω 0 ± µ will illuminate the ions, selectively exciting phonon modes as described in [3]. In this case, we choose the laser beams to propagate along the ±ê 3 direction, as defined in Fig. 1b, such that the laser beams are insensitive to the micromotion which is entirely radial. The phonon modes are excited in a spin-dependent way to generate effective spin-spin interactions which can be described by the Ising spin coupling matrix, J mn where σ z i is the Pauli spin matrix of ion i in theê 3 -direction and we have neglected the time-dependent terms of the spin couplings J mn . The explicit formula for J mn is [38] where the coefficient J 0 , depends on the carrier transition Rabi frequency (Ω), the difference in wavevector between the laser beams (δk), the ion mass (m), and the frequency of the center-of-mass mode, (ω CM ), and is given by It is expected that if one detunes, µ, to be larger than the center-of-mass mode frequency, then one can generate long-range spin-spin couplings that vary as a power law from 0 to 3, as they decay with distance. Hence, we fit the spin-spin couplings to a power law in distance as a function of detuning in Fig. 8 for N = 10 and 20. Note that our system is still rather small, so there are likely to be finite size effects that modify the simple power law behavior.

Quantum motional effects
The trap could also be used to examine different types of quantum motional effects of ions, similar to the recent work on the Aharonov-Bohm effect [21]. In order to examine such effects, one would need to cool the system to nearly the ground state. This can be accomplished by including Raman side-band cooling after Doppler cooling the system for all modes except the soft rotational mode, at least when the potential is large enough that the mode frequencies are sizeable. To cool the rotational mode, one would need to add a perturbation to the system that lifts the phonon mode frequency, side-band cool it, and then adiabatically reduce the frequency by removing the perturbation. This procedure will cool off that phonon mode, which can yield quite small quanta in it [21]. Once the system has been prepared in this state, then quantum tunneling effects, or coherent motional effects could be studied in the trap for a range of different ion configurations. It might also be interesting to extend these types of studies to cases where the ions no longer lie completely in one plane, but have deformed into a full three-dimensional structure (as long as the larger micromotion does not cause problems). Finally, many of these ideas would need to be used if one tried to examine time crystals, especially the cooling of the rotational mode to be able to see quantum effects.

Conclusion
In this work, we have studied 2D ion crystals in an oblate Paul trap for use in quantum simulations. With this system, one can trap a modest number of ions in 2D planar structures that are likely to be highly frustrated without needing a Penning trap, providing a controlled way to study the onset of frustration effects in quantum simulations. We calculated the equilibrium positions and the phonon frequencies for the proposed oblate Paul trap over its stable region. The equilibrium positions with N ≤ 5 form a single ring configuration and could potentially be used to study periodic boundary conditions and the Aharonov-Bohm effect when N = 4 or 5 (and possibly time crystals). Once N > 5, the equilibrium configurations have multiple rings that are nearly formed from triangular lattices. One can generate an effective Ising Hamiltonian by driving axial modes with a spin-dependent optical dipole force. When detuning is to the blue of the axial center-of-mass mode, the spin-spin coupling, J mn , has an approximate power law that is within the expected range of 0 to 3. In the future, as this trap is tested and performs simulations of spin models with ions, the work presented here will be critical to determining the parameters of the Hamiltonian and for selecting the appropriate configurations to use in the different simulations.